/* === General comments: Must use consistent units: e.g MN for force, Gg for mass, m for lengnth, s for time Misc Assumptions (that can readily be changed): * The calculated bending stresses sb are at the midsurface of the pipe wall. * The concrete coating stiffening effect (see "csf") is assumed to be as much for the axial rigidity as for the flexural rigidity === Key Nomenclature: _ appended to name for actual as opposed to normalized values E.g. n_ denotes the actual value of the effective axial force for which n denotes the normalized value. i appended to name for expression valid within span mnsh appended to name means maximum value near span shoulder mis appended to name means maximum in span o appended to name for expression valid in soil-supported zone p appended to name for derivative w.r.t. x, e.g. vp=dv/dx sa axial stress sb bending stress, positive for compression on the side of positive transverse displacements sb1 dynamic bending stress at amplitude odwc sc1 sum of absolute values of dynamic bending and axial stresses ua appended to name means value for unit amplitude of vibration v static displacement w modeshape displacement x axial coordinate xx used instead of x_ for actual value of x (x_ seems to have a special meaning in maxima) The inputs and writing results to a file have been deactivated. See "Inputs" and "write_data". DO NOT run this file as is. Intended for with from another file that first specifies the inputs that are deactivated in this file, then uses the command batch("C:/Users/Ralf Peek/Documents/w/Spans/spanL08.mac")$ and finally specifies where to print the results. Notes: {1} Dynamic axial stiffness ka is used only to calcuate the decay in dynamic axial force in the outer zone using characteristic length sqrt(ea/ka). */ fpi:float(%pi)$ /* Inputs Consistent Units used are: m for length gG for mass s for time MN for force */ n_:-0.315$ /* effective axial force */ cl_:38.5$ /* span length from zrb support to TDP */ k:1$ /* dynamic lateral soil spring stiffness */ ka:1$ /* dynamic axial soil stiffness {1} */ cka:10$ /* dynamic stiffness of zrb support */ frquig_:0.5$ /* initial guess for natural frequency when axially unconstrained */ od:20*0.0254$ /* outer diameter of steel pipe */ wt:0.018$ /* wall thickness of steel pipe */ tcc:0.003$ /* thickness of corrosion coating */ twc:0.040$ /* thickness of weight coating */ esteel:207000.$ /* Young's modulus of steel */ cggas:0.1$ /* Specific gravity of contents (e.g. gas or oil) */ cgsteel:7.85$ /* Specific gravity of steel */ cgcc:1.3$ /* Specific gravity of corrosion coating */ cgwc:2.4$ /* Specific gravity of weight coating */ cgsw:1.025$ /* Specific gravity of seawater */ rhofw:0.001$ /* Density of fresh water */ csf:0.$ /* Concrete stiffness factor */ cm:1.$ /* Added mass coefficient */ xiab:0.1$ /* axial seabed damping ratio */ xilb:0.1$ /* lateral seabed damping ratio */ xias:0.2$ /* zrb support axial damping ratio */ xiap:0.01$ /* pipe damping ratio for axial deformations */ xibp:0.01$ /* pipe damping ratio in bending */ /* to be phased out fpull:float(600*9.81e-6)$ /* pull force in pull-and-release test */ agr:9.81$ /* Gravity force per unit mass */ */ /* derived pipe section properties*/ omegaig_:2*fpi*frquig_$ id:od-2*wt$ /*diameters*/ odcc:od+2*tcc$ odwc:odcc+2*twc$ isteel:(od^4-id^4)*fpi/64$ /*2nd moment of area */ iconc:(odwc^4-odcc^4)*fpi/64$ agas:fpi*id^2/4$ /* cross sectional areas */ asteel:fpi*od^2/4$ acc:fpi*odcc^2/4$ awc:fpi*odwc^2/4$ atotal:awc$ awc:awc-acc$ acc:acc-asteel$ asteel:asteel-agas$ cgpipe:(cggas*agas+cgsteel*asteel+cgcc*acc+cgwc*awc)/atotal$ /* specific gravity */ meff: rhofw*atotal*(cgpipe+cgsw*cm)$ /*effective mass per unit length*/ mpipe: rhofw*atotal*cgpipe$ /*pipe & contents mass per unit length*/ qw:agr*rhofw*atotal*(cgpipe-cgsw)$ /*effective weight of pipe*/ ei:(1+csf)*esteel*isteel$ ea:(1+csf)*esteel*asteel$ rb:float(0.5*(od-wt))$ /* radius used to calculate the bending stress; yields stress at midsurface of pipe wall */ /* Other derived parameters */ xsh_:float(cl_)$ lca:sqrt(ea/ka)$ /* char length for dynamic axial force decay */ olea:xsh_+lca$ clea:olea+2*ea/cka$ clea:float(clea)$ /*Unit values for dynamics*/ x_1:(ei/k)^0.25$ n_1:sqrt(ei*k)$ omega_1:sqrt(k/meff)$ /*Calculation of normalized Inputs for dynamics*/ n:n_/n_1$ omegaig:omegaig_/omega_1$ cl:cl_/x_1$ /* Calculation of dimensionless dynamics results */ rl(omega):=block( gdi:n^2+4*omega^2, /*capital delta for inner region*/ a:sqrt((sqrt(gdi)+n)/2.0), /*real root of char eqn is + or - a*/ b:sqrt((sqrt(gdi)-n)/2.0), /*purely imaginary root is + or - ib*/ gdo:gdi-4, /*capital delta for outer region*/ gl2:(n+%i*sqrt(-gdo))/2.0, gl:sqrt(gl2), ga:abs(realpart(gl)), gb:abs(imagpart(gl)), /*General solution functions, f1 & f2 for inner region (symmetric), and f3 and f4 for outer region (decaying) */ xo:x-cl, f1:sinh(a*x), f2:sin(b*x), f3:exp(-ga*xo)*sin(gb*xo), f4:exp(-ga*xo)*cos(gb*xo), /*Boundary conditions (r1 to r4 should vanish at x=cl/2) */ r1:[f1,f2,-f3,-f4], r2:diff(r1,x), r3:diff(r2,x), r4:diff(r3,x), r:matrix(r1,r2,r3,r4), subst([x=cl],r))$ detrl(omega):=determinant(rl(omega))$ load("mnewton")$ omegaval:mnewton(detrl(omega),omega,omegaig)$ omegav:subst(omegaval,omega)$ /* Evaluate further details of solution */ rll:rl(omegav)$ f1:subst(omegaval,f1)$ f2:subst(omegaval,f2)$ f3:subst(omegaval,f3)$ f4:subst(omegaval,f4)$ c1:1.0$ ccval:linsolve( [rll[2,1]*c1+rll[2,2]*c2+rll[2,3]*c3+rll[2,4]*c4, rll[3,1]*c1+rll[3,2]*c2+rll[3,3]*c3+rll[3,4]*c4, rll[4,1]*c1+rll[4,2]*c2+rll[4,3]*c3+rll[4,4]*c4], [c2,c3,c4])$ ccval:float(ccval)$ c2:subst(ccval,c2)$ c3:subst(ccval,c3)$ c4:subst(ccval,c4)$ wi:f1*c1+f2*c2$ wo:f3*c3+f4*c4$ /*Normalize mode for unit maximum displacement */ wpi:diff(wi,x)$ xmval:mnewton(wpi,x,0.4*cl)$ w0:subst(xmval,wi)$ wi:float(wi/w0)$ wo:float(wo/w0)$ c1:float(c1/w0)$ c2:float(c2/w0)$ c3:float(c3/w0)$ c4:float(c4/w0)$ /*Differentiate Modeshape*/ wpi:float(diff(wi,x))$ wppi:float(diff(wpi,x))$ wpppi:float(diff(wppi,x))$ wpo:float(diff(wo,x))$ wppo:float(diff(wpo,x))$ wpppo:float(diff(wppo,x))$ /*Check continuity conditions at shoulder up to 3rd derivative of displ (optional)*/ xsh:cl$ /*x value at span shoulder*/ werr0:wi-wo$ werr1:diff(werr0,x)$ werr2:diff(werr1,x)$ werr3:diff(werr2,x)$ werr:subst([x=xsh],[werr0,werr1,werr2,werr3])$ /*check that this returns 4 zeroes*/ /*Plot normalized displacements*/ w:if x>xsh then wo else wi$ /*expression for modal displacement, not valid for x<-xsh*/ wpp:if x>xsh then wppo else wppi$ /*expression curvature, not valid for x<-xsh*/ /* plot2d([w,wpp],[x,0,2*xsh])$ */ /*Calculate modal strain energy, in bending, axially, and in soil Expected to be symmetric Therefore integrate over the positive half of the span only */ ewi:quad_qag(wi^2,x,0,xsh,3, 'epsrel=1d-9)$ ewo:quad_qagi(wo^2,x,xsh,inf, 'epsrel=1d-9)$ ewpi:quad_qag(wpi^2,x,0,xsh,3, 'epsrel=1d-9)$ ewpo:quad_qagi(wpo^2,x,xsh,inf, 'epsrel=1d-9)$ ewppi:quad_qag(wppi^2,x,0,xsh,3, 'epsrel=1d-9)$ ewppo:quad_qagi(wppo^2,x,xsh,inf, 'epsrel=1d-9)$ ew:ewi[1]+ewo[1]$ ewp:ewpi[1]+ewpo[1]$ ewpp:ewppi[1]+ewppo[1]$ usoil:ewo[1]$ /*modal strain energy in soil*/ ubend:ewpp$ /*modal strain energy in bending */ un:n*ewp$ /*modal potential energy of axial force */ utotal:ubend+un+usoil$ /*total modal strain energy*/ /* maximum moment near shoulder (mnsh) */ x2:xsh$ wppp2:subst(x=x2,wpppo)$ x1:x2$ wppp1:wppp2$ for i:1 while wppp1*wppp2 >= 0 do block( x1:x2,wppp1:wppp2,x2:x2+0.5,wppp2:subst(x=x2,wpppo) )$ xmnsh:find_root(wpppo,x,x1,x2)$ /* x value at max moment just beyond shoulder */ wppmnsh:subst([x=xmnsh],wppo)$ /* maximum moment in span (mis) */ x2:0$ wppp2:subst(x=x2,wpppi)$ x1:x2$ wppp1:wppp2$ for i:1 while wppp1*wppp2 >= 0 do block( x1:x2,wppp1:wppp2,x2:x2+0.5,wppp2:subst(x=x2,wpppi) )$ xmis:find_root(wpppi,x,x1,x2)$ /* x value at max moment in span */ wppmis:subst([x=xmis],wppi)$ /* other normalized mode shape parameters of interest */ wsh:subst([x=xsh],wi)$ /* summarize nondimensional dynamic results of interest */ yynd:[n,omegav,cl,xmnsh-xsh,wppmis,wppmnsh,wsh, ubend, un, usoil, utotal,ewi[1]]$ /* Calculation of dimensional dynamics results Note xx used instead of x_ for the dimensional x, since x_ seems to have a special meaning in Maxima */ omega_:omegav*omega_1$ frqu_:float(omega_/(2*%pi))$ xmnsh_:xmnsh*x_1$ xmis_:xmis*x_1$ sb1mis:(1+csf)*rb*esteel*odwc*wppmis/x_1^2$ sb1mnsh:(1+csf)*rb*esteel*odwc*wppmnsh/x_1^2$ a_:a/x_1$ b_:b/x_1$ ga_:ga/x_1$ gb_:gb/x_1$ dsp:[a_,b_,ga_,gb_,c1,c2,c3,c4]$ /* these are dimensional solution parameters in wi_=c1*sinh(a_*xx)+c2*sin(b_*xx) wo_=(c3*sin(gb_*xo_)+c4*cos(gb_*xo_))*exp(-ga*xo_) where xo_=xx-cl_ */ f1:sinh(a*x), f2:sin(b*x), f3:exp(-ga*xo)*sin(gb*xo), f4:exp(-ga*xo)*cos(gb*xo), wi_:float(subst([x=xx/x_1],wi))$ wo_:float(subst([x=xx/x_1],wo))$ w_:if xx0 side */ ewo_:ewo[1]*x_1$ ewp_:ewp/x_1$ ewpp_:ewpp/x_1^3$ cnf:cgg*ea/olea$ /* unit amplitude value of capital n sub f in notes */ cnd:cgg/(2/cka + olea/ea)$ /* unit amplitude value of capital n sub d in notes */ mmeff:2*ew_*meff$ /* effective modal mass */ ud0:2*cnd/cka$ /* unit-amplitude values for strain energy (both sides) */ cuab:0.5*cnd^2/sqrt(ea*ka)$ /* axial seabed springs */ cuap:cl_*cnd^2/ea+cuab$ /* axial in pipe */ cuas:2*cnd^2/cka$ frqa_:float(sqrt((cuap+cuab+cuas)/(meff*ew_)))/(2*fpi)$ frq_:sqrt(frqu_^2+frqa_^2)$ cubp:ei*ewpp_$ culb:k*ewo_$ cunp:n_*ewp_$ frqu2:float(sqrt((cubp+cunp+culb)/(meff*ew_)))/(2*fpi)$ /* as a check */ zchkf:frqu2-frqu_$ frqs_:[frqu_,frqa_,frq_,zchkf]$ cuk:(2*%pi*frq_)^2*meff*ew_$ /* unit-amplitude kinetic energy */ cuk:float(cuk)$ cutot:cuab+cuap+cuas+cubp+culb+cunp$ zchkcu:cuk-cutot$ zchks:[zchkg,zchkf,zchkcu]; cuall:[cuab,cuas,cuap,cubp,culb,cunp]; curatios:float(cuall/cutot); xiall:[xiab,xias,xiap,xibp,xilb,0]; zall:float(xiall*curatios); float(curatios.[1,1,1,1,1,1]); ztot:zall.[1,1,1,1,1,1]; /* contributions to modal damping */ zeta_ab:xiab*cuab/cutot; zeta_as:xias*cuas/cutot; zeta_ap:xiap*cuap/cutot; zeta_a:zeta_ab+zeta_as+zeta_ap; zeta_bp:xibp*cubp/cutot; zeta_lb:xilb*culb/cutot; zeta_u:zeta_bp+zeta_lb; zeta:zeta_a+zeta_u; sa1i:cnd*odwc/asteel$ /* dynamic axial stress amplitude at vertical vibration amplitude odwc */ sa1o:sa1i*exp(-abs(xx-xsh_)/lca)$ sa1:if xx