For verification purposes, I needed to come up with the fixed-end axial force for a beam subjected to uniform thermal expansion. It’s not rocket science, nowhere near an LPU, but a little pencil and paper and the free end axial deflection is 
Naturally, I tried this verification using OpenSees thermal capabilities, which are somewhat hit and miss. Not unexpectedly, things got interesting.
There is no thermal truss element, so I used the forceBeamColumnThermal element, which basically mirrors the standard forceBeamColumn element. Nor is there a thermal elastic section, so I used FiberSectionThermal with the ElasticThermal uniaxial material, which I added to the OpenSeesPy interpreter after a small bug fix (PR #1337).
The thermal loading is defined using the -beamThermal option, which requires temperature to be specified at the bottom and top of the section. For uniform heating through the section, the temperatures are equal. The load pattern is defined with a linear time series as for nonlinear material response, you’ll want to ramp up the temperature.
In the script below, use any reasonable values for element length, L, rectangular cross-section dimensions, b and d, elastic modulus, E, coefficient of thermal expansion, alpha, and temperature change, deltaT.
import openseespy.opensees as ops
ops.wipe()
ops.model('basic','-ndm',2,'-ndf',3)
ops.node(1,0,0); ops.fix(1,1,1,1)
ops.node(2,L,0)
ops.uniaxialMaterial('ElasticThermal',1,E,alpha)
ops.section('FiberThermal',1)
ops.patch('rect',1,2,2,-d/2,-b/2,d/2,b/2)
ops.beamIntegration('Lobatto',1,1,3)
ops.geomTransf('Linear',1)
ops.element('forceBeamColumnThermal',1,1,2,1,1)
ops.timeSeries('Linear',1)
ops.pattern('Plain',1,1)
ops.eleLoad('-ele',1,'-type','-beamThermal',deltaT,-d/2,deltaT,d/2)
Nsteps = 1
ops.integrator('LoadControl',1/Nsteps)
ops.analysis('Static','-noWarnings')
ops.analyze(Nsteps)
With the forceBeamColumnThermal element, I obtained the expected free end axial deflection of
But when I switched to dispBeamColumnThermal, the axial deflection changed based on the number of time steps, which should not happen because this analysis is linear. The results are shown below for increasing Nsteps in the script above.
Clearly, there is a bug in the dispBeamColumnThermal element. Digging a little further, analysis with this element formulation leads to an axial reaction at the fixed end, which should not happen.
I’ll keep working with this minimal example to find the issue. Should be something simple, but probably not. If you have any additional information, please let me know!
Portwood Digital 
It seems that thermal loads cannot be removed like other loads after Analyze. If we need to run other linear load cases after thermal loads, could you suggest what is the best approach to remove the thermal loads?
LikeLike
Maybe just set the thermal loads back to ambient temperature?
LikeLike
Thanks for the prompt response Prof. Scott. Here is a revised minimal code that shows there appears to be round off errors if thermal loads are simply reset to ambient temperature. Or could you suggest another way to reset the temperature loads? When the thermal loads are added to a multi-span continuous structure, the round-off errors results in large errors.
For a 100m long continuous structure with 4 spans, after setting thermal loads back to ambient temperature, I still see about 3mm displacement in the direction perpendicular to the thermal loading direction.
Thanks for helping!
L=100
E=27000*1000
alpha=1e-5
b=1
d=2.5
deltaT=50
import openseespy.opensees as ops
ops.wipe()
ops.model(‘basic’,’-ndm’,2,’-ndf’,3)
ops.node(1,0,0); ops.fix(1,1,1,1)
ops.node(2,L,0); ops.fix(2,0,0,1) # fixed for bending
ops.uniaxialMaterial(‘ElasticThermal’,1,E,alpha)
ops.section(‘FiberThermal’,1)
ops.patch(‘rect’,1,2,2,-d/2,-b/2,d/2,b/2)
ops.beamIntegration(‘Lobatto’,1,1,3)
ops.geomTransf(‘Linear’,1)
ops.element(‘forceBeamColumnThermal’,1,1,2,1,1)
ops.timeSeries(‘Linear’,1)
ops.pattern(‘Plain’,1,1)
ops.eleLoad(‘-ele’,1,’-type’,’-beamThermal’,deltaT,-d/2,0,d/2) # thermal load creates some bending
Nsteps = 1
ops.integrator(‘LoadControl’,1/Nsteps)
ops.analysis(‘Static’,’-noWarnings’)
ops.analyze(Nsteps)
disp1=ops.nodeDisp(2,1)
print (f’displacement with some thermal load: {disp1}’)
ops.wipeAnalysis()
ops.remove(‘loadPattern’, 1)
ops.pattern(“Plain”, 101, 1)
ops.eleLoad(‘-ele’,1,’-type’,’-beamThermal’,0,-d/2,0,d/2) # temperate is changed to 0 at all locations
Nsteps = 1
ops.integrator(‘LoadControl’,1/Nsteps)
ops.analysis(‘Static’,’-noWarnings’)
ops.analyze(Nsteps)
disp2=ops.nodeDisp(2,1)
print(f’displacement after temperature loading is reset to 0: {disp2}’)
LikeLike
I tried to run similar test using 3d beam-column finite elements, but fixing both ends of the beam. I divided the beam into multiple finite elements.
Unfortunately, I’m getting inconsistent results for the reaction forces at the beam ends.
For the ‘forceBeamColumnThermal’, the results are correct, but the ‘dispBeamColumnThermal’ gives me double the value of the expected reaction. This does not happen for 2d elements.
What can be the reason? I’m not sure I’m applying the thermal load for the 3d elements correctly…
Here is the python code:
import openseespy.opensees as ops
L = 1000
E = 210000
alpha = 1e-5
d = 50
b = 50
deltaT = 10
ne = 5 # Number of elements
ops.wipe()
ops.model(‘basic’,’-ndm’,3,’-ndf’,6)
ops.node(1,0,0,0)
ops.fix(1, 1, 1, 1, 1, 1, 1)
for n in range(2,ne+2):
ops.node(n,L/ne*(n-1),0,0)
ops.fix(n, 1, 1, 1, 1, 1, 1)
ops.uniaxialMaterial(‘ElasticThermal’,1,E,alpha)
ops.section(‘FiberThermal’,1, ‘-GJ’, 1.)
ops.patch(‘rect’,1,2,2,-d/2,-b/2,d/2,b/2)
ops.beamIntegration(‘Lobatto’,1,1,3)
vecxz = [0., 0., 1.]
ops.geomTransf(‘Linear’, 1, *vecxz)
for n in range(1,ne+1):
ops.element(‘dispBeamColumnThermal’,n,n,n+1,1,1)
# dispBeamColumnThermal # forceBeamColumnThermal
ops.timeSeries(‘Linear’,1)
ops.pattern(‘Plain’,1,1)
for n in range(1,ne+1):
ops.eleLoad(‘-ele’,n,’-type’,’-beamThermal’,deltaT,-d/2,deltaT,d/2)
Nsteps = 100
ops.integrator(‘LoadControl’,1/Nsteps)
ops.analysis(‘Static’,’-noWarnings’)
ops.reactions()
print(E*b*d*alpha*deltaT)
print(ops.nodeReaction(1))
print(ops.nodeReaction(ne+1))
LikeLike
I’m not sure. You can compare the 2d and 3d dispBeamColumnThernal source codes to look for differences.
LikeLike