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.

image.png

image-2.png

image-3.png

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

In [ ]:
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:

In [ ]:
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:

In [1]:
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

In [2]:
esize = round(1/0.2)
print('Sufficiently coverged element size: {} mm'.format(esize))
Sufficiently coverged element size: 5 mm

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:

In [ ]:
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:

In [3]:
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()
TOC Next Prev

Disclaimer:This site is designed for educational purposes only. There are most likely errors, mistakes, typos, and poorly crafted statements that are not detected yet... www.ntnu.edu/employees/nils.p.vedvik

Copyright 2024, All rights reserved