Chapmann模型的matlab数值模拟

%w为角频率，lamda，u为弹性参数，fai为孔隙度，e为裂隙密度，tao为时间尺度参数，en为流体粘滞系数，v,miu为泊松比 %r为裂隙纵横比，xiaok为渗透率，cv为裂隙体积,a为裂隙半径，kf为流体体积模量 clear all; vp=2790;vs=1463;rou=2.08;fai=0.3;tao=2*10^-5;e=0.1;miu=0.3104;v=miu;xiaok=1000; lamda=rou*vp*vp-2*rou*vs*vs;r=0.1; u=rou*vs*vs;kf=2E9; k=rou*(vp^2-4*v*vs/3); kc= (pi*u*r*(lamda + u))/(kf*(lamda + 2.0*u)); u=rou*vs*vs; %k=2000; en=0.01;a=0.1;l=2; f=5:1:104; %以下为输入的值； w=2*pi*f; cv=4*pi*a*a*a*r/3; sigma=tao*6*xiaok*l/(en*cv*(1+kc)); kf=(pi*u*u*r*r)/(2*(1-v)*sigma); %tao=4*en*a*a*a*(1-v)/(9*xiaok*l*u);%裂隙纵横比很小的时候 kp=4*u/(3*kf); y=(3*pi*(lamda+u)*(1+kp))/(4*(lamda+2*u)*(1+kc)); yyipie=y*(lamda+2*u)/((3*lamda+2*u)*(1+kp)); awfenzi=((ones(1,100)+i*w*y*tao)*(lamda+2*u)/(lamda+u))*(16*e/(27*fai*(1+kp))+(lamda+u)/(3*lamda+2*miu))+i*w*tao*(1/(3*(1+kc))-yyipie); awfenmu=1+i*w*tao+(1+i*w*y*tao)*(16*e*(1+kc)*(lamda+2*u))/(9*fai*(1+kp)*(lamda+u)); bwfenzi=1/(3*(1+kc))+9*(1+kp)*(lamda+u)/(16*(1+kc)*(3*lamda+2*miu))+i*w*tao*(yyipie-1/(3*(1+kc)))/(1+i*w*pi); bwfenmu=9*(1+kp)*(lamda+u)/(16*(1+kc)*(lamda+2*u))+(1+i*w*y*tao)/(1+i*w*tao); aw=awfenzi/awfenmu; bw=bwfenzi/bwfenmu; keff=k-e*(4*(3*lamda+2*miu)*(lamda+2*miu)*(1-aw)/(miu*(lamda+miu))-4*pi*aw*r)-fai*((3*lamda+2*miu)*((lamda+2*miu)/(3*lamda+2*miu)+bw)/4*miu-3*bw); ueff=u*ones(1,100)-((16/45)*e/(1+kc))*(u*(lamda+2*u)/(3*lamda+4*u))*(kc+1./(1+i*w*tao))-ones(1,100)*(32/45)*e*(u*(lamda+2*u)/(3*lamda+4*u))-ones(1,100)*(fai*15*u*(lamda+2*u)/(9*lamda+14*u)); eeff=keff+4*ueff/3; qp=(-real(eeff))./(imag(eeff)); qs=(-real(ueff))./(imag(ueff)); x=log(abs(f))/log(20); figure(1); qpni=1./qp; plot(x,qpni); qsni=1./qs; figure(2); plot(x,qsni);