非线性测井解释模型求解问题-以西门度公式为例

问题

对于含水饱和度等非线性方程求解问题，无法适用基础运算进行，因此在计算时可采用逐次搜索法求解的近似解。

即采用某采样间隔进行逐点计算，当相邻两点方程求解答案结果相乘小于零则可以认为方程求解结果为该两点之间的值，答案精度取决于步长，程序运行时间也取决于步长。

西门度公式

$\frac{1}{R_t} = \frac{V^c_{sh}}{R_{sh}} * S^{n/2}_w + \frac{\Phi^{m}}{a * R_w * (1-V^c_{sh})} * S^n_w$

要求解 $S_w$ ，解方程组可得

$\frac{R_t * \Phi^{m}}{a * R_w * (1-V^c_{sh})} * S^n_w + \frac{R_t * V^c_{sh}}{R_{sh}} * S^{n/2}_w = 1$

可以看做如下方程

$F(x)=A*x^n+B*x^{n/2}+C$ $A = \frac{R_t * \Phi^{m}}{a * R_w * (1-V^c_{sh})}$ $B = \frac{R_t * V^c_{sh}}{R_{sh}}$ $C=-1$

那么求解西门度方程的问题就变成了非线性方程 $F(x)=0$ 的解(0≤x≤1)

当 $F(x)*F(x+1)<=0$ 说明方程解为[x,x+1]

为节省运算内存，可以只关注F(x)和F(x+1)结果的正负值(-1,1)来进行F(x)*F(x+1)判断

程序输入

输入曲线
VSH：泥质含量曲线，%；
RT：电阻率曲线，Ω·m；
POR： $\Phi$ ，孔隙度曲线，%；
输入参数
RSH2：西门度公式泥质电阻率，Ω·m，(10)；
SOA2：西门度公式地层参数A，无量纲，(0.8578)；
SOM2：西门度公式胶结指数M，无量纲，(1.358)；
SON2：西门度公式饱和度指数N，无量纲，(3.502)；
SRW2：西门度公式地层水电阻率Rw，Ω·m，(0.035)；
SOC2：西门度公式泥质含量幂乘指数，无量纲(1)；
SW2D：搜索区间步长；
程序内部参数
SWFT：程序循环体判断结束标志符；
SW2A： $\frac{R_t * \Phi^{m}}{a*b*R_w (1-V_{sh})}$ ， $F(x)$ 参数A；
SW2B： $\frac{R_t * V^c_{sh}}{R_{sh}}$ ， $F(x)$ 参数B；

代码实现

     
C     SW默认值      
      SW=-9999.99
C     SW起始值
      SWQ=SW2D
      
C     FX 参数A
      SW2A=REAL(RT*((POR/100)**SOM2)/SOA2/SRW2/(1-(VSH/100)**SOC2))
C     FX 参数A
      SW2B=REAL(RT*((VSH/100)**SOC2)/RSH2)

C     F（Xn）
      FX1=REAL(SW2A*(SWQ**SON2)+SW2B*(SWQ**(SON2/2))-1)
      WRITE(*,*)'FX1=',FX1
C     F（Xn）判断变量，节省计算内存
      FX11=0
      IF(FX1.GE.0) THEN
      	FX11=1
      ELSE
      	FX11=-1
      END IF
      
C     F（Xn+1）      
      FX2=REAL(SW2A*((SWQ+SW2D)**SON2)+SW2B*((SWQ+SW2D)**(SON2/2))-1)
C     F（Xn+1）判断变量，节省计算内存
      FX22=0
      IF(FX2.GE.0) THEN
      	FX22=1
      ELSE
      	FX22=-1
      END IF
      
C     临界判断变量，当SWFT<0证明为取值区间
      SWFT=REAL(FX11*FX22)

C     执行循环体，按步长挪动区间
      DO WHILE(SWFT.GT.0) 
      	SWQ=SWQ+SW2D
      	
      	FX11=0
      	FX1=REAL(SW2A*(SWQ**SON2)+SW2B*(SWQ**(SON2/2))-1)
      	WRITE(*,*)'FX1=',FX1
      	IF(FX1.GE.0) THEN
      		FX11=1
      	ELSE
      		FX11=-1
      	END IF
      	
      	FX22=0
     	FX2=REAL(SW2A*((SWQ+SW2D)**SON2)+SW2B*((SWQ+SW2D)**(SON2/2))-1)
     	WRITE(*,*)'FX2=',FX2
      	IF(FX2.GE.0) THEN
      		FX22=1
      	ELSE
      		FX22=-1
      	END IF
   
      	SWFT=REAL(FX11*FX22)
C      	WRITE(*,*)'SWFT=',SWFT
      	
C     	SW范围是0-1	
      	IF(SWQ.GT.1) exit
      END DO
      
      SW=REAL(SWQ*100)