huge differences between static gravitational load and the dynamic one

[antonioAlvarez]

Hi,

I am calculating a wind turbine blade, it is cantilevered blade in which the centre of mass of the beam varies along their length, therefore I am connecting the different beams forming the blade among them using rigidLinks. When I run the static calculations I get the proper gravitational forces, but when I run the dynamic analysis the result given by the programme is much smaller than the static one, some differences can be expected, but unfortunatelly this is not the case. I am getting exactly the same results using both Transformation and Penalty constraints. I am using the following parameters for the dynamic calculations:

ops.numberer('RCM')
ops.system('BandGeneral')
ops.test('NormDispIncr', 1.0e-6, 100, 2)
ops.algorithm('KrylovNewton')
ops.integrator('Newmark', 0.5, 0.25)
ops.analysis('Transient')

I am not having any convergence issue.

If you have any insight in what can be going on in the analysis or if I am making any mistake, please do not hesitate in sharing your opinion or asking for more info if you think it will help to clarify the issue.

Thank you in advance and kind regards,

Antonio

[antonioAlvarez]

Hi,

Maybe I have found the issue. In order to avoid it the previous static analysis should not be wiped from the analysis. I am adding a code in which this can be checked.

#############################################################################
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import os
import openseespy.opensees as ops

from scipy import interpolate

import os, sys
from scipy.integrate import odeint

def createElasticSection(tag,YoungModulus,Area,InertiaZ,InertiaY,ShearModulus,TorsionalInertia):
#tag label for the section
ops.section('Elastic',tag,YoungModulus,Area,InertiaZ,InertiaY,ShearModulus,TorsionalInertia)

def createLinearTransform(tag,vecXZ,offsetIni,offsetFin):
#tag is the label of the transformation
#vecXZ is a list of three components, containing the z local axis of the element, referred to the global coordinates
#offsetIni is a list of three components, containing the offset of the initial node, referred to the global coordinates
#offsetFin is a list of three components, containing the offset of the final node, referred to the global coordinates

ops.geomTransf('Linear',tag,*vecXZ,'-jntOffset',*offsetIni,*offsetFin)

###############################################################################################################

###############################################################################################################

def createNonLinearBeam(tag,nodeIni,nodeFin,numIntPoints,secTag,transTag,maxIter,tol,distMass):
#tag label of the element
#nodeIni label of the initial node
#nodeFin label of the final node
#numIntPoints the number of sections in the beam
#secTag labels of the sections
#transTag label of the transformation
#maxIter maximum number of iterations
#tol tolerance
#distMass distributed Mass (mass per unit length)

#print("tag= ",tag)
#print("nodeIni= ",nodeIni)
#print("nodeFin= ",nodeFin)
#print("numIntPoints= ",numIntPoints)
#print("secTag= ",secTag)
#print("transTag= ",transTag)
#print("maxIter= ",maxIter)
#print("tol= ",tol)
#print("distMass= ",distMass)

ops.element('nonlinearBeamColumn',int(tag),int(nodeIni),int(nodeFin),numIntPoints,secTag,transTag,'-iter',maxIter,tol,'-mass',distMass)

#if isinstance(secTag,list) or isinstance(secTag,np.ndarray):
# if isinstance(secTag,np.ndarray):
# secTag=list(secTag)
# if isinstance(transTag,list) or isinstance(transTag,np.ndarray):
# if isinstance(transTag,np.ndarray):
# transTag=list(transTag)
# ops.element('nonlinearBeamColumn',tag,*[nodeIni,nodeFin],len(secTag),*[secTag],*[transTag],'-iter',maxIter,tol,'-mass',distMass)
# else:
# ops.element('nonlinearBeamColumn',tag,*[nodeIni,nodeFin],len(secTag),*[secTag],transTag,'-iter',maxIter,tol,'-mass',distMass)
#else:
# ops.element('nonlinearBeamColumn',tag,*[nodeIni,nodeFin],1,secTag,transTag,'-iter',maxIter,tol,'-mass',distMass)

###############################################################################################################
###############################################################################################################
###############################################################################################################

ops.wipe()
ops.model('basic', '-ndm', 3, '-ndf', 6) #ndm is the number of dimensions and ndf is the degrees of freedom

ops.timeSeries('Constant', 1)
ops.pattern('Plain', 1, 1) #(’Plain’, patternTag, timeSeriesTag, ’-fact’, fact)
ops.constraints('Transformation')
#ops.constraints('Penalty',9e12,9e12)#usualy 3 or 4 order of magnitude bigger than EI/L works well

numel=11
height=100
h=np.linspace(0,height,num=numel)

#section properties
tagSec=1
YoungModulus=25e9
Area=1*0.3
InertiaZ=(0.3**3*1)/12
InertiaY=(0.3*1)/12
ShearModulus=YoungModulus*0.2
TorsionalInertia=0.263*(1*0.3**3)
createElasticSection(tagSec,YoungModulus,Area,InertiaZ,InertiaY,ShearModulus,TorsionalInertia)

#linear transform
tagTrans=1
vecXZ=[1,0,0]
offsetIni=[0,0,0]
offsetFin=[0,0,0]
createLinearTransform(tagTrans,vecXZ,offsetIni,offsetFin)

#generate the node for the first tower
nodeTagsTw00=[]
tagCounter=0
for i in range(numel):
ops.node(tagCounter,0.0,0.0,h)
nodeTagsTw00.append(tagCounter)
tagCounter+=1

#generate the beam
numBeams=numel-1
beamCounter=0
numIntPoints=5
maxIter=10
tol=1e-8
distMass=2500*Area
beamTags00=[]
for i in range(numBeams):
createNonLinearBeam(beamCounter,nodeTagsTw00,nodeTagsTw00[i+1],numIntPoints,tagSec,tagTrans,maxIter,tol,distMass)
beamCounter+=1
beamTags00.append(beamCounter)

ops.fix(int(nodeTagsTw00[0]), 1, 1, 1, 1, 1, 1) #all dof fixed at the base, so it is a cantiliver


#beam joined with rigid links
nodeTagsTw01=[]
for i in range(numel):
ops.node(tagCounter,25.0,0.0,h)
nodeTagsTw01.append(tagCounter)
tagCounter+=1
if i>0 and i<numel-1:
ops.node(tagCounter,25.0,0.0,h)
nodeTagsTw01.append(tagCounter)
tagCounter+=1
print("len(nodeTagsTw01)= ",len(nodeTagsTw01))

#generate the beam
beamTags01=[]
for i in range(numBeams):
#print("2*i+1= ",2*i+1)
createNonLinearBeam(beamCounter,nodeTagsTw01[2*i],nodeTagsTw01[2*i+1],numIntPoints,tagSec,tagTrans,maxIter,tol,distMass)
beamCounter+=1
beamTags01.append(beamCounter)
if i>0:
ops.rigidLink('beam',nodeTagsTw01[2*i-1],nodeTagsTw01[2*i])


ops.fix(int(nodeTagsTw01[0]), 1, 1, 1, 1, 1, 1) #all dof fixed at the base, so it is a cantiliver
##############################################################################################
#modal analysis
ops.wipeAnalysis()

#then the number of modes is
numModes=10

#calculate the modes
eigenValues=ops.eigen('-fullGenLapack',numModes) #default solver is '-genBandArpack'

#as all the modes have the same damping
damping=[0.02]*numModes

#to apply the modal damping
ops.modalDamping(*damping)

#the frequencies of the modes is
freqRad=np.sqrt(np.array(eigenValues))
freqHz=freqRad/(2*np.pi)
##############################################################################################


#assign loads
paramsTw00=[]
paramCounter=1
l=height/numBeams
mass=l*2500*Area
for i in range(numBeams):
if i==numBeams-1:
ops.load(int(nodeTagsTw00[i+1]),*[1000,0,-mass*9.81,0,0,0])
else:
ops.load(int(nodeTagsTw00[i+1]),*[0,0,-mass*9.81,0,0,0])
ops.parameter(paramCounter,"loadPattern",1,'loadAtNode',int(nodeTagsTw00[i+1]),1)
paramsTw00.append(paramCounter)
paramCounter+=1


paramsTw01=[]
for i in range(numBeams):
if i==numBeams-1:
ops.load(int(nodeTagsTw01[-1]),*[1000,0,-mass*9.81,0,0,0])
else:
ops.load(int(nodeTagsTw01[2*i-1]),*[0,0,-mass*9.81,0,0,0])
ops.parameter(paramCounter,"loadPattern",1,'loadAtNode',int(nodeTagsTw01[2*i-1]),1)
paramsTw01.append(paramCounter)
paramCounter+=1

doStatic=True
if doStatic:
ops.constraints('Transformation')
ops.integrator('LoadControl', 1.0)
ops.algorithm('Linear')
ops.analysis('Static')
ops.analyze(1)
ops.reactions()

print("\n############################\nSTATIC\n")
nodeReactionTw00=ops.nodeReaction(int(nodeTagsTw00[0]),)
nodeReactionTw01=ops.nodeReaction(int(nodeTagsTw01[0]),)
print("nodeReactionTw00= ",nodeReactionTw00)
print("nodeReactionTw01= ",nodeReactionTw01)
baseDispTw00=ops.nodeDisp(int(nodeTagsTw00[0]),)
baseDispTw01=ops.nodeDisp(int(nodeTagsTw01[0]),)
print("baseDispTw00= ",baseDispTw00)
print("baseDispTw01= ",baseDispTw01)
tipDispTw00=ops.nodeDisp(int(nodeTagsTw00[-1]),)
tipDispTw01=ops.nodeDisp(int(nodeTagsTw01[-1]),)
print("tipDispTw00= ",tipDispTw00)
print("tipDispTw01= ",tipDispTw01)

ops.constraints('Transformation')
#ops.constraints('Penalty',9e15,9e15)#usualy 3 or 4 order of magnitude bigger than EI/L works well

ops.numberer('RCM')
ops.system('BandGeneral')
#ops.system('FullGeneral')
ops.test('NormDispIncr', 1.0e-6, 100, 2)
#ops.algorithm('Newton')
#ops.integrator('Newmark', 0.5, 0.25)

##ops.system('BandGeneral') # how to store and solve the system of equations in the analysis, BandGeneral=band matrix, BandSPD,ProfileSPD=sparse matrices, UmfPack=fast sparse
#ops.system('UmfPack') #with this one do not use NewtonRapsom for algorithm, the linear algorithm will not iterate to equilibrium
##ops.algorithm('Linear') # use Linear algorithm for linear analysis
ops.algorithm('KrylovNewton') #or BFGS
#ops.algorithm('Newton') #or BFGS
ops.integrator('Newmark', 0.5, 0.25) # determine the next time step for an analysis
#ops.integrator('CentralDifference')

ops.analysis('Transient')
ops.reactions()
T=1
deltaT=0.1
Nsteps = int(T/deltaT)
print("\n############################\nTRANSIENT\n")
for i in range(Nsteps):
ops.analyze(1,deltaT)
ops.reactions()
#ops.updateParameter(1,fuerza)

print("\n############################\nt=%f\n"%(ops.getTime()))
nodeReactionTw00=ops.nodeReaction(int(nodeTagsTw00[0]),)
nodeReactionTw01=ops.nodeReaction(int(nodeTagsTw01[0]),)
print("nodeReactionTw00= ",nodeReactionTw00)
print("nodeReactionTw01= ",nodeReactionTw01)
baseDispTw00=ops.nodeDisp(int(nodeTagsTw00[0]),)
baseDispTw01=ops.nodeDisp(int(nodeTagsTw01[0]),)
print("baseDispTw00= ",baseDispTw00)
print("baseDispTw01= ",baseDispTw01)
tipDispTw00=ops.nodeDisp(int(nodeTagsTw00[-1]),)
tipDispTw01=ops.nodeDisp(int(nodeTagsTw01[-1]),)
print("tipDispTw00= ",tipDispTw00)
print("tipDispTw01= ",tipDispTw01)