cesp.rar_Verruijt_圆孔扩张_孔扩张

% INTRODUCTION % This program works for the problem of Cavity Expansion in Semi-infinite Plate. % It was given by Verruijt (1998). % INPUT r = 2.08; h = 6.68;% m t = 0.2e6; E = 45e6;% Pa u = 0.3; xmax = 10; ymin = -20;% ranger % ANALYTICAL SOLUTION d = 0.5;% degree m = 2*ceil(xmax/d); n = ceil((0-ymin)/d); b = r/h; a = (1-sqrt(1-b^2))/b; c = -1i*h*(1-a^2)/(1+a^2); p = a^2*t*h/(1-a^2)/(1-a^4); rf = 0; for x=0:m for y=0:n z = (x-m/2)*d-1i*y*d; rz = abs(z+1i*h); if rz>=r s = (z-c)/(z+c); F1 = 4i*p*c/(z+c)^2*(1-(a/s)^2); G1 = 4i*p*c/(z+c)^2*(a^2+s-1/s^2-a^2/s^3); F2 = 4i*p*((2*c*((a^2*(c+z)^2)/(c-z)^2-1))/(c+z)^3-(c*((a^2*(2*c+2*z))/(c-z)^2+(2*a^2*(c+z)^2)/(c-z)^3))/(c+z)^2); TC = conj(z)*F2+G1; tx1 = real(2*F1-TC); txx(y+1,x+1) = tx1;% sigmax,horizontal normal stress ty1 = real(2*F1+TC); tyy(y+1,x+1) = ty1;% sigmay,vertical normal stress txy(y+1,x+1) = imag(TC);% tauxy,shear stress t1(y+1,x+1) = txx(y+1,x+1)/2+tyy(y+1,x+1)/2+sqrt((txx(y+1,x+1)/2-tyy(y+1,x+1)/2)^2+txy(y+1,x+1)^2); t2(y+1,x+1) = txx(y+1,x+1)/2+tyy(y+1,x+1)/2-sqrt((txx(y+1,x+1)/2-tyy(y+1,x+1)/2)^2+txy(y+1,x+1)^2); if t1(y+1,x+1)>0 if rf<rz rf = rz; end plot(real(z),imag(z),'g.'); hold on end end end end