simple supported beam subjected to moving sprung mass

[giacuni]

Hello, I would like to ask, if you have time, to help me to solve a simple model problem in ABAQUS.
I’m trying to model a simple supported beam subjected to moving sprung mass. This should represent the movement of a train travelling over a rail with constant speed v. The rail is modeled as a beam and the train is idealized as a lumped mass supported by the suspension systems, represented by a spring, which in turn is acting on the rail.

I have created two parts:
- A beam with a reference point, which represents the point to which the spring is attached
- A reference point with a mass, which represents the lumped mass
Then I’ve created a spring between the two reference points and I’ve applied the boundary conditions and the gravity load.

The fact is that I don't think it's the right way to model it because when I run it, it doesn't converge.
Do you know how to fix it?

Here below I've posted the script that I've done.


# %%%%%%%%%%%%%%%%%% INTRO %%%%%%%%%%%%%%%%%%%%%%%%%%%%
from abaqus import *
from abaqusConstants import *
# Create a model.
Mdb()
mdb.models.changeKey(fromName='Model-1', toName='Track')
myModel = mdb.models['Track']
# Create a new viewport
session.Viewport(name='2D Railtrack', origin=(0.0, 0.0), width=150, height=214)
session.viewports['2D Railtrack'].makeCurrent()
session.viewports['2D Railtrack'].maximize()

# %%%%%%%%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%
model_length = 25 # rail length
train_velocity = 100/3.6 # speed of the train
susp_stiffness = 1594E03 # stiffness of the suspended mass
rail_modulus = 2.87E09 # Beam's elastic modulus
rail_poisson = 0.2 # Poisson's ratio
rail_inertia = 2.90 # Moment of inertia
rail_massPerUnit = 2303 # kg/m of the rail
susp_mass = 5750 # mass of the suspended mass
num_elements = 50 # number of finite elements of the rail
mesh_size = 0.5 # mesh size of the rail
timePeriod = 0.001 # time integration step
rail_heigth = 0.159 # rail height
rail_width = 12 * 4.642e-5 / (rail_heigth*rail_heigth*rail_heigth) # calculate rail width which gives the same 2nd mom of area as a pair of rails
wheel_start = 1 # starting position from the beginning of the beam of the wheel
wheel_heigth = 0.5 # heigth of the wheel
num_steps=200
max_num_increments=10000

# %%%%%%%%%%%%%%%%%%%% MATERIALS OF PARTS %%%%%%%%%%%%%%%%%%%%
import material
## ------------------------- 1. RAIL MATERIAL -------------------------
myRailSteel = myModel.Material(name='railSteel')
elasticProperties = (rail_modulus, rail_poisson)
myRailSteel.Elastic(table=(elasticProperties, ))
equiv_steel_density = rail_massPerUnit/(rail_width*rail_heigth)
myRailSteel.Density(table=((equiv_steel_density, ), )) # true density of steel
myRailProfile = myModel.RectangularProfile(name='railProfile', a=rail_width, b=rail_heigth)

# %%%%%%%%%%%%%%% GEOMETRIES OF PARTS %%%%%%%%%%%%%%%%%%%%%%%%
import part
## --------------------------- 1. RAIL PART ---------------------------
mySketch = myModel.ConstrainedSketch(name='RailProfile',sheetSize=model_length)
mySketch.Line(point1=(-model_length*0.5,0), point2=(model_length*0.5,0))
myRail = myModel.Part(dimensionality=TWO_D_PLANAR, name='Rail', type=DEFORMABLE_BODY)
myModel.parts['Rail'].BaseWire(sketch=mySketch)
railEdge = myRail.edges.findAt(((0.0, 0.0, 0.0),))
myRail.Set(edges=railEdge, name='railEdge')

railMidPoint = (0.0, 0.0, 0.0)
myRail.DatumPointByCoordinate(coords=railMidPoint)
datumid = myRail.DatumPointByCoordinate(coords=railMidPoint).id
edges = myRail.edges.findAt(coordinates=railMidPoint)
points = myRail.datum[datumid]
myRail.PartitionEdgeByPoint(edge=edges, point=points)

springEndRP = myRail.ReferencePoint(point=(-model_length*0.5+wheel_start, 0.0, 0.0))
springEndSet = myRail.Set(referencePoints=(myRail.referencePoints[springEndRP.id],), name='springEndSet')


## ------------------------- 2. LUMPED MASS PART -------------------------
mySketch = myModel.ConstrainedSketch(name='lumpedMassProfile',sheetSize=model_length)
lumpedMass = myModel.Part(dimensionality=TWO_D_PLANAR, name='lumpedMass', type=DISCRETE_RIGID_SURFACE)
lumpedMassRP = lumpedMass.ReferencePoint(point=(-model_length*0.5+wheel_start, wheel_heigth, 0.0))
lumpedMassSet = lumpedMass.Set(referencePoints=(lumpedMass.referencePoints[lumpedMassRP.id],), name='lumpedMassSet')
lumpedMass.engineeringFeatures.PointMassInertia(name='mass', region=lumpedMassSet, mass=susp_mass, alpha=0.0, composite=0.0)

# %%%%%%%%%%%%%%%%%%%%%%%%%% MESH %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
import mesh
## --------------------------- 1. RAIL MESH ---------------------------
region = (railEdge,)
elemType = mesh.ElemType(elemCode=B21, elemLibrary=STANDARD)
myRail.setElementType(regions=region, elemTypes=(elemType,))
myRail.seedEdgeByBias(biasMethod=DOUBLE, centerEdges=railEdge,
minSize=mesh_size, maxSize=mesh_size, constraint=FINER)
myRail.generateMesh()

# %%%%%%%%%%%%%%%%%%%%%%%% SECTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
import section
import regionToolset

## --------------------------- 1. RAIL SECTION ---------------------------
myRailSection = myModel.BeamSection(name='railSection', integration=DURING_ANALYSIS, poissonRatio=rail_poisson, profile='railProfile',
material='railSteel',temperatureVar=LINEAR, consistentMassMatrix=False)
region = myRail.sets['railEdge']
myRail.SectionAssignment(region=region, sectionName='railSection', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
myRail.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(0.0, 0.0, -1.0))
myRail.MaterialOrientation(region=region, orientationType=GLOBAL, axis=AXIS_3,
additionalRotationType=ROTATION_NONE, localCsys=None, fieldName='', stackDirection=STACK_3)

# %%%%%%%%%%%%%%%%%%%% ASSEMBLY %%%%%%%%%%%%%%%%%%%%%%%%%%%
import assembly
myAssembly = myModel.rootAssembly # Create part instances.

## --------------------------- 1. RAIL INSTANCE ---------------------------
railInstance = myAssembly.Instance(name='railInstance', part=myRail, dependent=ON)

## ------------------------ 2. LUMPED MASS INSTANCE -----------------------
lumpedMassInstance = myAssembly.Instance(name='lumpedMassInstance', part=lumpedMass, dependent=ON)

## ------------------------ 3. SPRING CONNECTION -----------------------
rgn1pair0=myAssembly.instances['lumpedMassInstance'].sets['lumpedMassSet']
rgn2pair0=myAssembly.instances['railInstance'].sets['springEndSet']
region=((rgn1pair0, rgn2pair0), )
myAssembly.engineeringFeatures.TwoPointSpringDashpot(
name='Spring', regionPairs=region, axis=FIXED_DOF, dof1=2, dof2=2,
springBehavior=ON, springStiffness=susp_stiffness, dashpotBehavior=OFF,
dashpotCoefficient=0)

# %%%%%%%%%%%%%%%%%%%%%%%% STEPS %%%%%%%%%%%%%%%%%%%%%%%%%%%
import step
## --------------------------- 1. GRAVITY STEP ---------------------------
myModel.StaticStep(name='gravityLoad', previous='Initial', timePeriod=1.0, initialInc=0.1,
maxNumInc=1000, description='Load the top of the ground.')

## --------------------------- 2. DYNAMIC STEP ---------------------------
timePeriod = (model_length-(wheel_start*2))/train_velocity
timeInc = timePeriod/num_steps
myModel.ImplicitDynamicsStep(name='trainMove', previous='gravityLoad', description='Train move dynamic',
timePeriod=timePeriod, initialInc=timeInc, maxInc=timeInc, maxNumInc=max_num_increments)

# %%%%%%%%%%%%%%%%%%%% INTERACTIONS %%%%%%%%%%%%%%%%%%%%%%%
## ---------------------- 1. WHEEL-RAIL INTERACTION ----------------------
myModel.ContactProperty('WheelRailInteraction')
myModel.interactionProperties['WheelRailInteraction'].TangentialBehavior(formulation=FRICTIONLESS)
myModel.interactionProperties['WheelRailInteraction'].NormalBehavior(pressureOverclosure=HARD, allowSeparation=OFF, constraintEnforcementMethod=DEFAULT)

railRightEdge = (-0.1, 0.0, 0.0)
railLeftEdge = (0.1, 0.0, 0.0)
railRightEdge = myRail.edges.findAt((railRightEdge,) )
railLeftEdge = myRail.edges.findAt((railLeftEdge,) )
railEdge = (railLeftEdge, railRightEdge, )
myRail.Surface(side1Edges=railEdge, name='railSurface')

region1=myAssembly.instances['railInstance'].surfaces['railSurface']
region2=myAssembly.instances['railInstance'].sets['springEndSet']
myModel.SurfaceToSurfaceContactStd(name='WheelRailSurf',
createStepName='gravityLoad', master=region1, slave=region2,
sliding=FINITE, enforcement=NODE_TO_SURFACE, thickness=OFF,
interactionProperty='WheelRailInteraction', surfaceSmoothing=NONE,
adjustMethod=NONE, smooth=0.2, initialClearance=OMIT, datumAxis=None,
clearanceRegion=None)

# %%%%%%%%%%%%%%%%%%%%%%%% LOADS %%%%%%%%%%%%%%%%%%%%%%%%%

# ----------------------------- 1. GRAVITY LOAD ---------------------------
import load
myModel.Gravity(name='gravity', createStepName='gravityLoad',
comp2=-9.81, distributionType=UNIFORM, field='')

# %%%%%%%%%%%%%%%%%%% BOUNDARY CONDITIONS %%%%%%%%%%%%%%%%%%%%%
# ------------------------ 1. ENCASTRES THE LEFT END ---------------------
leftRailPoint = railInstance.vertices.findAt( ((-0.5*model_length, 0.0, 0.0), ))
leftRegion = (leftRailPoint,)
myModel.DisplacementBC(name='leftFixed', createStepName='gravityLoad',
region=leftRegion, u1=0.0, u2=0.0, ur3=0.0, amplitude=UNSET, fixed=OFF,
distributionType=UNIFORM, fieldName='', localCsys=None)

# ------------------------ 2. ENCASTRES THE RIGHT END ---------------------
rightRailPoint = railInstance.vertices.findAt( ((0.5*model_length, 0.0, 0.0), ))
rightRegion = (rightRailPoint,)
myModel.DisplacementBC(name='rightFixed', createStepName='gravityLoad',
region=rightRegion, u1=0.0, u2=0.0, ur3=0.0, amplitude=UNSET, fixed=OFF,
distributionType=UNIFORM, fieldName='', localCsys=None)

# ----------------------------- 3. TRAIN VELOCITY -------------------------
region=myAssembly.instances['lumpedMassInstance'].sets['lumpedMassSet']
myModel.EncastreBC(name='trainRefFixed', createStepName='gravityLoad',region=region, localCsys=None)
mdb.models['Track'].boundaryConditions['trainRefFixed'].deactivate('trainMove') # don't propagate this condition to the next step
myModel.VelocityBC(name='trainRefConstraint', createStepName='trainMove', region=region, v1=train_velocity, v2=0.0, v3=0.0, vr1=0.0, vr2=0.0, vr3=0.0,
amplitude=UNSET, localCsys=None, distributionType=UNIFORM, fieldName='')

region=myAssembly.instances['railInstance'].sets['springEndSet']
myModel.VelocityBC(name='railRefConstraint', createStepName='trainMove', region=region, v1=train_velocity, v2=0.0, v3=0.0, vr1=0.0, vr2=0.0, vr3=0.0,
amplitude=UNSET, localCsys=None, distributionType=UNIFORM, fieldName='')



# %%%%%%%%%%%%%%%%%%%% OUTPUT %%%%%%%%%%%%%%%%%%%%%%%%%%
myModel.historyOutputRequests['H-Output-1'].setValues(frequency=10)
myModel.fieldOutputRequests['F-Output-1'].setValues(frequency=10)

# %%%%%%%%%%%%%%%%%%%%%%% JOB %%%%%%%%%%%%%%%%%%%%%%%%%%%%
import job
jobName = 'train_run'
myJob = mdb.Job(name=jobName, model='Track',description='2D Rail track model')