MultiPhysicsVault/.raw/FiniteElementProcedures/FiniteElementProcedures_031.md

area

Metric	Value
(σ₁)ₘₐₓ	0.4940
δ	1.625

heatmap

(σ₁)max	Value
(σ₁)max	0.6000
δ	1.680

(a) Bands of maximum principal stress. Case of Poisson's ratio \nu = 0.30 . Sixteen and 16 \times 64 element mesh results
Figure 4.21 Analysis of cantilever bracket in plane strain conditions. Bracket is shown in Fig. 4.20(a). Same meshes as in Fig. 4.20 are used but with the nine-node mixed interpolated element (the 9/3 element). Compare the results shown with those given in Fig. 4.20.

ables, we need another equation to connect these two solution variables. This equation is provided by (4.131) written in integral form (see Example 4.31).

\boxed {\int_ {V} \left(\frac {p}{\kappa} + \epsilon_ {V}\right) \overline {{p}} d V = 0} \tag {4.143}

These basic equations can also be derived more formally from variational principles (see L. R. Herrmann [A] and S. W. Key [A]). We derive the basic equations in the following example from the Hu-Washizu functional.

area

Metric	Value
(σ₁)ₘₐₓ	0.4983
δ	1.349

surface_3d

Metric	Value
(σ₁)max	0.5998
δ	1.393

(b) Bands of maximum principal stress. Case of Poisson's ratio \nu = 0.499 . Sixteen and 16 \times 64 element mesh results
Figure 4.21 (continued)

EXAMPLE 4.31: Derive the u/p formulation from the Hu-Washizu variational principle.

The derivation is quite analogous to the presentation in Example 4.30 where we considered a mixed interpolation for a beam element.

We start by letting \tau = \mathbf{C}\epsilon in (4.114) to obtain the Hellinger-Reissner functional,

\Pi_ {\mathrm{HR}} ^ {*} (\mathbf {u}, \epsilon) = - \int_ {V} \frac {1}{2} \epsilon^ {T} \mathbf {C} \epsilon d V + \int_ {V} \epsilon^ {T} \mathbf {C} \partial_ {\epsilon} \mathbf {u} d V - \int_ {V} \mathbf {u} ^ {T} \mathbf {f} ^ {B} d V - \int_ {S _ {f}} \mathbf {u} ^ {S _ {f} T} \mathbf {f} ^ {S _ {f}} d S \tag {a}

where we assume that the displacement boundary conditions are satisfied exactly (hence, also the displacement variations will be zero on the surface of prescribed displacements).

Next we establish the deviatoric and volumetric contributions and postulate that the deviatoric contribution will be evaluated from the displacements. Hence, we can specialize (a) into

\tilde {\Pi} _ {\mathrm{HR}} ^ {*} (\mathbf {u}, p) = \int_ {V} \frac {1}{2} \epsilon^ {\prime T} \mathbf {C} ^ {\prime} \epsilon^ {\prime} d V - \int_ {V} \frac {1}{2} \frac {p ^ {2}}{\kappa} d V - \int_ {V} p \epsilon_ {V} d V - \int_ {V} \mathbf {u} ^ {T} \mathbf {f} ^ {B} d V - \int_ {S _ {f}} \mathbf {u} ^ {S _ {f} T} \mathbf {f} ^ {S _ {f}} d S (\mathrm{b})

where the prime denotes deviatoric quantities, \epsilon_{V} is the volumetric strain evaluated from the displacements, p is the pressure, and \kappa is the bulk modulus. Note that whereas in (a) the independent variables are u and \epsilon , in (b) the independent variables are u and p.

Invoking the stationarity of \Pi_{HR}^{*} with respect to the displacements and the pressure, we obtain

\int_ {V} \delta \epsilon^ {\prime T} \mathbf {C} ^ {\prime} \epsilon^ {\prime} d V - \int_ {V} p \delta \epsilon_ {V} d V = \mathcal {R}

and \int_{V}\left(\frac{p}{\kappa} +\epsilon_{V}\right)\delta p dV = 0

where R corresponds to the virtual work of the externally applied loading [see (4.7)].

It is interesting to note that we may also think of (b) as the total potential in terms of the displacements and the pressure plus a Lagrange multiplier term that enforces the constraint between the volumetric strains and the pressure,

\tilde {\Pi} _ {\mathrm{HR}} ^ {*} = - \int_ {V} \frac {1}{2} \epsilon^ {\prime T} \mathbf {C} ^ {\prime} \epsilon^ {\prime} d V + \int_ {V} \frac {1}{2} \frac {p ^ {2}}{\kappa} d V - \int_ {V} \mathbf {u} ^ {T} \mathbf {f} ^ {B} d V \tag {c}

- \int_ {S _ {f}} \mathbf {u} ^ {S _ {f} T} \mathbf {f} ^ {S _ {f}} d S - \int_ {V} \lambda \left(\epsilon_ {V} + \frac {p}{\kappa}\right) d V

In (c) the last integral represents the Lagrange multiplier constraint, and we find \lambda = p .

To arrive at the governing finite element equations, we can now use (4.140) and (4.143) as in Section 4.2.1, but in addition to interpolating the displacements we also interpolate the pressure p. The discussion in Section 4.2.1 showed that we need to consider the formulation of only a single element; the matrices of an assemblage of elements are then formed in a standard manner.

Using, as in Section 4.2.1,

\mathbf {u} = \mathbf {H} \hat {\mathbf {u}} \tag {4.144}

we can calculate

\epsilon^ {\prime} = \mathbf {B} _ {D} \hat {\mathbf {u}}; \quad \epsilon_ {V} = \mathbf {B} _ {V} \hat {\mathbf {u}} \tag {4.145}

The additional interpolation assumption is

p = \mathbf {H} _ {p} \hat {\mathbf {p}} \tag {4.146}

where the vector \hat{\mathbf{p}} lists the pressure variables [(see the discussion following (4.148)].

Substituting from (4.144) to (4.146) into (4.140) and (4.143), we obtain

\left[ \begin{array}{l l} \mathbf {K} _ {u u} & \mathbf {K} _ {u p} \\ \mathbf {K} _ {p u} & \mathbf {K} _ {p p} \end{array} \right] \left[ \begin{array}{l} \hat {\mathbf {u}} \\ \hat {\mathbf {p}} \end{array} \right] = \left[ \begin{array}{l} \mathbf {R} \\ \mathbf {0} \end{array} \right] \tag {4.147}

where \mathbf{K}_{uu} = \int_{V}\mathbf{B}_{D}^{T}\mathbf{C}^{\prime}\mathbf{B}_{D}dV

\mathbf {K} _ {u p} = \mathbf {K} _ {p u} ^ {T} = - \int_ {V} \mathbf {B} _ {V} ^ {T} \mathbf {H} _ {p} d V \tag {4.148}

\mathbf {K} _ {p p} = - \int_ {V} \mathbf {H} _ {p} ^ {T} \frac {1}{\kappa} \mathbf {H} _ {p} d V

and C' is the stress-strain matrix for the deviatoric stress and strain components.

The relations in (4.144) to (4.148) give the basic equations for formulating elements with displacements and pressure as variables. The key question for the formulation is now, What pressure and displacement interpolations should be used to arrive at effective elements? For example, if the pressure interpolation is of too high a degree compared to the displacement interpolation, the element may again behave as a displacement-based element and not be effective.

Considering for the moment only the pressure interpolation, the following two main possibilities exist and we label them differently.

The u/p formulation. In this formulation, the pressure variables pertain only to the specific element being considered. In the analysis of almost incompressible media (as so far discussed), the element pressure variables can be statically condensed out prior to the element assemblage. Continuity of pressure is not enforced between elements but will be a result of the finite element solution if the mesh used is fine enough.

The u/p-c formulation. The letter "c" denotes continuity in pressure.

In this formulation, the element pressure is defined by nodal pressure variables that pertain to adjacent elements in the assemblage. The pressure variables therefore cannot be statically condensed out prior to the element assemblage. Continuity of pressure between elements is directly enforced and will therefore always be a result of the solution irrespective of whether the mesh used is fine or coarse.

Consider the following two elements, one corresponding to each of the formulations.

EXAMPLE 4.32: For the four-node plane strain element shown, assume that the displacements are interpolated using the four nodes and assume a constant pressure. Evaluate the matrix expressions used for the u/p formulation.

text_image

y, v 2 1 2 x, u Young's modulus E Poisson's ratio v 3 4 2

Figure E4.32 A 4/1 element

This element is referred to as the u / p 4/1 element. In plane strain analysis we have

\boldsymbol {\epsilon} ^ {\prime} = \left[ \begin{array}{c} \boldsymbol {\epsilon} _ {x x} - \frac {1}{3} \left(\boldsymbol {\epsilon} _ {x x} + \boldsymbol {\epsilon} _ {y y}\right) \\ \boldsymbol {\epsilon} _ {y y} - \frac {1}{3} \left(\boldsymbol {\epsilon} _ {x x} + \boldsymbol {\epsilon} _ {y y}\right) \\ \gamma_ {x y} \\ - \frac {1}{3} \left(\boldsymbol {\epsilon} _ {x x} + \boldsymbol {\epsilon} _ {y y}\right) \end{array} \right] = \left[ \begin{array}{c} \frac {2}{3} \frac {\partial u}{\partial x} - \frac {1}{3} \frac {\partial v}{\partial y} \\ \frac {2}{3} \frac {\partial v}{\partial y} - \frac {1}{3} \frac {\partial u}{\partial x} \\ \frac {\partial u}{\partial y} + \frac {\partial v}{\partial x} \\ - \frac {1}{3} \left(\frac {\partial u}{\partial x} + \frac {\partial v}{\partial y}\right) \end{array} \right]; \quad \boldsymbol {\epsilon} _ {V} = \boldsymbol {\epsilon} _ {x x} + \boldsymbol {\epsilon} _ {y y} \tag {a}

and \mathbf{S} = \mathbf{C}'\epsilon' , where

\mathbf {C} ^ {\prime} = \left[ \begin{array}{c c c c} 2 G & 0 & 0 & 0 \\ 0 & 2 G & 0 & 0 \\ 0 & 0 & G & 0 \\ 0 & 0 & 0 & 2 G \end{array} \right]; \qquad G = \frac {E}{2 (1 + \nu)}

The displacement interpolation is as in Example 4.6,

\mathbf {u} = \mathbf {H} \hat {\mathbf {u}}

with \mathbf{u}(x,y) = \begin{bmatrix} u(x,y)\\ v(x,y) \end{bmatrix};\qquad \hat{\mathbf{u}}^T = \left[u_1\quad u_2\quad u_3\quad u_4\quad \vdots \quad v_1\quad v_2\quad v_3\quad v_4\right]

\mathbf {H} = \left[ \begin{array}{c c c c c c c c c} h _ {1} & h _ {2} & h _ {3} & h _ {4} & \vdots & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \vdots & h _ {1} & h _ {2} & h _ {3} & h _ {4} \end{array} \right] \tag {b}

h _ {1} = \frac {1}{4} (1 + x) (1 + y); \quad h _ {2} = \frac {1}{4} (1 - x) (1 + y)

h _ {3} = \frac {1}{4} (1 - x) (1 - y); \quad h _ {4} = \frac {1}{4} (1 + x) (1 - y)

Using (a) and (b), the strain-displacement interpolation matrices are

\mathbf {B} _ {D} = \left[ \begin{array}{c c c c c c c} \frac {2}{3} h _ {1, x} & \frac {2}{3} h _ {2, x} & \dots & \vdots & - \frac {1}{3} h _ {1, y} & - \frac {1}{3} h _ {2, y} & \dots \\ - \frac {1}{3} h _ {1, x} & - \frac {1}{3} h _ {2, x} & \dots & \vdots & \frac {2}{3} h _ {1, y} & \frac {2}{3} h _ {2, y} & \dots \\ h _ {1, y} & h _ {2, y} & \dots & \vdots & h _ {1, x} & h _ {2, x} & \dots \\ - \frac {1}{3} h _ {1, x} & - \frac {1}{3} h _ {2, x} & \dots & \vdots & - \frac {1}{3} h _ {1, y} & - \frac {1}{3} h _ {2, y} & \dots \end{array} \right]

and \mathbf{B}_V = [h_{1,x}\quad h_{2,x}\quad \ldots \quad \vdots \quad h_{1,y}\quad h_{2,y}\quad \ldots ]

For a constant pressure assumption we have

\mathbf {H} _ {p} = [ 1 ]; \quad \hat {\mathbf {p}} = [ p _ {0} ]

Since the degree of freedom \hat{p} = [p_{0}] pertains only to this element and not to the adjacent elements, we can use static condensation to obtain from (4.147) the element stiffness matrix corresponding to the nodal point displacement degrees of freedom only;

\mathbf {K} = \mathbf {K} _ {u u} - \mathbf {K} _ {u p} \mathbf {K} _ {p p} ^ {- 1} \mathbf {K} _ {p u}

The element is further discussed in Example 4.38.

EXAMPLE 4.33: Consider the nine-node plane strain element shown in Fig. E4.33. Assume that the displacements are interpolated using the nine nodes and that the pressure is interpolated using only the four corner nodes. Refer to the information given in Example 4.32 and discuss the additional considerations for the evaluation of the matrix expressions of this element.

text_image

y, v 2 5 1 2 6 9 8 x, u 3 7 4 2

• Displacement node
  ● Displacement and pressure node

Young's modulus E Poisson's ratio v

Figure E4.33 A 9/4-c element

This element was proposed by P. Hood and C. Taylor [A]. In the formulation the nodal pressures pertain to adjacent elements, and according to the above element nomination we refer to it as a u/p-c element (it is the 9/4-c element).

The deviatoric and volumetric strains are as given in (a) in Example 4.32. The displacement interpolation corresponds to the nine nodes of the element,

\left[ \begin{array}{l} u (x, y) \\ v (x, y) \end{array} \right] = \left[ \begin{array}{c c c c c c c} h _ {1} ^ {*} & \dots & h _ {9} ^ {*} & \vdots & 0 & \dots & 0 \\ 0 & \dots & 0 & \vdots & h _ {1} ^ {*} & \dots & h _ {9} ^ {*} \end{array} \right] \left[ \begin{array}{c} u _ {1} \\ \vdots \\ \vdots \\ u _ {9} \\ \dots \\ v _ {1} \\ \vdots \\ v _ {9} \end{array} \right] \tag {a}

where the interpolation functions h_i^* are constructed as explained in Section 4.2.3 (or see Section 5.3 and Fig. 5.4).

The deviatoric and volumetric strain-displacement matrices are obtained as in Example 4.32. The pressure interpolation is given by

p = \left[ \begin{array}{l l l l} h _ {1} & h _ {2} & h _ {3} & h _ {4} \end{array} \right] \left[ \begin{array}{l} p _ {1} \\ p _ {2} \\ p _ {3} \\ p _ {4} \end{array} \right]

where the h_i are those given in (b) in Example 4.32.

A main computational difference between this element and the four-node element discussed in Example 4.32 is that the pressure degrees of freedom cannot be statically condensed out on the element level because the variables p_{1}, \ldots, p_{4} pertain to the element we are considering here and to the adjacent elements, thus describing a continuous pressure field for the discretization.

Let us now return to the discussion of what pressure and displacement interpolations should be used in order to have an effective element.

For instance, in Example 4.32, we used four nodes to interpolate the displacements and assumed a constant pressure, and we may ask whether a constant pressure is the appropriate choice for the four-node element. Actually, for this element, it is a somewhat natural choice because the volumetric strain calculated from the displacements contains linear variations in x and y and our pressure assumption should be of lower order.

When higher-order displacement interpolations are used, the choice of the appropriate pressure interpolation is not obvious and indeed much more difficult: the pressure should not be interpolated at too low a degree because then the pressure prediction could be of higher order and hence be more accurate, but the pressure should also not be interpolated at too high a degree because then the element would behave like the displacement-based elements and lock. Hence, we want to use the highest degree of pressure interpolation that does not introduce locking into the element.

For example, considering the u/p formulation and biquadratic displacement interpolation (i.e., nine nodes for the description of the displacements), we may naturally try the following cases in plane strain analysis:

1. Constant pressure, p = p_0 (9/1 element)
2. Linear pressure, p = p_{0} + p_{1}x + p_{2}y (9/3 element)
3. Bilinear pressure, p = p_0 + p_1x + p_2y + p_3xy (9/4 element)

and so on, up to a quadratic pressure interpolation (which corresponds to the 9/9 element).

These elements have been analyzed theoretically and by use of numerical experiments. Studies of the elements show that the 9/1 element does not lock, but the rate of convergence of pressure (and hence stresses) as the mesh is refined is only of o(h) because a constant pressure is assumed in each nine-node element. The poor quality of the pressure prediction can of course also have a negative effect on the prediction of the displacements.

Studies also show that the 9/3 element is most attractive because it does not lock and the stress convergence is of o(h^{2}) . Hence, the predictive capability is optimal since if a biquadratic displacement expansion is used, no higher-order convergence in stresses can be expected. Also, the 9/3 element is effective for any Poisson's ratio up to 0.5 (but the static condensation of the pressure degrees of freedom is possible only for values of \nu < 0.5 ).

Hence, we may be tempted to always use the 9/3 element (instead of the displacement-based nine-node element). However, we find in practice that the 9/3 element is computationally slightly more expensive than the nine-node displacement-based element, and when \nu is less than 0.48, the additional terms in the pressure expansion of the displacement-based element allow a slightly better prediction of stresses.

The next u/p element of interest is the 9/4 element, and studies show that this element locks when v is close to 0.50; hence it cannot be recommended for almost incompressible analysis, in general.

In an analogous manner, other u/p elements can be constructed, and Table 4.6 summarizes some choices. Regarding these elements, we may note that the four-node two-dimensional and eight-node three-dimensional elements are extensively used in practice. However, the nine-node two-dimensional and 27-node three-dimensional elements are frequently more powerful.

As indicated in Table 4.6, the Q_{2} - P_{1} and P_{2}^{+} - P_{1} elements are the first members of two families of elements that may be used. That is, the quadrilateral elements Q_{n} - P_{n-1} , and the triangular elements P_{n}^{+} - P_{n-1}, n > 2 , are also effective elements.

In Table 4.6 we refer to the inf-sup condition, which we will discuss in Section 4.5.

From a computational point of view, the u/p elements are attractive because the element pressure degrees of freedom can be statically condensed out before the elements are assembled (assuming v < 0.5 but possibly very close to 0.5). Hence, the degrees of freedom for the assemblage of elements are only the same nodal point displacements that are also the degrees of freedom in the pure displacement-based solution.

However, the u/p-c formulation has the advantage that a continuous pressure field is always calculated. Table 4.7 lists some effective elements.

The Finite Element Solution of Totally Incompressible Conditions

If we want to consider the material to be totally incompressible, we can still use (4.140) and (4.143), but we then let \kappa \to \infty . For this reason, we refer to this case as the limit problem. Then (4.143) becomes

\int_ {V} \epsilon_ {v} \overline {{p}} d V = 0 \tag {4.149}

and (4.147) becomes, correspondingly,

\left[ \begin{array}{c c} \mathbf {K} _ {u u} & \mathbf {K} _ {u p} \\ \mathbf {K} _ {p u} & \mathbf {0} \end{array} \right] \left[ \begin{array}{l} \hat {\mathbf {u}} \\ \hat {\mathbf {p}} \end{array} \right] = \left[ \begin{array}{l} \mathbf {R} \\ \mathbf {0} \end{array} \right] \tag {4.150}

Hence, in the coefficient matrix, the diagonal elements corresponding to the pressure degrees of freedom are now zero. It follows that a static condensation of the element pressure degrees of freedom in the u/p formulation is no longer possible and that the solution of the equations of the complete assemblage of elements needs special considerations (beyond those required in the pure displacement-based solution) to avoid encountering a zero pivot element (see Section 8.2.5).

Suitable elements for solution are listed in Tables 4.6 and 4.7. These elements are effective (however — see notes in the Tables) because they have good predictive capability irrespective of how close the behavior of the medium is to a situation of total incompressibility (but the procedure for solving the governing finite element equations must take into account that the elements in K_{pp} become increasingly smaller as total incompressibility is approached).

As already noted earlier, we refer to the inf-sup condition in Tables 4.6 and 4.7. This condition is the basic mathematical criterion that determines whether a mixed finite element discretization is stable and convergent (and hence will yield a reliable solution). The condition was introduced as the fundamental test for mixed finite element formulations by I. Babuška [A] and F. Brezzi [A] and since then has been used extensively in the analysis of mixed finite element formulations. In addition to the inf-sup condition, there is also the ellipticity condition which has not received as much attention because frequently—as in the analysis of almost incompressible media—the ellipticity condition can be easily satisfied.

TABLE 4.6 Various effective u/p elements (displacements between elements are continuous and pressure variables pertain to individual elements) {}^{ \dagger }

Element	Nodal points		Remarks
	2-D element	3-D element
$Q_1 - P_0$ in 2-D: 4/1in 3-D: 8/1	$v_1$ $u_1$ Node 1 $p = p_0$	$w_1$ $v_1$ $u_1$ Node 1 $p = p_0$	The element predicts reasonably good displacements, but stresses may be inaccurate because of the constant pressure assumption and possible pressure fluctuations. The element does not satisfy the inf-sup condition (see discussion of element in Section 4.5.5).
$Q_2^r - P_0$ ; $Q_2^r = Q_2 \cap P_3$ in 2-D: 8/1in 3-D: 20/1	$p = p_0$	$p = p_0$	The element satisfies the inf-sup condition, but the constant pressure assumption may require fine discretizations for accurate stress prediction.*

†For the interpolation functions, see Figs. 4.13, 5.4, 5.5, 5.11, and 5.13. *The axisymmetric and 3-D elements do not satisfy the ellipticity condition and should be used with care.