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

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

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

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

问题

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

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

西门度公式

$$ \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)=Ax^n+Bx^{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}}{abR_w (1-V_{sh})} $$,$$ F(x) $$ 参数A;
    SW2B:$$ \frac{R_t * V^c_{sh}}{R_{sh}} $$,$$ F(x) $$ 参数B;

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
     
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)
文章作者: HibisciDai
文章链接: http://hibiscidai.com/2021/10/30/非线性测井解释模型求解问题-以西门度公式为例/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 HibisciDai