Merge pull request #22 from sandialabs/RM2

Updates for added mass

src/structs.jl

@@ -19,6 +19,7 @@ mutable struct FEAModel
     BC
     nodalTerms
     numModes
+    AddedMass_Coeff_Ca
 end
 
 # this way you can use defaults and pass in what is different, and it's mapped
@@ -54,7 +55,8 @@ end
         maxNumLoadSteps = 20,
         minLoadStepDelta = 0.0500,
         minLoadStep = 0.0500,
-        prescribedLoadStep = 0.0)
+        prescribedLoadStep = 0.0,
+        AddedMass_Coeff_Ca = 0.0)
 
 Model inputs for FEA analysis, struct
 
@@ -92,6 +94,7 @@ Model inputs for FEA analysis, struct
 * `minLoadStep`: used in static (steady state) analysis
 * `prescribedLoadStep`: used in static (steady state) analysis
 * `predef::Bool`: will update the elStorage array if passed into Unsteady() with the nonlinear strain stiffening, to be used for subsequent analyses
+* `AddedMass_Coeff_Ca::Float64`: added mass coefficient, also used as a flag depending on if it is 0.0 or not
 
 # Outputs:
 * `none`:
@@ -129,7 +132,8 @@ function FEAModel(;analysisType = "TNB",
     minLoadStepDelta = 0.0500,
     minLoadStep = 0.0500,
     prescribedLoadStep = 0.0,
-    predef = false)
+    predef = false,
+    AddedMass_Coeff_Ca = 0.0)
 
     if jointTransform==0.0
         jointTransform, reducedDOFList = OWENSFEA.createJointTransform(joint,numNodes,numDOFPerNode) #creates a joint transform to constrain model degrees of freedom (DOF) consistent with joint constraints
@@ -159,7 +163,7 @@ function FEAModel(;analysisType = "TNB",
 
     return FEAModel(analysisType,initCond,aeroElasticOn,guessFreq,airDensity,
     gravityOn,nlOn,spinUpOn,outFilename,RayleighAlpha,RayleighBeta,elementOrder,joint,
-    platformTurbineConnectionNodeNumber,jointTransform,reducedDOFList,nlParams,BC,nodalTerms,numModes)
+    platformTurbineConnectionNodeNumber,jointTransform,reducedDOFList,nlParams,BC,nodalTerms,numModes,AddedMass_Coeff_Ca)
 end
 
 """
@@ -567,8 +571,11 @@ mutable struct SectionPropsArray
     aeroCenterOffset
     xaf
     yaf
+    added_M22
+    added_M33
 end
-SectionPropsArray(ac,twist,rhoA,EIyy,EIzz,GJ,EA,rhoIyy,rhoIzz,rhoJ,zcm,ycm,a,EIyz,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,rhoIyz,b,a0,aeroCenterOffset) = SectionPropsArray(ac,twist,rhoA,EIyy,EIzz,GJ,EA,rhoIyy,rhoIzz,rhoJ,zcm,ycm,a,EIyz,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,rhoIyz,b,a0,aeroCenterOffset,nothing,nothing) #convenience function
+SectionPropsArray(ac,twist,rhoA,EIyy,EIzz,GJ,EA,rhoIyy,rhoIzz,rhoJ,zcm,ycm,a,EIyz,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,rhoIyz,b,a0,aeroCenterOffset) = SectionPropsArray(ac,twist,rhoA,EIyy,EIzz,GJ,EA,rhoIyy,rhoIzz,rhoJ,zcm,ycm,a,EIyz,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,rhoIyz,b,a0,aeroCenterOffset,nothing,nothing,zero(rhoA),zero(rhoA)) #convenience function
+SectionPropsArray(ac,twist,rhoA,EIyy,EIzz,GJ,EA,rhoIyy,rhoIzz,rhoJ,zcm,ycm,a,EIyz,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,rhoIyz,b,a0,aeroCenterOffset,xaf,yaf) = SectionPropsArray(ac,twist,rhoA,EIyy,EIzz,GJ,EA,rhoIyy,rhoIzz,rhoJ,zcm,ycm,a,EIyz,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6,rhoIyz,b,a0,aeroCenterOffset,xaf,yaf,zero(rhoA),zero(rhoA)) #convenience function
 
 """
 Internal, see ?Ort and ?SectionPropsArray

src/timoshenko.jl

@@ -147,6 +147,8 @@ function calculateTimoshenkoElementInitialRun(elementOrder,modalFlag,xloc,sectio
         rhoIzz = interpolateVal(sectionProps.rhoIzz,N)
         rhoJ   = interpolateVal(sectionProps.rhoJ,N)
         rhoIyz = interpolateVal(sectionProps.rhoIyz,N)
+        added_M22 = interpolateVal(sectionProps.added_M22,N) 
+        added_M33 = interpolateVal(sectionProps.added_M33,N)
 
         vprime = 0.0 #set to zero to deactivate nonlinearites from initial element calculations
         wprime = 0.0
@@ -188,9 +190,9 @@ function calculateTimoshenkoElementInitialRun(elementOrder,modalFlag,xloc,sectio
         calculateElement1!(rhoA,integrationFactor,N1,N1,M11)
         calculateElement1!(rhoA*zcm,integrationFactor,N1,N5,M15)
         calculateElement1!(-rhoA*ycm,integrationFactor,N1,N6,M16)
-        calculateElement1!(rhoA,integrationFactor,N2,N2,M22)
+        calculateElement1!(rhoA+added_M22,integrationFactor,N2,N2,M22)
         calculateElement1!(-rhoA*zcm,integrationFactor,N2,N4,M24)
-        calculateElement1!(rhoA,integrationFactor,N3,N3,M33)
+        calculateElement1!(rhoA+added_M33,integrationFactor,N3,N3,M33)
         calculateElement1!(rhoA*ycm,integrationFactor,N3,N4,M34)
         calculateElement1!(rhoJ,integrationFactor,N4,N4,M44)
         calculateElement1!(rhoIyy,integrationFactor,N5,N5,M55)
@@ -550,15 +552,15 @@ function calculateTimoshenkoElementNL(input,elStorage;predef=nothing)
         a_z_n = input.gravityOn[3]
     end
 
-    a_x = accelVec[1] #acceleration of body in hub frame (from platform rigid body motion)
-    a_y = accelVec[2]
-    a_z = accelVec[3]
+    a_x_body = accelVec[1] #acceleration of body in hub frame (from platform rigid body motion)
+    a_y_body = accelVec[2]
+    a_z_body = accelVec[3]
 
-    a_temp = CN2H*[a_x_n; a_y_n; a_z_n]
+    a_grav = CN2H*[a_x_n; a_y_n; a_z_n]
 
-    a_x = a_x + a_temp[1]
-    a_y = a_y + a_temp[2]
-    a_z = a_z + a_temp[3]
+    a_x = a_x_body + a_grav[1]
+    a_y = a_y_body + a_grav[2]
+    a_z = a_z_body + a_grav[3]
 
 
     #Integration loop
@@ -575,6 +577,8 @@ function calculateTimoshenkoElementNL(input,elStorage;predef=nothing)
 
         #..... interpolate for value at quad point .....
         rhoA = interpolateVal(sectionProps.rhoA,N) #struct mass terms
+        added_M22 = interpolateVal(sectionProps.added_M22,N) 
+        added_M33 = interpolateVal(sectionProps.added_M33,N)
 
         # Only used if (useDisp || preStress)
         EA = interpolateVal(sectionProps.EA,N)
@@ -631,8 +635,9 @@ function calculateTimoshenkoElementNL(input,elStorage;predef=nothing)
         zbarlocal = posLocal[3]
 
         fx = rhoA*a_x #let these loads be defined in the inertial frame
-        fy = rhoA*a_y
-        fz = rhoA*a_z
+        fy = rhoA*a_y + added_M22*a_y_body
+        fz = rhoA*a_z + added_M33*a_z_body
+
         rvec = [ 0; ycm; zcm]
 
         fi_hub = [fx;fy;fz]