diffusion equation in matlab

Question

Mohammadfarid ghasemi 2015-6-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/223253-diffusion-equation-in-matlab

评论： Walter Roberson 2015-6-10

Hi, I have a pressure diffusion equation on a quadratic boundary. I have write the following code to solve it, the pressure should increase with time as we have injection in one side, and constant pressure other side. top and bottom side have isolated. but the code works only when length of medium is so small(<1). can anybody tell me how can I solve it for large length?

%B*DP/Dt+A*(DP/Dx+DP/Dy)+f(t)=0
%f(t)=-2.2/(t-0.1)^0.5
%General input data
clear
t0=0.1;
Sh=20;
G=20000;
nu=0.25;
n=0.716;
E=2*G*(1+nu);
K=1.2;
q_inj=0.067;
h=100;
h_tip=2500;
h_end=2600;
t_max=100;
N_t=60;
N_x=20;
N_y=30;
t = linspace(0.1,t_max,N_t);
y=linspace(-50,50,N_y);
h_fracture=linspace(h_tip,h_end,N_y);
B_y=((1-nu)/G)*(h^2-4*y.^2).^0.5;
B_y_mean=mean(B_y);
%Wellbore pressure calculation
vis=0.85;
den_f=1900;
den_p=2700;
prop_con=0.35;
den_eff=(1-prop_con)*den_f+prop_con*den_p;
D_t=1.28;
Re=4*q_inj*den_f/(pi*vis*D_t);
Re_w=((1+3*n)/4*n)*Re;
X_P1=log10(Re_w)/log10(exp(1));
X_P=log10(X_P1)/log10(exp(1));
f=exp(28.135+(-29.379+(8.2405-0.86227*X_P)*X_P)*X_P);
P_Frac=(-2*f*vis^2/(D_t*den_f)+den_f*9.81*0.001)*h_fracture*0.001;
w0_t=B_y'*(P_Frac-Sh);
w0=mean(mean(w0_t))
P_ave=mean(P_Frac)
B_y_mean=w0/P_ave
E_prin=E/(1-nu^2);
delt_t=(t_max-t0)/N_t;
Lf0=50;
delt_t=(t_max/(N_t));
delt_y=h/N_y;
C_Leak=(9.84*10^-6);
%pde geometry
A0=-w0^3/(12*vis);
rect=[3,4,0,Lf0,Lf0,0,0,0,100,100]';
gd=[rect];
ns=char('rect');
ns=ns';
sf='(rect)';
g=decsg(gd,sf,ns);
model=createpde;
geometryFromEdges(model,g);
generateMesh(model,'jiggle','on')
p = model.Mesh.Nodes;
pdegplot(model,'EdgeLabels','on')
%initial condition
u0 = Sh;
pdegplot(model,'EdgeLabels','on')
tlist = linspace(t0,t_max,N_t);
c=-A0*den_eff;
d=den_eff*B_y_mean;
f='-2.2./((t).^0.5)';
a=0;
applyBoundaryCondition(model,'Edge',4,'g',q_inj*den_eff/1000,'q',0)    
applyBoundaryCondition(model,'Edge',[1,3],'g',0,'q',0)
applyBoundaryCondition(model,'Edge',2,'u',Sh)
u1 = parabolic(u0,tlist,model,c,a,f,d);
sizeu1=size(u1);
nx=sizeu1(1);
x_w=linspace(0,Lf0,nx);
y_w=linspace(-h/2,h/2,nx);
B_y_w=((1-nu)/G)*(h^2-4*(p(2,:)-h/2).^2).^0.5;
w11=(u1(:,1)'-Sh).*B_y_w;
w12=(u1(:,N_t)'-Sh).*B_y_w;
max1=max(max(w12));
min1=min(min(w11));
figure1=figure
for ii=1:N_t
    hold off
   t=tlist(ii);
      figure1=pdeplot(model,'xydata',u1(:,ii));
      colormap default
      minu1=min(min(u1));
      maxu1=max(max(u1));
      caxis([minu1,maxu1])
  title(['t= ',num2str(t),'  sec'])
      pause(0.01)
      if ii<N_t
      delete(figure1)
      end
  end
figure2=figure
for ii=1:N_t
    hold off
    t=tlist(ii);
    w1=B_y_w'.*(u1(:,ii)-Sh);
    figure2=pdeplot(model,'xydata',w1);
    colormap default
    title(['t= ',num2str(t),'  sec'])
    caxis([min1,max1])
    colorbar
    shading interp
    pause(0.01)
    if ii<N_t
    delete(figure2)
    end
end