SUBROUTINE DISP(U,KSTEP,KINC,TIME,NODE,NOEL,JDOF,COORDS) C INCLUDE 'ABA_PARAM.INC' C DIMENSION U(3),TIME(2),COORDS(3) C IF(KSTEP.EQ.1) THEN U(1)=45.*TIME(1)*1.E5 U(2)=0.0 U(3)=0.0 ELSE TMP1=1. TMP4=4. PI=ATAN(TMP1)*TMP4 TMP22=.22222222 TMP45=45. U(1)=2.*SIN(TMP22*PI*TIME(1))+TMP45 U(2)=(2.*TMP22*PI)*COS(TMP22*PI*TIME(1)) U(3)=-(2.*(TMP22*PI)**2)*SIN(TMP22*PI*TIME(1)) END IF RETURN END C User subroutines UWAVE used in C exa_riserdynamics_stokes_disp_uwave.inp. LAMBDA and MU (please C refer to Equation 6.2.3-7 and Equation 6.2.3-8 in "Stokes wave C theory," Section 6.2.3 of the Abaqus Theory Manual) used in the C subroutine can be obtained in the data file by running datacheck C on the input file with *PREPRINT, MODEL=YES subroutine uwave(v,a,pdyn,dpdyndz,surf,lpdyn, 1 lrecompute,luplocal,lupglobal, 1 lsurf,ndim,xcur,xintmed, 2 grav,density,elevb,elevs, 3 seed,nspectrum,wamp, 2 time,dtime,noel,npt,kstep,kinc) C include 'aba_param.inc' C dimension v(ndim),a(ndim),xcur(ndim),xintmed(ndim) dimension time(2),wamp(2,nspectrum) dimension e(6), rD(5), rE(5) C parameter(pi=3.14159265358979d0,two=2.0d0,abig=1.d36) parameter (term=50.d0,const2=2.d10,twopi=2.d0*pi ) parameter (half=0.5d0 ) C C tolexp added to check for overflow of exponential in Stokes waves tolexp = 0.9d0*log(abig) luplocal=0 lupglobal=0 if(lrecompute.ne.0) then C Only stochastic analysis with UWAVE can have lrecompute=1 C The user must set other flags accordingly; see User's manual. else C For regular Aqua analysis with UWAVE, lrecompute=0 always C C C Wave definition for stokes wave component: C Phase angle of waves: in radians phase=35.974d0*pi/180.d0 C C Wave travel direction: xdir=1.0d0 ydir=0.0d0 C C Period, wavelength, wave number, wave height, frequency: C period=9.d0 c the wavelength lambda and the parameter mu in the equations c (6.2.3-7) and (6.2.3-8) of the theory manual for Stokes wave. waveln=424.294026946382D0 rmu=0.146795734989733D0 wavenum=twopi/waveln c waveht=20.0d0 freq=1.0D0/period angvel=twopi*freq C depth=elevs-elevb tpdowln=twopi*depth/waveln if(tpdowln.lt.term)then s=sinh(tpdowln) c=cosh(tpdowln) else rLnTwo = log(two) argument = tpdowln - rLnTwo s=exp(argument) c=s end if rmu2=rmu*rmu rmu3=rmu2*rmu rmu4=rmu3*rmu rmu5=rmu4*rmu os=1.D0/s os2=os*os os3=os2*os os4=os3*os os5=os4*os os6=os5*os os7=os6*os oc=1.0D0/c oc2=oc*oc oc4=oc2*oc2 oc5=oc4*oc oc6=oc5*oc oc8=oc6*oc2 oc10=oc8*oc2 oc12=oc10*oc2 oc14=oc12*oc2 oc16=oc14*oc2 b55fac=8.0D0-11.d0*oc2+3.*oc4 c2fac=6.d0-oc2 bsub33=.046875d0*(8.0D0+oc6) b35top=88128.d0-208224.d0*oc2+70848.d0*oc4+54000.d0*oc6- 1 21816.d0*oc8+6264.d0*oc10-54.*oc12-81.d0*oc14 b35bot=12288.d0*c2fac bsub35=b35top/b35bot b55top=192000.d0-262720.d0*oc2+83680.d0*oc4+20160.d0*oc6- 1 7280.d0*oc8+7160.d0*oc10-1800.d0*oc12-1050.d0*oc14+ 2 225.d0*oc16 b55fac=8.0D0-11.d0*oc2+3.*oc4 b55bot=12288.*c2fac*b55fac bsub55=b55top/b55bot t=s/c t2=t*t t3=t2*t t4=t2*t2 t6=t4*t2 t8=t6*t2 t9=t8*t t10=t6*t4 t12=t10*t2 a11=os a13=-.125d0*(5.d0+oc2)*os/t4 a15=-(1184.d0-1440.d0*oc2-1992.d0*oc4+2641.d0*oc6- 1 249.d0*oc8+18.d0*oc10)*os/(1536.d0*t10) a22=.375d0*os4 a24=(192.d0-424.d0*oc2-312.d0*oc4+480.d0*oc6-17.d0*oc8)* 1 os2/(768.d0*t8) a33=.015625d0*(13.d0*oc2-4.0d0)*os5/t2 a35=(512.d0+4224.d0*oc2-6800.d0*oc4-12808.d0*oc6+16704.d0 1 *oc8-3154.*oc10+107.d0*oc12)*os3/(4096.d0*t10*c2fac) a44=(80.d0-816.d0*oc2+1338.d0*oc4-197.d0*oc6)*os6/ 1 (1536.d0*t4*c2fac) a55=-(2880.d0-72480.d0*oc2+324000.d0*oc4-432000.d0*oc6+ 1 163470.d0*oc8-16245.d0*oc10)*os7/(61440.d0*t4*c2fac*b55fac) b22=0.25d0*(two+oc2)/t3 b24=(272.d0-504.*oc2-192.d0*oc4+322.d0*oc6+ 1 21.d0*oc8)/(384.d0*t9) b33=bsub33/t6 b35=bsub35/t12 b44=(768.d0-448.d0*oc2-48.d0*oc4+48.d0*oc6+106.d0*oc8- 1 21.d0*oc10)/(384.d0*t9*c2fac) b55=bsub55/t10 rD(1)=rmu*a11+rmu3*a13+rmu5*a15 rD(2)=-rmu2*a22-rmu4*a24 rD(3)=rmu3*a33+rmu5*a35 rD(4)=-rmu4*a44 rD(5)=rmu5*a55 rE(1)=-rmu rE(2)=rmu2*b22+rmu4*b24 rE(3)=-rmu3*b33-rmu5*b35 rE(4)=rmu4*b44 rE(5)=-rmu5*b55 if (lsurf.eq.1) then C Calculate the instantaneous water surface only, no C wave kinematics are required: sn=xdir*xcur(1) if (ndim.eq.3) sn=sn+ydir*xcur(2) termt=wavenum*sn-angvel*time(2)+phase eta=0.0D0 do k1=1,5 eta=eta+rE(k1)*cos(dble(k1)*termt) end do eta=eta/wavenum surf=elevs+eta else C Calculate wave kinematics: C C lpdyn = 0 - velocity and acceleration for C drag or inertia loads C lpdyn = 1 - dynamic pressure and gradient of C dynamic pressure for buoyancy loads C Calculate some terms needed in the computations: wavespeed=waveln/period w11=twopi*wavespeed/period C STOKES' 5th order wave theory: C C Calculate velocity and acceleration due to C Stokes' waves using current stuctural point. C Also, calculate the dynamic pressure and derivatives. sn=xdir*xcur(1) if(ndim.eq.3)sn=sn+ydir*xcur(2) termz=wavenum*(xcur(ndim)-elevb) termz=max(termz,0.0D0) termt=wavenum*sn-angvel*time(2)+phase do k1=1,6 e(k1)=0.0D0 end do do k1=1,5 fk1=dble(k1) tz=fk1*termz tt=fk1*termt cost=cos(tt) sint=sin(tt) t1=fk1*rD(k1) if(tz.lt.50.0d0) then et1=exp(tz) etm1=exp(-tz) coshz=half*(et1+etm1) sinhz=half*(et1-etm1) e(1)=e(1)+t1*coshz*cost e(2)=e(2)+t1*sinhz*sint t1=t1*fk1 e(3)=e(3)+t1*coshz*sint e(4)=e(4)+t1*sinhz*cost if (lpdyn.eq.1) then e(5)=e(5)+fk1*t1*sinhz*sint e(6)=e(6)+fk1*t1*coshz*cost end if else if(t1.ne.0.0D0) then t1sign=1.0D0 if(t1.lt.0.0D0) t1sign=-1.0D0 C check overflow argument = log(abs(t1))+tz-log(two) if (argument.gt.tolexp) argument=tolexp t1cosh=t1sign*exp(argument) t1sinh=t1cosh else t1cosh=0.0D0 t1sinh=0.0D0 end if e(1)=e(1)+t1cosh*cost e(2)=e(2)+t1sinh*sint e(3)=e(3)+fk1*t1cosh*sint e(4)=e(4)+fk1*t1sinh*cost if (lpdyn.eq.1) then e(5)=e(5)+fk1*fk1*t1sinh*sint e(6)=e(6)+fk1*fk1*t1cosh*cost end if end if end do term11=-wavespeed*e(1) v(1)=v(1)+term11*xdir if(ndim.ne.2) v(2)=v(2)+term11*ydir v(ndim)=v(ndim)-wavespeed*e(2) term11=w11*(e(2)*e(4)-(1.0D0+e(1))*e(3)) a(1)=a(1)+term11*xdir if(ndim.ne.2) a(2)=a(2)+term11*ydir a(ndim)=a(ndim)+w11*((1.0D0+e(1))*e(4)+e(2)* 1 e(3)) if (lpdyn.eq.1) then C Dynamic pressure = rho * ( d phi / d t - v.v/2 ) C - lambda * cbar / tau termdyn = - waveln*wavespeed/period termdyn = termdyn * e(1) C Norm of the velocity squared: vnorm2 = v(1)*v(1) + v(2)*v(2) if(ndim.eq.3) vnorm2 = vnorm2 + v(3)*v(3) pDyn = density * ( termdyn - half*vnorm2 ) C Dynamic pressure gradient -- C d pDyn/d z = - rho*2*pi*cbar/tau*( e4 + e1*e4 + e2*e3 ) termgrad = - density * w11 DpdynDz = termgrad*( e(4) + e(1)*e(4) + e(2)*e(3) ) C Second derivative of the dynamic pressure, C d2 pDyn/d z2 = - (rho*2*pi*cbar/tau) * (2pi/lambda) * C ( e6 + e1*e6 + e4*e4 + e2*e5 + e3*e3 ) D2pDz2 = termgrad*twopi/waveln 1 * ( e(6) + e(1)*e(6) + e(4)*e(4) + e(2)*e(5) + e(3)*e(3) ) end if end if end if C return end