Skip to content

Commit ec9992f

Browse files
andrew-plattandrew-platt
authored
Merge pull request #2941 from luwang00/b/HD_Hst
HD: Switched to an equivalent but numerically more reliable expression for the hydrostatic moment on (partially wetted) circular endplates
2 parents 47f153c + 19994dc commit ec9992f

File tree

1 file changed

+81
-76
lines changed

1 file changed

+81
-76
lines changed

‎modules/hydrodyn/src/Morison.f90‎

Lines changed: 81 additions & 76 deletions
Original file line numberDiff line numberDiff line change
@@ -4759,27 +4759,27 @@ SUBROUTINE GetSectionFreeSurfaceIntersects_Cyl( pos0, FSPt, k_hat, y_hat, z_hat,
47594759
! Note: if d2=0, the current section is perfectly parallel to the free surface
47604760
IF ( (d2 >= c*c) .and. (.not. EqualRealNos(d2,0.0_DbKi)) ) THEN ! Has intersection
47614761
d = SQRT(d2)
4762-
IF (b>=0.0) THEN
4762+
IF (b>=0.0_DbKi) THEN
47634763
alpha = ACOS(a/d)
47644764
ELSE
47654765
alpha = -ACOS(a/d)
47664766
END IF
47674767
beta = ACOS(c/d)
47684768
theta1 = alpha - beta
47694769
theta2 = alpha + beta
4770-
IF ( dot_product( (cos(theta2)-cos(theta1))*z_hat-(sin(theta2)-sin(theta1))*y_hat, n_hat) < 0.0 ) THEN
4770+
IF ( dot_product( (cos(theta2)-cos(theta1))*z_hat-(sin(theta2)-sin(theta1))*y_hat, n_hat) < 0.0_DbKi ) THEN
47714771
tmp = theta1
47724772
theta1 = theta2
4773-
theta2 = tmp + 2.0*PI_D
4773+
theta2 = tmp + 2.0_DbKi*PI_D
47744774
END IF
47754775
secStat = 1;
4776-
ELSE IF (c > 0.0) THEN ! Section is fully submerged
4777-
theta1 = -1.5*PI_D
4778-
theta2 = 0.5*PI_D
4776+
ELSE IF (c > 0.0_DbKi) THEN ! Section is fully submerged
4777+
theta1 = -1.5_DbKi*PI_D
4778+
theta2 = 0.5_DbKi*PI_D
47794779
secStat = 2;
47804780
ELSE ! Section is completely dry
4781-
theta1 = -0.5*PI_D
4782-
theta2 = -0.5*PI_D
4781+
theta1 = -0.5_DbKi*PI_D
4782+
theta2 = -0.5_DbKi*PI_D
47834783
secStat = 0;
47844784
END IF
47854785

@@ -4809,8 +4809,8 @@ SUBROUTINE GetSectionHstLds_Cyl( origin, pos0, k_hat, y_hat, z_hat, R, dRdl, the
48094809
cosPhi = SQRT(k_hat(1)**2+k_hat(2)**2)
48104810

48114811
C0 = Z0*dTheta + R*cosPhi*(cosTheta1 -cosTheta2)
4812-
C1 = Z0*(sinTheta2-sinTheta1) + 0.5*R*cosPhi*(cosTheta2**2-cosTheta1**2)
4813-
C2 = Z0*(cosTheta1-cosTheta2) + 0.5*R*cosPhi*(dTheta-sinTheta2*cosTheta2+sinTheta1*cosTheta1)
4812+
C1 = Z0*(sinTheta2-sinTheta1) + 0.5_DbKi*R*cosPhi*(cosTheta2**2-cosTheta1**2)
4813+
C2 = Z0*(cosTheta1-cosTheta2) + 0.5_DbKi*R*cosPhi*(dTheta-sinTheta2*cosTheta2+sinTheta1*cosTheta1)
48144814

48154815
dFdl(1:3) = -R *dRdl*C0*k_hat + R*C1*y_hat + R*C2*z_hat
48164816
dFdl(4:6) = -R**2*dRdl*C2*y_hat + R**2*dRdl*C1*z_hat + CROSS_PRODUCT((pos0-origin),dFdl(1:3))
@@ -4844,21 +4844,21 @@ SUBROUTINE GetSectionHstLds_Rec( origin, pos0, k_hat, x_hat, y_hat, Sa, Sb, dSad
48444844
z0 = origin(3)
48454845

48464846
! Global coordinates of the four vertices of the section
4847-
rv(:,1) = pos0 + x_hat * (-0.5*Sa) + y_hat * (-0.5*Sb)
4848-
rv(:,2) = pos0 + x_hat * ( 0.5*Sa) + y_hat * (-0.5*Sb)
4849-
rv(:,3) = pos0 + x_hat * ( 0.5*Sa) + y_hat * ( 0.5*Sb)
4850-
rv(:,4) = pos0 + x_hat * (-0.5*Sa) + y_hat * ( 0.5*Sb)
4847+
rv(:,1) = pos0 + x_hat * (-0.5_DbKi*Sa) + y_hat * (-0.5_DbKi*Sb)
4848+
rv(:,2) = pos0 + x_hat * ( 0.5_DbKi*Sa) + y_hat * (-0.5_DbKi*Sb)
4849+
rv(:,3) = pos0 + x_hat * ( 0.5_DbKi*Sa) + y_hat * ( 0.5_DbKi*Sb)
4850+
rv(:,4) = pos0 + x_hat * (-0.5_DbKi*Sa) + y_hat * ( 0.5_DbKi*Sb)
48514851

48524852
! Unit normal vector of side walls
4853-
n(:,1) = y_hat + k_hat * 0.5 * dSbdl
4854-
n(:,2) = -x_hat + k_hat * 0.5 * dSadl
4855-
n(:,3) = -y_hat + k_hat * 0.5 * dSbdl
4856-
n(:,4) = x_hat + k_hat * 0.5 * dSadl
4853+
n(:,1) = y_hat + k_hat * 0.5_DbKi * dSbdl
4854+
n(:,2) = -x_hat + k_hat * 0.5_DbKi * dSadl
4855+
n(:,3) = -y_hat + k_hat * 0.5_DbKi * dSbdl
4856+
n(:,4) = x_hat + k_hat * 0.5_DbKi * dSadl
48574857

48584858
! Check for and count vertices in water
48594859
numVInWtr = 0
48604860
do i = 1,4
4861-
vInWtr(i) = ( Dot_Product(rv(:,i)-rFS,nFS) <= 0.0 )
4861+
vInWtr(i) = ( Dot_Product(rv(:,i)-rFS,nFS) <= 0.0_DbKi )
48624862
if ( vInWtr(i) ) then
48634863
numVInWtr = numVInWtr + 1
48644864
end if
@@ -4874,7 +4874,7 @@ SUBROUTINE GetSectionHstLds_Rec( origin, pos0, k_hat, x_hat, y_hat, Sa, Sb, dSad
48744874
end if
48754875

48764876
! Compute the hydrostatic force and moment on each of the 4 sides
4877-
dFdl = 0.0
4877+
dFdl = 0.0_DbKi
48784878
do i = 1,4
48794879
v1 = i
48804880
if (i == 4) then
@@ -4896,29 +4896,29 @@ SUBROUTINE GetSectionHstLds_Rec( origin, pos0, k_hat, x_hat, y_hat, Sa, Sb, dSad
48964896
if (EqualRealNos(s,0.0_DbKi)) cycle
48974897
if ( vInWtr(v1) .and. vInWtr(v2) ) then
48984898
! Side fully submerged
4899-
s1 = 0.0
4899+
s1 = 0.0_DbKi
49004900
s2 = s;
49014901
else if ( vInWtr(v1) .or. vInWtr(v2) ) then
49024902
! Side partially wetted
49034903
sInt = s * DOT_PRODUCT(rFS-rv(:,v1),nFS) / DOT_PRODUCT(rv(:,v2)-rv(:,v1),nFS)
49044904
if ( vInWtr(v1) ) then
4905-
s1 = 0.0
4905+
s1 = 0.0_DbKi
49064906
s2 = sInt
49074907
else
49084908
s1 = sInt
49094909
s2 = s
49104910
end if
49114911
else
49124912
! Side fully out of water
4913-
s1 = 0.0;
4914-
s2 = 0.0;
4913+
s1 = 0.0_DbKi;
4914+
s2 = 0.0_DbKi;
49154915
end if
49164916

49174917
dFdl(1:3) = dFdl(1:3) + &
4918-
-n(:,i) * ( z1*(s2-s1) + 0.5*(z2-z1)/s*(s2*s2-s1*s1) )
4919-
C(1) = (z2-z1)*(x2-x1)/3.0/(s*s)*(s2**3-s1**3) + 0.5*((z2-z1)*(x1-x0)/s+(x2-x1)*z1/s)*(s2*s2-s1*s1) + z1*(x1-x0)*(s2-s1)
4920-
C(2) = (z2-z1)*(y2-y1)/3.0/(s*s)*(s2**3-s1**3) + 0.5*((z2-z1)*(y1-y0)/s+(y2-y1)*z1/s)*(s2*s2-s1*s1) + z1*(y1-y0)*(s2-s1)
4921-
C(3) = (z2-z1)*(z2-z1)/3.0/(s*s)*(s2**3-s1**3) + 0.5*((z2-z1)*(z1-z0)/s+(z2-z1)*z1/s)*(s2*s2-s1*s1) + z1*(z1-z0)*(s2-s1)
4918+
-n(:,i) * ( z1*(s2-s1) + 0.5_DbKi*(z2-z1)/s*(s2*s2-s1*s1) )
4919+
C(1) = (z2-z1)*(x2-x1)/3.0_DbKi/(s*s)*(s2**3-s1**3) + 0.5_DbKi*((z2-z1)*(x1-x0)/s+(x2-x1)*z1/s)*(s2*s2-s1*s1) + z1*(x1-x0)*(s2-s1)
4920+
C(2) = (z2-z1)*(y2-y1)/3.0_DbKi/(s*s)*(s2**3-s1**3) + 0.5_DbKi*((z2-z1)*(y1-y0)/s+(y2-y1)*z1/s)*(s2*s2-s1*s1) + z1*(y1-y0)*(s2-s1)
4921+
C(3) = (z2-z1)*(z2-z1)/3.0_DbKi/(s*s)*(s2**3-s1**3) + 0.5_DbKi*((z2-z1)*(z1-z0)/s+(z2-z1)*z1/s)*(s2*s2-s1*s1) + z1*(z1-z0)*(s2-s1)
49224922
dFdl(4:6) = dFdl(4:6) - CROSS_PRODUCT(C,n(:,i))
49234923

49244924
end do
@@ -4951,16 +4951,16 @@ SUBROUTINE getElementHstLds_Mod2_Cyl( pos1In, pos2In, FSPtIn, k_hatIn, y_hatIn,
49514951
INTEGER(IntKi) :: errStat2
49524952
CHARACTER(ErrMsgLen) :: errMsg2
49534953
ErrStat = ErrID_None
4954-
ErrMsg = ""
4954+
ErrMsg = ""
49554955

49564956
pos1 = REAL(pos1In,DbKi)
49574957
pos2 = REAL(pos2In,DbKi)
49584958
r1 = REAL(r1In,DbKi)
49594959
r2 = REAL(r2In,DbKi)
49604960
dl = REAL(dlIn,DbKi)
49614961
dRdl = (r2-r1)/dl
4962-
rMid = 0.5*( r1+ r2)
4963-
posMid = 0.5*(pos1In+pos2In)
4962+
rMid = 0.5_DbKi*( r1+ r2)
4963+
posMid = 0.5_DbKi*(pos1In+pos2In)
49644964
FSPt = REAL(FSPtIn,DbKi)
49654965
k_hat = REAL(k_hatIn,DbKi)
49664966
y_hat = REAL(y_hatIn,DbKi)
@@ -4997,11 +4997,11 @@ SUBROUTINE getElementHstLds_Mod2_Cyl( pos1In, pos2In, FSPtIn, k_hatIn, y_hatIn,
49974997
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
49984998

49994999
! Distribute the hydrostatic load to the two end nodes
5000-
F_B1(1:3) = 0.5 * F_B(1:3)
5001-
F_B2(1:3) = 0.5 * F_B(1:3)
5000+
F_B1(1:3) = 0.5_DbKi * F_B(1:3)
5001+
F_B2(1:3) = 0.5_DbKi * F_B(1:3)
50025002
F_B(4:6) = F_B(4:6) - CROSS_PRODUCT(k_hat*dl,F_B2(1:3))
5003-
F_B1(4:6) = 0.5 * F_B(4:6)
5004-
F_B2(4:6) = 0.5 * F_B(4:6)
5003+
F_B1(4:6) = 0.5_DbKi * F_B(4:6)
5004+
F_B2(4:6) = 0.5_DbKi * F_B(4:6)
50055005

50065006
END SUBROUTINE getElementHstLds_Mod2_Cyl
50075007

@@ -5040,9 +5040,9 @@ SUBROUTINE getElementHstLds_Mod2_Rec( pos1In, pos2In, FSPtIn, k_hatIn, x_hatIn,
50405040
dl = REAL(dlIn,DbKi)
50415041
dSadl = (Sa2-Sa1)/dl
50425042
dSbdl = (Sb2-Sb1)/dl
5043-
SaMid = 0.5*( Sa1+ Sa2)
5044-
SbMid = 0.5*( Sb1+ Sb2)
5045-
posMid = 0.5*(pos1In+pos2In)
5043+
SaMid = 0.5_DbKi*( Sa1+ Sa2)
5044+
SbMid = 0.5_DbKi*( Sb1+ Sb2)
5045+
posMid = 0.5_DbKi*(pos1In+pos2In)
50465046
FSPt = REAL(FSPtIn,DbKi)
50475047
k_hat = REAL(k_hatIn,DbKi)
50485048
x_hat = REAL(x_hatIn,DbKi)
@@ -5077,11 +5077,11 @@ SUBROUTINE getElementHstLds_Mod2_Rec( pos1In, pos2In, FSPtIn, k_hatIn, x_hatIn,
50775077
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
50785078

50795079
! Distribute the hydrostatic load to the two end nodes
5080-
F_B1(1:3) = 0.5 * F_B(1:3)
5081-
F_B2(1:3) = 0.5 * F_B(1:3)
5080+
F_B1(1:3) = 0.5_DbKi * F_B(1:3)
5081+
F_B2(1:3) = 0.5_DbKi * F_B(1:3)
50825082
F_B(4:6) = F_B(4:6) - CROSS_PRODUCT(k_hat*dl,F_B2(1:3))
5083-
F_B1(4:6) = 0.5 * F_B(4:6)
5084-
F_B2(4:6) = 0.5 * F_B(4:6)
5083+
F_B1(4:6) = 0.5_DbKi * F_B(4:6)
5084+
F_B2(4:6) = 0.5_DbKi * F_B(4:6)
50855085

50865086
END SUBROUTINE getElementHstLds_Mod2_Rec
50875087

@@ -5128,10 +5128,10 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Cyl( origin, pos1, posMid, pos2, FSPt,
51285128
ErrStat = ErrID_None
51295129
ErrMsg = ""
51305130

5131-
posMidL = 0.5*(pos1+posMid)
5132-
posMidR = 0.5*(posMid+pos2)
5133-
rMidL = 0.5*(r1+rMid)
5134-
rMidR = 0.5*(rMid+r2)
5131+
posMidL = 0.5_DbKi*(pos1+posMid)
5132+
posMidR = 0.5_DbKi*(posMid+pos2)
5133+
rMidL = 0.5_DbKi*(r1+rMid)
5134+
rMidR = 0.5_DbKi*(rMid+r2)
51355135

51365136
! Avoid sections coincident with the SWL
51375137
IF ( ABS(k_hat(3)) > 0.999999_ReKi ) THEN ! Vertical member
@@ -5144,7 +5144,7 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Cyl( origin, pos1, posMid, pos2, FSPt,
51445144
END IF
51455145

51465146
! Total hydrostatic load on the element (Simpsons Rule)
5147-
F_B_3pt = (dFdl1 + 4.0*dFdlMid + dFdl2) * dl/6.0
5147+
F_B_3pt = (dFdl1 + 4.0_DbKi*dFdlMid + dFdl2) * dl/6.0_DbKi
51485148

51495149
! Mid point of left section
51505150
CALL GetSectionFreeSurfaceIntersects_Cyl( posMidL, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), rMidL, theta1, theta2, secStatMidL)
@@ -5154,7 +5154,7 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Cyl( origin, pos1, posMid, pos2, FSPt,
51545154
CALL GetSectionFreeSurfaceIntersects_Cyl( posMidR, FSPt, REAL(k_hat,ReKi), REAL(y_hat,ReKi), REAL(z_hat,ReKi), REAL(n_hat,ReKi), rMidR, theta1, theta2, secStatMidR)
51555155
CALL GetSectionHstLds_Cyl( origin, posMidR, k_hat, y_hat, z_hat, rMidR, dRdl, theta1, theta2, dFdlMidR)
51565156

5157-
F_B_5pt = (dFdl1 + 4.0*dFdlMidL + 2.0*dFdlMid + 4.0*dFdlMidR + dFdl2) * dl/12.0
5157+
F_B_5pt = (dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*dFdlMid + 4.0_DbKi*dFdlMidR + dFdl2) * dl/12.0_DbKi
51585158

51595159
error = ABS(F_B_3pt - F_B_5pt)
51605160
tolMet = .TRUE.
@@ -5177,8 +5177,8 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Cyl( origin, pos1, posMid, pos2, FSPt,
51775177
END IF
51785178

51795179
IF (refine) THEN ! Recursively refine the load integration if tolerance not met
5180-
CALL RefineElementHstLds_Cyl(origin,pos1,posMidL,posMid,FSPt,r1,rMidL,rMid,0.5*dl,dRdl,secStat1,secStatMidL,secStatMid,k_hat,y_hat,z_hat,n_hat,dFdl1,dFdlMidL,dFdlMid, recurLvl+1, tmp, ErrStat, ErrMsg)
5181-
CALL RefineElementHstLds_Cyl(origin,posMid,posMidR,pos2,FSPt,rMid,rMidR,r2,0.5*dl,dRdl,secStatMid,secStatMidR,secStat2,k_hat,y_hat,z_hat,n_hat,dFdlMid,dFdlMidR,dFdl2, recurLvl+1, F_B_5pt, ErrStat, ErrMsg)
5180+
CALL RefineElementHstLds_Cyl(origin,pos1,posMidL,posMid,FSPt,r1,rMidL,rMid,0.5_DbKi*dl,dRdl,secStat1,secStatMidL,secStatMid,k_hat,y_hat,z_hat,n_hat,dFdl1,dFdlMidL,dFdlMid, recurLvl+1, tmp, ErrStat, ErrMsg)
5181+
CALL RefineElementHstLds_Cyl(origin,posMid,posMidR,pos2,FSPt,rMid,rMidR,r2,0.5_DbKi*dl,dRdl,secStatMid,secStatMidR,secStat2,k_hat,y_hat,z_hat,n_hat,dFdlMid,dFdlMidR,dFdl2, recurLvl+1, F_B_5pt, ErrStat, ErrMsg)
51825182
F_B_5pt = F_B_5pt + tmp
51835183
END IF
51845184

@@ -5227,12 +5227,12 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec( origin, pos1, posMid, pos2, FSPt,
52275227
ErrStat = ErrID_None
52285228
ErrMsg = ""
52295229

5230-
posMidL = 0.5*(pos1+posMid)
5231-
posMidR = 0.5*(posMid+pos2)
5232-
SaMidL = 0.5*(Sa1+SaMid)
5233-
SbMidL = 0.5*(Sb1+SbMid)
5234-
SaMidR = 0.5*(SaMid+Sa2)
5235-
SbMidR = 0.5*(SbMid+Sb2)
5230+
posMidL = 0.5_DbKi*(pos1+posMid)
5231+
posMidR = 0.5_DbKi*(posMid+pos2)
5232+
SaMidL = 0.5_DbKi*(Sa1+SaMid)
5233+
SbMidL = 0.5_DbKi*(Sb1+SbMid)
5234+
SaMidR = 0.5_DbKi*(SaMid+Sa2)
5235+
SbMidR = 0.5_DbKi*(SbMid+Sb2)
52365236

52375237
! Avoid sections coincident with the SWL
52385238
IF ( ABS(k_hat(3)) > 0.999999_ReKi ) THEN ! Vertical member
@@ -5245,15 +5245,15 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec( origin, pos1, posMid, pos2, FSPt,
52455245
END IF
52465246

52475247
! Total hydrostatic load on the element (Simpsons Rule)
5248-
F_B_3pt = (dFdl1 + 4.0*dFdlMid + dFdl2) * dl/6.0
5248+
F_B_3pt = (dFdl1 + 4.0_DbKi*dFdlMid + dFdl2) * dl/6.0_DbKi
52495249

52505250
! Mid point of left section
52515251
CALL GetSectionHstLds_Rec( origin, posMidL, k_hat, x_hat, y_hat, SaMidL, SbMidL, dSadl, dSbdl, FSPt, n_hat, dFdlMidL, secStatMidL)
52525252

52535253
! Mid point of right section
52545254
CALL GetSectionHstLds_Rec( origin, posMidR, k_hat, x_hat, y_hat, SaMidR, SbMidR, dSadl, dSbdl, FSPt, n_hat, dFdlMidR, secStatMidR)
52555255

5256-
F_B_5pt = (dFdl1 + 4.0*dFdlMidL + 2.0*dFdlMid + 4.0*dFdlMidR + dFdl2) * dl/12.0
5256+
F_B_5pt = (dFdl1 + 4.0_DbKi*dFdlMidL + 2.0_DbKi*dFdlMid + 4.0_DbKi*dFdlMidR + dFdl2) * dl/12.0_DbKi
52575257

52585258
error = ABS(F_B_3pt - F_B_5pt)
52595259
tolMet = .TRUE.
@@ -5276,9 +5276,9 @@ RECURSIVE SUBROUTINE RefineElementHstLds_Rec( origin, pos1, posMid, pos2, FSPt,
52765276
END IF
52775277

52785278
IF (refine) THEN ! Recursively refine the load integration if tolerance not met
5279-
CALL RefineElementHstLds_Rec(origin,pos1,posMidL,posMid,FSPt,Sa1,SaMidL,SaMid,Sb1,SbMidL,SbMid,0.5*dl,dSadl,dSbdl, &
5279+
CALL RefineElementHstLds_Rec(origin,pos1,posMidL,posMid,FSPt,Sa1,SaMidL,SaMid,Sb1,SbMidL,SbMid,0.5_DbKi*dl,dSadl,dSbdl, &
52805280
secStat1,secStatMidL,secStatMid,k_hat,x_hat,y_hat,n_hat,dFdl1,dFdlMidL,dFdlMid,recurLvl+1,tmp,ErrStat,ErrMsg)
5281-
CALL RefineElementHstLds_Rec(origin,posMid,posMidR,pos2,FSPt,SaMid,SaMidR,Sa2,SbMid,SbMidR,Sb2,0.5*dl,dSadl,dSbdl, &
5281+
CALL RefineElementHstLds_Rec(origin,posMid,posMidR,pos2,FSPt,SaMid,SaMidR,Sa2,SbMid,SbMidR,Sb2,0.5_DbKi*dl,dSadl,dSbdl, &
52825282
secStatMid,secStatMidR,secStat2,k_hat,x_hat,y_hat,n_hat,dFdlMid,dFdlMidR,dFdl2,recurLvl+1,F_B_5pt,ErrStat,ErrMsg)
52835283
F_B_5pt = F_B_5pt + tmp
52845284
END IF
@@ -5321,32 +5321,37 @@ SUBROUTINE GetEndPlateHstLds_Cyl(pos0, k_hat, y_hat, z_hat, R, theta1, theta2, F
53215321
dy = y2-y1
53225322
sz = z1+z2
53235323
dy_3 = y2*y2*y2-y1*y1*y1
5324+
dz = z2-z1
53245325
dz_2 = z2_2-z1_2
53255326
dz_3 = z2_3-z1_3
53265327
dz_4 = z2_4-z1_4
53275328
tmp1 = y1*z2-y2*z1
53285329
tmp2 = z1_2+z1*z2+z2_2
53295330

53305331
! End plate force in the k_hat direction
5331-
Fk = -0.5*Z0*(R_2*dTheta-tmp1) + cosPhi/6.0*( 2.0*dy_3 - z1*z2*dy - z1_2*(y2+2.0*y1) + z2_2*(y1+2.0*y2) )
5332+
Fk = -0.5_DbKi*Z0*(R_2*dTheta-tmp1) + cosPhi/6.0_DbKi*( 2.0_DbKi*dy_3 - z1*z2*dy - z1_2*(y2+2.0_DbKi*y1) + z2_2*(y1+2.0_DbKi*y2) )
53325333
F(1:3) = p%WaveField%WtrDens * g * Fk * k_hat
53335334

53345335
! End plate moment in the y_hat and z_hat direction
5335-
My = Z0/6.0*( 2.0*dy_3 + 2.0*dy*tmp2 + 3.0*tmp1*sz ) & ! y_hat component
5336-
+ cosPhi/24.0*( -3.0*R_4*dTheta + 3.0*y1*z1*(2.0*z1_2-R_2) - 3.0*y2*z2*(2.0*z2_2-R_2) &
5337-
+ 6.0*dy*sz*(z1_2+z2_2) + 8.0*tmp1*tmp2 )
5338-
IF (EqualRealNos(z1, z2)) THEN ! z_hat component (Nonzero only when z1 /= z2)
5339-
Mz = 0.0
5340-
ELSE
5341-
dz = z2-z1
5342-
a = dy/dz
5343-
b = tmp1/dz
5344-
tmp1 = a*a+1.0
5345-
tmp2 = a*b
5346-
tmp3 = b*b-R_2
5347-
Mz = -Z0/ 6.0*( tmp1*dz_3 + 3.0*tmp2*dz_2 + 3.0*tmp3*dz ) &
5348-
-cosPhi/24.0*(3.0*tmp1*dz_4 + 8.0*tmp2*dz_3 + 6.0*tmp3*dz_2)
5349-
END IF
5336+
My = Z0/6.0_DbKi*( 2.0_DbKi*dy_3 + 2.0_DbKi*dy*tmp2 + 3.0_DbKi*tmp1*sz ) & ! y_hat component
5337+
+ cosPhi/24.0_DbKi*( -3.0_DbKi*R_4*dTheta + 3.0_DbKi*y1*z1*(2.0_DbKi*z1_2-R_2) - 3.0_DbKi*y2*z2*(2.0_DbKi*z2_2-R_2) &
5338+
+ 6.0_DbKi*dy*sz*(z1_2+z2_2) + 8.0_DbKi*tmp1*tmp2 )
5339+
! IF (EqualRealNos(z1, z2)) THEN ! z_hat component (Nonzero only when z1 /= z2)
5340+
! Mz = 0.0
5341+
! ELSE
5342+
! dz = z2-z1
5343+
! a = dy/dz
5344+
! b = tmp1/dz
5345+
! tmp1 = a*a+1.0
5346+
! tmp2 = a*b
5347+
! tmp3 = b*b-R_2
5348+
! Mz = -Z0/ 6.0*( tmp1*dz_3 + 3.0*tmp2*dz_2 + 3.0*tmp3*dz ) &
5349+
! -cosPhi/24.0*(3.0*tmp1*dz_4 + 8.0*tmp2*dz_3 + 6.0*tmp3*dz_2)
5350+
! END IF
5351+
5352+
Mz = -Z0/ 6.0_DbKi*( dz*( y2*y2 + y2*y1 + y1*y1 + tmp2 - 3.0_DbKi*R_2 ) ) &
5353+
-cosPhi/24.0_DbKi*( dz_2*(3.0_DbKi*z2_2+3.0_DbKi*z1_2+2.0_DbKi*y1*y2-6.0_DbKi*R_2) + dz*dz*(y1*y1-y2*y2) + 4.0_DbKi*dz*(y1*y1*z1+y2*y2*z2) )
5354+
53505355
F(4:6) = p%WaveField%WtrDens * g * (My*y_hat + Mz*z_hat)
53515356

53525357
END SUBROUTINE GetEndPlateHstLds_Cyl

0 commit comments

Comments
 (0)