Introduction

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 notebook shows the necessary steps for everything except the CAD model generation.

The modeling process as such results in a system of partial differential equations (PDEs) that can be solved with NDSolve, ParametricNDSolve and NDEigensystem.

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.

Overview example and analysis types

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

Creating a solid mechanics model always comprises the same steps:

• Geometry creation

• Material selection

• Boundary condition specification

• Mesh generation

• Solving the partial differential equations

• Post processing the solution

The dependent variables , , and represent the displacement in the , , and -directions, respectively.

The output has been abbreviated since it is large and hard to read, but this is why using the PDEModel components is convenient.

In many cases like this there is no need to generate a mesh, passing the geometry to the solver is sufficient.

This gives an overview of the workflow. The next sections describe the steps in more detail on a more elaborate geometric model.

Geometry

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.

Sometimes it is useful to check the imported boundary mesh for defects. If defects are found a legend with the defects will be displayed. In this case no defects are found.

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.

In this case the bracket geometry does not have a screw thread. This will be taken care of when specifying the boundary conditions.

Generally speaking, one should start with an as simple geometric model as possible and only add details to see if and how they affect the overall solution.

More information on generating or importing 3D geometric models can be found in the Using OpenCascadeLink tutorial.

Material parameters

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.

A more convenient way to do material parameter setup is to select a material from an Entity. A convenient way to get to an Entity is it make use of the free-form input for WolframAlpha.

Material parameters where the values are given as Quantity objects can be used. 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.

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.

Units

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. If not explicitly stated otherwise, examples in this notebook use the default "SIBase" units.

Boundary conditions

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.

A way to find the positions where the boundary conditions apply is to visualize them together with the outline of the geometry defined previously in the geometry section.

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.

Next, we construct the constraint for the screws.

We simply use two large enough cylindrical regions that cover the screw holes and chamfer.

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.

Mesh generation

To perform a finite element analysis, the boundary mesh representation of the geometric model needs to be discretized into a mesh.

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.

Stationary analysis

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.

Post processing

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 displacements are often termed , and and are the computed displacements in the -, - and -direction, respectively. Sometimes a displacement vector is used.

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.

Sometimes it is useful to know the value of the automatic scaling factor computed.

Inspecting the total displacement is a way to get an overview of the displacement of a constrained object under load. The total displacement is given by:

where , and are the computed displacements in the -, - and -direction, respectively.

Strain is a quantity that describes the amount of deformation within a body [2] 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.

Loosely speaking strain is the gradient of the displacement. Strain plots show where the gradient of displacement is large not where there is much deformation.

Note, in this example the absolute values of the strain in the -direction are smaller than absolute values in the -direction, even though the displacement in the negative -direction is much larger than in the positive -direction.

Sometime it is useful to visualize the strain results with the same color scale. For this we find the minimum and maximum of all strain values and use those to set the bounds for the color scale.

Sometimes one would like to identify the regions where strains are larger than a specific value.

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.

The equivalent strain computes:

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.

Sometimes it is more convenient to plot the strain over the deformed body.

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.

Principal strain values are ordered according to their size.

At each evaluation coordinate there are three principal strains. Each of those principal strains can be interpreted as a vector and visualized.

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

where the the normal stresses are denoted by and the shear stresses by . The first index tells you the direction of the force, while the second index specifies the direction of surface the force is acting on.

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.

The von Mises stress computes:

Sometimes it is more convenient to plot the stresses over the deformed body.

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.

Evaluated at a coordinate returns a list of three values: the first, second and third principal stress value.

A legend is not meaningful here, as an eigenvector multiplied by a scalar is still an eigenvector.

The section Boundary Load: Tension has an example that shows the use of PrincipalEigensystem to find stresses in a non-axial loaded cylinder.

Since PrincipalEigensystem sorts the eigenvalues from largest to smallest, it cannot be used for a symbolic computation of principal eigenvalues as sorting of symbolic expressions is not possible.

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.

We compute the principal eigenvalues of stressed body to transform the stresses into load direction independent stresses, which then can be used to compute the von Mises stress.

The equivalent von Mises stress can be computed with the von Mises stress functions if the shear stresses are set to 0.

For ductile material plastic deformation starts when the stress reaches the yield strength (for brittle materials when stress reaches the ultimate strength).

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.

Every object that is constrained and at the same time under load will have reaction forces that balance the sum of the loads. These reaction force will be at the constrained parts of the body.

The function PDEModels`ReactionForce returns the displacements just as a normal call to NDSolve would and a second set of InterpolatingFunction that are the reaction forces.

The sum of the reaction forces matches the values of the boundary loads.

The shown reaction forces counteract the (not shown) downward force that acts on the bracket.

In the current version, the reaction forces can only be computed for linear stationary solid mechanics models.

Time dependent analysis

A time dependent analysis will have to be used if some of the PDE model components have time dependent behaviour. A common time dependent component are time dependent boundary conditions.

An important numerical aspect of dynamically loading an object is that the load is not applied instantaneously. For one such a model is unphysical and will result in long simulation times.

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.

To monitor a the movement of a specific point within the body over time a query point is set up. The query point can then be traced through time.

Note, that the system is not damped and the amplitude of the displacement remains the same over the cycles.

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].

Eigenmode analysis

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.

There are multiple reasons to compute eigenfrequencies and eigenmodes among them are:

• Objects in resonance are a common cause of failure. Knowing these frequencies helps avoid the shapes being exposed to them.

• The resonant frequencies are a basis for a frequency response analysis.

The generalized equations of motion for solid mechanics for a linear elastic material is given as:

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 [7] and leaves us with:

In an eigenfrequency analysis constraints are considered but loads are not considered thus .

For an eigenfrequency analysis one is interested in the solution of an equation of this form

Here, is the solid mechanics PDE component operator that provides and . is the eigenvalue .

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 eigen values are related to the angular frequencies by:

The angular frequency is related to the frequency measured in by:

The natural frequencies can be computed from the eigen values by:

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.

Note, that the first 6 eigenvalues are relatively close to 0 and represent the rigid body modes. Also note that rigid body modes can have eigen values with small imaginary parts.

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.

If a modal analysis reveals rigid body modes then the object is not constraint enough. An eigen analysis is a good check to verify that the body is sufficiently constrained.

Parametric analysis

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.

Note that the results remain qualitatively the same, but the magnitude of the total displacement changes linearly with the force applied.

A parametric analysis can be use to create force-displacement plots. For a specific query point the displacement is tracked while the force increases.

In this case we see a linear behaviour between force and displacement, which is expected for the linear material model that is used by default.

It is possible to have parametric material laws. For example, the Poisson ratio, or a heat transfer coefficient could be parametric. This cases works like any other parametric case. The only thing to be aware of is that material law parameter values need to be respecified when computing the strain and stress.

In this case the strain computation does not depend on the Poisson ratio parameter and can be done in the usual way.

The stress computation, however, does depend on the Poisson ratio parameter . In this case the parametric parameters need to have the actual parameter value for the computation. This can be done by replacing the symbolic parameter value with the actual value.

Frequency response analysis

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.

In the solid mechanics domain a frequency response analysis in effect solves:

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 this case we consider, like above, the undamped case such that .

To perform a frequency response analysis a frequency dependent load or constraint needs to be set up. This is in contrast to an eigenmode analysis where boundary loads are not relevant.

The remaining boundary conditions are the stationary boundary conditions.

Choose a frequency range between from the minimal eigenfrequency to the largest eigenfrequency computed.

To get a better understanding of the result we compare frequency responses with the stationary solution. For the frequency response analysis we used a frequency dependent load

of 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.

We see that close to the resonance frequency the maximal displacement of the query point is about one order of magnitude larger than the displacement of that same point in the stationary analysis.

In the previous section the eigenvalues of the constrained bookshelf were computed.

Overall this means that when the body is in resonance the amount of deformation can be much larger than what is caused by a comparable stationary force.

Note that evaluating the parametric harmonic function at the resonant frequencies computed during the eigenmode analysis may not be possible.

Equations

Overview

Solid mechanics is the analysis of bodies under load and constraints. The entire topic of solid mechanics hinges around four components:

• Applied Load

• Stress

• Strain

• Displacement

These components are related in the following manner:

Before going into too much detail we start with a broad overview of the different solid mechanics components.

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

where a function relates the strain and other field quantities to the stress. The most important relation is Hooke's law where the relation between stress and strain is linear and given by:

where the proportionality constant is the Young's modulus. Most metals and alloys can be modeled as an elastic material if the strains remain small.

If we can not assume a linear relation between stress and strain we have a nonlinear constitutive equation. This is also called a nonlinear material law.

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 current version allows for elastic material models that be linear or nonlinear. Plasticity is reserved for a future version.

Equilibrium equations

No material properties (constitutive equations) are used in the derivation of the equilibrium equations and as such they are applicable to all materials. The equilibrium equation is given by:

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.

The units of the solid mechanics model terms are force density in .

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.

Kinematic equations

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.

A quantity related to the displacement is the displacement gradient often denoted by :

For a rigid body translation the displacement gradient is 0.

Deformation describes the relative movement of particles next to each other while displacement describes absolute movement. Deformation captures changes in size and shape while displacement does not.

Strain is a quantity that describes the amount of deformation or distortion within a body [2] 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.

The infinitesimal strain measure is also called the engineering strain measure or small strain measure.

A simple strain definition for a 1 dimensional rod is given by

where is the length in of an undeformed object, is the length in of the deformed object and is the change in length in due to an applied load.

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 same 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.

First, we consider the changes in length. Consider a line segment parallel to the axis in the original configuration. A displacement moves this line to in the final configuration.

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:

Since is parallel to the axis the increments are functions of only and we can write

Substituting back we get

The normal strain in the -direction is then defined as

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

Thus the normal strain is approximated with

The same reasoning can be applied for the remaining normal strains in the - and -direction such that we get

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

Along and along such that

Using this we get a shear stress of

Proceeding in a similar way we get

Now, one could go ahead and group the components in a matrix like this

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 strains are defined as engineering shear strains

The strain tensor components are symmetric and the complete infinitesimal strain tensor looks like this

or written out as

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.

Using the tensorial strain we can succinctly write the tensor as

where the vector is the displacement gradient.

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 function maps to :

We express the infinitesimal line segment :

Linearize:

After reinserting we get:

where is the identity matrix and is the deformation gradient with components

where is the Kronecker delta. In matrix form can be written as

Or, written succinctly:

The deformation gradient can also be expressed in term of the spatial and material coordinates. We note that:

from with follows that:

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 deformation gradient is at the heart of nonlinear solid mechanics.

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

where is the length in of an undeformed object and is the length in of the deformed object.

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.

In the case of the Green-Lagrange strain the normal and shear strains are zero, as expected for rigid body motions.

The Green-Lagrange strain measure is implemented by computing the deformation gradient

where is the identity matrix. The Green-Lagrange tensorial strain is given by:

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:

where is the length in of an undeformed object and is the length in of the deformed object.

The Euler-Almansi tensorial strain is given by:

The actually used strain measure can be inspected.

Poisson's ratio gives the relation between the lateral contraction and longitudinal expansion when a body is pulled in the longitudinal direction

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:

Because the material is isotropic and the relation of to and is Poisson's ratio:

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.

Constitutive equations

The constitutive equations relate the equilibrium equations with the kinematic equations via the material model. Before discussing these we will introduce the concept of stress.

Stress is a quantity the describes the distribution of internal forced within a body [2] and is measured in or

where is the force in and is the initial area in the force acts on.

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.

At a minimum the data needs to contain the load applied by the machine and the measured strain.

For the stress computation the diameter of the specimen is needed.

Typically, the last data point is a recording of the rapture of the material and needs to be removed.

The above material data is from a cold rolled steel and as such is a ductile material. A more brittle material, for example Aluminum, would have a much less distinct flat top part.

The areas under the curves correspond to the absorbed energy.

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.

The following points of interest are marked:

• 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.

• E: Fracture marks the point where the material breaks.

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 [14].

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.

Consider an undeformed cuboid with length and area shown in here as a 2D cross section. The specimen is subject to forces under which it deforms and now has length and area .

The true stress is given as:

The engineering stress is given as:

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

The engineering 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.

Assuming a constant volume:

Under the constant volume assumption the following conversions can be established:

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 Wolfram Language uses engineering strain.

Linear elastic material models

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 :

Sometimes the elasticity matrix is also called the constitutive matrix or the material stiffness matrix.

In the most general case the elasticity matrix of the linear elastic constitutive equation is given as:

The are the shear stresses. These are used to indicate that we used engineering shear strains .

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. We typically distinguish between, isotropic and anisotropic materials. In an isotropic material all material properties are the same in all directions. Most metals are isotropic material. All materials that are not isotropic are called anisotropic. In the most general case, anisotropic means that a material behaves different in all possible directions.

Because the general anisotropic case is indeed very general, special named sub-cases exist. For example, in an orthotropic material the material properties are different in each of the object's axis directions. Wood is considered an orthotropic material where the material properties depend on the direction of the wood grain.

Generally speaking, anisotropic material are typically compound material or biological tissue where the material properties vary in some or all directions. Also fibre reenforced material where the weave effects the material properties in all directions are considered anisotropic.

For anisotropic material, one thing to keep in mind is that the orientation of your geometric model must match the orientation of the material properties.

A linear elastic material model is good for the vast majority engineering design calculations, where components cannot exceed yield [11, c.1.1.5].

The subsequent sections will introduce these linear material models.

In the case where both the constitutive equation and the kinematic equation are linear the equations can be simplified into a form like:

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].

The elasticity matrix for an isotropic material is given by:

The Young's modulus is and Poisson's ratio are the only elastic coefficients needed to specify the elasticity matrix.

This is the linear variant of the term .

For the linear elastic isotropic case the commonly used Young modulus and Poisson ratio can be expressed in terms other elastic moduli.

Currently, the functionality that the various moduli can be expressed in terms of each other is available in 3D only.

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.

The general form for orthotropic materials is given as:

To keep things simple the coefficients of elasticity matrix are typically not given explicitly but are defined through the compliance matrix . Here we have :

Note that the elasticity matrix can be found by inverting the compliance matrix such that .

This is the linear variant of the term .

Transverse isotropy is a special type of orthotropy. Orthotropic materials are symmetric with respect to three orthogonal planes. Transversely isotropic materials are symmetric with respect to one orthogonal plane.

Transverse isotropic materials have the same physical properties in an infinitesimal thin layer. Sedimentary rocks are an example. These rocks are made up from various layers, say in the - and -direction. Material properties, like stiffness, are the same in the - and -directions. In the direction of the layers, the -direction, however, the physical properties change and the material is orthotropic with respect to that direction.

In other words, transversely isotropic materials have one axis of symmetry around which the material is symmetric. This means that if we rotate the material about the axis it behaves like an isotropic material in the plane to which the axis is normal. Transversely isotropic material is a material class that sits between isotropic and orthotropic material.

Another example of a transverse isotropic object would be a biological membrane. In this case the material properties of the membrane will be the same in the plane of the membrane but different from those in the perpendicular direction.

The material properties of sedimentary rocks are the same in the - and -direction, which form the plane of isotropy. In the -direction the material properties differ.

A second type of materials falling into the category of transversely isotropic materials are materials that are reinforced by one family of fibers running in the same direction. In this case the fibers take the role of the axis around which is the material symmetric. Fiber-reinforced materials are usually much stiffer in the direction of the fibers than in the transverse direction.

A fiber-reinforced composite with a single type of a fiber embedded in a material matrix.

A difference between fiber-reinforces materials and sedimentary rocks is the structure of their cross section. A sedimentary rock is homogeneous in the plane of isotropy, while a fiber-reinforced material is not.

The left view shows the fiber-reinforced material from a front view, while the right illustration shows a 2D cross section. This cross section is the plane of isotropy but shows that, strictly speaking, the plane of isotropy is not homogeneous.

So how do we obtain material parameters? Material parameters of some composites can be found in relevant literature, for example for fiberglass see [25]. If the parameters of our material are nowhere to be found, then the best way is to do some experiments. This option is however not always feasible. If material parameters of the matrix and fibers are known than not everything is lost. One approach would be to use a homogenisation simulation to obtain the material parameters of the composite. For the process see for example [26], a survey of the homogenization theory and related techniques can be found in [27]. Another possibility is to just simply use the material parameters of the matrix and fibers. The material parameters of the matrix can be used in the transverse direction and the parameters of the fibers can be used in the axial direction.

The transverse isotropic case is only suitable for fibers running in a single direction, it is not applicable to woven fibres.

A general overview of the mechanics of transversely isotropic materials can be found in [16]. An extensive overview with focus on fiber-reinforced materials can be found in [24].

Since transverse isotropy is a special type of orthotropy, we can start from the general form of the orthotropic material, which has 9 independent components:

Now, if we consider a transversely isotropic material with the axis of symmetry in the direction, we obtain the following:

The number of independent components is now reduced to 5. The compliance matrix is commonly written as:

where is the Young modulus, is the Poisson ration and is the shear modulus with the condition that . This relation is similar to the one that holds for isotropic materials.

Since and are material parameters in the direction of the axis, we can denote them with a subscript . Material parameters and are associated with the transverse direction, therefore we can denote them with a subscript . This gives the following formulation, in which the material parameters are split into two groups, one for the transverse direction and one for the axial direction:

once again with the condition .

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].

To set up an anisotropic material one provides the full elasticity matrix with the parameter "ElasticityMatrix".

The fully anisotropic elasticity matrix needs to be specified in Voigt notation with the ordering .

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.

An alternative is to specify the full elasticity tensor.

This is the linear variant of the term .

An example of how to set up a multi-material compliance matrix can be found on the SolidMechanicsPDEComponent reference page. The same procedure can be applied for a multi-material elasticity tensor.

So far we have considered only orthotropic materials where the axes of symmetry were aligned with the -, - and -axis. The same goes for the transversely isotropic materials where we assumed that the axis of symmetry is in the -direction. It's not always possible to meet these assumptions: Wood grains or steel beams in concrete, for example, are not always straight. Additionally, the set up of the model may be complicated: Consider, for example, a cubic block of wood with diagonally aligned grains. We will show, that instead of aligning the geometry with the -direction, it is much easier to align the material model to grain direction.

For simplicity, we will consider transversely isotropic material as an example, but the approach can be used for orthotropic materials or even fully anisotropic materials.

The equation of linear elasticity can be written as

where is stress, is strain and is elasticity matrix. In our case we cannot simply use this equation as the orientation of our material varies from point to point. We can, however, use it locally in an appropriately oriented coordinate system.

At each point we can choose a coordinate system in such a way, that the orientation of our material model is in the -direction. We then transform the -direction aligned coordinate system to some fixed coordinate system using an appropriate rotation matrix. We can then set up our model in our fixed global coordinate system and the transverse isotropy of the material will be satisfied as it will be locally enforced.

In the case of a non-axis aligned material symmetry, the rotation matrix rotates the elasticity matrix from the fixed basis coordinate system in such a way that the orientation of the elasticity matrix is in the -direction at each point in the non-aligned object. The orientation of the model is given by the axis of symmetry in the case of a transversely isotropic material.

Let us take some fixed coordinate system and denote the direction of the axis of symmetry in the reference configuration by . Then we have a rotation matrix , which rotates our fixed coordinate system such that the axis of symmetry is in the -direction and therefore:

This rotated coordinate system is sometimes called privileged, laboratory or anisotropic and we will denote quantities in this coordinate system with superscript . We can thus transform the infinitesimal strain tensor and the stress tensor from the fixed basis into the privileged basis back and forth:

where is the infinitesimal strain tensor in the privileged basis and is the stress tensor in the privileged basis. Notice that the rotation matrix does not need to be constant, it can be a function of space. This allows us to deal with materials in which the orientation varies.

Now, with and we can use the linear elastic constitutive equation (1), because in the privileged basis the material is locally properly oriented. This gives us the following relation:

In case of transversely isotropic material, this translates to

Let us now take a look at how we can use this approach in an actual computation. We will present an example of a transversely isotropic material with a variable axis of symmetry. We will consider a unit cube, which will be held at one side, while pressure will be applied on the other side.

Note, the the orientation matrix is a function of the spatial coordinates.

Note, that we also set the "EngineeringStrain" to False. In this version an orientation matrix can not be used with the engineering strain. This is worth while mentioning, as the engineering strain is the default strain for linear elastic material but not in this case.

The actually used linear material model can be inspected.

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:

The next sections describe the various components.

Initial strains can be used to model offsets in the strain caused by physical phenomena like temperature. The initial strain is subtracted from the computed strain and thus:

Scroll to the right to see the strain contribution.

For multi material models the initial strain can be set up with branching construct like an If statement.

The generation and usage of ElementMarker is explained in the ElementMesh generation tutorial.

Thermally induced elasticity is modeled as an additional strain. The thermal strain is a common form of an initial strain and given by:

Temperature changes induce thermal strains. Thermal strains can result in thermal stresses if

• thermal gradients exist

• multiple material with different thermal expansions are present

• constraints restrict free movement

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

where is the specific heat capacity, is the mass density and thermal conductivity . For a constant heat flux and a time scale that is larger than a steady state analysis is sufficient.

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.

To verify the result the axial reaction force can be computed on one of the constrained surfaces. This can be done with

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.

Next, the reaction force on one of the end caps of the cylinder is computed.

Next, for an orthotropic material the thermal expansion coefficient has to be specified in the respective directions:

Scroll to the right to see the thermal strain contribution.

Next, for an anisotropic material the thermal expansion coefficient has to be specified in the respective directions:

Scroll to the right to see the thermal strain contribution.

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.

Here the maximal displacement is in the negative -direction and has thus a negative sign.

A similar examples of coupling heat transfer with solid mechanics can be found in the application examples: Thermal and Structural Analysis of a Disc Brake and Thermal Load on a Beam.

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:

• Use the thermal strain function itself

• Specifying the secant coefficient of thermal expansion and a reference temperature

• Specify the tangent coefficient of thermal expansion

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.

The following illustration helps understand the difference between the tangent and the secant coefficient of thermal expansion.

The tangent coefficient of thermal expansion is defined as:

It follows that the thermal strain is

If we define an average coefficient of thermal expansion as

The thermal strain is then

The Wolfram language PDE models require that is given, not .

If the strain temperature over a given region is linear then the secant and tangent coefficient are the same.

Initial stresses are modeled by an addition of an initial stress term to the constitutive equation.

Scroll to the right to see the initial stress contribution in the output.

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].

For multi material models the initial strain can be set up with branching construct like an If statement.

The generation and usage of ElementMarker is explained in the ElementMesh generation tutorial.

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 a 2D simplification 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.

The derivation of the 2D PDE models starts from the 3D isotropic elastic model. Orthotropic and anisotropic 2D models are currently not supported.

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; in fact that is not the case. These points should just 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 the -direction compared to the -plane expansion for the plane strain model to be valid. Analysis of dams and shafts fall into this category.

The starting assumption for the plane strain model is that there is no displacement in the -direction. This means that

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

This simplifies to

This can be written as

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.

If there is a thermally induced strain with then the strain vector becomes

and the same analysis as above is done. This leads to

and thermal strain

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 quarter pipe structure exploits a symmetry condition -direction on the left.

Inside we have a pressure of 20 units in the outward direction.

The remaining sides are free to move.

As discussed above, in the plane strain case there is no -directional strain component but there is a -directional stress component . You can compute these manually, but for convenience there are functions that compute these components. The functions take the computed displacement, strain and stress, which are structured arrays of dimension 2 by 2 and convert them to structured arrays of dimension 3 by 3 and the component is then found in the {3,3} position.

For a plane strain set up the strain component in the -direction is 0.

For a plane strain set up the stress component in the -direction is given by:

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.

The plane stress derivation starts from the inverse strain - stress relation

This can be written as

Expressed as a stress strain relation this becomes

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.

If there is a thermally induced strain with then the strain vector becomes

The following example explains the usage of the "PlaneStress" model form with multiple materials and thermal expansion.

The generation and usage of ElementMarker is explained in the ElementMesh generation tutorial.

In the plane stress case there is a -directional strain component but no -directional stress component . These values can be computed manually, but for convenience there are functions that compute them. The functions take the computed displacement, strain and stress, which are structured arrays of dimension 2 by 2 and convert them to structured arrays of dimension 3 by 3 and the component is then found in the {3,3} position.

For a plane stress the strain component in the -direction is given by:

For a plane stress set up the stress component in the -direction is 0.

An axisymmetric simulation can be performed when a 3D solid has an axis of revolution. An axisymmetric model uses a truncated cylindrical coordinate system. Compared to a Cartesian coordinate system, the radial axis direction takes the place of the direction, the angular direction the place of and becomes the symmetry axis. The order of the independent variables is fixed and can not be changed.

The illustration shows a body that has an axis of revolution around the dashed -axis. The body can be represented in a axisymmetric setting by making use of a 2D area that is rotated around the -axis, following the angular direction . In the case that the material properties and loading are also symmetric about the -axis a 2D axisymmetric model can be used, which is depicted below as the 2 dimensional area embedded in 3D.

The derivation of the 2D axisymmetric solid mechanics equations starts from the 3D Cartesian version. The displacements are still called , and . Now, is the displacement in the radial direction , the displacement in the angular direction and in the axial direction . In cylindrical coordinates the infinitesimal strain is defined as:

Remembering that the tensorial shear strains are defined as engineering shear strains :

We arrive at:

The simplifying assumption for the axisymmetric model is that there is no displacement in the -direction. This means that . Also the radial displacement and the axial displacement are independent of :

Next, we use the fact that there is no displacement in the -direction. The simplified axisymmetric strain tensor is given as:

The fact that there is no displacement in the -direction also means that the strains and . Note that the strain is not zero and the same is true for the stress, that is also non zero.

Note, that to compute the strain tensor above we directly used the dependent variable specification .

The stress-strain relation is coordinate system independent and does not change. The linear elastic isotropic relation is given by:

Using the axisymmetric shear strain and stress specifications the axisymmetric stress-strain relation is now given by:

As we have seen above, the fact that there is no displacement in the -direction also means that the strains and . From this follows that the stresses and .

The reduced axisymmetric strain tensor can equivalently be given as:

Compared to the plane strain and plane stress case, the axisymmetric additionally requires that the equilibrium equations be adapted to the cylindrical coordinate system too. The equilibrium equations in cylindrical coordinates are given by:

In the above only the first and third equation contribute to the equilibrium equation.

Let's look at an example.

One important point to realize with respect to the boundary conditions is that if the geometry starts from that mean that the rotated body in 3D does not have hole in around the -axis, and you setting boundary conditions on does not make physical sense.

In the axisymmetric case the stress computation needs the computed displacement as an argument.

Very much like in the plane strain and plane stress case there are strain and stress components in the direction. For the extended strain computation the function takes the computed strain and stress, which are structured arrays of dimension 2 by 2 and convert them to structured arrays of dimension 3 by 3 and the components are then found in the {2,2} position.

For an axisymmetric setup the strain component in the -direction is given by:

For an axisymmetric set up the stress component in the -direction is given by:

In the axisymmetric case there is another alternative way to immediately get the strain and stress in the direction. To make use of this the dependent variable is set to 0 in the dependent variable set up.

With setting the axisymmetric case is treated like a 3D set up. The actual computation, however, will still be done in 2D. This has several advantages. First, the retuned displacement will also have set and, secondly, both stress and strain will have the components set. Vector valued parameters are now specified as 3 component vectors.

Note that the second, the second component of the displacement .

A time dependent axisymmetric equation can be set up in exactly the same manner as in the 3D case.

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.

• Strains may not be small and large deformations need to be considered.

• Plastic, irreversible, deformation may take place.

The question of when an elastic model is applicable is defined in terms of stress:

The function is referred to as a yield function. A material is said to yield at . Two common yield functions are the Tresca and the von Mises criterion.

The Tresca criterion is given by:

where are the principal stresses (eigenvalues) of and is a material parameter.

The von Mises criterion is given by:

with a material parameter.

The material parameters and define the condition after which the linear elastic model no longer holds and inelastic deformations take place.

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

where is a 4 tensor and the displacement vector. Let's illustrate what the issue is.

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

Let's take a linear elastic material law where we have:

This form of the material law is a short hand. When we convert the Voigt notation to the full 4-tensor notation we get the coefficients for .

Nonlinear elastic material models - Hypoelastic models

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 can be 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.

• The usage of hypoelastic material is not with controversy.

Hypoelastic material models are used to model behavioral phenomena.

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.

Conceptually speaking, these models are derived by defining a strain energy density function and finding the derivative with respect to the strain:

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 [11] but the hypoelastic constitutive equation is given as:

with parameters:

where is the slope of the uniaxial stress-strain curve at . The model parameter , and are parameters that need to be estimated from experimental data. This will be shown later.

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.

To try this custom material model out we use a simple setup. A rectangular region is constrained on the left and bottom and a pressure is applied in the positive -direction.

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.

A good way to evaluate the same model for various pressure inputs and an updated InitialSeeding is to make use of ParametricNDSolveValue.

In the following sampling we use the previous solution to start the nonlinear solve for the current solution:

The above output shows a highly nonlinear relation between the applied force (through setting a pressure) and the displacement the object shows.

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.

We want to visualize the data. The last data point, however, is the data point where the specimen ruptures and thus we remove that.

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.

The next step is to estimate the model parameter .

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.

Hyperelasticity

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.

Failure theory

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 [12].

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.

The Tresca failure theory is also called the maximum shear stress theory.

Yielding occurs when the maximum shear stress is equal to the shear stress at yielding in an uniaxial tensile test [12].

In a uniaxial test is equal to the applied stress, since that is the uniaxial direction in which the stress is applied. Correspondingly the and stresses are 0. Therefore or

The Tresca factor of safety can be computed with:

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 consistent with the observation that hydrostatic stresses do not contribute to yielding in ductile materials.

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.

The Tresca theory is a special case of the Coulomb-Mohr theory.

The following example demonstrates computing the factor of safety for a stress tensor.

See also the safety factor section in the overview section for an example.

The von Mises failure theory is also called the maximum distortion energy theory.

Yielding occurs when the maximum distortion energy is equal to the distortion energy at yielding in an uniaxial tensile test.

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.

The Coulomb-Mohr theory is used for brittle material. The Tresca theory is a special case of the Coulomb-Mohr theory.

The modified Mohr theory extends the Coulomb-Mohr theory to make it fit better with experimental data.

Multiple materials

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.

While not strictly necessary, it is better to have the geometric model represent the material boundary by an internal boundary.

Note the division in the middle of the beam. There will be one material to the left and one to the right of the division.

Less manual and more general methods to create geometries with internal boundaries is described in the section ElementMesh with Subregions. In 3D using the OpenCascadeLink is a good option.

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.

Should the materials to be used not be available as an Entity then the material parameters can be given as conditional statements, like an If statement.

Note that the maximal total displacement is to the right of the material boundary. Lead is a much softer material than Titanium and so this shift is expected.

Subsequent analysis are done in exactly the same way as in the single materials case.

For more involved geometries it might be easier to make use of ElementMarker to specify sub-regions, boundary parts and boundary nodes. The generation and usage of ElementMarker is explained in the ElementMesh generation tutorial. That section also explains how to evaluate functions that contain ElementMarker on the mesh.

Damping

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.

The generalized equations of motion for solid mechanics for a linear elastic material in matrix form is given by:

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:

and transforms the equations into a system of a generalized wave equation, that additionally to the second order time derivative now also has a first order time derivative.

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.

In order to monitor the progress during the time integration an EvaluationMonitor has been added. Also note that a MaxStepSize has been added to make sure no solution feature will be stepped over.

For completeness, we also compute the stationary solution.

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.

Getting good values for the Rayleigh damping parameters and can be tricky. It can be shown that the damping ratio for different natural cyclic frequences can be calculated as:

The following plot visualizes various damping types by plotting the cyclic frequency against the damping ratio .

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:

With the additional constrain:

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.

Now, we assume damping ratios of () and () for the first and second frequency respectively. If values for the damping ratios are not known, we show a procedure for measuring them further down.

The mode based damping is now a combination of a damping ratio of 1% mass damping at frequency and ratio of 2% stiffness damping at frequency .

If values for the damping ratios are not available they can be measured with an impulse test. For this the logarithmic decrement, defined as

is used. Here is the number of periods, the amplitude of a peak and the amplitude of a peak periods away.

It is important but difficult to excite just single mode. The damping ratio for that single specific mode can then be computed as:

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.

The goal is to find the damping ratio for this system ring-down. First, we find the times at two peaks of the decaying solution:

With this procedure, one damping ratio coefficient was found, the experiment is repeated to find the second damping ratio coefficient.

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.

We've shown and used a fictitious ring-down data set since obtaining and using actual data is quite complicated.

In the finite element programming tutorial there is a section that shows how Rayleigh damping can be implemented on a system matrix level.

Solid mechanics boundary conditions

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.

To explain various boundary conditions a solid mechanics PDE component with default units is set up.

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.

Displacement constraints

To explain the various constraint cases a simple beam is modeled. The beam is fixed at the left.

Note how the right hand face has been shifted downward by 1 unit in the negative -direction. The displacement in the and -direction is zero, as specified.

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.

When prescribed displacement of 1 unit in only the negative -direction is specified, note that there is also a displacement in the -direction at the right end.

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.

Generally speaking a roller constraint would constraint the dependent variables along the normal in the following manner:

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 has been specified. For good reason: Not specifying a DirichletCondition will make the system of equations singular and the solution can only be found, at 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 its 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 the degrees of freedom that constrain the rigid body movement need to be constrained. The continuum finite elements used 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.

Let’s look at a 2D example:

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.

Other choices are possible for example:

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.

To illustrate that the actual choice of the minimal constraints is arbitrary, we consider a rectangular geometry with an unsymmetrically placed hole inside.

The example use a plane stress model form. All outer boundaries are under pressure.

We solve this model twice, each time imposing the minimal number of constraints at different positions.

In the next step we solve the problem again, this time with the constraints at different positions.

Clearly the displacements are different. In both cases the node at {0,0} has constraints on u and v. In the first case, on the left, the second constraint on v is set at the lower right. In the second case the second constraint on u is set at the top left. While the displacements are different for the different choices of the constrains, the strain and stress are not.

We can see that the vonMises stress is the same in both cases.

Boundary load

When considering boundary load, we often speak of traction. Traction in units of []is a vector that is defined as:e

In other words traction is a force per area, just like stress. Written out in component form we have

When the stress tensor is symmetric we can also say:

In the absence of body forces, when the body is in equilibrium the integral over the traction along the boundary should be zero:

Since the Gauß theorem also holds for tensors we have:

The left hand side boundary integral is a NeumannValue as explained on the NeumannValue reference page or in the derivation of the Neumann boundary value in Partial Differential Equations and Boundary Conditions. The right hand side is the equilibrium equation . So one can use NeumannValue to specify traction vector components.

A load is either a pressure or a force acting on a surface. Both can be specified with 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 .

SolidBoundaryLoadValue will evaluate to a list of NeumannValue.

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.

To illustrate the various loading scenarios a simple cylinder geometry is created.

In all subsequent examples the rod is fixed to a wall on the left.

Boundary loads need not be in axial directions. For example consider a pressure field acting normal to the body of the cylinder.

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.

Note that the stress corresponds with the applied surface pressure.

Note, that there is a stress concentration at the left, the fixed end - that is expected.

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.

To keep things simple we use an arbitrary rotation matrix to rotate the cylinder and also rotate the boundary load to be normal to the right end of the cylinder.

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.

We see that the largest eigenvalue, the first principal stress, now corresponds to the stress value if the cylinder and load had been axis aligned.

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.

In this example a rod is fixed at the left and a torque of is applied on the right end.

The torque needs to be converted into a surface pressure. Starting from

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.

A point rotated by will end up at position . Now, with the small angle assumption the point will end up at . This is what is picked up by the element mesh deformation plot.

The shear angle is a little less than one degree for a moment of torque of .

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.

We are barely above the yield strength of Titanium with this amount of torque.

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.

The small angle requirement is no longer satisfied.

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.

Next, we plot the deformation.

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.

Symmetry conditions

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.

One each of the cut surfaces the symmetry constraint fixes the object in the plane of the cut and at the same time allows for movement in the other directions.

Appendix

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.

Model parameters

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.

Boundary condition predicates

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 markers 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.

The found ElementMarker ID can then be used to specify a surface load value.

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.

The exact details of the boundary condition specification will vary between different analysis types and are shown in the appropriate sections.

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.

The exact details of the boundary condition specification will vary between different analysis types and are shown in the appropriate sections.

Stress singularities

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.

To get a better feeling for the density of the mesh, it is instructive to look at the wireframe of the generated mesh.

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.

The stress value has further increased and it is even more evident that the stress is maximal along the inward corner. To avoid these stresses, typically, corners are filleted.

Note that now, even after a refinement the stress value is about the same. Using fillets is also physically more realistic than the infinitely sharp corner.

Verification

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.

Large scale finite element models

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.

• Solving Memory-Intensive PDEs

• Linear Solver options

Nomenclature

References

1.  Backstrom, G; Simple Deformation and Vibration; GB Publishing; 2006

2.  The Efficient Engineer, https://www.youtube.com/watch?v=aQf6Q8t1FQE; retrieved 2021/4/27

3.  Cook, R; et. al.; Concepts and Applications of Finite Element Analysis; Wiley 2002; 4th Edition

4.  Bhatti, M.; Fundamental Finite Element Analysis and Application; Wiley 2005

5.  Bhatti, M.; Advanced Topics in Finite Element Analysis of Structures; Wiley 2006

6.  Brand, M.; Grundlagen FEM mit Solidworks; Vieweg+Teuber 2011; ISBN: 978-3-8348-1306-0

7.  https://en.wikipedia.org/wiki/Modal_analysis_using_FEM, retrieved 01/06/2021

8.  Nasdala, L.; FEM-Formelsammlung Statik und Dynamik; Vieweg+Teuber 2010

9.  Alawadhi, E. M.; Finite Element Simulations Using ANSYS; CRC Press 2010

10.  Constantinescu A., Korsunsky A.; Elasticity with Mathematica; Cambridge University Press 2007

11.  Bower, A.; Applied Mechanics of Solids; CRC Press 2010; ISBN-13: 978-1439802472; (http://solidmechanics.org/contents.php)

12.  The Efficient Engineer, https://www.youtube.com/watch?v=xkbQnBAOFEg; retrieved 2021/9/13

13.  Kelly, P.A.; Mechanics Lecture Notes: An introduction to Solid Mechanics. (http://homepages.engineering.auckland.ac.nz/~pkel015/SolidMechanicsBooks/index.html)

14.  Moosbrugger, C. (Editor); Atlas of Stress-strain Curves; ASM International 2002; ISBN 0-87170-739-X (https://books.google.com/books?id=up5KS9fd_pkC)

15.  Sifakis, E., Barbič, J.; Finite Element Method Simulation of 3D Deformable Solids; M&C Publishers 2016; ISBN 9781627054423

16.  Holzapfel, A.; Nonlinear Solid Mechanics; Wiley 2020; ISBN 978-0-471-82319-3

17.  Bonet, J., Wood, R.; Nonlinear Continuum Mechanics for Finite Element Analysis; Cambridge University Press 1997; ISBN 0-521-57272-X

18.  Mooney, M. Large elastic deformations of isotropic materials. Journal of Applied Physics 1940. 11: 582–592. https://doi.org/10.1063/1.1712836

19.  Rivlin, R. S. Large elastic deformations of isotropic materials. Philosophical Transactions of the Royal Society of London. Mathematical and Physical Sciences 1948. 241: pg 379–397. https://doi.org/10.1098/rsta.1948.0024

20.  Arruda, M. Boyce, M. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids. 1993. Volume 41, Issue 2. Pages 389-412,. https://doi.org/10.1016/0022-5096(93)90013-6

21.  Gent, A. N. A New Constitutive Relation for Rubber. Rubber Chemistry and Technology. 1996. 69 (1): 59–61. https://doi.org/10.5254/1.3538357

22.  Treolar, L.R.G.; Stress-strain data for vulcanised rubber under various types of deformation; Transactions of the Faraday Society, 1943

23.  Ning, X., Zhu, Q., Lanir, Y., Margulies, S. S.: A Transversely Isotropic Viscoelastic Constitutive Equation for Brainstem Undergoing Finite Deformation. Journal of Biomechanical Engineering, 2006. 128: pg 925-933. DOI: 10.1115/1.2354208

24.  Spencer, A. J. M.: Deformations of fibre-reinforced materials, 1972

25.  Nekliudova, E. A., Semenov, A. S., Melnikov, B. E., Semenov, S. G.: Experimental research and finite element analysis of elastic and strength properties of fiberglass composite material. Magazine of Civil Engineering 3 (47), 2014

26.  Otero F., Oller S., Martinez X., Salomón O.: Numerical homogenization for composite materials analysis. Comparison with other micro mechanical formulations. Composite Structures. 2015 Apr 1;122:405-16.

27.  Charalambakis, N.: Homogenization techniques and micromechanics: A survey and perspectives. 2010: 030803.

28.  Pascon, J. P. (2019). Large deformation analysis of plane-stress hyperelastic problems via triangular membrane finite elements. International Journal of Advanced Structural Engineering, vol 3, pp. 331-350, https://doi.org/10.1007/s40091-019-00234-w.

29.  Yosibash, Z. Weiss, D. Hartmann, S. (2014). High-order FEMs for thermo-hyperelasticity at finite strains. Computers and Mathematics with Applications, vol 67, pp. 477-496, https://doi.org/10.1016/j.camwa.2013.11.003

30.  Wriggers, P. Nonlinear finite element methods; Springer 2008

31.  Gasser, T. C. Vascular Biomechanics: Concepts, Models and applications; Springer 2021 https://doi.org/10.1007/978-3-030-70966-2

Related Tech Notes

▪
• Solid Mechanics Verification Tests
▪
• Solid Mechanics Model Collection
▪
• Hyperelasticity
▪
• PDEModels Overview