Drive shaft¶
The theory and objectives underlaying this example is described in the course Lightweight Design.
The purpose of the model is to find the critical load for torsional buckling of a drive shaft shaped as a thin-walled tube.
The drive shaft is modeled by shell elements, where reference points at both ends are connected to respective edges through a multipoint constraint (MPC) type Beam. Boundary conditions are imposed on the reference points, as well as the applied torque.
The buckling load is equal to the applied torque multiplied with the eigenvalue. Note that eigenvalues can be negative, indicating the buckling load if all loads where reversed.
Main function¶
from abaqus import *
from abaqusConstants import *
def driveShaft(modelname,r, L, t, T, esize, E, v, rho):
mod = mdb.Model(name=modelname, modelType=STANDARD_EXPLICIT)
# Part
ske = mod.ConstrainedSketch(name='__profile__', sheetSize=200.0)
ske.CircleByCenterPerimeter(center=(0.0, 0.0), point1=(0.0, r))
prt = mod.Part(name='Tube', dimensionality=THREE_D, type=DEFORMABLE_BODY)
prt.BaseShellExtrude(sketch=ske, depth=L)
del mod.sketches['__profile__']
# Sets of edges
for edge in prt.edges:
if edge.pointOn[0][2] == 0:
prt.Set(name='Edge-1', edges=prt.edges[edge.index:edge.index+1])
if edge.pointOn[0][2] == L:
prt.Set(name='Edge-2', edges=prt.edges[edge.index:edge.index+1])
# Material and section
mod.Material(name='Mat1')
mod.materials['Mat1'].Density(table=((rho, ), ))
mod.materials['Mat1'].Elastic(table=((E, v), ))
mod.HomogeneousShellSection(name='Sec-1', material='Mat1', thickness=t)
region = prt.Set(faces=prt.faces, name='Set-faces')
prt.SectionAssignment(region=region, sectionName='Sec-1', offsetType=MIDDLE_SURFACE)
prt.seedPart(size=esize, deviationFactor=0.1, minSizeFactor=0.1)
prt.generateMesh()
# Assembly
ass = mod.rootAssembly
ass.DatumCsysByDefault(CARTESIAN)
ins = ass.Instance(name='TubeIns', part=prt, dependent=ON)
# MPCs
rp1=ass.ReferencePoint(point=(0.0, 0.0, 0.0))
rreg1=ass.Set(referencePoints=(ass.referencePoints[rp1.id], ), name='RP1')
ereg1=ins.sets['Edge-1']
mod.MultipointConstraint(name='MPC-1', controlPoint=rreg1, surface=ereg1, mpcType=BEAM_MPC)
rp2=ass.ReferencePoint(point=(0.0, 0.0, L))
rreg2=ass.Set(referencePoints=(ass.referencePoints[rp2.id], ), name='RP2')
ereg2=ins.sets['Edge-2']
mod.MultipointConstraint(name='MPC-2', controlPoint=rreg2, surface=ereg2, mpcType=BEAM_MPC)
# Initial Bondary conditions
mod.DisplacementBC(name='BC-1', createStepName='Initial', region=rreg1,
u1=SET, u2=SET, u3=SET, ur1=SET, ur2=SET, ur3=SET)
mod.DisplacementBC(name='BC-2', createStepName='Initial', region=rreg2,
u1=SET, u2=SET, u3=SET, ur1=SET, ur2=SET)
# Step and applied moment
mod.BuckleStep(name='Step-Buckl', previous='Initial', numEigen=2, vectors=4, maxIterations=500)
mod.Moment(name='Torque-Buckl', createStepName='Step-Buckl', region=rreg2, cm3=T)
# Job
job = mdb.Job(name=modelname, model=modelname)
job.submit(consistencyChecking=OFF)
job.waitForCompletion()
# Result
odb = session.openOdb(name=modelname+'.odb')
try:
ev = abs(float(odb.steps['Step-Buckl'].frames[1].description.split('=')[1].strip()))
return [r, t, ev]
except:
return [0.0, 0.0, 0.0]
# Example:
res = driveShaft(modelname='M1',r=25.0, L=1.0E3, t=1.0, T=1.0E6, esize=12, E=2.0E5, v=0.31, rho=7850.0E-12)
print res
The script can easily be extended to a parametric study. For example, buckling can be very sensitive to the element size. Let everyting be fixed except the parameter ezize, and replace the two last line in the script by:
results = []
for es in [12, 10, 8, 6, 5, 4, 3, 2]:
name = 'M-ez-'+str(es)
res = driveShaft(modelname=name, r=25.0, L=1000.0, t=1.0, T=1.0E6, esize=es, E=2.0E6, v=0.31, rho=7850.0E-12)
results.append(res)
print results
Taking the print-out from the script and assign it to a variable results to be run in this notebook:
results = [[25.0, 1.0, 2.0532],
[25.0, 1.0, 1.958],
[25.0, 1.0, 1.8962],
[25.0, 1.0, 1.8537],
[25.0, 1.0, 1.836],
[25.0, 1.0, 1.8208],
[25.0, 1.0, 1.8096],
[25.0, 1.0, 1.8014]]
es = [12, 10, 8, 6, 5, 4, 3, 2] # element size
ev = [res[2] for res in results] # eigenvalues
inverse_es = [1/float(e) for e in es]
import matplotlib.pyplot as plt
plt.plot(inverse_es, ev)
plt.grid()
plt.show()
We may judge the change above 0.2 as insignificant and keep going with an element size of
esize = round(1/0.2)
print('Sufficiently coverged element size: {} mm'.format(esize))
In the lightweight desing case study, the required thicknesses for corresponding radii based on yield strength were found.
Replacing the previous loop with the following:
t = [1.1027, 0.9113, 0.7657, 0.6525, 0.5626, 0.4901, 0.4307, 0.3815, 0.3403, 0.3054, 0.2757] # from TMM4175
r = [25.0, 27.5, 30.0, 32.5, 35.0, 37.5, 40.0, 42.5, 45.0, 47.5, 50.0] # from TMM4175
results = []
i = 1
for ti, ri in zip(t,r):
name = 'M-yield-'+str(i)
res = driveShaft(modelname=name, r=ri, L=1.0E3, t=ti, T=1.0E6, esize=5.0, E=2.0E5, v=0.31, rho=7850.0E-12)
results.append(res)
i = i + 1
print results
Using the data in this notebook:
results = [[25.0, 1.1027, 2.3234],
[27.5, 0.9113, 1.5923],
[30.0, 0.7657, 1.1678],
[32.5, 0.6525, 0.91573],
[35.0, 0.5626, 0.76811],
[37.5, 0.4901, 0.63284],
[40.0, 0.4307, 0.49342],
[42.5, 0.3815, 0.39757],
[45.0, 0.3403, 0.33086],
[47.5, 0.3054, 0.28474],
[50.0, 0.2757, 0.24565]]
plt.plot([res[0] for res in results],[res[2] for res in results], label = 'Eigenvalue' )
plt.plot([res[0] for res in results],[res[1] for res in results], label = 'Thickness' )
plt.xlabel('Radius [mm]')
plt.ylabel('Eigenvalue / Thickness')
plt.legend()
plt.grid()
plt.show()
