【程序杂谈】岩石物理模型算法程序

以前的作业。

DEM模型

function[kd,ud]=dem(km,um,kf,uf,alpha,phy)
dp=0.0001;
j=0;
k(1)=km;
u(1)=um;
 for p=0:dp:phy
      j=j+1;
   [P,Q]=pqvalue(k(j),u(j),alpha,kf,uf);
   k(j+1)= k(j) + (kf-k(j))*dp*P/( 1-p);
   u(j+1)= u(j)+(uf-u(j))*dp*Q/(1-p);
 end
 kd=k(j+1);
 ud=u(j+1);

KT模型

function [k,u]=kt(km,um,kf,uf,alpha,phy)
%KT模型的实现程序，[k,u]=Kuster(km,um,kf,uf,a,phy,vm)
%km,um -----背景介质
%kf,uf---------包含物介质
%phy------------包含物体积分数
%a------------孔隙纵横比
%vm---------泊松比
%k ,u -----返回干岩石的体模量和剪切模量
[P,Q]=pqvalue(km,um,alpha,kf,uf);
sm=um./6.*(9.*km+8.*um)./(km+2.*um);
k=(4.*km.*um + 3.*km.^2 + 4.*P.*kf.*um.*phy - 4.*P.*km.*um.*phy)./(3.*km + 4.*um - 3.*P.*kf.*phy + 3.*P.*km.*phy);
u=(sm.*um + um.^2 + Q.*sm.*uf.*phy - Q.*sm.*um.*phy)./(sm + um - Q.*uf.*phy + Q.*um.*phy);

SC模型

function [k,u]=sc(km,um,ki,ui,alpha,phy)
% km,um -----背景介质
% ki,ui---------包含物介质
% phy------------体积分数
% alpha------------孔隙纵横比
[P,Q]=pqvalue(km,um,alpha,ki,ui);
k=km+phy*(ki-km)*P;
u=um+phy*(ui-um)*Q;
岩石骨架模量
function [km,um]=vrh(k1,u1,v1,k2,u2,v2)
%  Hill平均估算模量
kv=v1*k1+v2*k2;
uv=v1*u1+v2*u2;
kr=1/(v1/k1+v2/k2);
ur=1/(v1/u1+v2/u2);
km=0.5*(kv+kr);
um=0.5*(uv+ur);

包含物任意宽高比 

function [p,q]=pqvalue(km,um,e,ki,ui)
%求扁圆球体的P，Q值
%function [P, Q]=ellipsess(Km,um,Ki,ui,arf)
%参数说明：1：基质体积模量，2：基质剪切模量。它们是由VRH模量得到。
%         3：每种矿物体积模量，4：剪切模量。它们是由实验得到或者经验值。
%         5：椭球形的扁度，
%         P，Q为返回的方程参数。
A=ui/um-1;
B=1/3*(ki/km-ui/um);
R=(3*um)/(3*km+4*um);
if e==1
    e=0.999;
end
if e>1
    o=e/(e^2-1)^(3/2)*(e*(e^2-1)^(1/2)-acosh(e));
else
   o=e/(1-e^2)^(3/2)*(acos(e)-e*(1-e^2)^(1/2));
end
f=e^2/(1-e^2)*(3*o-2);
F1=1+A*(3/2*(f+o)-R*(3/2*f+5/2*o-4/3));
F2=1+A*(1+3/2*(f+o)-(R/2)*(3*f+5*o))+B*(3-4*R)+A/2*(A+3*B)*(3-4*R)*(f+o-R*(f-o+2*o^2));
F3=1+A*(1-(f+3/2*o)+R*(f+o));
F4=1+A/4*(f+3*o-R*(f-o));
F5=A*(-f+R*(f+o-4/3))+B*o*(3-4*R);
F6=1+A*(1+f-R*(f+o))+B*(1-o)*(3-4*R);
F7=2+A/4*(3*f+9*o-R*(3*f+5*o))+B*o*(3-4*R);
F8=A*(1-2*R+f/2*(R-1)+o/2*(5*R-3))+B*(1-o)*(3-4*R);
F9=A*((R-1)*f-R*o)+B*o*(3-4*R);
T1=3*F1/F2;
T2=2/F3+1/F4+(F4*F5+F6*F7-F8*F9)/(F2*F4);
p=1/3*T1;
q=1/5*T2;