📄 poromodel.m
字号:
function poromodel
t0=0;
g=30;
t=2000;
b=0;
p11=2167;
p22=191;
p12=-83;
P=2.0332*10^10;
Q=0.0953*10^10;
RR=0.0331*10^10;
N=0.684*10^10;
q=0.1;
vf=3210;
vs=1180;
vr=1790;
arfa=-0.337*10^(-10);
bata=0.719*10^(-11);
af=1.1;
as=-16.2;
R11=p11/(p11*p22-p12^2);
R12=p12/(p11*p22-p12^2);
R22=p22/(p11*p22-p12^2);
dt=0.1*10^(-3);dx=2.5;m=85;
mup=80;mdown=80;mleft=80;mright=80;
ipsi=2*pi^2*g^2;
tt0=1/g;
ep=0.8;
aa=max([mleft,mright,mup,mdown,1]);
A=(0.015*ep)^2/aa^2;
for i=1:m+mleft+mright
for j=1:m+mup+mdown
if i<=mleft||i>=(m+mleft)||j<=mup||j>=(m+mup)
x1=max([mleft-i,i-m-mleft,mup-j,j-m-mup]);
tp(i,j)=exp(-A*x1^4);
else
tp(i,j)=1.0;
end
end
end
for i=1:m+mleft+mright
for j=1:m+mup+mdown
H(i,j)=0;
Vs(i,j)=0;
Vf(i,j)=0;
Us(i,j)=0;
Uf(i,j)=0;
R(i,j)=0;
T(i,j)=0;
S(i,j)=0;
end
end
for k=1:t
for i=3:m-1+mleft+mright
for j=3:m-1+mup+mdown
Us(i,j)=Us(i,j)+dt/dx*R22*(9/8*(R(i,j)-R(i-1,j))-1/24*(R(i+1,j)-R(i-2,j)))...
-dt/dx*R12*(9/8*(S(i,j)-S(i-1,j))-1/24*(S(i+1,j)-S(i-2,j)))...
+dt/dx*R22*(9/8*(H(i,j)-H(i,j-1))-1/24*(H(i,j+1)-H(i,j-2)));
Uf(i,j)=Uf(i,j)-dt/dx*R12*(9/8*(R(i,j)-R(i-1,j))-1/24*...
(R(i+1,j)-R(i-2,j)))+dt/dx*R11*(9/8*(S(i,j)-S(i-1,j))-1/24*(S(i+1,j)-S(i-2,j)))...
-dt/dx*R12*(9/8*(H(i,j)-H(i,j-1))-1/24*(H(i,j+1)-H(i,j-2)));
end
end
for i=2:m-2+mleft+mright
for j=2:m-2+mup+mdown
Vs(i,j)=Vs(i,j)+dt/dx*R22*(9/8*(H(i+1,j)-H(i,j))...
-1/24*(H(i+2,j)-H(i-1,j)))+dt/dx*R22*(9/8*(T(i,j+1)-T(i,j))-...
1/24*(T(i,j+2)-T(i,j-1)))-dt/dx*R12*(9/8*(S(i,j+1)-S(i,j))-...
1/24*(S(i,j+2)-S(i,j-1)));
Vf(i,j)=Vf(i,j)-dt/dx*R12*(9/8*(H(i+1,j)-H(i,j))...
-1/24*(H(i+2,j)-H(i-1,j)))-dt/dx*R12*(9/8*(T(i,j+1)-T(i,j))-...
1/24*(T(i,j+2)-T(i,j-1)))+dt/dx*R11*(9/8*(S(i,j+1)-S(i,j))-...
1/24*(S(i,j+2)-S(i,j-1)));
end
end
for i=3:m-1+mleft+mright
for j=2:m-2+mup+mdown
H(i,j)=H(i,j)+dt/dx*N*(9/8*(Us(i,j+1)-Us(i,j))...
-1/24*(Us(i,j+2)-Us(i,j-1))+9/8*(Vs(i,j)-Vs(i-1,j))...
-1/24*(Vs(i+1,j)-Vs(i-2,j)));
end
end
for i=2:m-2+mleft+mright
for j=3:m-1+mup+mdown
R(i,j)=R(i,j)+dt/dx*P*(9/8*(Us(i+1,j)-Us(i,j))-1/24*...
(Us(i+2,j)-Us(i-1,j)))+dt/dx*(P-2*N)*(9/8*(Vs(i,j)-...
Vs(i,j-1))-1/24*(Vs(i,j+1)-Vs(i,j-2)))+...
dt/dx*Q*(9/8*(Uf(i+1,j)-Uf(i,j))-1/24*(Uf(i+2,j)-Uf(i-1,j)))...
+dt/dx*Q*(9/8*(Vf(i,j)-Vf(i,j-1))-1/24*(Vf(i,j+1)-Vf(i,j-2)));
T(i,j)=T(i,j)+dt/dx*(P-2*N)*(9/8*(Us(i+1,j)-Us(i,j))-1/24*...
(Us(i+2,j)-Us(i-1,j)))+dt/dx*P*(9/8*(Vs(i,j)-...
Vs(i,j-1))-1/24*(Vs(i,j+1)-Vs(i,j-2)))+...
dt/dx*Q*(9/8*(Uf(i+1,j)-Uf(i,j))-1/24*(Uf(i+2,j)-Uf(i-1,j)))...
+dt/dx*Q*(9/8*(Vf(i,j)-Vf(i,j-1))-1/24*(Vf(i,j+1)-Vf(i,j-2)));
S(i,j)=S(i,j)+dt/dx*Q*(9/8*(Us(i+1,j)-Us(i,j))-1/24*(Us(i+2,j)-Us(i-1,j))...
+9/8*(Vs(i,j)-Vs(i,j-1))-1/24*(Vs(i,j+1)-Vs(i,j-2)))...
+dt/dx*RR*(9/8*(Uf(i+1,j)-Uf(i,j))-1/24*(Uf(i+2,j)-Uf(i-1,j))...
+9/8*(Vf(i,j)-Vf(i,j-1))-1/24*(Vf(i,j+1)-Vf(i,j-2)));
end
end
%if k<61
% R(mleft+m/2,mup+m/2)=R(mleft+m/2,mup+m/2)+(1-2*(pi*g*(k-40)/1000)^2)*exp(-(pi*g*(k-40)/1000)^2)...
% -(1-2*(pi*g*(k-40-1)/1000)^2)*exp(-(pi*g*(k-40-1)/1000)^2);
% end
% R(mleft+m/2,mup+m/2)=R(mleft+m/2,mup+m/2)-2*ipsi*(exp(1/(2*ipsi)))^.1/2 *exp(-ipsi*(k-tt0)^2)...
% *(k-tt0);
% R(mleft+m/2,mup+m/2)=R(mleft+m/2,mup+m/2)+2*(ipsi/10^6)*(exp(1/(2*ipsi/10^6)))^0.5*exp(-ipsi/10^6*(k-tt0).^2).*(k-tt0);
% R(mleft+m/2,mup+m/2)=R(mleft+m/2,mup+m/2)-2*ipsi*(exp(1/(2*ipsi)))^0.5*exp(-ipsi*((k+1/2)/1000-tt0)^2)*((k+1/2)/1000-tt0);
% R(mleft+m/2,mup+m/2)=R(mleft+m/2,mup+m/2)-2*ipsi*(exp(1/(2*ipsi)))^0.5*exp(-ipsi*((k+1/2)/1000-tt0)^2)*((k+1/2)/1000-tt0);
% Ricker sourse
% R(mleft+1,mup+1)=R(mleft+1,mup+1)-2*ipsi*(exp(1/(2*ipsi)))^0.5*exp(-ipsi*((k+1/2)/10000-tt0)^2)*((k+1/2)/10000-tt0);
% T(mleft+1,mup+1)=T(mleft+1,mup+1)-2*ipsi*(exp(1/(2*ipsi)))^0.5*exp(-ipsi*((k+1/2)/10000-tt0)^2)*((k+1/2)/10000-tt0);
% S(mleft+1,mup+1)=S(mleft+1,mup+1)-2*ipsi*(exp(1/(2*ipsi)))^0.5*exp(-ipsi*((k+1/2)/10000-tt0)^2)*((k+1/2)/10000-tt0);
% Gaussian sourse
% R(mleft+3,mup+3)=R(mleft+3,mup+3)+exp(-ipsi*(k/10000-tt0)^2);
% T(mleft+3,mup+3)=T(mleft+3,mup+3)+exp(-ipsi*(k/10000-tt0)^2);
R(mleft+43,mup+43)=R(mleft+43,mup+43)-(1-q)*(exp(-ipsi*(k/10000-tt0)^2)-exp(-ipsi*((k-1)/10000-tt0)^2));
T(mleft+43,mup+43)=T(mleft+43,mup+43)-(1-q)*(exp(-ipsi*(k/10000-tt0)^2)-exp(-ipsi*((k-1)/10000-tt0)^2));
S(mleft+43,mup+43)=S(mleft+43,mup+43)-q*(exp(-ipsi*(k/10000-tt0)^2)-exp(-ipsi*((k-1)/10000-tt0)^2));
% H(mleft+3,mup+3)=-(1-q)*exp(-ipsi*(k/10000-tt0)^2);
% S(mleft+3,mup+3)=-exp(-ipsi*(k/10000-tt0)^2);
% Gaussian sourse acting on replacement
% Uf(mleft+3,mup+3)=(1-q)*exp(-ipsi*(k/10000-tt0)^2);
% Vf(mleft+3,mup+3)=(1-q)*exp(-ipsi*(k/10000-tt0)^2);
% Us(mleft+3,mup+3)=q*exp(-ipsi*(k/10000-tt0)^2);
% Vs(mleft+3,mup+3)=q*exp(-ipsi*(k/10000-tt0)^2);
for i=1:m+mleft+mright
for j=1:m+mup+mdown
Uf(i,j)=Uf(i,j)*tp(i,j);
Vf(i,j)=Vf(i,j)*tp(i,j);
Us(i,j)=Us(i,j)*tp(i,j);
Vs(i,j)=Vs(i,j)*tp(i,j);
H(i,j)=H(i,j)*tp(i,j);
T(i,j)=T(i,j)*tp(i,j);
R(i,j)=R(i,j)*tp(i,j);
S(i,j)=S(i,j)*tp(i,j);
end
end
for i=1:m+mleft+mright
uft45(i,k)=Uf(i,mup+45);
vft45(i,k)=Vf(i,mup+45);
rt45(i,k)=R(i,mup+45);
ht45(i,k)=H(i,mup+45);
tt45(i,k)=T(i,mup+45);
ust45(i,k)=Us(i,mup+45);
vst45(i,k)=Vs(i,mup+45);
st45(i,k)=S(i,mup+45);
end
for i=1:m+mleft+mright
uft43(i,k)=Uf(i,mup+43);
vft43(i,k)=Vf(i,mup+43);
rt43(i,k)=R(i,mup+43);
ht43(i,k)=H(i,mup+43);
tt43(i,k)=T(i,mup+43);
ust43(i,k)=Us(i,mup+43);
vst43(i,k)=Vs(i,mup+43);
st43(i,k)=S(i,mup+43);
end
if k==200
for i=1:m+mleft+mright
for j=1:m+mup+mdown
uf200(i,j)=Uf(i,j);
vf200(i,j)=Vf(i,j);
r200(i,j)=R(i,j);
h200(i,j)=H(i,j);
t200(i,j)=T(i,j);
us200(i,j)=Us(i,j);
vs200(i,j)=Vs(i,j);
s200(i,j)=S(i,j);
end
end
end
if k==400
for i=1:m+mleft+mright
for j=1:m+mup+mdown
uf400(i,j)=Uf(i,j);
vf400(i,j)=Vf(i,j);
r400(i,j)=R(i,j);
h400(i,j)=H(i,j);
t400(i,j)=T(i,j);
us400(i,j)=Us(i,j);
vs400(i,j)=Vs(i,j);
s400(i,j)=S(i,j);
end
end
end
if k==600
for i=1:m+mleft+mright
for j=1:m+mup+mdown
uf600(i,j)=Uf(i,j);
vf600(i,j)=Vf(i,j);
r600(i,j)=R(i,j);
h600(i,j)=H(i,j);
t600(i,j)=T(i,j);
us600(i,j)=Us(i,j);
vs600(i,j)=Vs(i,j);
s600(i,j)=S(i,j);
end
end
end
if k==700
for i=1:m+mleft+mright
for j=1:m+mup+mdown
uf700(i,j)=Uf(i,j);
vf700(i,j)=Vf(i,j);
r700(i,j)=R(i,j);
h700(i,j)=H(i,j);
t700(i,j)=T(i,j);
us700(i,j)=Us(i,j);
vs700(i,j)=Vs(i,j);
s700(i,j)=S(i,j);
end
end
end
wt(k)=Vs(mleft+83,mup+83);
ut(k)=Us(mleft+83,mup+83);
wtf(k)=Vf(mleft+83,mup+83);
utf(k)=Uf(mleft+83,mup+83);
rt(k)=R(mleft+83,mup+83);
ht(k)=H(mleft+83,mup+83);
tt(k)=T(mleft+83,mup+83);
st(k)=S(mleft+83,mup+83);
end
zuida=max(abs(wt));
wt=wt/zuida*0.9;
a=1:t;
plot(a*10^(-4),wt)
% figure
% plot(a*10^(-4),ut)
% figure
% plot(a*10^(-4),wtf)
% figure
% plot(a*10^(-4),utf)
% [fp1,msg] = fopen('D:\MATLAB6p5\work2\liangceng.txt','w');
% fprintf(fp1,'%f ',dx);
% fprintf(fp1,'%f ',m);
% fprintf(fp1,'%f ',dt*10^3);
% fprintf(fp1,'%f ',t);
% fprintf(fp1,'%f ',g);
% for i=1:m
% for j=1:t
% fprintf(fp1,'%E ',ut(mleft+i,j));
% end
% end
% fclose(fp1);
save poromodeldx25changesourcequan Uf Vf Us Vs H R T S wt ut wtf utf st rt ht tt m dx mleft mdown mup mright ...
uft43 vft43 rt43 ht43 tt43 ust43 vst43 st43 uft45 vft45 rt45 ht45 tt45 ust45 vst45 st45...
uf200 vf200 r200 h200 t200 us200 vs200 s200 uf400 vf400 r400 h400 t400 us400 vs400 s400...
uf600 vf600 r600 h600 t600 us600 vs600 s600 uf700 vf700 r700 h700 t700 us700 vs700 s700
beep
warndlg('work is complete!')⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -