Deformation

In any structural problem the analyst describes the initial configuration of the structure and is interested in its deformation throughout the history of loading. The material particle initially located at some position $\mathbf{X}$ in space will move to a new position $\mathbf{x}$: since we assume material cannot appear or disappear, there will be a one-to-one correspondence between $\mathbf{x}$ and $\mathbf{X}$, so we can always write the history of the location of a particle as

and this relationship can be inverted—we know $\mathbf{X}$ when we know $\mathbf{x}$ and t. Now consider two neighboring particles, located at $\mathbf{X}$ and at $\mathbf{X}+d\mathbf{X}$ in the initial configuration. In the current configuration we must have

using the “mapping” Equation 1.

The matrix

is called the deformation gradient matrix, and Equation 2 is written

As the material behavior depends on the straining of the material and not on its rigid body motion, those parts of the motion in the vicinity of a material point must be distinguished. Looking at an infinitesimal gauge length $d\mathbf{X}$ emanating from the particle initially at $\mathbf{X}$, we can measure its initial and current lengths as

so the “stretch ratio” of this gauge length is

If $\lambda =1$, there is no strain of this infinitesimal gauge length—it has undergone rigid body motion only. Now using Equation 4,

so that, from Equation 5,

where $\mathbf{N}$ is a unit vector in the direction of the gauge length $d\mathbf{X}$.

Equation 6 shows how to measure the stretch ratio associated with any direction, $\mathbf{N}$, at any material point defined by $\mathbf{X}$ (or by $\mathbf{x}$). Useful results are obtained when we vary the direction defined by $\mathbf{N}$ at a particular material point and look for stationary values of the stretch ratio, $\lambda$. Since $\mathbf{N}$ must always be a unit vector, stationary values of ${\lambda }^2$ are obtained by solving the constrained variational equation

where $p$ is a Lagrange multiplier, introduced to retain the constraint

Taking the variation gives back the constraint (conjugate to $\delta p$) and, conjugate to $\delta \mathbf{N}$, gives

Taking the dot product of the left-hand side of this equation with $\mathbf{N}$ and comparing with Equation 6 identifies $p={\lambda }^2$, so Equation 7 is

This problem is an eigenvalue one that can be solved for the three extreme values of ${\lambda }^2$. Since $\lambda$ is always real and positive (and nonzero), ${\lambda }^2>0$, and hence ${\mathbf{F}}^T\cdot \mathbf{F}$ must be positive definite. Equation 8 thus gives three real, positive eigenvalues, ${\lambda }_I$, ${\lambda }_{II}$, ${\lambda }_{III}$, the “principal stretches,” with three corresponding eigenvectors, ${\mathbf{N}}_I$, ${\mathbf{N}}_{II}$, ${\mathbf{N}}_{III}$, which will be orthogonal, by Equation 8, if the corresponding eigenvalues are different and can be orthogonalized otherwise. The ${\mathbf{N}}_I$ are the principal directions of strain.

Now let ${\mathbf{n}}_I$, ${\mathbf{n}}_{II}$, ${\mathbf{n}}_{III}$ be unit vectors corresponding to ${\mathbf{N}}_I$, ${\mathbf{N}}_{II}$, ${\mathbf{N}}_{III}$, but in the current configuration, so that, using Equation 4,

Then

by the orthogonality results just mentioned. Thus, ${\mathbf{n}}_I$, ${\mathbf{n}}_{II}$, and ${\mathbf{n}}_{III}$ will also be an orthogonal set. Since each is a unit vector,

where $\mathbf{R}$ is the same pure rigid body rotation matrix in each of these equations. A pure rigid body motion matrix has the property that its inverse is its transpose: ${\mathbf{R}}^T={\mathbf{R}}^{-1}$. Comparing the principal stretch directions in the current and original configurations, therefore, isolates the rigid body rotation and the stretch. Finding the principal stretch ratios and their directions thus provides one solution to the problem of isolating straining motion and rigid body motion in the vicinity of a material point.

Now consider a gauge length in the reference configuration, d${\mathbf{X}}_I$, directed along ${\mathbf{N}}_I$. The same infinitesimal material line in the current configuration will be along ${\mathbf{n}}_I$ and will be stretched by ${\lambda }_I$, so that

Similarly, along the other principal directions,

and

Since (${\mathbf{N}}_I$, ${\mathbf{N}}_{II}$, ${\mathbf{N}}_{III}$) is an orthonormal set of base vectors in the reference configuration, any infinitesimal material line (gauge length) $d\mathbf{X}$ at $\mathbf{X}$ can be written in terms of its components in this basis:

where

Each of the vectors $d{\mathbf{X}}_I$ moves and stretches to the corresponding $d{\mathbf{x}}_I$, as defined above. Thus, the current gauge length, $d\mathbf{x}$, is

which we write as

where

is the “left stretch” matrix, which is the sum of three dyadic products.

Comparison with the definition of the deformation gradient, Equation 4, shows that

which is the polar decomposition theorem—that any motion can be represented as a pure rigid body rotation, followed by a pure stretch of three orthogonal directions. The polar decomposition theorem is important because it allows us to distinguish the straining part of the motion from the rigid body rotation. Specifically, $\mathbf{F}$ completely defines the relative motions of material particles in the infinitesimal neighborhood of the material particle that was at $\mathbf{X}$ in the reference configuration; and the left stretch matrix, $\mathbf{V}$, completely defines the deformation of the material particles at $\mathbf{X}$. The rotation matrix $\mathbf{R}$ defines the rigid body rotation of the principal directions of strain (${\mathbf{N}}_I$ in the reference configuration; ${\mathbf{n}}_I$ in the current configuration). $\mathbf{R}$ represents only the rigid body rotation of the material at the point under consideration in some average sense: in a general motion, each infinitesimal gauge length emanating from a material particle has a different amount of rotation. This distinction between the rotation of the principal directions of strain, $\mathbf{R}$, and the rotations of individual directions in the material becomes significant when we must discuss large deformations of nonisotropic materials. Nevertheless, we have established an important result: if $\mathbf{F}=\mathbf{R}$ only, we know there is no deformation of the material in the immediate neighborhood of the point originally at $\mathbf{X}$ and currently at $\mathbf{x}$, since in this case $\mathbf{V}=\mathbf{I}$ and so ${\lambda }_I={\lambda }_{II}={\lambda }_{III}=1$. $\mathbf{V}-\mathbf{I}$ must be nonzero for there to be any deformation of the material at the point in question: in this sense $\mathbf{V}-\mathbf{I}$ (and, hence, $\mathbf{V}$ itself) is sufficient to define the deforming part of the motion (it contains complete information about all except pure rigid body rotation of the point). For this reason—so that, later in the development, we will be able to link the kinematics to the stressing of the material—we will need to be able to isolate $\mathbf{V}$ from $\mathbf{F}$. It is easy to obtain $\mathbf{V}\cdot \mathbf{V}$, for

since ${\mathbf{R}}^T={\mathbf{R}}^{-1}$ and $\mathbf{V}$ is symmetric.

Since we originally defined $\mathbf{V}$ from the principal stretches and their principal directions in the current configuration as

then

We see that the eigenvalues of $\mathbf{F}\cdot {\mathbf{F}}^T$, are $\lambda _I{}^2$, $\lambda _{II}{}^2$, and $\lambda _{III}{}^2$, and the corresponding eigenvectors are ${\mathbf{n}}_I$, ${\mathbf{n}}_{II}$, and ${\mathbf{n}}_{III}$. We can then construct $\mathbf{V}$. The deformation at the point is, thus, readily obtained by multiplying a $3\times 3$ matrix with its transpose ($\mathbf{F}\cdot {\mathbf{F}}^T$) and solving the real eigenproblem for the resulting (symmetric) matrix. We can then obtain the rotation $\mathbf{R}$ as

Since $\mathbf{V}$ has been constructed from its eigenvalues and eigenvectors, its inverse is immediately available:

So far we have written the results quite generally, without reference to any particular coordinate system. To perform computations we must choose a basis system to express these results as arrays of individual numbers. We now do so with some generality with respect to the choice of basis system. The justification for retaining generality at this stage is twofold: as an exercise, to provide a little more familiarity in the notation system we have chosen to use in this guide, and because we do need some—but, as it turns out, not all—of the generality when we have to deal with shell elements, where it is undesirable to use the rectangular Cartesian base vectors of the global, spatial system because the natural orientation of the shell reference surface causes us to prefer to choose two of the base vectors to be tangent to the shell's reference surface and the other to be normal to this surface. This preference causes us to need two basis systems: one associated with the body in its current configuration, when the point in question is at $\mathbf{x}$, and one associated with the body in its reference configuration, when the same point was at $\mathbf{X}$, because the orientation of the shell's reference surface—which determines our choice of basis vectors—will be quite different in these two configurations. We will write ${\mathbf{e}}_{\alpha }$, $\alpha =1,2,3$, as the basis vectors chosen to write components associated with the current configuration (so that any vector $\mathbf{a}$ associated with the current configuration is written as $\mathbf{a}=a^{\alpha }{\mathbf{e}}_{\alpha }$) and ${\mathbf{E}}_{\alpha }$, $\alpha =1,2,3$, as the basis at the same material point but in the reference configuration. (Since we assume that both of these basis systems are adequate to express any vector-valued function by its components in the basis system—that is, the basis vectors are not linearly dependent—either would, in principal, serve for both configurations. We introduce two distinct systems by preference, because each is chosen as particularly suitable for a particular configuration.) Since we do not yet impose any particular restrictions on the ${\mathbf{e}}_{\alpha }$ or the ${\mathbf{E}}_{\alpha }$ (except for the requirement that the vectors must not be linearly dependent), we cannot assume that they will be orthogonal or of unit length: we will, therefore, need to use the corresponding contravariant vectors defined by

and the contravariant metric tensors

We can express the deformation gradient, $\mathbf{F}$, numerically by projecting it onto the bases:

Recall the definition of $\mathbf{F}$:

Since the components of $d\mathbf{x}$ along ${\mathbf{e}}_{\alpha }$ are $d{x}_{\alpha }=d\mathbf{x}\cdot {\mathbf{e}}_{\alpha }$ and we can write $d\mathbf{X}=d{X}^{\beta }{\mathbf{E}}_{\beta }$,

Thus, writing $d{x}_{\alpha }=F_{\alpha \beta }d{X}^{\beta }$ defines

We must continue to bear in mind that the first index of $F_{\alpha \beta }$ is associated with a component of $\mathbf{F}$ along a base vector in the current configuration (${\mathbf{e}}_{\alpha }$ in this case), while its second index is associated with a component of $\mathbf{F}$ along a base vector in the reference configuration (${\mathbf{E}}_{\beta }$).

From Equation 13 we can write

where $G^{\alpha \beta }$ is the contravariant metric of the basis system that we have chosen in the reference configuration.

The eigenproblem for the squared principal stretch ratios and their directions is solved by finding the eigenvalues of the matrix of numbers $F_{\alpha \beta }{G}^{\beta \gamma }{F}_{\delta \gamma }$. The eigenvectors will appear as the components ${\left(n_I\right)}_{\alpha }$ along the ${\mathbf{e}}^{\alpha }$ base vectors in the current configuration. Since we have defined the left stretch on the current configuration as

we will write its components on the basis in the current configuration as

and, since

The polar decomposition gives

so

where $g^{\gamma \delta }$ is the contravariant metric tensor of the basis system we have chosen to use in the current configuration and—as with $F_{\alpha \beta }$—we see that the first index of $R_{\alpha \beta }$ is associated with the contravariant base vector ${\mathbf{e}}^{\alpha }$ in the current configuration, while the second index is associated with the contravariant base vector ${\mathbf{E}}^{\beta }$ in the reference configuration.

We should take care to understand the distinction between the direct matrix expression of the rigid body rotation of the principal directions of strain of the material, $\mathbf{R}$, and the components of $\mathbf{R}$ expressed on a particular basis. Suppose, for example, that the rigid body rotation at a point is zero (that is, $\mathbf{R}=\mathbf{I}$) but we, nevertheless, have chosen different basis systems ${\mathbf{e}}_{\alpha }$ and ${\mathbf{E}}_{\alpha }$. In this case $R_{\alpha \beta }={\mathbf{e}}_{\alpha }\cdot \mathbf{R}\cdot {\mathbf{E}}_{\beta }={\mathbf{e}}_{\alpha }\cdot \mathbf{I}\cdot {\mathbf{E}}_{\beta }={\mathbf{e}}_{\alpha }\cdot {\mathbf{E}}_{\beta }$. This implies that, even though $\mathbf{R}$ is a unit matrix (in the sense that operating on any vector with this matrix makes no change in that vector), the numerical values we have chosen to store the matrix—the $R_{\alpha \beta }$—do not form a unit matrix of numbers unless the ${\mathbf{e}}_{\alpha }$ and the ${\mathbf{E}}_{\alpha }$ are coincident and orthonormal. Thus, our choice of quite general basis systems that are not the same in the current and reference configurations (introduced as being “natural” for writing results for shells) somewhat complicates the interpretation of the numbers we store.

In the previous few paragraphs we have chosen to explore the expression of the basic results we have derived so far for the kinematics of the total motion in terms of quite general basis systems, ${\mathbf{e}}_{\alpha }$ and ${\mathbf{E}}_{\alpha }$. In Abaqus we wish to express results as simply and directly as possible, and we can do so by choosing particular sets of basis vectors that offer the most convenience for our purposes. First, we take the ${\mathbf{e}}_{\alpha }$ (and, by extension, the ${\mathbf{E}}_{\alpha }$, since these are just the ${\mathbf{e}}_{\alpha }$ at the beginning of the motion) to be a local, orthonormal system at each point. Although it is not possible to construct a Cartesian system with orthonormal base vectors over a general shell surface, we can always project the general results onto such a system when that system is chosen specifically at each point where we need to make the projection—typically at the integration points of the elements. The choice of which system is used as this local orthonormal basis is made in Abaqus at two levels: we distinguish continuum (“solid”) elements from structural (shell and beam) elements, and we distinguish the default choice of directions from the particular choice of directions (orientation) specified by the user. For continuum elements the default ${\mathbf{E}}_{\alpha }$ are unit vectors along the axes of the global Cartesian system chosen for the problem. At points where the orientation is defined by the user, the specified ${\mathbf{E}}_{\alpha }$ are used. For shells (and membranes) we take ${\mathbf{E}}_1$ and ${\mathbf{E}}_2$ tangent to the shell's reference surface and ${\mathbf{E}}_3$ normal to that surface at the point under consideration. By default, ${\mathbf{E}}_1$ is the projection of the global x-axis onto the reference surface or, if the global x-axis is almost normal to that surface at the point, ${\mathbf{E}}_1$ is the projection of the global z-axis onto the surface. If the orientation is defined by the user, ${\mathbf{E}}_1$ and ${\mathbf{E}}_2$ are the projections of the two specified axes onto the reference surface at the point. In all cases ${\mathbf{E}}_3$ is normal to the shell's reference surface. For beams ${\mathbf{E}}_1$ is along the beam axis, with ${\mathbf{E}}_2$ and ${\mathbf{E}}_3$ defined from the beam section definition option and beam normals given as part of the nodal coordinate definition. For continuum elements the same schemes are applied by default to define the basis system in the current configuration. For continuum elements with the orientation specified by the user and in all cases for shells, beams, and membranes, the ${\mathbf{e}}_{\alpha }$ are defined by

These schemes all have the same property: at any point in time the ${\mathbf{e}}_{\alpha }$ are orthonormal vectors: ${\mathbf{e}}_{\alpha }\cdot {\mathbf{e}}_{\beta }={\delta }_{\alpha \beta }$, so ${\mathbf{e}}^{\alpha }={\mathbf{e}}_{\alpha }$ and, thus, $g_{\alpha \beta }=g^{\alpha \beta }=g_{\alpha }^{\beta }={\delta }_{\alpha \beta }$, and—in particular—${\mathbf{E}}_{\alpha }\cdot {\mathbf{E}}_{\beta }={\delta }_{\alpha \beta }$ and, thus, $G^{\alpha \beta }={\delta }^{\alpha \beta }$. This simplifies the understanding of all the quantities we write, since the components of any tensor $T_{\alpha \beta \mathrm{\dots }}$ are always the physical projections of that tensor-valued quantity on the local orthogonal basis system ${\mathbf{e}}_{\alpha }$ and we need not distinguish covariant and contravariant components as we did in the general development above. In practical terms the only price we must pay for this simplicity is in shells when we have to use a separate basis system at each point under study, since we cannot construct a single system with the orthonormal property on a general curved surface. (In an axisymmetric system we also have to use $d{x}^3=rd\theta$ to ensure that the ${\mathbf{e}}_3$ base vector is a unit vector, but this is a minor point.) The simplifications are valuable and, from our perspective of studying finite element formulations, they are bought at modest cost, since we generally only consider a single integration point at a time. Throughout the rest of this guide, whenever we need to write down particular components of a tensor, we shall assume that the basis on which they are written has the orthonormal property ${\mathbf{e}}_{\alpha }\cdot {\mathbf{e}}_{\beta }={\delta }_{\alpha \beta }$.

The material also undergoes rigid body translation, but this is not important in the development since we need consider only relative motion of neighboring points because we are interested in the deformation of the material to link the kinematics of the motion to the material's constitutive behavior. Numerically, rigid body translation is significant only for two reasons. One is that the spatial discretization must allow rigid body translation without giving strain, which is important in choosing interpolation functions for the finite elements. The other is that care must be exercised to ensure that the strain and rotation are calculated accurately when the rigid body motion is large, since then the strain and rotation depend on the difference between two very large motions.