|Overview example and Analysis types||Appendix|
|Solid Mechanics Boundary Conditions||References|
The analysis and behaviour of solids under loads and constraints is of fundamental importance in mechanics. Solid mechanics deals with the mechanics of solid bodies in three dimensions, while the topic of structural mechanics encompasses a wider range of objects, such as thin shells or beams, for example.
This tutorial gives an introduction to modeling solid mechanics with partial differential equations. Equations and boundary conditions that are relevant for performing solid mechanics analysis are derived and explained.
Modeling solid mechanics with partial differential equations (PDEs) is not the only way to model solid mechanics. Other techniques include setting up ordinary differential equations (ODEs). This approach is followed by the Wolfram System Modeler. Roughly speaking the system modeler approach is more suitable for large systems of solid bodies interacting, while the partial differential equation approach is more suitable for a fine grained analysis of a specific body. In some cases it is beneficial to use a combination of the two approaches.
The approach taken here is that in an introductory section a single solid body, a bookshelf bracket, is used to introduce various solid mechanics analysis types and the functionality available. This will be followed by a more theoretical explanation of the underlying ideas and concepts. The theoretical background is much easier to understand once an intuition for the various analysis types exist. After that the available boundary conditions are discussed.
Solid mechanics is typically considering three dimensional solid objects. Special cases exist to deal with two dimensional simplifications. These simplifications, however, have some pitfalls that are avoidable if the understanding for the three dimensional scenario is correct. For this reason the initial examples will be three dimensional examples. At a later stage plane (two dimensional) stress and strain and their limitations will be introduced.
The goal of a solid mechanics analysis is to find the deformation of a body under load. A subsequent step then finds strains and stresses within the deformed body. The analysis and interpretation of these physical quantities are useful to create a better quality engineering design of the body under consideration. For example, structurally weak parts of an object can be identified and improved upon.
The solid mechanics analysis process is typically done in stages. First, for the body to be analyzed a geometric model needs to be created. The geometric model is typically created within a computer aided design (CAD) process. CAD models can either be imported or created in product. To import geometries common file formats like DXF, STL or STEP are supported. These geometries can be imported with Import. The alternative is to create the geometrical models in product with, for example, by using OpenCascadeLink. Once the geometric model is made available some thought needs to be put into what what type of analysis is to be performed. Currently supported analysis types are static analysis, time dependent analysis, eigenmode analysis and frequency response analysis. The next step is the setup of boundary conditions and constraints. Materials to be used further specify the PDE model. Once the PDE model is fully specified, the subsequent finite element analysis will then compute the desired quantities of the body under investigation. These quantities are then post processed. Either by visualizing them or some derived quantities are computed. This notebooks shows the necessary steps for everything except the CAD model generation.
The accuracy and the effectiveness of the solid mechanics PDE model is validated in the separate notebook entitled Solid Mechanics Model Verification Tests.
Many of the animations of the simulation results shown in this notebook are generated with a call to Rasterize. This is to reduce the disk space this notebook requires. The downside is that the visual quality of the animations will not be as crisp as without it. To obtain high quality graphics remove or comment out the call to Rasterize.
The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.
To illustrate the usage of the finite element method in solid mechanics it is instructive to present a simple example and give an overview of the setup, various analysis types and post processing steps possible.
The first example is introduce the workflow of setting up a simple solid mechanics PDE model. Subsequent sections will then use more sophisticated and complete models. For now, consider the following setup where we have a long beam that is fixed at the left and has a downward force applied
For the following overview of available solid mechanics analysis types a predefined boundary element mesh is loaded. For more information on how this boundary element mesh was generated refer to the Bookshelf Bracket tutorial that explains the creation of this specific geometric model.
One thing to keep in mind is the scale used in the geometric model. If the length of the boundary mesh is, for example, in units of meters then the material parameters will need to be specified in consistent units. In this specific case the units of the boundary mesh are in meters so the bracket is 0.2 meters long.
The amount of detail of the geometry will effect the number of elements a finite element mesh will need to represent that geometry. The amount of elements in a mesh have a direct influence on the CPU time and memory needed to solve a particular problem. To reduce the number of elements one can consider that geometric details, such as the screw holes, often only have an influence in their closer surroundings [11, c.1]. Saint Venant's principle suggests that a detail with a characteristic length will affect the surroundings in a distance of about . If one is interested in the mechanical performance of the object in a specific region that is far enough away from a particular detail, that detail might be left out of the geometric model and thus save computational effort. Also for vibration analysis geometric details that are smaller than about 10% of the geometric cross section can usually be neglected.
More information on generating or importing 3D geometric models can be found in the Using OpenCascadeLink tutorial.
The next step is to assign material parameters. Generally, all parameters for a solid mechanics model are collected in an Association that include the necessary parameter values.
Should a material not be available or different units be needed then the material properties can be added be specifying parameters like YoungModulus and PoissonRatio directly. The exact property names needed can be found on the reference page of SolidMechanicsPDEComponent.
Material parameters where the values are given as Quantity objects can be used. Should the units of the geometry be different from the material units, then the material units can be scaled.
Internally, all material data units are converted to "SIBase" units. As a consequence the default unit of length is "Meters". If the units of the geometry are also in meters then nothing needs to be changed. If the units of the geometry are not in meters then either the PDE and material properties need to be scaled to the units the geometry or the geometry needs to be scaled to "Meters". To scale the units of the PDE and material parameters the parameter "ScaleUnits" can be given.
The default material model is a linear elastic isotropic material. As a rule of thumb, a linear elastic material models is applicable until a maximal stretching of 5% is reached [8, p. 159].
Geometries consisting of several materials can also be used and an example of such a use case is presented in the section on multiple materials.
Since solid mechanics is about the deformation of objects under load and constraints, boundary conditions are an essential component of solid mechanics modeling. A boundary load is a force or pressure that is applied on the surface of an object. For a load to be applicable to an object that object must also be constrained in some way, for example screwed to a wall, as otherwise the object would not pose a resistance to the load. Loads and constraints are set up by specifying boundary conditions. Various boundary conditions are can be used and will be discussed in more detail in the boundary condition section. For the purpose of this overview a boundary surface load and the constraints introduced by a wall and screws will be sufficient.
The purpose of this section is to establish the positions where the boundary conditions are to be applied. The positions of boundary conditions apply will remain the same, regardless of the analysis type preformed. The exact specification of the values of the boundary conditions, however, depend on the different analysis types and will follow in the respective sections.
In a next step we would like to specify a constraint that models that the bracket is mounted to a wall. To fix the bracket to the wall two constraints play a role. First, the back of the bracket cannot penetrate the wall and secondly screws fix the bracket to the wall.
Starting with the constraints the screws impose, we recall that the geometry did not provide the screw threads. We can remedy that by completely constraining the movement in all directions of the surfaces that make the screw hole. This means at the points of contact of the screw and the bracket no movement is possible at all.
The bracket cannot penetrate the wall and thus we should constraint the movement of the back of the bracket in the negative -direction. Because the two screws press the bracket to the wall and the bracket cannot bend into the wall, a reasonable approach is to also limit the movement in the positive (out from the wall) direction. This is a common simplification. If this simplification is not justified then a contact problem arises.
For more complicated geometries a different technique to specify boundary condition predicates may be appropriate and is given in the appendix in the section boundary condition predicates.
More information about the mesh generation process can be found in the ElementMesh generation tutorial. Another option is to import a mesh. Some common mesh file formats can be imported with the help of the FEMAddOns.
A solid mechanics analysis will seek the displacement of an object which is a consequence of applied forces and constraints. The dependent variables for the displacement are called , and and represent the displacements in the -, - and -axis direction, respectively.
The default setup generates a model for a linear elastic isotropic material with a small deformation assumption. This model does not include the self-weight of the object. This can easily be added with the parameter "BodyLoadValue". Adding that is explained in the section Body Load.
The result is a list of three InterpolatingFunction objects. Each InterpolatingFunction gives the displacement for its respective dependent variable. The three InterpolatingFunction objects make up the displacement vector.
More information about the solution process and its options can be found in the NDSolve Options for Finite Elements tutorial.
The primary solution of a solid mechanics PDE model is the displacement that results due to the acting forces. This displacement can be used to visualize how the body deforms under the load and constraints. Also of interest are, however, the stresses and strains in the object. The strains in the body are recovered from the displacement. The stresses are then recovered from the strains. So, both the strain and stress are secondary unknowns. The recovered strains and stresses can further more be combined in overview concepts like the equivalent strain or the von Mises stress. The terminology post processing is an umbrella term for all computations and visualizations done after the displacements have been computed.
The deformation of an object under load is simply the computed displacement added to the coordinates of the original object. A point in an object originally at is moved to when the object is under load.
The deformation shown is scaled by the minimal size of the geometries' extension and the maximal displacement. This can be switched off or changed by specifying a "ScalingFactor" for ElementMeshDeformation. For more information see the ElementMesh Visualization tutorial.
ElementMeshDeformation by default will make the deformation visible, regardless how small it is, as long as it is nonzero; however, the visualization does not depict the true deformation. The true deformation can be shown by setting "ScalingFactor"->None.
In most cases deformations are small and so "ScalingFactor"->None will not give an interesting plot.
Another consequence of the automatic scaling factor computation is that if the body load is increased the deformation plot will remain the same as the values will be scaled back again. The value of the deformation plot is in seeing how the body deforms. This can give an idea if the model works as expected.
Strain is a quantity that describes the amount of deformation within a body  and is a ratio and unit-less. Strain is not to be confused with the amount of displacement. Strain is related to the displacement by the gradient of a given displacement. The concept of strain will be explained in more detail in the theory section about strain.
The SolidMechanicsStrain function computes the normal and shear strains for each independent direction. The function returns a SymmetrizedArray that contains an expression representing the strain components in the various directions. The SymmetrizedArray enforces symmetry when present and is compact. The strains returned have the following order:
where the the normal strains are denoted by and the shear strains by . The software uses engineering strains by default and . It is possible to switch off the usage of engineering strains by setting "EngineeringStrain"->False in the parameters . This may be of interest when comparing linear material laws with nonlinear material laws, where true stain is used. More information can be found in the reference page of SolidMechanicsStrain.
Note, how in this example the absolute values of the strain in the -direction are small then absolute values in the -direction, even though the displacement in the negative -direction is much larger then in the positive -direction.
Inspecting all strain components can be cumbersome. The equivalent strain is a simplification that combines all strain components into a scalar and as such is easier to grasp. However, the equivalent strain does not include the complete picture of the strains.
Some implementations of the equivalent strain use a factor of , where is Poisson's ratio. We have chosen not to include this factor to make the equivalent strain compatible with the vonMises stress, which is the the stress analog of the equivalent strain.
Principal strain values and direction vectors
Principal strain values are the main strain values. They are the computed eigenvalues of the strain matrix and correspond to a local coordinate system where the shear strains are 0. The principal strain values are the in the local coordinates and give the load direction independent strains.
The section Boundary Load: Tension has an example that shows the use of PrincipalEigensystem to find stresses in a non-axial loaded cylinder.
Stress is a quantity the describes the distribution of internal forces within a body  and is measured in or
The Stress function computes the normal stress and shear stress for each independent direction. The function returns a SymmetrizedArray that contains an expression representing the stress component in that direction. The stresses returned have the following order:
von Mises Stress
Inspecting all stress variables can be a cumbersome. A simplification that combines all stress components in a single expression is the von Mises stress. The von Mises stress is a scalar and as such easier to grasp. However, the von Mises stress does not include the complete picture of the stresses present within a body. The von Mises stress is to the stresses the same as the equivalent strain is for the strains.
Principal stress values and direction vectors
Principal stress values are the main stress values. They are the computed eigenvalues of the stress matrix and correspond to the diagonal stress values where the shear stresses are 0. The principal stress values give stresses independent of load direction.
The section Boundary Load: Tension has an example that shows the use of PrincipalEigensystem to find stresses in a non-axial loaded cylinder.
The factor of safety is a number telling us by how large a factor the maximal resulting stresses are below a stress limit. There are multiple definitions of the factor of safety, abbreviated as FOS. Because the material at hand is a metal, which is ductile and not brittle, the von Mises failure criterion is a reasonable choice. For ductile materials failure is considered to start at the onset of plastic deformation (while for brittle materials failure is considered at fracture). For more information see the section Failure Theory.
A safety factor below one is problematic. In this case we assumed a linear stress strain relation; the maximum equivalent stress may be different if a nonlinear stress-strain relation is used. When a linear material model is used outside its realm of validity, the stresses computed are typically higher than the actual stresses. Nevertheless, finding stresses that are higher than the yield stress is something to be very cautious about.
Since the time dependent solid mechanics equations are second order time dependent, the initial conditions comprise an initial displacement and an initial velocity field. Note that the initial conditions affect the whole body. For example, specifying a non zero initial displacement means the entire body is displaced by the specified amount at the initial time.
A time dependent analysis can take a long time and require a lot of memory. For this reason it is important to be able to estimate if an analysis really needs to be transient. The time it takes for a shear wave to travel though an elastic solid with a characteristic length is approximately [11, c.1]
where is the mass density and is the shear modulus. Stress values decay to their stationary values in about . Similar considerations can be made for stresses induced by acceleration in rotating objects [11, c.1].
An eigenmode analysis computes the resonant frequencies of an object. These frequencies are also called natural frequencies. At each of these eigenfrequencies the object under investigation deforms into a distinct shape called eigenmode. The result of an eigenmode analysis is a list of frequencies and their corresponding modes.
- Objects in resonance are a common cause of failure. Knowing these frequencies helps avoid the shapes being exposed to them.
where is the mass matrix, the damping matrix, the stiffness matrix and the load vector. is the displacement vector of dependent variables , and . is the derivative of the displacement vector, the velocity vector and the second derivative of the displacement, the acceleration. For vibrational analysis the damping is generally ignored  and leaves us with:
It is worthwhile noting that contrary to other analysis types the eigenmode analysis needs to be specified as a parameter. In the other cases SolidMechanicsPDEComponent can deduce the analysis type from the variable specification. The eigenmode analysis is very similar to the time dependent setup and can not be deduced from the variable input alone.
First, we look at a contained eigenmode analysis. This means the body is constrained by boundary conditions. The only conditions relevant for a constrained eigenmode analysis are boundary conditions that result in DirichletCondition where a dependent variable is set to 0. Boundary forces and non zero SolidDisplacementCondition are not relevant.
The amount of the deformation in the deformed shapes are not actual deformations. Keep in mind that the amplitude for the deformations shown are arbitrary. Any constant times an eigenmode is still the same eigenmode. An eigenmode analysis computes the eigenfrequency and the qualitative shape of the modes at the respective frequencies.
A special situation occurs when there are no constraints on the object. Then the first modes will be zero and are called rigid body modes. We speak of a rigid body mode when a body translates or rotates without deformations. Eigen values and modes capture the fundamental shape of deformations a body is capable. A zero eigen mode accounts for the fact that the body could translate or rotate. The fact that a body can translate or rotate is no news and thus these modes are of no interest. In three dimensions there are 6 rigid body modes, 3 for the translation in each direction and 3 for a rotation around each axis.
These are the rigid body modes. A rigid body that is translated or rotated has a zero eigenvalue. All in all there are three possible translations in the corresponding -, - and -axis directions and three rotations around the same axis. The eigenmodes found by the eigen solver are not the pure rigid body modes but are a superposition of combinations of the fundamental rigid body movements.
Sometimes one would like to vary a parameter of a PDE model and repetitively solve the same PDE for a variety of parameters. A convenient way to do so is a parametric analysis. In the bookshelf example, it might be of interest how the bookshelf deforms under various surface loads. In this case the boundary load force is made a parameter. The simulation is set up in exactly the same way as in a non parametric analysis, only using the ParametricNDSolve family of functions and specifying the name of the parameter in the model.
A frequency response analysis gives information of how a specific point in a body reacts to a sweep through a frequency range. The result of a frequency response analysis is a frequency response plot showing the relation between real displacement at the range of frequencies. The frequency response analysis is also called harmonic analysis.
Before performing a frequency response analysis it can be useful to perform an eigenmode and a static analysis. The eigenmode analysis will find the critical frequencies and their related modes that will play a part in the frequency response analysis. The static analysis will provide the maximum displacement without any frequency component.
Here is the angular frequency, the imaginary init, the resulting displacement. , and are the mass, damping and stiffness of the solid mechanics PDE. These are provided by SolidMechanicsPDEComponent.in the negative -direction. In the stationary analysis a comparable boundary load of in the negative -direction was used. It is interesting to find the displacement of the query point from the stationary solution. This allows us to put the frequency response displacements in relation to the displacement of the query point in the stationary solution.
The equilibrium equations relate a force density and the stress . Loosely speaking, stress is the resistance of an internal point to the applied load. The time dependent equilibrium equation is given by:
Where is the mass density and the displacement vector. A load can either be acting on the entire body or the boundary. To simplify things we ignore the time dependent term for now. The component form of the time independent equilibrium equation in three dimensional space for the stress components and body forces are given by:
The kinematic equations relate the displacement to the strain . Strain describes the relative displacement between points in the body. There are various strain measures. Here we distinguish between two approaches. In the infinitesimal deformation theory, it is assumed that the displacements and strains are small. This is the default method in SolidMechanicsPDEComponent. The infinitesimal strain measure is given by:
For the finite deformation theory no assumptions are made. If the strain-displacement relation is nonlinear then we speak of a geometric nonlinearity. The finite deformation theory is such a geometric nonlinearity. Geometric nonlinearity is also called a nonlinear kinematic equation. A geometric nonlinearity accounts for the fact that the geometry is evolving during the loading.
The relation of the equilibrium equations - the forces - and the kinematic equations - the description of deformation - is done through the constitutive equations. The constitutive equations, also known as the material model, relate stress to strain. Generally speaking this is a relationship of the form
Besides linearity there are other properties of interest in a material law. The deformation of an object is called elastic if after removing a force the objects returns to its original configuration. We speak of a plastic deformation if the object does not return to its original configuration when loads are removed. Energy is lost whenever an object experiences plastic deformation. Plastic deformation is permanent. Plasticity is a special form of material nonlinearity.
The presented solid mechanics formulation is based on displacements as primary unknowns. This is a common approach [4, p.480], as the finite element solver will compute the displacements. The secondary unknowns such as strain and stress will be recovered from the displacements. While it is conceivable that equations are set up in such a manner that all unknowns are solved for, the resulting system of equations will become prohibitively large to solve in a reasonable time and with limited computer memory available.
Body loads are forces that act on the entire volume of the object and arise from external force fields. These forces are also called volume forces and are specified as a force density in . Body loads are often also called body forces, but really they are a force per unit volume. Body loads can be specified with the parameter "BodyLoad" and are specified as a vector field. Gravity is an example of a constant body load or a position dependent centrifugal forces in rotating objects is another example.
To model, for example, the self-weight of a body one would have to specify the product of the mass density with the gravitational acceleration such that the body load becomes . Because this is a common case an additional parameter "BodyLoadValue" can be specified. When "BodyLoadValue" is specified the mass density is multiplied automatically and yields a body load.
Kinematics is the mathematical formulation of the movement and deformation of objects. The formulations found will be independent of the forces causing these deformations. The displacement and strain are introduced and defined in the following sections.
The solution of the solid mechanics equations gives a set of three displacement functions , and which are are the displacements in the -, - and -directions, respectively. An arbitrary point in an object originally at is moved to when the object is under load. The displacements are often collected in the displacement vector .
If all points in a body experience the same displacement, there is no deformation. In other words, a body only deforms when there is nonuniform displacement. If there is only displacement but no deformation we speak of a rigid body motion of the object. A rotation or translation of a body are rigid body motions.
Strain is a quantity that describes the amount of deformation or distortion within a body  and is a ratio and unit-less.
A body that undergoes a rigid body motion like rotation or translation does not experience deformation; the name rigid body motion expresses this. A rigid body motion is pure displacement. Strain, however, is a consequence of a change in size or a change in shape. A good strain measure captures this requirement: strain should be zero for rigid body movements since it should only measures deformation.
Strain is not a physical quantity like temperature and there are various strain measures. Almost all strain measures in use today [13, c. 4.1, p. 90] are constructed in such a way that they give more or less the same results when the deformation is small. For small deformations it does not matter which strain measure is used.
A strain can be non-zero even when the displacement is zero. Consider a rod that is fixed at the left and pulled on the right in the positive -direction. At the fixed end the displacement will be zero, yet the strain is non-zero.
Infinitesimal strain and engineering strain
The above definition of strain is not what is actually used in the Wolfram Language but is useful for building an intuition. For the derivation of the definition that is actually in use we follow [4, p.475]. To describe the deformation of a body we consider a point in an original configuration and that some point in a final, deformed, configuration. We consider an infinitesimal cube where we characterize the change in size and shape. The change in size is determined by the changes of length of the cube. The change in shape is determined by the change in angle from the originally between the sides.
Let's say point is at . Then point is at since it is parallel to the x axis. The point is a displaced by such that point is at . The point is displaced by such that its position in the final configuration is . From these coordinates the original and deformed length of segments and can be computed as follows:
So far this is exact. Now comes a crucial simplification that specifically leads to the infinitesimal strain measure. It is assumed that the displacements are small and that they are much smaller compared to 1. This means that
What remains to be done is to quantify the changes in the angles which are initially at . This will lead to an expression for the shear strains. The shear strains quantify the change in angle. Consider a line segment to line segment in the original configuration. In the final configuration these are displaced to and , respectively.
The point originally at is displaced to . Since we already have made the assumption that displacements are small and since we only want to describe the change in angle and thus ignore the change in length the point moves up by relative to point . Point moves to the right by relative to . Again, changes in length are not considered when we want to describe the change in angle. With the
This, however, would have some disadvantages. For one the object above does not behave like a tensor and rules out other succinct ways of representing the object. Engineers measured shear strains long before tensors were invented. To remedy this situation the tensorial shear strain are defined as engineering shear strains
This is what is called the infinitesimal strain measure. It is very important to realize that this strain definition can only be used if the deformation and rotations are small. If that is not the case then other strain measures have to be used.
It is possible to switch off the usage of engineering strains by setting "EngineeringStrain"->False in the parameters . This may be of interest when comparing linear material laws with nonlinear material laws, where true stain is used. More information can be found in the reference page of SolidMechanicsStrain.
The infinitesimal strain measure is the default strain measure that the Wolfram Language uses, as the default material model is a linear elastic material model. The choice was made as small deformations are the most common scenario that is modeled. Furthermore, the infinitesimal strain measure is linear and does not per-se result in a system of nonlinear equations which take longer to solve. The infinitesimal strain measure is useful for modeling small deformations of concrete, stiff plastics, metals, linear viscoelastic materials such as polymeric materials, porous media such as soils and clays at moderate loads; in fact almost any material can be modeled with the infinitesimal strain measure if the load is not too high. The infinitesimal strain measure is inadequate for rubber materials, soft tissue or large deformations in general [13, c. 4.1, p. 95].
The infinitesimal strain measure has, however, limitations that one needs to be are of. The infinitesimal strain does not produce a zero strain for a rigid body rotation. This is best illustrated by an example.
Note that the normal strains are not zero. The infinitesimal strain measure is only valid for small rotations. If rotations are large other strain measures, such as the Green-Lagrange measure, are a better choice.
Unfortunately, it is not as simple as to just replace the infinitesimal strain with a, say, a Green-Lagrange strain as the Green-Lagrange strain is not compatible with the Cauchy stress. This will be explained in much more detail in the section on hyperelasticity.
Before we consider strain measures for large deformations the concept of a deformation gradient is established. For small deformations the difference between the object before and after the deformation is small compared to the size of the object and is neglected. In the case of large deformations that is no longer the case and the deformation needs to be accounted for.
We begin by distinguishing between the undeformed object and the deformed object. In this section the displacement vector will not be shown in bold to make the notation more consistent with the notation used in literature.
The undeformed object is placed in a coordinate system with basis vectors . This domain is referred to as the reference configuration. We follow the convection that capital letters are used to refer to entities from this domain. When this object undergoes deformation every material point is displaced to a material point the deformed object. We use lower case letters to denote all entities from the deformed domain. To keep things simple both the object in the reference configuration and the object in deformed configuration make use of the same coordinate system . The reference configuration is sometimes also called the material configuration or the initial configuration or the Lagrangian description. While the deformed configuration is sometimes called the spatial configuration or the current or final configuration or the Eulerian description.
The deformation gradient matrix is the Jacobian matrix of the deformation map . The deformation gradient tensor allows to describe the relative position of two neighboring particles in the deformed configuration in terms of their relative particles position in the reference configuration [17, p. 81].
The Green-Lagrange strain measure is a nonlinear strain measure that does not suffer from the small displacement and rotation limitations the infinitesimal strain measure is limited by. Conceptually the Green-Lagrange strain measure models
To show that the Green-Lagrange strain measure does not suffer from the small deformation limit we consider the same example as in the infinitesimal strain section but make use of the Green-Lagrange strain measure.
The Green-Lagrange strain measure is a material strain tensor that describes the strain in the reference configuration. This is in contrast to, for example, the Euler-Almansi strain measure that describes the strain in the deformed configuration.
For completeness the Euler-Almansi strain is also given. The Euler-Almansi strain measure is a spatial strain tensor that describes the strain in the deformed configuration. Conceptually the Euler-Almansi strain models:
An isotropic material is a material that behave the same in all directions. We assume that the -direction is the longitudinal direction and the - and -direction are the lateral directions. In the elastic region we can say that:
For a isotropic linear elastic material Poisson's ration is a value between . For most metals the Poisson's ratio is about 1/3. Many materials have a Poisson's ration of 0.2-0.3. On the extreme ends, Cork, for example has a Poisson's ratio of 0. Rubber, for example, can have a Poisson's ratio of close to . Poisson ration values of mean that the material is incompressible and pose a problem for numerical simulation. Artificial material, such as auxcetic material, can have a negative Poisson's ratio. If an auxcetic material is stretched in one direction it becomes thicker in the other direction.
Stress is a quantity the describes the distribution of internal forced within a body  and is measured in or
Finding the stresses in an object is an important task as it allows to predict when the object will fail. Say a material can withstand a maximum stress of . A material will fail if . The maximum force that can be applied before the material fails is
For more information see the section Failure Theory.
The relation between stress and strain is a material property and can be visualized in a stress-strain diagram. The data for these diagrams come from tensile tests. For these tests a material specimen is clamped into a machine that has two clamps and pulls the specimen apart. The force applied and the strain produced are recorded until a fracture occurs.
The stress-strain curve relates the measured strain to the force and hence stress applied. An archetypical stress strain curve has several sections worth discussing. The following illustrates an exaggerated stress-strain curve for a ductile material. Several points of interest have been marked. Especially the points A, B and C may be very close together, or the connection between A and C may be a straight line. What is happening physically at these points of the stress-strain curve is that the material slips along internal crystal boundaries. This happens very fast and is hard to detect.
- A: The proportionality limit. Up until this point the stress-strain curve is linear. Hooke's law applies here.
- B: The elastic limit. This point is beyond the linear stress-strain relation and marks the end of the nonlinear elastic region. This means that a load that is applied up until this point will not cause a permanent deformation. Up until this point the material is still fully elastic and after removal of the load the material will return to its original form. Beyond this point we have a permanent plastic deformation.
- C: The yield point. Beyond this point the strain will increase rapidly. If the yield point of a material is unknown, it is typically estimated with the 0.2% offset method. For this, a line parallel to the linear part of the curve is drawn that starts at strain of 0.002 (or 0.2%). The point where that parallel line intersects the stress-strain curve is taken as the yield point. Some industries may have a lower offset value. The elastic limit and the yield point are typically very close.
- D: The ultimate strength is the maximum stress value of a material. The ultimate strength is also called the tensile strength. Beyond this point the material will exhibit a phenomena called necking. Here the cross-sectional area of the test specimen will start become thinner. Much like the neck of a droplet before breaking off.
The initial section is the linear section. Here the relation between stress and strain is linear and known as Hooke's law. A load applied in this region is fully reversible. This region is also called the elastic region. In this region we have
where is Young's modulus and a material property. Young's modulus is also called modulus of elasticity. Young's modulus has the same units as stress . The material property Young's modulus is measured with a uniaxial tension test.
More in depth information and material data for the the stress-strain relation can be found in .
Stress strain curves are obtained from tensile tests. In these tests a specimen is pulled apart by a force, which applies a stress and the strain is recorded. During the test the specimen will deform and as a consequence the cross sectional area will change. The cross sectional area, however, enters the equation for the stress as:
It is difficult to measure the instantaneous cross sectional area during testing. So only the initial form of the specimen is recorded. The true stress and strain, however, take the change of form into account.
The difference comes during measurement of stresses in a test specimen. Typically only the original area of the specimen is considered. The same holds true for the strains. The true strain is given as
Next, we assume that the volume of the specimen stays constant. This assumption is valid in the elastic region because volume changes in the elastic region will be small. The assumption is also valid in the first part of the plastic region because materials are considered incompressible during this part of plastic deformation. The assumption, however, is no longer valid in the second part of the plastic region after the process of necking has started.
Also note that for small strain values, when we are in the elastic region the difference between the true and engineering stresses and strains is small compared to the difference in the plastic region. In cases where there is significant plastic deformation the true stress strain curve should be used.
Using the data from above section Stress-Strain relation we convert the engineering stress and strain into true stress and strain.
The constitutive equation is also called the material model. The constitutive equation describes how stress and strain are related. In the linear case stress and strain are related by a generalization of Hooke's law through the elasticity matrix :
The number of the elastic constants can be reduced. This simplification is achieved by making several assumptions [4, p.478]. It is almost universally assumed that the elasticity matrix is symmetric. This assumption reduces the 36 elastic constants to 21.
Further reductions are possible if the material exhibits symmetries about some planes. The exact form of the elasticity matrix depends on these symmetries and materials are thus classified as isotropic, orthotropic or anisotropic. In the simplest case the material properties are the same in all directions, this is the isotropic case. Most metals are an example of an isotropic material.
If the material is only symmetric about three orthogonal planes then we have an orthotropic material. Wood is considered an orthotropic material where the material properties depend on the direction of the wood grain.
The most general case are anisotropic material. Anisotropic material are typically compound material or biological tissue where the material properties vary in all directions. Also fibre reenforced material where the weave effects the material properties in all directions are considered anisotropic.
A linear elastic material model is good for the vast majority engineering design calculations, where components cannot exceed yield [11, c.1.1.5].
where and transforms the elasticity matrix into its full tensor form . The full tensor form is what is seen from the output of SolidMechanicsPDEComponent.
This is explained in more detail in the section The output of SolidMechanicsPDEComponent.
The isotropic linear elastic material model is the default material model used in the Wolfram language. We speak of a material as being isotropic if the material the behaves the same in all directions. In other words their material properties are direction independent. An isotropic linear elastic material model is good for polycrystalline metals, ceramics, glass and polymers undergoing small deformations and low loads [11, c.1].
Alternative Material Parameters
Any elastic modulus can be expressed in terms of any other two moduli. With this is possible to use any of the common moduli and SolidMechanicsPDEComponent will find the moduli it needs for its operation.
Orthotropic material have different properties in the three dimensions. Wood is an example that can be modeled as an orthotropic material. Consider a wood log. The wood then is most stiff along the grain, somewhat stiff in the circumferential direction and least stiff in the radial direction.
For anisotropic material the material properties are different in all of the material directions. The anisotropic linear elastic material model is good for reinforced composites, wood, single crystals of metals and ceramics [11, c.1].
If only the inverse of the elasticity matrix - the compliance matrix - is available then that can be specified with the parameter "ComplianceMatrix" and the inversion to get the elasticity matrix will be done automatically. Specifying an elasticity matrix will overwrite the compliance matrix if specified.
The linear elastic equations can be generalized in several ways. Initial strains or thermal induced strains can be considered or initial stresses can be modeled. Adding these, the constitutive equation changes in the following way:
A change in temperature of a body has a change of volume, and thus length, as a response. This temperature induced change of length can be modeled by an change of strain, the thermal strain . The thermal strain is then subtracted from the strain such that then is the temperature free strain. The default model for the thermal strain is given as the linear relation between temperature and strain
where is the constant coefficient of linear thermal expansion (CLTE), measured in , is the temperature of the body in and is the reference temperature, also in at which there is no thermal strain in the body.
To estimate if a thermal strain is relevant for a simulation [11, c.1] one can use
Here and are the maximal and minimal expected product of in the object. If is a significant fraction of the stress in the object then a thermal strain should be considered. Similar estimates can be made to decide if a transient heat conduction analysis is needed or not [11, c.1]. A solid with characteristic length will reach steady state temperature in time
One has to be aware that itself is temperature dependent. The default model is suitable as long as and the linear coefficient of thermal expansion does not vary too much over a temperature range. The exact specifics of what "does not vary too much over a temperature range" means depends on the accuracy needed in a specific application. If a high accuracy is needed then a temperature dependent will be needed and doing so is explained further down. To keep things simple in the first example, a setup is chosen where the coefficient of thermal expansion does not depend on temperature. No thermal analysis is needed in this case.
As an example a cylinder of length along the -axis and radius is constraint at both ends. When the cylinder is thermally loaded it cannot move in the -direction and only expand in the - and -directions. For this setup an analytical solution for the reaction force acting on the end caps is available and will be used to verify the result.
One thing to keep in mind is that if only a temperature difference is to be specified , then one should use "ThermalStrainTemperature" for setting that up since the "ThermalStrainReferenceTemperature" is set to zero when not given.
Fixing the movement in the -direction is not sufficient to fully constrain the cylinder as that only constrains the displacement represented by dependent variable . To also constrain the displacement dependent variables and the points and will be fixed such that no movement at all is possible.
where, is Young's modulus, the coefficient of linear thermal expansion, is the temperature of the body and is the reference temperature at which there is no thermal strain in the body. is the relevant surface area. Because the geometry is discretized and the verification of the result should not consider the discretization error the area of the discretized surface is used.
Coupled heat transfer equation
In the previous example the temperature distribution was constant throughout the domain. The next example considers a non constant temperature field. The temperature field is coupled to the solid mechanics PDE model by specifying the "ThermalStrainTemperature" either as an interpolating function from a previous thermal simulation or by generating a fully coupled solid mechanics thermal model. Making use of an interpolating function is the preferred approach if the deformed body does not influence the temperature distribution and the coupling is only in the direction from the thermal field to the solid mechanics. Using an InterpolatingFunction object as a "ThermalStrainTemperature" source also has the advantage that the maximal memory requirement to solve the fields sequentially will be less than a fully coupled PDE. The fully coupled PDE is more general, though and will be shown here, although not strictly necessary.
The following example models a bimetallic strip used as a thermocouple [9, p. 296]. The bimetallic strips consist of two metals attached on top of each other. The metals have a different coefficient of linear thermal expansion (CLTE). At the left hand side the bimetallic strip is held at 100 . The top and bottom of the strip is exposed to HeatTransferValue. At the left the structure is fixed to the wall. Due to the different CLTE the bimetallic strip is expected to bend. The maximal deflection is sought.
More information on setting up multi material region can be found in the section Multiple materials.
Note, that the "ThermalStrainTemperature" is now the dependent thermal variable . The value of is computed by the coupled heat transfer model. For more information about heat transfer modeling please consult the Heat Transfer monograph. The "ThermalStrainReferenceTemperature" is the temperature at which no thermal strains are expected.
Temperature dependent coefficient of thermal expansion
Up until now the coefficient of thermal expansion was linear and constant (CLTE). The linear approach is applicable only in a material dependent temperature range. If a larger range or a more accurate result is sought then assuming a linear coefficient is not sufficient. Typically a material supplier will have collected this information in a data sheet. At this point a different scheme must be used to calculate the thermal strain at a given temperature. Several choices are available to do so:
Using the thermal strain function itself can be done by specifying an "InitialStrain" in model parameters section. Setting the "InitialStrain" is explained in the section Initial Strain and is not further explained here.
There are two types of CTE (coefficient of thermal expansion) definitions, the tangent coefficient and the secant coefficient. The PDE models in the Wolfram Language use the secant coefficient of thermal expansion. This section discusses different definitions of the CTE, how to convert between them and how to make use of them in solid mechanics PDE models.
Initial stresses or residual stresses, as they are also called, can come up when, for example, during the manufacturing process the object under consideration was exposed to large strains or high temperatures. The body can be in a pre stressed state which can be modeled by specifying this initial stress [10, p. 77].
In this section we concern ourself with simplifications of solid mechanics PDE models that result in two space dimensional models. By this we mean that we deal with idealized two dimensional solids where both the region dimension and the embedding dimension is two. This is in contrast to the three dimensional solids where both the region and embedding dimension was three, or shell elements where the region dimension is three but the embedding dimension is two.
The main reason to use 2D simplifications is time. Typically 2D models are faster to set up and need less resources to solve. However, 2D models are simplifications and cannot provide all the details a full 3D model provides.
In this section we will introduce plane strain and plane stress PDE models. These are simplification of the full 3D solid mechanics PDE model that can be used in 2D. Caution: even though these are 2D models they still have components in the third dimension. 2D modeling forces you to make assumptions on the boundaries of the 3D object. It is important to understand where these models come from and what their limits are. Typical cases for plane strain models are cross sections of very thick objects and typical cases for plane stress models are thin sheets.
A plane strain situation can be exploited if the movement of a 3D object is constrained in the -direction, like a rod clamped between a wall, and all forces act in the -plane. Note that this does not mean that the -direction end points of the 3D object are also constrained in the -plane; that is not the case they just should be constrained in the -direction. Should that not be the case and the -direction end points are also fully constrained in the -plane, then, and only then, the object has to be very long in then -direction the compared to the -plane expansion of the object for the plane strain model to be valid. Analysis of dams and shafts fall into this category.
Starting from the constitutive law, relating stress and strain by . We insert the assumption that there is no displacement in the -direction in the strain measure . A direct consequence is that the strains in the -direction are 0 such that
A very important point to realize here is that even though the strains in the -direction are set to zero the normal stress in the -direction is not 0. In the current version SolidMechanicsStress does not recover the -direction stress in the plane strain case. This needs to be done manually.
The following example demonstrates the use of the "PlaneStrain" model form. The model domain is a quarter cross section through a pipe. With matching units of an inner radius an outer radius and a thickness . At the left boundary we have a symmetry constraint such that the pipe can move up and down and at the right bottom we have a second symmetry constraint such that the pipe can move left and right. A pressure of is acting within the pipe in an outward direction. The remaining boundaries are free to move. Young's modulus is given as and Poisson's ratio is .
The derivation of the plane stress case is similar to the plane strain case but models very different cases, explained in a second. In the plane stress case we assume that all -direction stresses are 0. This is an approximation and the plane stress model becomes more accurate the closer the thickness of the object approaches zero. Setting all -direction stresses to 0 is equivalent to a 0 boundary load on these surfaces and the closer these boundaries are the more accurate the condition is inside the plane.
Note that the normal strain in the -direction is not zero. In the current version SolidMechanicsStrain does not recover the -direction strain in the plane stress case. This needs to be done manually.
Here are a few points that indicate the limitations of linear elasticity [10, p. 78].
- Although strains remain small the linearity of material behaviour is not given. Some materials (polymers) may undergo shape changes that are not linear but reversible.
For more information see the section Failure Theory.
This section explains the relation of the output of SolidMechanicsPDEComponent and the equilibrium equation. The equilibrium equation for the stationary case when there is no external load is is given as
When we look at the SolidMechanicsPDEComponent operator it is given by
However, if we look at the output of SolidMechanicsPDEComponent that is more like
The default output equation deviates from the equilibrium form given in the equations. The question is, why is that? Before we answer that question let us assume we want to use or own strain and stress functions to overwrite the system functions. This can be done with specifying "StrainFunction" and "StressFunction". When these functions are specified SolidMechanicsPDEComponent assumes that these function are nonlinear. For simplicity we just use SolidMechanicsStrain and SolidMechanicsStress.
Now, we have the Div form expected. When the solver sees this Div operator will need to match to either a DiffusionPDETerm, a ConvectionPDETerm or to the DerivativePDETerm, simply because these are the only terms available (and the only terms needed). Because we do not have a function of but of the derivatives of the only choice is to parse it as a DerivativePDETerm. This works well and is the most general form possible. However, this has a disadvantage as it will trigger the nonlinear solver. We solve for not for the derivatives of . That in it self is not a problem, but the nonlinear solver will have to make at least two evaluations of the function to converge so this general approach is a bit slower than it could be. After all a linear equation should be solvable in one step. How do we get to a linear formulation when the equations are linear? We do this by rewriting the equations into an equivalent form that is making use of the DiffusionPDETerm and uses that results in a form like
A special class of nonlinear elastic material models with small deformations are called hypoelastic material models and are not to be confused with hyperelastic material models where deformations are large. A hypoelastic material model is an elastic material model with small deformations, but the constitutive equation is no longer linear. The hypoelastic material model presented here is based on [11, c.3.3 and c.8.3]. Hypoelastic material models are not commonly used because they are not an accurate model of any real material. Yet, there are some reasons why to include them:
- Hypoelastic material models sometime are used to model measured stress strain data that is beyond the limits of the linear stress strain relationship.
- Hypoelastic material models sometime are used to model material in the plastic regime while not using a full plastic model. A necessary prerequisite for doing so is that the object can not be unloaded. This means once a load has been applied to an object it can not be removed because a true plastic deformation would behave differently. In this case the hypoelastic material model is used to mimic the loading of a plastic deformation.
- A hypoelastic material model is included here because it makes a good example for how to write a custom material model, in this case with a nonlinear stress strain relationship.
A hypoelastic material shows a nonlinear stress strain relationship even at small strains but at the same time is fully reversible. Fully reversible means that if an applied load is removed the object returns exactly to its initial state. The material is fully elastic and there is no plastic deformation involved. For this model strains and rotations are assumed to be small and the infinitesimal strain measure is used. As a stress model the Cauchy stress is used.
The concept of finding the derivative of an energy density function to express the stress is explored in more detail in the section about hyperelastic materials.
The models specific strain energy density function is not given in  but the hypoelastic constitutive equation is given as:
We write a function that specifies this material model. The arguments to the function are the variables vars, the parameters pars and data data, that contains data such as the default strain measure. The material model function given below contains some model parameters hard coded into the function to not obscure the essence of the model. These parameters could very well be placed in the parameters pars and parsed within the function.
Sometimes, the nonlinearity may be so strong that it is not possible to get to the solution directly. In these cases the nonlinear solver will issue a message that it can not converge. Then not all hope is lost, one can try to get to the final solution in steps. This is done by starting with a much smaller than the final pressure that we want to achieve. In a next step we slightly increase the pressure but at the same time use the solution of the last, lower pressure result as an InitialSeeding for the nonlinear solver to start over. We will see an example of that in a minute.
To verify that the behaviour of our model is indeed a model of a nonlinear stress-strain relation we construct a plot of force versus the displacement. To create the plot we compute how much displacement we have at the point a in the a direction for a given pressure. We use pressure as a proxy for force.
What remains to be done is find the model parameters , and . For this we load an example of a measured stress-strain curve and create fits to the data to estimate the model parameters. and will be estimated from the linear part of the measured stress-strain data and will be fitted to the hypoelastic model.
Now, we extract the linear part of the stress-strain curve. For this we truncate the data. It might not be obvious where the linear part of the stress-strain curve ends and where to stop; this is in the discretion of the analyst. The noise of the measurement can make this difficult.
The point of this is to find model parameters that create a good fit. Found by trial an error, we evaluate the fitted model just beyond the linear section of the curve. In this case this gives a good overall fit to the measured stress strain curve.
We know have a all model parameters to adjust the model to our data. The last step is to verify that these parameters create a model of the measurement data. For this we evaluate the model at various steps and visually compare the evaluated model to the stress-strain curve measurement data.
Materials like rubber or foam can be exposed to large deformations and still remain fully elastic. This means that if the load is removed, the deformation is fully reversible with no plastic deformation. These are called hyperelastic materials. Beyond rubber and foam, some biological tissue or polymers, which can have rubbery regimes, can fall into the hyperelastic material category. In contrast to hypoelastic materials, hyperelastic materials can be subjected to large deformations and rotations.
For more information see the notebook on hyperelasticity and hyperelastic material laws.
Different failure theories are available for different material types. For example for ductile material there are the von Mises and Tresca failure theory while for brittle materials there are the Coulomb-Mohr and Modified Mohr theories, to name a few.
Most failure theories compare the principal stresses to material properties. For ductile materials the principal stress is compare to the yield strength. For brittle materials the principal stress is compared to the ultimate strength or the point of fracture. The principal stresses are computed to make the stress values comparable to the material data like yield and ultimate strength which come from uni-axial tensile tests. The section Boundary Load: Tension has an example that shows the use of PrincipalEigenvalue to find stresses in a non-axial loaded cylinder.
Failure theories need to account for different observed behaviour. For example, brittle material tend to have larger compressive strength than tensile strength. Also, hydrostatic stresses do not cause yielding in ductile materials. Hydrostatic stresses are stresses that cause a change in volume as opposed to deviatoric stresses that cause a change in shape.
Next, we discuss different failure theories .
The Rankine failure theory considers an object failed if the absolute value of the maximum or minimum value of the principal stresses reach the yield strength for ductile material or the ultimate strength for brittle materials. The Rankine failure theory is also called the maximal principal stress theory. The Rankine theory is a simple theory but not a particularly accurate theory, especially for ductile materials. For ductile materials it does not respect the observation that hydrostatic stresses do not cause yielding in ductile materials. Hence, if the Rankine failure theory is used at all then it is used for brittle materials.
Yielding occurs when the maximum shear stress is equal to the shear stress at yielding in an uniaxial tensile test .
A safety factor below one is problematic. Typically objects are engineered so that the FOS always stays well above 1 in all usage scenarios. Finding stresses that are higher than the yield stress is something to be cautious about.
The Tresca theory is easier to apply and more conservative than the von Mises failure theory and the maximum difference between the two theories can be calculated to be 15.5%. The von Mises theory agrees better with experimental data.
See also the safety factor section in the overview section for an example.
The Tresca theory is easier to apply and more conservative than the von Mises failure theory and the maximum difference between the two theories can be calculated to be 15.5%. The von Mises theory agrees better with experimental data.
A composite is a body made up of multiple materials. To illustrate the set up of a multi material region a simple two material bar is subjected to a surface force while at the same time being constrained at both ends.
The internal material boundary has been maintained during the creation of the full mesh. To better illustrate the affect of a force acting on the two materials, two extreme elements are chosen: On the left the bar is made from titanium and right hand part is made of lead.
The time dependent solid mechanics models we have seen to far are undamped. This means vibrating objects will continue the vibrate indefinitely. This is not what we experience in real life. For example the amplitude of a struck tuning fork will decay over time; the amplitude will dampen out. With damping, energy dissipates over time. There are various models for damping; the model we are considering here is viscous damping where the vibrating object is dissipating energy because it is surrounded by a viscous medium that counteracts the vibration.
Modeling damping requires time dependent PDE models. For an introduction to time integration see the section Time dependent analysis.
Rayleigh damping introduces a method to construct the damping matrix . The term represents viscous damping and effects the velocity of the solid. Rayleigh damping assumes that is a linear combination of the mass matrix and the stiffness matrix
The mass damping parameter has units of and the stiffness damping parameter has units of . It should be noted that this linear combination assumption does not have sufficient physical basis, but works reasonably well to be widely used. Because Rayleigh damping is a phenomenological model and not a physical model, identifying the values for the parameters and can be difficult task. We will discuss some aspects of it below.
The general form of the equilibrium equation with the Rayleigh damping parameter is given by:
As an introductory example we set up a rectangular region of length meters, height and thickness with a plane stress model form. The beam is fixed at the left end and at the right hand side an downward load is acting. In this example the downward force is a constant, time independent value. This essentially means that the load is applied as a unit step function from the beginning of the time integration. We make an undamped and a damped model. The damped model has a mass and stiffness parameter defined. The values for these either come from the literature or are computed for resonance frequencies, which is discussed further below.
Next, we set up an auxiliary function that will take a set of model parameters and compute the displacement. This makes sure the same setup is used for both cases we want to simulate. The auxiliary function includes the PDE model, the boundary and initial conditions, the region and the time range.
With the damping parameters specified the damped PDE model shows a decay of the amplitude of the displacement of the query point while the undamped model does not show a minimal decay. This minimal decay is due to numerical error and can be reduced by setting a smaller MaxStepSize. Both oscillate around the value of the stationary solution.
When we have stiffness damping with . Stiffness damping is a linear relation. Looking at the plot we see that stiffness damping, and thus the stiffness damping parameter dampens out high frequencies. For we have mass damping with and we see that the mass damping parameter dampens low frequencies, mass damping is like 'underwater' damping. Rayleigh damping is the linear combination of mass and stiffness damping.
When we manipulate the previous equation we get an alternative to directly defining and . The alternative is to specifying two resonance frequencies and in and the two corresponding damping ratios and which gives the following for the damping parameters:
For most engineering metals values of < 5% can be assumed for the damping ratio, many are even below 2%. For elastomeric materials like rubber values will be larger. For now assume we have values for , then the next question is how to obtain the frequencies and . For this we perform an eigenvalue analysis of the undamped system, as shown in the section Eigenmode analysis. This will give us the eigenfrequencies and their corresponding modes. We then choose the two lowest frequency modes that relate most to the motion we wish to dampen out. The related eigenfrequencies are then the frequencies and .
Since we want to dampen the downward motion we select the first two modes. The third mode is not relevant, as that is in the -direction. The higher modes are left out. The idea is that you choose the lowest, most relevant modes.
The following example demonstrates the procedure. For this a single specific mode of the system is excited. The amplitude of the displacement versus time diagram, which is sometimes called the system ring-down, is recorded. To illustrate the procedure, we first generate a fictitious ring-down data set. Typically this would come from an experiment, however.
To verify that the computed damping ratio is correct, we plot the envelope of the simulated ring-down. Note, that this can only be done because in our generation of the fictitious experimental data we also specified the natural frequency .
Specifying lower will result in larger amplitudes and larger stresses, so specifying lower damping ratio is more conservative. If in doubt a lower damping ratio should be specified. The damping ratio is a function of frequency. This and the fact that we only using two modes results in the disadvantage of Rayleigh damping: it rarely matches the necessary damping over a large frequency range. The tendency is that for too large a frequency range, there is too little damping in the mid frequency domain and too much damping in the low and high frequency range.
Boundary conditions for solid mechanics applications fall into one of two categories. One are boundary constraints and the other are boundary loads. The constraints restrict or prescribe the displacement of an object in one or several directions. An example is an object that is fixed to a wall. At the fixation points no displacement of the object is possible. These type of boundary constraints are realised by Dirichlet conditions. Thus their names include the the term condition. Typically, at least one condition type boundary condition must be specified to make the differential equation solvable. On the other hand we have boundary loads also called traction. These are force, pressures or moments acting on surfaces. An example is an external force, like the weight of a book on a bookshelf, acting on a surface. These boundary loads are realised with Neumann value boundary conditions and their names include the term value. If no boundary condition is specified on some part of the boundary then there is no traction on that part and the object is free to move.
For boundary loads it is typically more convenient to make use of ElementMarker as a predicate for boundary conditions. For constraint conditions either geometric predicates should be used or a mesh generated with PointMarkers → BoundaryDeduced in conjunction with SelectPointMarkerFromBoundaryMarker should be used.
None can be used to exclude a direction from a prescribed displacement. This can be used, for example, to prescribe a displacement in the negative -direction while at the same time the object is then still free to move in the other directions.
A roller constraint is used to constrain the displacement of an object normal to the face the constraint is applied to. This means the object cannot move normal to the face but the face is free to move in other directions. This is as if the face was resting on a roller. Note that the face is also restricted in the negative normal direction. In other words the face cannot move along the normal direction, but can move in tangential directions perpendicular to the normal.
Here , and are the dependent variables and is the normal in the various directions. Currently, at the top level, roller constraints are only supported when the normal is in an axial direction and thus two of the normal components are zero. In that case the roller condition can be realized with a DirichletCondition. For the more general approach the low level finite element functions need to be used.
The following graphics illustrates the workings of a roller constraint. The beam is fixed at the lower left edge, shown in black. On the top a downward facing pressure is active, indicated by the blue arrows. In the area highlighted in red, a roller constraint is active. In this area the beam cannot move normal to the red surface.
Not applying any DirichletCondition type boundary conditions and solving the system of equations will lead NDSolve to emit a warning message that no DirichletCondition or Robin type NeumannValue have been specified. For good reason: Not specifying a DirichletCondition will make the system of equations singular and the solution can only be found, at the best, up to a constant. This is behaviour is independent of the solid mechanics application but generally true and is shown in the reference page of DirichletCondition and NeumannValue. How many constraints need to be specified to sufficiently constrain an object? The goal is to eliminate rigid body movement. The reason is that it should not matter to the overall performance of a body and it’s loads if the body and the loads are translated or rotated. In 3D any body has 6 degrees of freedom. It can translate in 3 directions and rotate about all 3 axis. In 2D the body can translate in two directions and rotate in the plane. So in 2D we have 3 degrees of freedom.
To create a fully constrained PDE model that is solvable at least these degrees of freedom need to be constrained. The continuum elements use in the Wolfram Language doe not provide rotational degrees of freedom and consequently those can not directly be set. In this case the choice of restricting translational degrees of freedom need to be such that a rotational movement is restricted.
Above, the node at is constrained in the - and -direction, which corresponds to a 0 DirichletCondition in the - and -displacements. For the rotational restriction another constraint is added at the lower right. Here we have a single constraint in the -direction. This means the body is free to move in the -direction. Note, that this second constraint prohibits movement normal to the direction between the two constraints.
In 3D an analogous procedure us used. Here we add a constraint that restricts -, - and -translation, a second constraint that restricts two directions, say, the - and -direction and a third constrain that restricts a single direction, say, the -direction.
A load is either a pressure or a force acting on a surface. Both can be specified in SolidBoundaryLoadValue. In the case of a force acting on a surface that force will be automatically converted into a pressure by taking the surface area the force is acting on into account. In other words a surface force will automatically be converted in a surface pressure.
Pressure or force have always to be given in units of force per area; also in two dimensions. In the two dimensional case the pressure is internally multiplied with the value of the "Thickness" model parameter to get to a force per length unit. In the case of a length unit of meters this then converts the pressure unit to a unit of .
The force components are divided by the area on which the boundary load is active. Any of the pressure or force components can be 0. There is, however, no None specification like there is for SolidDisplacementCondition.
To illustrate various concepts of boundary loads some standard loading types will be shown. The standard loading types are shown in the below illustration and will be explained in the subsequent sections.
Axial direction compression
Non-axial direction compression
A more general approach is to make use of the BoundaryUnitNormal expression and extract the normal component needed. Computing the BoundaryUnitNormal is computationally more expensive than specifying the direction vector explicitly. In cases where the direction vector is readily available it is thus beneficial to make use of it. In other cases, for example, when the geometry is complicated using the BoundaryUnitNormal is the better option.
See also this verification example.
As a quick test we can also proceed to compute the strains and stresses and using that the force over the area is the stress in in -direction: . We specified the surface pressure which should be recoverable.
This stress concentration is referred to as Poisson's effect. It can be prevented by changing the boundary conditions that fix the cylinder on the left. This will be shown just now. Before that, we look at the volume where the stress in the -direction is larger than 1% of the given boundary pressure.
As mentioned above the stress concentration at the left boundary is induced by the boundary condition. This is expected and acceptable. There is a slightly different way to set the boundary condition up; that however, is also a slightly different way to model a wall boundary. This is done by adding a constraint for the -direction and leave all other directions free to move, except at one point where the and -direction are also constrained.
Note that the over all stress distribution is now close to the surface pressure specified. The remaining disparity is numerical noise. This approach works well if there is a symmetry in the object that can be exploited, like here where we put the fixed condition at . In nonsymmetric geometries it may not always be obvious where to place that condition.
In the previous example we had a cylinder that was aligned to the axis and the boundary load was also aligned to the same axis. Because of that the stress component reflected the surface pressure. If an object is not loaded parallel to an axis the stress (or strain) components will also not be axes aligned.
Now, the normal stress is no longer the same as the surface load. This is expected as the load is no longer axis aligned. It is worthwhile mentioning that the von Mises stress still gives the surface pressure because the von Mises stress combines the various normal and shear stresses.
It is still possible to find the maximal normal stress as if the cylinder and load had been axis aligned. This is done by computing the PrincipalEigensystem as indicated in the section Principal stress values and direction vectors.
See also this verification example.
Ultimately all boundary loads need to be converted to pressures acting on surfaces. Converting a boundary force to a pressure is a matter of computing the area the surface force is acting on and then dividing the force by the computed area. A similar process needs to be done for torsion.
where is the shear stress (a pressure), the radius and the second moment of area in . After rearranging we get
The shaft appears to be growing radially at the right hand end. This is physically incorrect and caused by two issues. One issue is that the view is exaggerated by the mesh deformation plot. The element mesh deformation plot computes a fixed amount of enlargement to better visualize the effect of a deformation. So while there is a small spurious radial growth in the solution it is further exaggerated by the deformation plot. The fact that there is a small radial deformation at all is due to the linear theory used in this case.
The linear elastic model and its associated strain measure make a small strain and small rotation assumption. Care has to be taken that the shear angle remains small. The rotational angle should be small such that the trigonometric simplifications and hold. In cases where this approximation does not hold large deformations have to be taken into account. Consider the graphics below.
The shear angle is small enough for the small strain and small rotation assumption to be valid. Even for the small shear angle there is a slight growth in the radial direction, which is picked up and exaggerated by the deformation plot.
As a side note, hollow beams are more efficient to carry a torsional load because the central part of the beam only resists a small part of the torsional load so removing that material is not having much effect on the performance of a beam under a torsion load.
See also this verification example.
In this section we will pursue a somewhat academic exercise. Let's ignore the yield strength and apply an even larger moment that will result in a large rotation. We know from the section on strain measures that infinitesimal strain measure is not well suited for large rotations.
Because we use a linear elastic material model the default strain measure is the infinitesimal strain measure. We can recompute the same PDE model but replace the infinitesimal strain with the nonlinear Green-Lagrange strain, that does require the small rotation limit.
Here we see that element mesh deformation actually indicates a compression. The element mesh deformation plot tries to emphasize deformations such that they are visible. In this case we know that effect shown by the deformation plot is purely from the visualization and not from the infinitesimal strain small rotation limitation.
Since the deformation is small we can plot the displacement with a scaling factor of 1 that can be specified with "ScalingFactor"1 or "ScalingFactor"None.
We would like to emphasize that the amount of torsion is well beyond the yield strength of the material. The large torque was used to emphasize the difference in the infinitesimal strain versus the Green-Lagrange strain.
The section on Roller constraints has an example of this kind of boundary load.
Next we study a uniform compressive load on a spherical shell such as that on a spherical bathyscaphe deep in the ocean. The geometry for this model is a hollow ball with a boundary load acting normal on all surfaces. In other words there are only SolidBoundaryLoadValue values acting on the surface. This poses a problem as the NeumannValue the SolidBoundaryLoadValue evaluate to are not sufficient to solve the PDE uniquely. More information on why that is the case can be found in the NeumannValue reference page. The solution is to use reduced geometry that exploits the regions symmetry much like in the 2D plane strain example. So instead of modeling the entire hollow ball, only an 1/8 segment will be used. This allows us to introduce symmetry constrains which are modeled with SolidDisplacementCondition which evaluate to the necessary DirichletCondition.
This section contains a variety of useful information. It contains utility functions that may be useful to better understand how the solid mechanics PDE models work. A different workflow for boundary condition predicate specification and a possible issues section about stress singularities.
The model parameters given are enriched with additional information. For example if a "Material" is specified then the appropriate "YoungModulus" and "PoissonRatio" are extracted from the material and stored in the model parameters. It is possible to inspect these enriched model parameters.
To perform a finite element analysis, the boundary mesh representation of the geometric model needs to meshed. To ease the set up of boundary conditions, the full mesh is generated before the boundary conditions are specified. In essence there are two types of boundary conditions. One are boundary conditions that operate on surfaces and essentially are of NeumannValue type. All solid mechanics boundary conditions names that end with a term "Value" are of this type. The second type of boundary conditions are of type DirichletCondition and operate on surface nodes of mesh. All solid mechanics boundary conditions names that end with a term "Condition" are of this type.
The position where a boundary condition is active is called a predicate. Various forms to specify these predicates exist. For one it is possible specify predicates as equations like in . Another approach is to select boundary surfaces and boundary points with ElementMarker. The generation and usage of ElementMarker is explained in the ElementMesh generation tutorial. The element markers used for boundary values in NeumannValue and boundary conditions in DirichletCondition are distinct. Boundary values use makers from boundary elements, while conditions use markers from point element. This allows for great flexibility. For complex geometries it may be desirable to set up both boundary values and conditions based on the boundary surfaces elements. This relates the point marker to be boundary markers, which is normally not necessarily the case. To do so the option "PointMarkers""BoundaryDeduced" needs to be specified.
Solid mechanics concerns itself with the computation of the deformation of objects under load and constraints. A load is a force or pressure that is applied on the surface of an object. For a load to be applicable to an object that object must also be constrained in some way, for example screwed to a wall, as otherwise the object would not pose a resistance to the load. Loads and constraints are set up by specifying boundary conditions. Various boundary conditions are available and are discussed in more detail in the boundary condition section. For the purpose of this example a boundary surface load and constraints introduced by a wall and screws will be sufficient.
In essence these conditions can be modeled with DirichletCondition and NeumannValue. An easy way to find the positions where boundary conditions apply on the surface is to browse through the boundaries of the geometric model.
A surface load is to be applied on the top of the books shelf bracket in a downward direction. Browsing through the surfaces by pressing + in the above manipulate the boundary marker identification (ID) number that corresponds to the top surface is found.
In a next step we would like to specify a constraint that models that the bracket is mounted to the wall. To fix the bracket to the wall two constraints play a role. For one, the back of the bracket cannot penetrate into the wall. Thus that boundary will be constrained to not be able to move in the -direction at all. It is obvious that the bracket cannot penetrate the wall and a constraint in the negative -direction is thus warranted.
One constraint we have is that the bracket is screwed to the wall. This means at the points of contact of the screw and the bracket no movement is possible at all. A second constraint is given in that bracket cannot penetrate the wall and a constraint in the negative -direction is thus warranted. Since the two screws press the bracket to the wall a reasonable approach is to also limit the movement in the positive -direction. This is a common simplification. If this simplification is not used then a contact problem arises because the bracket could in principal detach in the positive -direction from the wall.
Next, we construct the point marker for the screws. For this we select the relevant boundary element markers that make up the bold holes. We exclude the nodes that are contributed from the back plane.
In this section it will be illustrated that stresses are singular in inward facing corners. Following [6, p. 21], a simple bracket is fixed at the top and a force is acting downward on the right front face, illustrated by the blue arrows. The red line indicates where the stress singularities are expected. The exact dimensions and material parameters are not relevant only that the length scale is in millimeters. Such that we can scale the material parameters to millimeters. Once the region is created and material parameters are set up a sequence of refined meshes is used to show how a stress singularity is developing at the inward facing corner.
There is a small helper function that hides the computation of the von Mises stress and the total displacement. This function is defined in the Helper functions Appendix section. The left plot visualizes the von Mises stress with an indication of where the maximum stress is located. This maximum will be shown to be singular. At the same time the right plot visualizes the total deformation and while the stresses are singular the deformation is constant throughout the simulations.
The stress is highest in the inward facing corner. To better resolve the stresses the mesh is refined in the next step. As a first refinement approach the entire mesh is refined. This increases the number of elements.
Note, how the value of the computed von Mises stress has increased. At the same time the total deformation of the bracket remains the same. One could come to the conclusion that a further refinement would lead to a convergence of the stress value. One idea is to refine the mesh specifically in the section of interest.
See the SolidMechanics verification tests. It is important to realize that the verification models verify the mathematical model. Whether the mathematical model represents the actual situation at hand or not is a different matter. It is essential that limitations of the models used are understood. For example a low error in a plane stress model may not mean much if the model is not applicable in a specific scenario in the first place.
Setting up and solving solid mechanics PDE models can result in large scale PDE models that need a long time and a large amount of memory to solve. Since this is a challenge for all large scale PDE models various solution methods are presented in a different section of the finite element method documentation.