program f90_AMS_org !********************************************* ! Buckling pressure of steel pipe under external ground water pressure ! (Amstutz original 1970) !********************************************* implicit none integer::nnn,iii real(8)::Es !Elastic modulus of steel pipe real(8)::pos !Poisson's ratio real(8)::As !Thermal expansion coefficient real(8)::DD !Internal diameter of steel pipe real(8)::t0 !Design thickness of steel pipe real(8)::eps !Margin thickness of steel pipe real(8)::Sigy !Yield stress of steel pipe real(8)::Siga !Allowable stress of steel pipe real(8)::Jc !Joint efficiency real(8)::Temp !Temperature change real(8)::Bg !Plastic deformation modulus of rock real(8)::Ess !Modified elastic modulus real(8)::ny !Coefficient of support effect real(8)::SigF !Modified yield stress real(8)::t !pipe thickness without margin real(8)::rm !rm=(DD+t0)/2 real(8)::r0d !r0d=(DD+2*t0)/2 real(8)::k0 !Gap size between outer surface of steel pipe and inner surface of backfilled concrete real(8)::SigN !Axial stress of pipe at the moment of buckling real(8)::Pk !Buckling pressure of steel pipe real(8)::Phi,Psi,Omega,theta !Coefficients for Amstutz's formula real(8)::offset !Offset at welded joint real(8)::DeltaD !Initial deflection real(8)::m !Parameter for offset real(8)::f1,f2,fi,x1,x2,xi !Variables for bisection method real(8)::pkr(1:5) !Memory of buckling pressure of steel pipe character::linebuf*1000,fmt1*200 fmt1="(f10.3,',',f5.0,',',f10.3,5(',',f10.3))" Es=206000.0D0 pos=0.3D0 As=1.2D-5 t0=30.0D0 eps=1.5D0 Jc=1.00D0 Temp=20.0D0 Bg=1.0D0 write(6,*) '#D0,t0,rm/t0,Pk1,Pk2,Pk3,Pk4,Pk5' do nnn=35,140 DD=2.0D0*dble(nnn)*t0 do iii=1,5 select case(iii) case(1) !HT100 Sigy=885.0D0 Siga=400.0D0 case(2) !SHY685 Sigy=685.0D0 Siga=330.0D0 case(3) !SM570 Sigy=450.0D0 Siga=240.0D0 case(4) !SM490 Sigy=315.0D0 Siga=175.0D0 case(5) !SM400 Sigy=235.0D0 Siga=130.0D0 end select !DeltaD=0.0D0 !offset=0.0D0 DeltaD=DD/300.0D0 offset=0.05D0*t0 Ess=Es/(1.0D0-pos*pos) ny=1.5D0-0.5D0*1.0D0/(1.0D0+0.002D0*Es/Sigy)/(1.0D0+0.002D0*Es/Sigy) SigF=ny*Sigy/sqrt(1.0D0-pos+pos*pos) t=t0-eps rm=(DD+t0)/2.0D0 r0d=(DD+2.0D0*t0)/2.0D0 k0=(As*Temp+Bg*Siga*Jc/Es)*r0d/(1.0D0+Bg) !k0=0.4e-3*rm !Initial deflection rm=rm*(1.0D0+1.52240775D0*DeltaD/DD) r0d=r0d*(1.0D0+1.52240775D0*DeltaD/DD) !Offset m=1.0D0+3.0D0*offset/t !Bisection method x1=Ess*(3.0D0*3.0D0-1.0D0)/12.0D0*t*t/rm/rm x2=SigF do xi=0.5D0*(x1+x2) SigN=x1 call CALC(Ess,rm,t,SigN,Phi,Psi,Omega,theta) f1=FUNC(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi) SigN=x2 call CALC(Ess,rm,t,SigN,Phi,Psi,Omega,theta) f2=FUNC(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi) SigN=xi call CALC(Ess,rm,t,SigN,Phi,Psi,Omega,theta) fi=FUNC(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi) if(f1*fi<0.0D0)x2=xi if(fi*f2<0.0D0)x1=xi if(f1*f2>0.0D0)then write(6,*) 'Imposible',x1,x2,f1,f2 stop end if if(abs(x2-x1)<1.0D-6)exit end do Pk=SigN/(rm/t)/(1.0D0+2.0D0*Omega*rm/t*(SigF-m*SigN)/Ess) pkr(iii)=Pk end do write(linebuf,FMT=fmt1) DD,t0,0.5*DD/t0,(pkr(iii),iii=1,5) call del_spaces(linebuf) write (6,'(a)') trim(linebuf) end do stop contains subroutine CALC(Ess,rm,t,SigN,Phi,Psi,Omega,theta) real(8),parameter::pi=3.14159265358979323846D0 real(8),intent(in)::Ess,rm,t,SigN real(8),intent(out)::Phi,Psi,Omega,theta real(8)::ang,cc,ee,alpha,delta,beta,cota,gamma ! ang=90->0 degree ! ee*ang=270->180 degree ee=sqrt(1.0D0+12.0D0*rm*rm/t/t*SigN/Ess) ang=270.0D0/ee do ang=ang-0.0001D0 alpha=ang*pi/180.0D0 cc=ee*tan(alpha)-tan(ee*alpha) if(0.0D0<=cc)exit end do delta=(ee*ee-1.0D0)*(1.0D0-cos(ee*alpha)) beta=(ee-1.0D0/ee)*(ee*alpha*cos(ee*alpha)-sin(ee*alpha)) if(ang==90.0D0)then cota=0.0D0 else cota=1.0D0/tan(alpha) end if gamma=-ee*(ee*alpha-sin(ee*alpha)*cos(ee*alpha)& +ee*alpha*sin(ee*alpha)*sin(ee*alpha)/sin(alpha)/sin(alpha)& -ee*sin(ee*alpha)*sin(ee*alpha)*cota) Phi=ee*ee*ee*beta/pi/delta Psi=-gamma/4.0D0/beta/delta Omega=-cos(ee*alpha)/(1.0D0-cos(ee*alpha)) theta=alpha/pi*180.0D0 end subroutine CALC real(8) function FUNC(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi) real(8),intent(in)::Ess,rm,t,k0,m,SigF,SigN,Phi,Psi real(8)::a,b a=(k0/rm+SigN/Ess)*(1.0+12.0*rm*rm/t/t*SigN/Ess)**1.5 b=2.0*Phi*rm/t*(SigF-m*SigN)/Ess*(1.0-2.0*Psi*rm/t*(SigF-m*SigN)/Ess) FUNC=a-b end function FUNC subroutine del_spaces(s) character (*), intent (inout) :: s character (len=len(s)) tmp integer i, j j = 1 do i = 1, len(s) if (s(i:i)==' ') cycle tmp(j:j) = s(i:i) j = j + 1 end do s = tmp(1:j-1) end subroutine del_spaces end program f90_AMS_org