3. Weak Forms

3.1. Introduction

A key feature of the finite element method is the form of the boundary value problem (or initial/boundary value problem in the case of dynamics) that is discretized. More specifically, the finite element method is one of a large number of variational methods that rely on the approximation of integral forms of the governing equations. In this chapter we briefly examine how such integral (alternatively, weak or variational) forms are constructed for the linear elastic system we introduced in Chapter Section 2.

3.2. Quasistatic Case

Consider the quasistatic case first, we recall (2.27) – (2.30) and explore an alternative manner in which the conditions can be stated. We consider a collection of vector-valued functions \(\mathbf{w}\), which we call weighting functions for reasons that will soon be clear. We require that these functions \(\mathbf{w} : \bar \Omega \to \mathbb{R}^3\) satisfy

(3.1)\[\mathbf{w} = \mathbf{0} \text{ on } \Gamma _{u} .\]

Furthermore it is assumed that these functions are sufficiently smooth that all necessary partial derivatives can be computed. Suppose we have the solution \(\mathbf{u}\) of (2.27) and (2.28). We can then take any smooth function \(\mathbf{w}\) satisfying (3.1) and compute its dot product with (2.27), which must produce

(3.2)\[\mathbf{w} \cdot ( \nabla \cdot \mathbf{T} + \mathbf{f} ) = 0 \text{ on } \Omega\]

at each time \(t \in (0,T)\). We can then integrate (3.2) over \(\Omega\) to obtain

(3.3)\[\int_{\Omega} \mathbf{w} \cdot ( \nabla \cdot \mathbf{T} + \mathbf{f} ) \mathrm{d} \Omega = 0 .\]

(3.3) can be manipulated further by noting that

(3.4)\[\mathbf{w} \cdot ( \nabla \cdot \mathbf{T} ) = \nabla \cdot (\mathbf{T} \mathbf{w}) - ( \nabla \mathbf{w} ): \mathbf{T}\]

(product rule of differentiation), and by also taking advantage of the divergence theorem from multivariate calculus:

(3.5)\[\int_{\Omega} \nabla \cdot ( \mathbf{T} \mathbf{w} ) \, \mathrm{d} \Omega = \int_{\partial \Omega} ( \mathbf{n} \cdot \mathbf{T} \mathbf{w} ) \, \mathrm{d} \Gamma.\]

Note that \(\mathbf{n}\) is the outward normal directed normal on \(\partial \Omega\) and \(\mathrm{d} \Gamma\) is a differential area of this surface. Use of (3.4) and (3.5) in (3.3) and rearranging gives

(3.6)\[\int_{\Omega} (\nabla \mathbf{w} ): \mathbf{T} \, \mathrm{d} \Omega = \int_{\Omega} \mathbf{w} \cdot \mathbf{f} \, \mathrm{d} \Omega + \int_{\partial \Omega} (\mathbf{n} \cdot \mathbf{T} \mathbf{w} ) \, \mathrm{d} \Gamma .\]

Now, taking advantage of the symmetry of \(\mathbf{T}\) and noting, from (2.16), that the surface traction \(\mathbf{t}\) equals \(\mathbf{Tn}\), we can write

(3.7)\[\int_{\partial \Omega} ( \mathbf{n} \cdot \mathbf{T} \mathbf{w} ) \, \mathrm{d} \Gamma = \int_{\partial \Omega} ( \mathbf{w} \cdot \mathbf{T} \mathbf{n} ) \, \mathrm{d} \Gamma = \int_{\partial \Omega} \mathbf{w} \cdot \mathbf{t} \, \mathrm{d} \Gamma.\]

We now recall the restrictions in (2.1), which tell us that \(\partial \Omega\) is the union of \(\Gamma _u\) and \(\Gamma _{\sigma}\). Since by definition \(\mathbf{w} = 0\) on \(\Gamma _u\), we can write

(3.8)\[\int_{\partial \Omega} \mathbf{w} \cdot \mathbf{t} \, \mathrm{d} \Gamma = \int_{\Gamma _u} \mathbf{w} \cdot \mathbf{t} \, \mathrm{d} \Gamma + \int_{\Gamma _{\sigma}} \mathbf{w} \cdot \mathbf{t} \, \mathrm{d} \Gamma = \int_{\Gamma _{\sigma}} \mathbf{w} \cdot \bar{\mathbf{t}} \, \mathrm{d} \Gamma\]

where the last equality incorporates the boundary condition \(\mathbf{t} = \bar{\mathbf{t}}\) on \(\Gamma _{\sigma}\).

We collect these calculations to conclude that

(3.9)\[\int_{\Omega} (\nabla \mathbf{w} ): \mathbf{T} \, \mathrm{d} \Omega = \int_{\Omega} \mathbf{w} \cdot \mathbf{f} \, \mathrm{d} \Omega + \int_{\Gamma _{\sigma}} \mathbf{w} \cdot \bar{\mathbf{t}} \, \mathrm{d} \Gamma ,\]

which must hold for the solution \(\mathbf{u}\) of (2.27) – (2.30) for any \(\mathbf{w}\) satisfying condition (3.1).

To complete our alternative statement of the boundary value problem, the concepts of solution and variational spaces need to be introduced. We define the solution space \(\mathbb{S}_t\) corresponding to time \(t\) via

(3.10)\[\mathbb{S}_t = \left\{ \mathbf{u} \,|\, \mathbf{u} = \bar{\mathbf{u}}(t) \quad \mbox{on} \, \Gamma_u , \quad \mathbf{u} \, \mbox{is smooth} \right\}\]

and the weighting space \(\mathbb{W}\) as

(3.11)\[\mathbb{W} = \left\{ \mathbf{w} \,|\, \mathbf{w} = 0 \quad \mbox{on} \, \Gamma_u , \quad \mathbf{w} \, \mbox{is smooth} \right\}.\]

With these two collections of functions in hand, we consider the following alternative statement of the boundary value problem summarized by (2.27) – (2.30):

Given the boundary conditions \(\bar{\mathbf{t}}\) on \(\Gamma_{\sigma} \times (0,T)\), \(\bar{**u**}\) on \(\Gamma_{u} \times (0,T)\) and the distributed body force \(\mathbf{f}\) on \(\Omega \times (0,T)\), find the \(\mathbf{u} \in \mathbb{S}_t\) for each time \(t \in (0,T)\) such that

(3.12)\[\int_{\Omega} (\nabla \mathbf{w} ): \mathbf{T} \, \mathrm{d} \Omega = \int_{\Omega} \mathbf{w} \cdot \mathbf{f} \, \mathrm{d} \Omega + \int_{\Gamma _{\sigma}} \mathbf{w} \cdot \bar{\mathbf{t}} \, \mathrm{d} \Gamma\]

for all \(\mathbf{w} \in \mathbb{W}\), where \(\mathbb{S}_t\) is as defined in (3.10), \(\mathbb{W}\) is as defined in (3.11), and the Cauchy stress, \(\mathbf{T}\), is given by

(3.13)\[\mathbf{T} = \mathbf{C} : ( \nabla _s \mathbf{u} ).\]

This statement of the boundary value problem is often referred to as a weak formulation, since it explicitly requires only a weighted integral of the governing partial differential equations to be zero, rather than the differential equation itself.

It should be clear, based upon the above derivation of the weak form, that the solution \(\mathbf{u}\) of (2.27) – (2.30), sometimes referred to as the strong form, will satisfy our alternative statement summarized by (3.12) and (3.13). Less clear is the fact that solutions of the weak form will satisfy the strong form whenever this formulation admits a solution. Since the continuity requirements for existence of a strong solution are more stringent than for the analogous weak formulation (hence the adjective strong), equivalence between these two forms is restricted to the case when both exist, i.e., whenever a solution of the strong form of the boundary value problem exists, then a weak solution also exists, and these solutions are identical.

It is important to note that the existence of a solution to the weak form of the boundary value problem does not necessarily imply existence of a solution to the strong form. The strong form’s constraints upon solution smoothness imply that for some problems (e.g., point sources that induce jumps in derivative terms), a weak form might exist, but no strong form can be constructed without substantially revising some basic principles of differential calculus. So the existence of a weak solution does not necessarily imply that an identical strong solution exists: only that if a strong solution can be found, it will be identical to the weak solution.

In practice, the existence of a weak solution in these cases turns out to be one of the most important advantages of finite element techniques, because the integral formulations that form the mathematical foundation of finite element approximations permit accurate simulation of important problems that are not readily solved via competing differential techniques derived from strong formulations. Many of the most important problems of computational mechanics (e.g., contact, material discontinuity, structural failure) often admit only weak solutions, and that is one of the main reasons why weak formulations are important in practice.

So the equivalence between strong and weak forms is restricted to those cases where strong solutions exist, and in that case, the strong solution is identical to the analogous weak solution. Although not shown here this equivalence can be rigorously established; the interested reader should consult Reference [[1]] at the end of this chapter for details. We simply remark in the present discussion that the equivalent argument depends crucially on the satisfaction of (3.12) for all \(\mathbf{w} \in \mathbb{W}\), with the arbitrariness of \(\mathbf{w}\) rendering the two statements equivalent whenever the strong solution exists.

Given the requirement of efficient numerical implementation, we can also remark that approximate methods will in effect narrow our definitions of the solution and weighting spaces to finite-dimensional subspaces. Simply stated, this means that rather than including an infinite number of smooth \(\mathbf{u}\) and \(\mathbf{w}\) satisfying the requisite boundary conditions in our problem definition, we will restrict our attention to some finite number of functions comprising subsets of \(\mathbb{S}_t\) and \(\mathbb{W}\).

In so doing we introduce a difference between the solution of our (now approximate) weak form and the strong form, where the degree of approximation is directly related to the difference between the full solution and weighting spaces and the subsets of them used in the numerical procedure. In fact it is this difference that is at the heart of solution verification, an important activity to ensure that an appropriate subset of spaces (i.e., discretization or mesh refinement) is chosen. Solution verification as part of the broader question of verification is discussed in the Solid Mechanics Verification Manual.

Finally, it is worthwhile at this point to make a connection to so-called virtual work methods which may be more familiar to those versed in linear structural mechanics. In this derivation we will work in index notation so that the meaning of the direction notation used above can be reinforced. Accordingly, for a possible solution \(u_i\) of the governing equations, we write the expression for the total potential energy of the system,

(3.14)\[P(u_{i}) = \frac{1}{2} \int _{\Omega} u_{(i,j)} C_{ijkl} u_{(k,l)} \mathrm {d} \Omega - \left[ \int _{\Omega} u_i f_i \, \mathrm{d} \Omega - \int _{\Gamma _{\sigma}} u_i \bar{t_i} \, \mathrm{d} \Gamma \right].\]

Note that the first term on the right hand side represents the total strain energy associated with \(u_i\) and the last two terms represent the potential energy of the applied loadings \(f_i\) and \(\bar{t_i}\). A virtual work principle for this system simply states that the potential energy defined in (3.14) should be minimized by the equilibrium solution. Accordingly, let \(u_i\) now represent the actual equilibrium solution. We can represent any other kinematically admissible displacement field via \(u_i + \epsilon w_i\), where \(\epsilon\) is a scalar parameter (not necessarily small) and \(w_i\) is a so-called virtual displacement, which we assume to obey the boundary conditions outlined in (3.1). This restriction on the \(w_i\) causes \(u_i + \epsilon w_i\) to satisfy the Dirichlet boundary conditions (hence the term kinematically admissible) because the solution \(u_i\) does. We can write the total energy associated with any of these possible solutions via

(3.15)\[\begin{split}P(u_{i} + \epsilon w_i) = \frac{1}{2} \int_{\Omega} \left( u_{(i,j)} + \epsilon w_{(i,j)} \right) C_{ijkl} \left( u_{(k,l)} + \epsilon w_{(k,l)}\right) \mathrm {d} \Omega \\ - \int_{\Omega} \left( u_i + \epsilon w_i \right) f_i \, \mathrm{d} \Omega - \int _{\Gamma _{\sigma}} \left( u_i +\epsilon w_i \right) \bar{t_i} \, \mathrm{d} \Gamma .\end{split}\]

Note that if the potential energy associated with \(u_i\) is to be lower that that of any other possible solution \(u_i + \epsilon w_i\), then the derivative of \(P(u_i + \epsilon w_i )\) with respect to \(\epsilon\) at \(\epsilon =0\) (i.e., at the solution \(u_i\)) should be zero for any \(w_i\) satisfying the conditions in (3.1), since \(u_i\) is an extremum point of the function \(P\). Computing this derivative of (3.15), and setting the result equal to zero, yields

(3.16)\[\left. \frac{d}{d \epsilon} \right|_{\epsilon = 0} P(u_{i}+ \epsilon w_i) = \int _{\Omega} w_{(i,j)} C_{ijkl} u_{(k,l)} \mathrm {d} \Omega - \int _{\Omega} w_i f_i \, \mathrm{d} \Omega - \int _{\Gamma _{\sigma}} w_i \bar{t_i} \, \mathrm{d} \Gamma = 0\]

which must hold for all \(w_i\) satisfying the boundary condition on \(\Gamma _u\). (3.16) can be manipulated further by noting that

(3.17)\[w_{(i,j)} C_{ijkl} u_{(k,l)} = w_{(i,j)} C_{ijkl} E_{kl} = w_{(i,j)} T_{ij} = w_{i,j} T_{ij} .\]

The last equality in (3.17), while perhaps not intuitively obvious, holds because of the symmetry of \(T_{ij}\):

(3.18)\[w_{(i,j)} T_{ij} = \frac{1}{2} (w_{i,j} +w_{j,i}) T_{ij} = \frac{1}{2} (w_{i,j} T_{ij} + w_{j,i} T_{ji} ) = w_{i,j} T_{ij} .\]

Use of (3.17) in (3.16) yields

(3.19)\[\int _{\Omega} w_{i,j} T_{ij} \mathrm {d} \Omega - \int _{\Omega} w_i f_i \, \mathrm{d} \Omega - \int _{\Gamma _{\sigma}} w_i \bar{t_i} \, \mathrm{d} \Gamma = 0 ,\]

which is simply the index notation counterpart of (2.27). Summarizing, we see that the weak or integral form of the governing equations developed previously can be interpreted as a statement of the principle of minimum potential energy. This alternative viewpoint is the reason that the weighting functions \(w_i\) are sometimes called variations or virtual displacements, with the terminology used often depending upon the mathematical and physical arguments used to develop the weak form.

Despite the usefulness of this physical interpretation, it should be noted that the presence of an energy principle is somewhat specific to the case at hand and may be difficult or impossible to deduce for many of the nonlinear systems to be considered in our later study. For example, many systems are not conservative, including those featuring inelasticity, so at best our thermodynamic understanding must be expanded if we insist on formulating such problems in terms of energy principles. Thus, while the energy interpretation is enlightening for many systems, including those featuring elastic continuum and/or structural response, insistence on this approach for more general applications of variational methods can be quite limiting. Conversely, the derivation given in (3.2) – (3.9) does not depend on the system being conservative, nor even upon the form of the constitutive equation used. We will exploit the generality of this weighted residual derivation as we increase the level of nonlinearity and complexity in the chapters to come.

3.3. Fully Dynamic Case

Another advantage of the weighted residual approach is that it can be straightforwardly applied to dynamic problems. Before examining the dynamic case in detail, whose development parallels that of quasistatic problems, it is worthwhile to emphasize again the definitions of the weighting and solution spaces and to highlight the differences between them. Examining the definition of \(\mathbb{S}_t\) in (3.10) and that of \(\mathbb{W}\) in (3.11), we see that \(\mathbb{S}_t\) depends on \(t\) through the boundary conditions on \(\Gamma _u\), while \(\mathbb{W}\) is independent of time. We retain these definitions in the current case and pose the following problem corresponding to the quasistatic system posed previously:

Given the boundary conditions \(\bar{\mathbf{t}}\) on \(\Gamma_{\sigma} \times (0,T)\) and \(\bar{\mathbf{u}}\) on \(\Gamma_{u} \times (0,T)\), the initial conditions \(\mathbf{u}_0\) and \(\mathbf{v}_0\) on \(\Omega\), and the distributed body force \(\mathbf{f}\) on \(\Omega \times (0,T)\), find the \(\mathbf{u} \in \mathbb{S}_t\) for each time \(t \in (0,T)\) such that

(3.20)\[\int_{\Omega} \rho \mathbf{w} \cdot \frac{\partial ^2 \mathbf{u}}{\partial t^2} \, \mathrm{d} \Omega + \int_{\Omega} (\nabla \mathbf{w} ): \mathbf{T} \, \mathrm{d} \Omega = \int_{\Omega} \mathbf{w} \cdot \mathbf{f} \, \mathrm{d} \Omega + \int_{\Gamma _{\sigma}} \mathbf{w} \cdot \bar{\mathbf{t}} \, \mathrm{d} \Gamma\]

for all \(\mathbf{w} \in \mathbb{W}\), where \(\mathbb{S}_t\) is as defined in (3.10), \(\mathbb{W}\) is as defined in (3.11), and the Cauchy stress, \(\mathbf{T}\), is given by

(3.21)\[\mathbf{T} = \mathbf{C} : ( \nabla _s \mathbf{u} ).\]

In addition, the solution \(\mathbf{u}\) is subject to the following conditions at \(t=0\):

(3.22)\[\int_{\Gamma} \mathbf{w} \cdot \left( \mathbf{u}(0) - \mathbf{u}_0 \right) \, \mathrm{d} \Omega = 0\]

and

(3.23)\[\int_{\Gamma} \mathbf{w} \cdot \left( \frac{\partial \mathbf{u}}{\partial t}(0) - \mathbf{v}_0 \right) \, \mathrm{d} \Omega = 0,\]

both of which must hold for all \(\mathbf{w} \in \mathbb{W}\).

The integral form of the dynamic equations given in (3.20) is obtained, just as in the quasistatic case, by taking the dynamic governing partial differential equation, (2.21), multiplying it by a weighting function, integrating over the body, and applying integration by parts to the stress divergence term. The new ingredients in the current specification are the initial conditions summarized by (3.22) and (3.23), which are simple weighted residual expressions of the strong form of the initial conditions given in (2.25).

Before leaving this section, we reemphasize the fact that the weighting functions are time independent while the solution spaces remain time dependent. This fact will have important consequences later when numerical algorithms are discussed, because we wish to use the same classes of functions in our discrete representations of \(\mathbb{W}\) and \(\mathbb{S}_t\). These discretizations will involve spatial approximation, which in the case of \(\mathbb{S}_t\) will leave the time variable continuous in the discrete unknowns of the system to be solved.

This semi-discrete approach to transient problems is pervasive in computational mechanics and has its origin in the difference between the weighting and solution spaces.

The reference for this chapter is [[1]].