CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COLUMN ANALYSIS PROGRAM C - Combining Flexure and Shear Responses C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C input units are N and mm (exception Po in kN). C output units are kN and m. C C Variables for Input Data for Concrete real fpco, ept, rcc C Variables for Input Data for Steel real fy_L, esh_L, rsh_L, fsu_L, esu_L, fy_T, esh_T, rsh_T, + fsu_T, esu_T, db_L, db_T, db_tie C Variables for Input Data for Section real bw, h, cover, s integer section, nbar, nlayer, nx_tie, ny_tie C Variables for Input Data for Geometry real Po, asp, pdelta, L integer end_fixity, NGauss, exec C Variables for Parameter Calculation for Concrete real Ec, epco, fpt, Esec, Esec_ten, fplx, fply, fpl, rK, nc, + rco, fc2epco, bracket, ecspall, fpcc, epcc, ncc, beta C Variables for Parameter Calculation for Steel real Es, ns, ey_L, ey_T, AsL1, Ab, Ab_tie, Ast, AshR, rhot, + rhov, AsL(25), rhocc, rhovx, rhovy, rhos, pL, pT, eystar C Variables for Parameter Calculation for Section real dc_T, dc_L, jd, d, Ag, Av, y(25), bc, dc, sprime, wix, + wiy, Ai, Ae, ke, bw_eff, Acc, eo, phi C Variables for Parameter Calculation for Geometry real xi(10,10), w(10,10), P, M, Lf, Ls C Variables for Flexural Analysis real Pcom, Mcom, Pten, Mten, PP(51), tolP, delP, Mn, My, + phipy, phiy, eoold, eonew, phinew, Pold, Mold, Pnew, Mnew, + eobal, phibal, Pbal, Mbal, phiold, deo, dphi, Meff, Peff, + esy, deltapy, the, Lpy, Mcr, EIun, EIcr, phicr integer ctype C Variables for Shear Mechanism real crack, alpha, gamma, jeta, rf, x1, AshT, r, T_th, T_th_s C Variables for Combination of Flexure and Shear Responses real Kfun, tolVf, tolrf, Vuold, Vfold, Vunew, Vfnew, Vtotal, + tDISP, sDISP, fDISP, Mmax, dgamma, dVudgamma, delV, rf1, + rf2, tolVu, Kphi, phiplastic, drift, fdrift, sdrift, Kfcr, + Vs, Vc, Vp, Vu, Vumax, fDISP_peak, D_Mmax integer flag C Variables for Functions real cot, fco, fcc, f1avg, fsL, fsT C Character Variables character*80 dummy, pjtitle C Statement for Common Arguments common /group01/ epco, nc, rco, ecspall, fc2epco, bracket common /group02/ epcc, ncc, fpcc, rcc common /group03/ fy_L, esh_L, fsu_L, esu_L, pL common /group04/ bw, Ast, AsL, dc_T, Ag, y, xi, w, beta common /group05/ Lf, phipy, My, Lpy, Mcr, EIun, EIcr common /group06/ AshT, jd, s, Av, x1 common /group07/ eystar, fsu_T, esu_T, pT common /group08/ rhov, ns, T_th_s common /group09/ T_th, Esec_ten common /group10/ fpt common /group11/ Es common /group12/ h common /group13/ alpha common /group14/ crack common /group15/ Ec common /group16/ ept, r common /group17/ tolesT, tole1, tolbw common /group18/ fpco common /groupi1/ section, nlayer, NGauss common /groupi2/ end_fixity C C Open Files C open (unit=1, file='input.dat', status='old') open (unit=2, file='para.out', status='new') C C READ DATA FILE C read (1,901) dummy read (1,901) pjtitle read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy C Concrete Property read (1,*) fpco, ept, rcc, rK read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy C Longitudinal Steel Property read (1,*) fy_L, esh_L, rsh_L, fsu_L, esu_L read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy C Transverse Steel Property read (1,*) fy_T, esh_T, rsh_T, fsu_T, esu_T read (1,901) dummy read (1,901) dummy C Section type read (1,*) section read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy C Section Geometry for both Rectangular and Circular Columns read (1,*) h, nbar, cover, db_L, db_T, s if (section.eq.1) then read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy C Section Geometry for Rectangular Column only read (1,*) bw, nlayer, db_tie, nx_tie, ny_tie else if (section.eq.2) then read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy endif read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy read (1,901) dummy C General Geometry read (1,*) end_fixity, L, Po, asp, pdelta, NGauss read (1,901) dummy read (1,901) dummy C Execution Code read (1,*) exec read (1,901) dummy C C PARAMETER CALCULATION FOR FLEXURE C C Constants pi=3.141592 C Zeros and Weights for Gauss-Legendre Quadrature xi(2,1)=-0.5773502692 xi(2,2)=0.5773502692 xi(3,1)=-0.7745966692 xi(3,2)=0.0 xi(3,3)=0.7745966692 xi(4,1)=-0.8611363116 xi(4,2)=-0.3399810436 xi(4,3)=0.3399810436 xi(4,4)=0.8611363116 xi(5,1)=-0.9061798459 xi(5,2)=-0.5384693101 xi(5,3)=0.0 xi(5,4)=0.5384693101 xi(5,5)=0.9061798459 xi(6,1)=-0.9324695142 xi(6,2)=-0.6612093865 xi(6,3)=-0.2386191861 xi(6,4)=0.2386191861 xi(6,5)=0.6612093865 xi(6,6)=0.9324695142 w(2,1)=1.0 w(2,2)=1.0 w(3,1)=5./9. w(3,2)=8./9. w(3,3)=5./9. w(4,1)=0.3478548451 w(4,2)=0.6521451549 w(4,3)=0.6521451549 w(4,4)=0.3478548451 w(5,1)=0.2369268850 w(5,2)=0.4786286705 w(5,3)=0.5688888889 w(5,4)=0.4786286705 w(5,5)=0.2369268850 w(6,1)=0.1713244924 w(6,2)=0.3607615730 w(6,3)=0.4679139346 w(6,4)=0.4679139346 w(6,5)=0.3607615730 w(6,6)=0.1713244924 C beta parameter for ACI stress block approach if (fpco.le.27.6) then beta=0.85 else if (fpco.gt.55.2) then beta=0.65 else beta=0.85-(fpco-27.6)/6.9*0.05 endif endif C Unconfined Concrete in Compression for Tsai's Equation C Ec=8200.*(fpco)**(3./8.) Ec=4700.*sqrt(fpco) epco=(fpco)**(1./4.)/1153. Esec=fpco/epco nc=Ec/Esec rco=fpco/5.2-1.9 fc2epco=nc*2.*fpco/(1.+(nc-rco/(rco-1.))*2.+2.**rco/(rco-1.)) bracket=(1.+2.*(nc-rco/(rco-1.))+2.**rco/(rco-1.))**2 ecspall=2.*epco+epco*fc2epco/fpco/nc/abs(1.-2.**rco)*bracket C Concrete in Tension for Popovics Equation fpt=sqrt(fpco)/3. Esec_ten=fpt/ept C Steel Es=200000.0 ns=Es/Ec C Longitudinal Steel ey_L=fy_L/Es pL=Es*rsh_L*(esu_L-esh_L)/(fsu_L-fy_L) C Transverse Steel ey_T=fy_T/Es pT=Es*rsh_T*(esu_T-esh_T)/(fsu_T-fy_T) eystar=(fy_T-Es*rsh_T*esh_T)/(Es-Es*rsh_T) C General Section dc_T=cover+0.5*db_T dc_L=cover+db_T+0.5*db_L jd=h-2.0*dc_L d=h-dc_L AsL1=pi/4.0*db_L**2.0 Ab=pi/4.0*db_T**2.0 Ast=AsL1*nbar IF (section.eq.1) THEN C Parameters for Rectangular Section Ab_tie=pi/4.0*db_tie**2.0 AshR=2.0*Ab+Ab_tie*ny_tie Ag=bw*h rhot=Ast/Ag Av=bw*jd rhov=AshR/bw/s do 101 k=1,nlayer y(k)=jd*(-0.5+(k-1.)/(nlayer-1.)) if (k.eq.1.or.k.eq.nlayer) then AsL(k)=AsL1*0.5*(nbar-2.0*(nlayer-2)) else AsL(k)=2.0*AsL1 endif 101 continue C Confinement Effectiveness for Rectangular Column bc=bw-2.0*dc_T dc=h-2.0*dc_T sprime=s-db_T rhocc=Ast/bc/dc wix=(bw-2.0*dc_L)/(AsL(1)/AsL1-1)-db_L wiy=jd/(nlayer-1)-db_L Ai=2.0*(AsL(1)/AsL1-1)*wix**2.0/6.0+2.0*(nlayer-1)*wiy**2.0/6.0 Ae=(bc*dc-Ai)*(1.0-sprime/2.0/bc)*(1.0-sprime/2.0/dc) ke=Ae/bc/dc/(1.0-rhocc) rhovx=(2.0*Ab+Ab_tie*nx_tie)/s/dc rhovy=(2.0*Ab+Ab_tie*ny_tie)/s/bc fplx=ke*rhovx*fy_T/fpco fply=ke*rhovy*fy_T/fpco ELSE IF (section.eq.2) THEN C Parameters for Circular Section Ag=pi/4.0*h**2 rhot=Ast/Ag Av=pi/4.0*(h-2.0*dc_T)**2.0 bw_eff=Av/jd rhov=2.0*Ab/s/(h-2.0*dc_T) the=2.0*pi/nbar if (int(nbar/2.0).ne.(nbar/2.0)) then nlayer=1+(nbar-1)/2 else nlayer=1+nbar/2 endif do 102 k=1,nlayer y(k)=-jd/2.0*cos((k-1)*the) if (int(nbar/2.0).ne.(nbar/2.0)) then if (k.eq.1) then AsL(k)=AsL1 else AsL(k)=2.0*AsL1 endif else if ((k.eq.1).or.(k.eq.nlayer)) then AsL(k)=AsL1 else AsL(k)=2.0*AsL1 endif endif 102 continue C Confinement Effectiveness for Circular Column sprime=s-db_T rhocc=nbar*AsL1/Av Ae=Av*(1.0-sprime/2.0/(h-2.0*dc_T))**2.0 Acc=Av*(1.0-rhocc) ke=Ae/Acc rhos=2.0*rhov fpl=0.5*ke*rhos*fy_T rK=-1.254+2.254*sqrt(1+7.94*fpl/fpco)-2.0*fpl/fpco END IF C Confined Concrete in Compression for Tsai's Equation fpcc=rK*fpco epcc=epco*(1.+5.*(rK-1.)) ncc=Ec*epcc/fpcc C C PARAMETER CALCULATION FOR SHEAR C C Parameters Pertaining to Two-Point Gauss-Legendre Quatrature x1=0.5+xi(2,1)/2. if (end_fixity.eq.1) then jeta=x1+(1.-x1)**2.*(1.-2.*x1) else jeta=2.-3.*x1+5.*x1**2.-2.*x1**3. endif C Corner-to-Corner Diagonal Angle alpha=atan(jd/L) C Crack Angle crack1=atan(((rhov*ns+jeta*(rhov/rhot)*(Av/Ag)) + /(1.+rhov*ns))**0.25) if (end_fixity.eq.1) then crack=crack1 else crack=(crack1+alpha)/2. endif C Parameter r in Popovics Equation for Conc Avg Principal Stress r=Ec/(Ec-Esec_ten) C Parameter T_th for Principal Strain T_th=((1.-(tan(crack))**2.)**2.*x1**2.+(tan(crack))**2.)**2.+ + ((1.-(tan(crack))**2.)**2.*(1.-x1)**2.+(tan(crack))**2.)**2. C Parameter T_th_s for transverse tie strain T_th_s=(1.+x1**2.*(cot(crack))**2.)**2.+ + (1.+(1.-x1)**2.*(cot(crack))**2.)**2. C Effective Section Area of Shear Steel in Columns if (section.eq.1) then AshT=AshR else AshT=2.*s/jd*tan(crack)*sin(pi/2./(1.+s/jd*tan(crack)))/ + sin(pi/2./(1.+jd/s*cot(crack)))*Ab endif C C COLUMN UNIT LENGTH FOR FLEXURE AND SHEAR C if (end_fixity.eq.1) then Lf=L/2. else Lf=L endif Ls=jd*cot(crack) C Length of Yield Penetration Lpy=32.*sqrt(db_L) C C PRINTING INFORMATION C write (2,901) pjtitle write (2,*) write (2,*) write (2,*) 'Parameter Calculation' write (2,*) write (2,*) write (2,*) 'Concrete' write (2,*) write (2,*) ' Modulus of Elasticity Ec', Ec, ' MPa' write (2,*) ' Tensile Strength f''t', fpt, ' MPa' if (section.eq.1) then print *, ' Ratio f''lx/f''co for rectangular column', fplx write (2,*) ' Ratio f''lx/f''co', fplx print *, ' Ratio f''ly/f''co for rectangular column', fply write (2,*) ' Ratio f''ly/f''co', fply elseif (section.eq.2) then write (2,*) ' Confinement Effectiveness K', rK endif write (2,*) write (2,*) write (2,*) 'Steel' write (2,*) write (2,*) ' Yield Strain of Longitudinal Steel ey_L', ey_L write (2,*) ' Yield Strain of Transverse Steel ey_T', ey_T if (section.eq.1) then continue elseif (section.eq.2) then write (2,*) ' Number of longitudinal bar layers', nlayer end if write (2,*) ' Section Area of Longitudinal Steel Ast', + Ast, ' mm^2' write (2,*) ' Volumetric Ratio of Longitudinal Steel rhot', rhot write (2,*) ' Transverse Steel Ratio rhov', rhov write (2,*) write (2,*) write (2,*) 'Column Section Geometry' write (2,*) if (section.eq.1) then write (2,*) ' Section Type: Rectangular' elseif (section.eq.2) then write (2,*) ' Section Type: Circular' end if write (2,*) ' Ratio of Modulus Es/Ec', ns write (2,*) ' Gross Section Area Ag', Ag, ' mm^2' write (2,*) ' Shear Area Av', Av, ' mm^2' write (2,*) ' Lever Arm jd', jd, ' mm' if (section.eq.1) then continue elseif (section.eq.2) then write (2,*) ' Effective Width of Circular Section', + bw_eff, ' mm' endif write (2,*) write (2,*) write (2,*) 'Shear Property' write (2,*) write (2,*) ' Crack Angle', crack1/pi*180.0, ' degree' write (2,*) ' Shear Steel Section Area', AshT, ' mm^2' C Screen for Execution if (exec.eq.1) then go to 999 else open (unit=3, file='pcom.prn', status='new') open (unit=4, file='pten.prn', status='new') open (unit=5, file='plon.prn', status='new') open (unit=6, file='ptra.prn', status='new') open (unit=7, file='pm.prn', status='new') open (unit=8, file='phim.prn', status='new') open (unit=9, file='pflex.prn', status='new') open (unit=10, file='pvs.prn', status='new') open (unit=11, file='pvc.prn', status='new') open (unit=12, file='pvp.prn', status='new') open (unit=13, file='pvu.prn', status='new') open (unit=14, file='pfd.prn', status='new') endif C C Plots for Material Properties C C Concrete in Compression - Unconfined Cover dx=0.1*epco x=0.0 write (3,*) x, x do 81 k=1, 40 x=x+dx z=fco(x) write (3,*) x, z 81 continue C Concrete in Compression - Confined Core dx=0.1*epcc x=0.0 write (3,*) write (3,*) x, x do 82 k=1, 40 x=x+dx z=fcc(x) write (3,*) x, z 82 continue C Concrete in Tension dx=0.1*ept x=0.0 write (4,*) x, x do 83 k=1, 200 x=x+dx z=f1avg(x) write (4,*) x, z 83 continue C Longitudinal Steel dx=0.2*ey_L x=0.0 write (5,*) x, x 84 x=x+dx z=fsL(x) write (5,*) x, z if (x.gt.(esu_L+esh_L)) then continue else go to 84 endif C Transverse Steel dx=0.2*ey_T x=0.0 write (6,*) x, x 85 x=x+dx z=fsT(x) write (6,*) x, z if (x.gt.(esu_T+esh_T)) then continue else go to 85 endif C C P-M Interaction Diagram C C Compression fiber strain is assumed as -0.003. C C Tolerance in Axial Load (kN) tolP=0.1 ctype=1 C Compression Capacity phi=0.0 eo=-0.003 call pm(ctype,eo,phi,Pcom,Mcom) C Tension Capacity eo=ey_L call pm(ctype,eo,phi,Pten,Mten) C Effective Maximum Compression Capacity phieff=0.003/(h/beta) eoeff=-0.003*beta/2. call pm(ctype,eoeff,phieff,Peff,Meff) C Balanced Failure Point phibal=(ey_L+0.003)/d eobal=phibal*h/2.-0.003 call pm(ctype,eobal,phibal,Pbal,Mbal) C Data Points on PM Interaction Diagram Do 94 k=1,50 PP(k)=(1.-k/50.)*Peff 94 continue PP(51)=0.5*Pten C Generation of Data Points by Iteration write (7,*) Mcom, Pcom write (7,*) Meff, Peff eoold=eoeff Pold=Peff do 103 k=1,51 deo=0.000001 95 eonew=eoold+deo phinew=(eonew+0.003)/(h/2.) if (phinew.lt.0.0) then phinew=phieff else continue endif call pm(ctype,eonew,phinew,Pnew,Mnew) delP=PP(k)-Pnew if (abs(delP).le.tolP) then write (7,*) Mnew, Pnew Pold=Pnew eoold=eonew else if (abs(delP*deo).gt.(0.01*abs(Pnew-Pold))) then deo=abs(deo)*sgn(delP)/(sgn(Pnew-Pold)/sgn(deo)) else deo=delP/((Pnew-Pold)/deo) endif Pold=Pnew eoold=eonew go to 95 end if 103 continue write (7,*) Mten, Pten C C Nominal Moment Capacity of Studied Case C Pold=Pbal eoold=eobal deo=0.000001 104 eonew=eoold+deo phinew=(eonew+0.003)/(h/2.) call pm(ctype,eonew,phinew,Pnew,Mnew) delP=-Po-Pnew if (abs(delP).le.tolP) then Mn=Mnew write (2,*) write (2,*) write (2,*) 'Calculated Parameters' write (2,*) write (2,*) ' Column Axial Load in Studied Case', Po, ' kN' write (2,*) ' Nominal Moment Strength', Mn, ' kN-m' else deo=delP/((Pnew-Pold)/deo) Pold=Pnew eoold=eonew go to 104 endif C C First Yield Moment and Curvature C ctype=2 eoold=0.00001 phiold=(ey_L-eoold)/(d-h/2.) call pm(ctype,eoold,phiold,Pold,Mold) deo=0.000001 105 eonew=eoold+deo phinew=(ey_L-eonew)/(d-h/2.) call pm(ctype,eonew,phinew,Pnew,Mnew) delP=-Po-Pnew if (abs(delP).le.tolP) then My=Mnew phipy=phinew write (2,*) ' First Yield Moment', My, ' kN-m' write (2,*) ' Curvature at First Yield Moment', phipy, ' rad' else deo=delP/((Pnew-Pold)/deo) Pold=Pnew eoold=eonew go to 105 endif C C Cracking Moment C eoold=-0.00001 phiold=(ept-eoold)/(h/2.) call pm(ctype,eoold,phiold,Pold,Mold) deo=-0.000001 71 eonew=eoold+deo phinew=(ept-eonew)/(h/2.) call pm(ctype,eonew,phinew,Pnew,Mnew) delP=-Po-Pnew if (abs(delP).le.tolP) then Mcr=Mnew phicr=phinew write (2,*) ' Moment at Cracking', Mcr, ' kN-m' else deo=delP/((Pnew-Pold)/deo) Pold=Pnew eoold=eonew go to 71 endif C C Calculation of Effective Flexural Stiffness C deltapy=0.0 phiold=0.0 Mold=0.0 deo=0.000001 eoold=0.0 Do 106 k=1,5 esy=0.2*k*ey_L 107 eonew=eoold+deo phinew=(esy-eonew)/(d-h/2.) call pm(ctype,eonew,phinew,Pnew,Mnew) delP=-Po-Pnew if (abs(delP).le.tolP) then deltapy=deltapy+Lf**2.*(Mnew-Mold)*(phinew*(2.*Mnew+Mold)+ + phiold*(Mnew+2.*Mold))/6./My**2. phiold=phinew Pold=Pnew Mold=Mnew else deo=delP/((Pnew-Pold)/deo) Pold=Pnew eoold=eonew go to 107 endif 106 continue C Effective Flexural Rigidity and Stiffness EIun=My*1000.*Lf**2./3./deltapy if (end_fixity.eq.1) then Kfun=3.*EIun/2./Lf**3. EIcr=Es*Ast*L**2.*(tan(alpha))**2./12./jeta/1000. else Kfun=3.*EIun/Lf**3. EIcr=Es*Ast*L**2.*(tan(alpha))**2./3./jeta/1000. endif write (2,*) ' Flexural Elastic Stiffness before Cracking', + Kfun, ' kN/mm' Kfcr=Es*Ast/L*(tan(alpha))**2./jeta/1000. write (2,*) ' Flexural Elastic Stiffness after Cracking', + Kfcr, ' kN/mm' D_Mmax=(Mcr/Kfun+(My-Mcr)/Kfcr)*1000./Lf C C Moment-Curvature & Flexural Force-Displacement Analysis C phi=0.0 phip=0.0 fdrift=0.0 M=0.0 Vfold=0.0 Vfnew=0.0 tolVf=0.1 Mmax=My fDISP_peak=D_Mmax write (8,*) phi, M write (9,*) fdrift, Vfnew phiy=phipy*Mn/My 108 if (phi.le.(2.*phiy)) then dphi=0.1*phiy elseif ((phi.gt.(2.*phiy)).and.(phi.le.4.*phiy)) then dphi=0.25*phiy elseif (phi.gt.(4.*phiy)) then dphi=0.5*phiy endif phi=phi+dphi 109 P=-Po+asp*Vfold call fd(phi,phip,P,M,Vfnew,fDISP,fDISP_peak,Mmax,rf) if ((Vfnew-Vfold).le.tolVf) then continue else Vfold=Vfnew go to 109 endif write (8,*) phi, M fdrift=fDISP/L write (9,*) fdrift, Vfnew if (fdrift.lt.0.04) then if (M.gt.Mmax) then Mmax=M fDISP_peak=fDISP else continue endif phiplastic=(phi-phicr)-(phipy-phicr)*(M-Mcr)/(My-Mcr) if (phiplastic.gt.phip) then phip=phiplastic else continue endif Vfold=Vfnew go to 108 else continue endif C C Shear Only Analysis C tolesT=ey_T/100. tole1=ept/100. tolbw=0.001 gamma=0.0 Vu=1.0 Vumax=1.0 tolVu=0.1 sdrift=0.0 Vs=0.0 Vc=0.0 Vp=0.0 Vunew=0.0 P=-Po write (10,*) sdrift, Vs write (11,*) sdrift, Vc write (12,*) sdrift, Vp write (13,*) sdrift, Vunew 110 if (gamma.le.0.005) then dgamma=0.0001 elseif ((gamma.gt.0.005).and.(gamma.le.0.01)) then dgamma=0.0005 elseif (gamma.gt.0.01) then dgamma=0.001 endif gamma=gamma+dgamma 111 rf=My*1000./Vumax/Lf if (rf.gt.1.0) then rf=1.0 else continue endif call shear(gamma,rf,P,Vs,Vc,Vp,Vunew) if (abs(Vu-Vunew).le.tolVu) then sDISP=gamma*L sdrift=sDISP/L write (10,*) sdrift, Vs write (11,*) sdrift, Vc write (12,*) sdrift, Vp write (13,*) sdrift, Vunew Vu=Vunew P=-Po+asp*Vunew if (Vu.gt.Vumax) then Vumax=Vu else continue endif if (gamma.gt.0.03) then go to 116 else go to 110 endif else Vu=Vunew go to 111 endif C C Combination of Flexure and Shear Responses C 116 tolrf=0.001 Mmax=My fDISP_peak=D_Mmax Vuold=0.0 Vfold=0.0 tDISP=0.0 Vtotal=0.0 gamma=0.0 phi=0.0 phip=0.0 Kphi=My*1000./Lf/phiy write (2,*) ' Flexural Elastic Stiffness for Curvature', + Kphi, ' kN-mm/rad' flag=0 write (14,*) tDISP, Vtotal C Increment in Shear Deformation 120 if (gamma.le.0.005) then dgamma=0.0001 elseif ((gamma.gt.0.005).and.(gamma.le.0.01)) then dgamma=0.0005 elseif (gamma.gt.0.01) then dgamma=0.001 endif if (flag.eq.0) then continue else P=-Po+asp*Vtotal go to 119 endif C Increment in Flexural Deformation if (phi.le.(2.*phiy)) then dphi=0.05*phiy elseif ((phi.gt.(2.*phiy)).and.(phi.le.4.*phiy)) then dphi=0.25*phiy elseif (phi.gt.(4.*phiy)) then dphi=1.25*phiy endif C Flexure phi=phi+dphi 121 P=-Po+asp*Vtotal call fd(phi,phip,P,M,Vfnew,fDISP,fDISP_peak,Mmax,rf) C Shear 119 gamma=gamma+dgamma 122 call shear(gamma,rf,P,Vs,Vc,Vp,Vunew) rf1=rf dVudgamma=(Vunew-Vuold)/dgamma C Combination IF ((Vunew.gt.Vuold).or. + ((Vunew.lt.Vuold).and.(dVudgamma.ge.0.0))) THEN if (Vunew.gt.Vfnew) then tolVu=0.0025*Vfnew kk=0 123 delV=Vfnew-Vunew if (flag.eq.0) then continue else Vfnew=Vunew phi=phi+delV/Kphi fDISP=fDISP+delV/Kfcr go to 118 endif if (abs(delV).le.tolVu) then continue else if (kk.gt.20) then go to 120 else dgamma=delV/((Vunew-Vuold)/dgamma) if (-dgamma.gt.gamma) then dgamma=-gamma/2. else continue endif gamma=gamma+dgamma Vuold=Vunew call shear(gamma,rf,P,Vs,Vc,Vp,Vunew) kk=kk+1 go to 123 endif endif 118 Vtotal=Vfnew elseif (Vunew.lt.Vfnew) then tolVf=0.1 124 delV=Vunew-Vfnew if (flag.eq.0) then continue else Vfnew=Vunew phi=phi+delV/Kphi fDISP=fDISP+delV/Kfcr go to 117 endif if (abs(delV).le.tolVf) then rf2=rf else dphi=delV/((Vfnew-Vfold)/dphi) if (-dphi.gt.phi) then dphi=-phi/2. else continue endif phi=phi+dphi Vfold=Vfnew call fd(phi,phip,P,M,Vfnew,fDISP,fDISP_peak, + Mmax,rf) go to 124 endif 117 if (abs(rf1-rf2).le.tolrf) then Vtotal=Vunew else go to 122 endif endif ELSEIF ((Vunew.lt.Vuold).and.(dVudgamma.lt.0.0)) THEN Vfnew=Vunew Vtotal=Vunew phi=phi+(Vfnew-Vfold)/Kphi fDISP=fDISP+(Vfnew-Vfold)/Kfcr flag=1 ENDIF if ((flag.eq.1).and.(Vunew.gt.Vuold)) then flag=0 else continue endif if (M.gt.Mmax) then Mmax=M fDISP_peak=fDISP else continue endif phiplastic=(phi-phicr)-(phipy-phicr)*(M-Mcr)/(My-Mcr) if (phiplastic.gt.phip) then phip=phiplastic else continue endif Vuold=Vtotal Vfold=Vtotal sDISP=gamma*L tDISP=fDISP+sDISP drift=tDISP/L write (14,*) drift, Vtotal if (drift.lt.0.04) then go to 120 else go to 999 endif C C Format Statement C 901 format(a80) 999 stop end C C Flexural Section Analysis C SUBROUTINE pm(ctype,eo,phi,P,M) C Variables for Input Data for Concrete real fpco, rcc C Variables for Input Data for Steel real fy_L, esh_L, fsu_L, esu_L C Variables for Input Data for Section real bw, h integer section, nlayer, NGauss C Variables for Parameter Calculation for Concrete real epco, nc, rco, fc2epco, bracket, ecspall, fpcc, epcc, ncc, + c, ecy1, ecy2, yeff1, yeff2, xcov1, xcov2, xcore, Ccov1, + Ccov2, Ccore, ecov1, ecov2, ecore, Fconc, Mconc, Cc, Mc, + beta C Variables for Parameter Calculation for Steel real Es, Ast, AsL(25), pL, esy, fstress, Fsteel, Msteel, Fs, Ms C Variables for Parameter Calculation for Section real dc_T, Ag, y(25), eo, phi C Variables for Parameter Calculation for Geometry real xi(10,10), w(10,10), P, M C Variables for Flexural Analysis integer ctype C Variables for Functions real sgn, epsilon, fco, fcc, fsL C Statement for Common Arguments common /group01/ epco, nc, rco, ecspall, fc2epco, bracket common /group02/ epcc, ncc, fpcc, rcc common /group03/ fy_L, esh_L, fsu_L, esu_L, pL common /group04/ bw, Ast, AsL, dc_T, Ag, y, xi, w, beta common /group11/ Es common /group12/ h common /group18/ fpco common /groupi1/ section, nlayer, NGauss if (phi.eq.0.0) then if (eo.lt.0.0) then if (ctype.eq.1) then Cc=0.85*fpco*Ag else Cc=fco(abs(eo))*Ag endif else Cc=0.0 endif Mc=0.0 go to 202 else continue endif ecy1=epsilon(eo,phi,-h/2.) ecy2=epsilon(eo,phi,-h/2.+dc_T) if (ecy1.ge.0.0) then Cc=0.0 Mc=0.0 go to 202 else c=abs(ecy1)/phi endif C ACI Code based P-M Interaction Diagram Construction if (ctype.eq.1) then if (beta*c.le.h) then if (section.eq.1) then Cc=0.85*fpco*beta*c*bw Mc=Cc*(h-beta*c)/2. else the=acos(1.-2.*beta*c/h) ybar=0.5*h*(1.-2./3.*(sin(the))**3./ + (the-sin(the)*cos(the))) Cc=0.85*fpco*(h/2.)**2.*(the-sin(the)*cos(the)) Mc=Cc*(h/2.-ybar) endif else Cc=0.85*fpco*Ag Mc=0.0 endif go to 202 else continue endif C Contribution of Concrete if (abs(ecy1).lt.ecspall) then yeff1=h/2. else yeff1=h/2.-c*(1.-ecspall/abs(ecy1)) endif if (abs(ecy2).lt.ecspall) then yeff2=h/2.-dc_T else yeff2=yeff1 endif if (c.gt.h) then c=h else continue endif Cc=0.0 Mc=0.0 do 201 k=1,NGauss xcov1=0.5*(-h/2.+c-yeff1)+0.5*(-h/2.+c+yeff1)*xi(NGauss,k) xcov2=0.5*(-h/2.+c-yeff2)+0.5*(-h/2.+c+yeff2)*xi(NGauss,k) xcore=0.5*(-h+c+dc_T)+0.5*(c-dc_T)*xi(NGauss,k) ecov1=abs(epsilon(eo,phi,xcov1)) ecov2=abs(epsilon(eo,phi,xcov2)) ecore=abs(epsilon(eo,phi,xcore)) if (c.lt.dc_T) then if (section.eq.1) then Ccov1=0.5*(-h/2.+c+yeff1)*w(NGauss,k)*fco(ecov1)*bw Ccov2=0. Ccore=0. else Ccov1=0.5*(-h/2.+c+yeff1)*w(NGauss,k)*fco(ecov1)*h* + sqrt(1.-(2.*xcov1/h)**2.) Ccov2=0. Ccore=0. endif else if (section.eq.1) then Ccov1=0.5*(-h/2.+c+yeff1)*w(NGauss,k)*fco(ecov1)*bw Ccov2=0.5*(-h/2.+c+yeff2)*w(NGauss,k)*fco(ecov2)* + (bw-2.*dc_T) Ccore=0.5*(c-dc_T)*w(NGauss,k)*fcc(ecore)*(bw-2.*dc_T) else Ccov1=0.5*(-h/2.+c+yeff1)*w(NGauss,k)*fco(ecov1)*h* + sqrt(1.-(2.*xcov1/h)**2.) Ccov2=0.5*(-h/2.+c+yeff2)*w(NGauss,k)* + fco(ecov2)*(h-2.*dc_T)* + sqrt(1.-(2.*xcov2/(h-2.*dc_T))**2.) Ccore=0.5*(c-dc_T)*w(NGauss,k)*fcc(ecore)*(h-2*dc_T)* + sqrt(1.-(2.*xcore/(h-2.*dc_T))**2.) endif endif Fconc=Ccov1-Ccov2+Ccore Mconc=Ccov1*abs(xcov1)-Ccov2*abs(xcov2)+Ccore*abs(xcore) Cc=Cc+Fconc Mc=Mc+Mconc 201 continue C Contribution of Steel 202 if (phi.eq.0.0) then if (ctype.eq.1) then Fs=sgn(eo)*fy_L*Ast else Fs=sgn(eo)*fsL(abs(eo))*Ast endif Ms=0.0 else Fs=0.0 Ms=0.0 do 203 k=1,nlayer esy=epsilon(eo,phi,y(k)) if (ctype.eq.1) then fstress=esy*Es if (abs(fstress).gt.fy_L) then fstress=sgn(esy)*fy_L else continue endif else fstress=sgn(esy)*fsL(abs(esy)) endif Fsteel=fstress*AsL(k) Msteel=Fsteel*y(k) Fs=Fs+Fsteel Ms=Ms+Msteel 203 continue endif C Total (units are kN, m) P=(Cc-Fs)/1000. M=(Mc+Ms)/1000000. return end C C Moment-Curvature Relationship C SUBROUTINE mphi(phi,P,M) real phi, P, M, deo, eonew, eoold, Pnew, Pold, delP, tolP, Mnew integer ctype ctype=2 tolP=0.1 eoold=0.0 Pold=0.0 deo=0.000001 204 eonew=eoold+deo call pm(ctype,eonew,phi,Pnew,Mnew) delP=P-Pnew if (abs(delP).le.tolP) then M=Mnew go to 205 else if (abs(Pnew-Pold).lt.0.01) then tolP=tolP+0.01 else continue endif if (abs(delP*deo).gt.(0.01*abs(Pnew-Pold))) then deo=abs(deo)*sgn(delP)/(sgn(Pnew-Pold)/sgn(deo)) else deo=delP/((Pnew-Pold)/deo) endif Pold=Pnew eoold=eonew go to 204 end if 205 return end C C Flexural Force-Displacement Relation C SUBROUTINE fd(phi,phip,P,M,Vf,fDISP,fDISP_peak,Mmax,rf) real phi, P, Vf, fDISP, M, Lf, phip, phipy, My, THETA_p, EIun, + Lpc, Lpy, fDISP_e, fDISP_p, Mmax, rf, alpha, Mcr, fD_e, EIcr, + rf_m, rf_d, fDISP_peak integer end_fixity common /group05/ Lf, phipy, My, Lpy, Mcr, EIun, EIcr common /group13/ alpha common /groupi2/ end_fixity call mphi(phi,P,M) if (M.le.Mcr) then fD_e=M*1000.*Lf**2./3./EIun else fD_e=M*1000.*Lf**2./3./EIun+ + Lf**2./6./M**2.*(M-Mcr)**2.* + (Mcr+2.*M)*(1./EIcr-1./EIun)*1000. endif Lpc=(1.-My/Mmax)*Lf THETA_p=phip*(Lpc/3.+Lpy) if (end_fixity.eq.1) then fDISP_e=2.*fD_e fDISP_p=2.*THETA_p*(Lf-0.25*Lpc) else fDISP_e=fD_e fDISP_p=THETA_p*(Lf-0.25*Lpc) endif Vf=(M*1000.)/Lf fDISP=fDISP_e+fDISP_p rf_m=My/Mmax rf_d=fDISP_peak/fDISP if (rf_m.gt.1.) then rf_m=1. else continue endif if (rf_d.gt.1.) then rf_d=1. else continue endif rf=rf_m*rf_d**0.75 return end C C Shear Analysis C SUBROUTINE shear(gamma,rf,P,Vs,Vc,Vp,Vu) real gamma, rf, Vu, Vs, Vc, Vp, bwsbw, bwcbw, bwpbw, esTold, + esTnew, fsTnew, AshT, jd, s, crack, e1old, e1new, f1new, + Av, Kpth, Ec, h, alpha, THpy, P, Qp, tolesT, tole1, tolbw, + Vs_cd2, Vc_cd2, Vp_cd2, fpco, x1 integer end_fixity, k C Variables for Functions real cot, esT, fsT, e1avg, f1avg C Statement for Common Arguments common /group06/ AshT, jd, s, Av, x1 common /group12/ h common /group13/ alpha common /group14/ crack common /group15/ Ec common /group17/ tolesT, tole1, tolbw common /group18/ fpco common /groupi2/ end_fixity k=1 bwsbw=1./3. bwcbw=1./3. bwpbw=1./3. esTold=0.000001 e1old=0.0000001 C Steel Mechanism 206 esTnew=esT(gamma,bwsbw,esTold) if (abs(esTnew-esTold).le.tolesT) then fsTnew=fsT(esTnew) Vs=AshT*fsTnew*jd/s*cot(crack)/1000. else esTold=esTnew go to 206 endif C Concrete Mechanism 207 e1new=e1avg(gamma,bwcbw,e1old) if (abs(e1new-e1old).le.tole1) then f1new=f1avg(e1new) Vc=f1new*Av*cot(crack)*(1.-2.*(sin(crack))**2.)/1000.*rf else e1old=e1new go to 207 endif C Axial Load Transfer Mechanism Kpth=0.5*Ec*Av*bwpbw*(1.5*h/jd-1.)*(sin(alpha))**2. if (end_fixity.eq.1) then THpy=tan(alpha)/(1.+Kpth/(P*1000.)) else THpy=tan(alpha)/(1.+2.*Kpth/(P*1000.)) endif Qp=-THpy/(tan(alpha)-THpy) Vp=Kpth*gamma*(Qp+(1.-Qp)/(1.+(gamma/THpy)**20.)**0.05)/1000. C Total Shear Strength Vu=Vs+Vc+Vp if (k.gt.100) then go to 208 else continue endif C Effective Width of A Section for Each Mechanism if((abs(bwsbw-Vs/Vu).gt.tolbw).or.(abs(bwcbw-Vc/Vu).gt.tolbw) + .or.(abs(bwpbw-Vp/Vu).gt.tolbw)) then bwsbw=Vs/Vu bwcbw=Vc/Vu bwpbw=Vp/Vu k=k+1 go to 206 else continue endif C Limited by Concrete Strut Strength 208 Vs_cd2=fpco*Av*bwsbw*cot(crack)/2./(0.8+170.*e1new)/ + (1.+(1.-x1)**2.*(cot(crack))**2.)/1000. if (Vs.gt.Vs_cd2) then Vs=Vs_cd2 else continue endif Vc_cd2=fpco*Av*bwcbw*tan(crack)*(1.-(tan(crack))**2.)/ + 2./(0.8+170.*e1new)/ + ((1.-(tan(crack))**2.)**2.*(1.-x1)**2.+(tan(crack))**2.)/ + 1000. if (Vc.gt.Vc_cd2) then Vc=Vc_cd2 else continue endif Vp_cd2=fpco*Av*bwpbw*(1.5*h/jd-1.)/2./(0.8+170.*e1new)/ + cot(alpha)/1000. if (Vp.gt.Vp_cd2) then Vp=Vp_cd2 else continue endif C Total Shear Strength Vu=Vs+Vc+Vp return end CCC CCC Functions CCC C sign function function sgn(x) real sgn, x if (x.lt.0.) then sgn=-1. else sgn=1. endif return end C calculation of plane strain function epsilon(eo,phi,y) real epsilon, eo, phi, y epsilon=eo+phi*y return end C cotangent of trigonometry function cot(x) real cot, x cot=1./tan(x) return end C C Unconfined Cover Concrete C C Tsai's Equation Modified in Descending Branch C function fco(ecy) real fco, ecy, epco, nc, fpco, rco, ecspall, fc2epco, bracket common /group01/ epco, nc, rco, ecspall, fc2epco, bracket common /group18/ fpco x=ecy/epco if (ecy.lt.(2.*epco)) then fco=nc*x*fpco/(1.+(nc-rco/(rco-1.))*x+x**rco/(rco-1.)) elseif ((ecy.ge.(2.*epco)).and.(ecy.lt.ecspall)) then fco=fc2epco-nc*abs(1.-2.**rco)*(x-2.)*fpco/bracket elseif (ecy.ge.ecspall) then fco=0. endif return end C C Confined Core Concrete C C Tsai's Equation Modified in Descending Branch C function fcc(ecy) real fcc, ecy, epcc, ncc, fpcc, rcc common /group02/ epcc, ncc, fpcc, rcc x=ecy/epcc fcc=ncc*x*fpcc/(1.+(ncc-rcc/(rcc-1.))*x+x**rcc/(rcc-1.)) return end C C Longitudinal Steel C C Menegotto-Pinto Equation Modified by Chang and Mander C function fsL(esL) real fsL, esL, Es, fy_L, esh_L, fsu_L, esu_L, pL common /group03/ fy_L, esh_L, fsu_L, esu_L, pL common /group11/ Es fsL=Es*esL/(1.+(Es*esL/fy_L)**20.)**0.05+(1.+sgn(esL-esh_L))/2. + *(fsu_L-fy_L)*(1.-abs((esu_L-esL)/(esu_L-esh_L))**pL) return end C C Transverse Steel C C Transverse Tie Strain C function esT(gamma,bwsbw,esTold) real cot, esT, gamma, crack, rhov, ns, T_th_s, bwsbw, eystar, Es, + fsu_T, esu_T, pT common /group07/ eystar, fsu_T, esu_T, pT common /group08/ rhov, ns, T_th_s common /group11/ Es common /group14/ crack esT=gamma*cot(crack)/(1.+2.*rhov*ns*T_th_s/bwsbw*((1.+(esTold/ + eystar)**20.)**(-0.05)+((1.+sgn(esTold-eystar))/2./Es/esTold)* + (fsu_T-Es*eystar)*(1.-(abs((esu_T-esTold)/ + (esu_T-eystar)))**pT))) return end C C Menegotto-Pinto Equation using Average Stress C function fsT(esT) real fsT, esT, Es, eystar, fsu_T, esu_T, pT common /group07/ eystar, fsu_T, esu_T, pT common /group11/ Es fsT=Es*esT/(1.+(esT/eystar)**20.)**0.05+(1.+sgn(esT-eystar))/2. + *(fsu_T-Es*eystar)*(1.-(abs((esu_T-esT)/(esu_T-eystar)))**pT) return end C C Concrete in Tension C C Concrete Average Principal Strain C function e1avg(gamma,bwcbw,e1old) real cot, e1avg, gamma, crack, T_th, bwcbw, Ec, Esec_ten, + e1old, ept, r common /group09/ T_th, Esec_ten common /group14/ crack common /group15/ Ec common /group16/ ept, r e1avg=gamma*cot(crack)*(cos(crack))**2./ + (1.+2.*T_th*(cos(crack))**4.*(cot(crack))**4./bwcbw/ + (1.+(Ec/Esec_ten-1.)*(e1old/ept)**r)) return end C C Concrete Average Principal Stress Using Popovics Equation C function f1avg(e1) real f1avg, e1, fpt, ept, r common /group10/ fpt common /group16/ ept, r f1avg=fpt*r*(e1/ept)/(r-1.+(e1/ept)**r) return end