Chapter 04

Exercise 04.01 - Heat transfer in a thin rod

Task

Let us consider a thin rod of length \(L\), having temperature \(t_0\) at the endpoint \(x=0\) and insulated at the other endpoint \(x=L\). Let us suppose that the cross-section of the rod has constant area equal to \(A\) and that the perimeter of \(A\) is \(p\). The temperature \(t\) of the rod at a generic point \(x \in (0,L)\) then satisfies the following mixed boundary-value problem:

(1)\[\begin{split} \left\{ \begin{array}{ll} - k A t'' + \sigma p t = 0 & x \in (0,L), \\ t(0) = t_0, & t'(L) = 0, \end{array} \right. \end{split}\]

having denoted by \(k\) the thermal conducitivity coefficient and by \(\sigma\) the convective transfer coefficient.

Verify that the exact solution of this problem is

\[ \hat{t}(x) = t_0 \frac{\cosh \left[ m \left( L - x \right) \right]}{\cosh \left( mL \right)}, \]

with \(m = \sqrt{\sigma p / k A}\). Write the weak formulation of (1), then its Galerkin-finite element approximation. Show how the approximation error in the \(H^1_0 (0, L)\)-norm depends on the parameters \(k\), \(\sigma\), \(p\), and \(t_0\).

Finally, solve this problem using linear and quadratic finite elements on uniform grids, then evaluate the approximation error.

Proof. Let us evaluate the derivatives of \(\hat{t}\) to check if is the solution of (1).

\[ \hat{t}' (x) = - m t_0 \frac{\sinh \left[ m \left( L - x \right) \right]}{\cosh \left( mL \right)} \]

\[ \hat{t}'' (x) = m^2 t_0 \frac{\cosh \left[ m \left( L - x \right) \right]}{\cosh \left( mL \right)} \]

Then

\[\begin{split} &- k A \hat{t}'' + \sigma p \hat{t} \\ &= - k A m^2 t_0 \frac{\cosh \left[ m \left( L - x \right) \right]}{\cosh \left( mL \right)} + \sigma p t_0 \frac{\cosh \left[ m \left( L - x \right) \right]}{\cosh \left( mL \right)} \\ &= - \cancel{k} \cancel{A} \frac{\sigma p }{\cancel{k} \cancel{A}} t_0 \frac{\cosh \left[ m \left( L - x \right) \right]}{\cosh \left( mL \right)} + \sigma p t_0 \frac{\cosh \left[ m \left( L - x \right) \right]}{\cosh \left( mL \right)} = 0 \end{split}\]

and the boundary conditions are satisfied, so \(\hat{t}\) is a solution.

To derive the weak formulation, we first need to apply a lifting (see [Quarteroni, 2017, Section 3.2.2]). Let \(\Omega := (0, L)\), and \(y = t - R\), where the constant lifting function \(R\) is simply defined as \(R(x) = t_0\) in \(\Omega\). Then, since \(t = y + R\), we have

\[ - k A t'' + \sigma p t = - k A y'' + \sigma p y + \sigma p t_0 \qquad \text{for } x \in (0, L). \]

Hence, solving (1) is equivalent solving to the homogeneous problem

(2)\[\begin{split} \left\{ \begin{array}{ll} - k A y'' + \sigma p y = - \sigma p t_0 & x \in (0,L), \\ y(0) = 0, & y'(L) = 0, \end{array} \right. \end{split}\]

Now, let \(V = \{ v \in H^1 (\Omega) : v(0) = 0 \}\). We obtain the weak formulation of (2) after multiplying the differential equation for a test function \(\varphi \in V\) and integrating in \(\Omega\). In particular, we observe that

\[ - k A \int_\Omega y'' \varphi \dx + \sigma p \int_\Omega y \varphi \dx = \cancel{- k A \left[ y' \varphi \right]_0^L} + k A \int_\Omega y' \varphi' \dx + \sigma p \int_\Omega y \varphi \dx, \]

since \(y'\varphi\) is zero both in \(x=0\) (because \(\varphi(0)=0\)) and in \(x=L\) (because \(y'(L)=0\)).

Putting together the weak formulation of left and right hand side of (2), we get

\[ k A \int_\Omega y' \varphi' \dx + \sigma p \int_\Omega y \varphi \dx = - \sigma p t_0 \int_\Omega \varphi \dx \quad \forall \varphi \in V. \]

The Galerkin approximation is obtained by defining an approximate space \(V_h\) and a basis \(\{ \varphi_j \}_{j=1}^{N_h}\) of \(V_h\). The space \(V_h\) is a Hilbert space with the norm of \(V\). The approximate problem is now

\[ a(y_h, \varphi) = f(\varphi) \qquad \forall \varphi \in V_h, \]

or equivalently

\[ a(y_h, \varphi_j) = f(\varphi_j) \qquad \text{for } j = 1, \dots, N_h. \]

where

\[ a(y_h, \varphi_j) = k A \int_\Omega y_h' \varphi_j' \dx + \sigma p \int_\Omega y_h \varphi_j \dx \]

and

\[ f(\varphi_j) = - \sigma p t_0 \int_\Omega \varphi_j \dx \]

Notice that a solution to the Galerkin approximation exist for the Lax-Milgram lemma (see [Quarteroni, 2017, Lemma 3.1]), since:

• \(a\) is bilinear (trivial to prove)

• \(a\) is continuous:

  \[ \left| a(u,v) \right| = \left| k A \int_\Omega u' v' \dx + \sigma p \int_\Omega u v \dx \right| \le \max\{ kA, \sigma p\} \| u \|_V \| v \|_V \]

• \(a\) is coercive:

  \[ a(u,u) = k A \left| y \right|_{H^1(\Omega)} + \sigma p \| y \|_{L^2(\Omega)} \ge \min\{kA, \sigma p\} \| u \|_V \]

• \(f\) is linear (trivial to prove)

• \(f\) is bounded (so continuous)

  \[ \left| f(v) \right| \le \sigma p t_0 \| v \|_{L^2(\Omega)}. \]

Using the Céa Lemma (see [Quarteroni, 2017, (4.10)]), we get

\[ \| y - y_h \|_V \le \frac{M}{\alpha} \inf_{w_h \in V_h} \| y - w_h \|_V, \]

where \(M\) is the continuity constant of \(a\), i.e. \(M = \max\{ kA, \sigma p\}\) and \(\alpha\) is the coercivity constant, i.e. \(\alpha = \min\{kA, \sigma p\}\). The rest of the estimate depends on the choice of the approximate space \(V_h\).

■

Exercise 04.02 - Temperature of a fluid between two parallel plates.

Task

We consider a viscous fluid located between two horizontal parallel plates, at a distance of \(2H\). Suppose that the upper plate, which has temperature \(t_{\mathsf{sup}}\), moves at a relative speed of \(U\) with respect to the lower one, having temperature \(t_{\mathsf{inf}}\). In such case the temperature \(t: (0, 2H) \to \mathbb{R}\) of the fluid satisfies the following Dirichlet problem:

(3)\[\begin{split} \left\{ \begin{array}{ll}\displaystyle - t''(x) = \alpha (H - x)^2, & x \in (0, 2H), \\ t(0) = t_{\mathsf{inf}}, & t(2H) = t_{\mathsf{sup}}, \end{array} \right. \end{split}\]

where \(\alpha = \frac{4 U^2 \mu}{H^4 k}\), \(k\) being the thermal conductivity coefficient and \(\mu\) the viscosity of the fluid. Find the exact solution \(t(x)\), then write the weak formulation and the Galerkin finite element formulation.

Proof. By integrating twice, we get

\[\begin{split} t(x) &= \int^x \int^z t''(x') \, \dx' \, \mathrm{d}z \\ &= - \alpha \int^x \int^z (H - x')^2 \, \dx' \, \mathrm{d}z \\ &= \alpha \int^x \left[ \frac{1}{3} (H - z)^3 + c_1 \right] \, \mathrm{d}z \\ &= - \frac{\alpha}{12} (H - x)^4 + \alpha \, c_1 x + c_2. \end{split}\]

By applying the boundary values, we get

\[\begin{split} \left\{ \begin{array}{l} \displaystyle t_{\mathsf{inf}} = t(0) = - \frac{\alpha}{12} H^4 + c_2 \\ \displaystyle t_{\mathsf{sup}} = t(2H) = - \frac{\alpha}{12} H^4 + 2 \, \alpha H c_1 + c_2 \end{array} \right. \end{split}\]

The system gives

\[ c_2 = t_{\mathsf{inf}} + \frac{\alpha H^4}{12} \]

and

\[ c_1 = \frac{1}{2 \, \alpha H} \left( t_{\mathsf{sup}} + \frac{\alpha H^4}{12} - c_2 \right) = \frac{1}{2 \, \alpha H} \left( t_{\mathsf{sup}} + \frac{\alpha H^4}{12} - t_{\mathsf{inf}} - \frac{\alpha H^4}{12} \right) = \frac{t_{\mathsf{sup}} - t_{\mathsf{inf}}}{2 \, \alpha H}. \]

So the solution reads

\[ t(x) = - \frac{\alpha}{12} (H - x)^4 + \frac{t_{\mathsf{sup}} - t_{\mathsf{inf}}}{2 \, H} x + t_{\mathsf{inf}} + \frac{\alpha H^4}{12} . \]

To derive the weak formulation, we first need to apply a lifting (see [Quarteroni, 2017, Section 3.2.2]). Let \(\Omega := (0, 2H)\), and \(y = t - R\), where the linear lifting function \(R\) is defined as

\[ R(x) = t_\mathsf{inf} \frac{2H - x}{2H} + t_\mathsf{sup} \frac{x}{2H}, \qquad x \in \Omega. \]

Then, since \(t = y + R\), we have \(t'' = y''\) for all \(x \in \Omega\).

Therefore, to solve (3) we can first solve the homogeneous problem

(4)\[\begin{split} \left\{ \begin{array}{ll}\displaystyle - y''(x) = \alpha (H - x)^2, & x \in (0, 2H), \\ y(0) = 0, & y(2H) = 0, \end{array} \right. \end{split}\]

and then define \(t = y + R\).

Now, let us define \(V = \{ v \in H^1(\Omega) : v(0) = v(2H) = 0 \}\) and \(g(x) = \alpha (H - x)^2\) in \(\Omega\). We obtain the weak formulation of (4) after multiplying the differential equation for a test function \(\varphi \in V\) and integrating in \(\Omega\). By the divergence theorem and by using that \(\varphi\) is zero on the boundary, we get

\[ - \int_\Omega y'' \varphi \dx = - \cancel{\left[ y' \varphi \right]_0^{2H}} + \int_\Omega y' \varphi' \dx \qquad \forall \varphi \in V. \]

So the weak formulation of (4) reads

\[ \int_\Omega y' \varphi' \dx = \int_\Omega g \varphi \dx \qquad \forall \varphi \in V. \]

Similarly to Exercise 1, the Galerkin formulation is

\[ a(y_h, \varphi) = f(\varphi) \qquad \forall \varphi \in V_h. \]

■

Exercise 04.03 - Deformation of a rope

Task

Let us consider a rope with tension \(T\) and unit length, fixed at the endpoints. The function \(u(x)\), measuring the vertical displacement of the rope when subject to a transversal charge of intensity \(w\), satisfies the following Dirichlet problem:

(5)\[\begin{split} \left\{ \begin{array}{ll}\displaystyle - u'' + \frac{k}{T}u = \frac{w}{T}, & x \in (0, 1 ), \\ u(0) = 0, & t(1) = 0, \end{array} \right. \end{split}\]

having indicated with \(k\) the elasticity coefficient of the rope. Write the weak formulation and the Galerkin-finite element formulation.

Proof. Let \(\Omega := (0, 1)\), \(V = \{ v \in H^1(\Omega) : v(0) = v(2H) = 0 \}\), and \(g(x) = \frac{w}{T}\) in \(\Omega\). As in the previous exercises, we obtain the weak formulation of (5) after multiplying the differential equation for a test function \(\varphi \in V\) and integrating in \(\Omega\). By the divergence theorem and by using that \(\varphi\) is zero on the boundary, we get

\[ - \int_\Omega u'' \varphi \dx = - \cancel{\left[ u' \varphi \right]_0^{2H}} + \int_\Omega u' \varphi' \dx \qquad \forall \varphi \in V. \]

So the weak formulation of (5) reads

\[ \int_\Omega u' \varphi' \dx + \frac{k}{T} \int_\Omega u \varphi \dx = \int_\Omega g \varphi \dx \qquad \forall \varphi \in V. \]

Similarly to Exercise 1, the Galerkin formulation is

\[ a(y_h, \varphi) = f(\varphi) \qquad \forall \varphi \in V_h. \]

■

Exercise 04.04 - Symmetry of stiffness matrix corresponds to symmetry of bilinear form

Task

Prove that the matrix \(A\) associated with a bilinear form \(a(\cdot, \cdot)\) is symmetric if and only if the form is symmetric, i.e.

\[ A \text{ is symmetric} \iff a \text{ is symmetric}. \]

Proof. Let us prove \((\Rightarrow)\) first. Let \(u, v \in V\) be arbitrary functions. Then

\[\begin{split} &a(u, v) \\ &= a \left( \sum_{j=1}^{N_h} u_j \varphi_j , \sum_{i=1}^{N_h} v_i \varphi_i \right) \qquad (\textit{decomposition of $u$ and $v$ into bases}) \\ &= \sum_{j=1}^{N_h} \sum_{i=1}^{N_h} u_j v_i a(\varphi_j, \varphi_i) \qquad (\textit{using bilinearity of $a$}) \\ &= \sum_{j=1}^{N_h} \sum_{i=1}^{N_h} u_j v_i A_{i,j} \qquad (\textit{Definition of $A$}) \\ &= \sum_{j=1}^{N_h} \sum_{i=1}^{N_h} u_j v_i A^\top_{i,j} \qquad (\textit{Symmetry of $A$}) \\ &= \sum_{j=1}^{N_h} \sum_{i=1}^{N_h} u_j v_i A_{j,i} \qquad (\textit{Definition of transpose}) \\ &= \sum_{j=1}^{N_h} \sum_{i=1}^{N_h} u_j v_i a(\varphi_i, \varphi_j) \qquad (\textit{Definition of $A$}) \\ &= a \left( \sum_{i=1}^{N_h} v_i \varphi_i , \sum_{j=1}^{N_h} u_j \varphi_j \right) \qquad (\textit{bilinearity of $a$}) \\ &= a(v, u). \end{split}\]

The direction \((\Leftarrow)\) is even easier. Assuming \(a\) to be symmetric, for any \(i,j\) the following chain holds:

\[ A_{i,j} = a (\varphi_j, \varphi_i) = a(\varphi_i, \varphi_j) = A_{j,i}. \]

■

Exercise 04.05 - Céa Lemma for symmetric, continuous, coercive, bilinear form

Task

Prove that for the solution \(u\) of an elliptic problem

\[ a(u, v) = F(v) \qquad \forall v \in V \]

with \(a\) symmetric, continuous, coercive, bilinear form, and assuming \(V_h\) is an approximate subspace of \(V\), and for \(u_h \in V_h\) the approximate solution, the following inequality holds:

\[ \| u - u_h \|_V \le \sqrt{\frac{M}{\alpha}} \inf_{w_h \in V_h} \| u - w_h \|_V, \]

where \(M\) is the continuity constant and \(\alpha\) the coercivity constant.

Proof. The coercivity makes the bilinear form positive definite. A symmetric, positive definite bilinear is an inner product with induced norm \(\|\cdot\|_a := \sqrt{a(\cdot, \cdot)}\) (the so-called energy norm). And for inner products the Cauchy-Schwarz inequality holds:

\[ | a( u, v ) | \le \| u \|_a \| v \|_a. \]

Let us recall the Galerkin orthonormality of the solution \(u_h\):

\[ a( u - u_h , v_h ) = 0 \qquad \forall v_h \in V_h. \]

Now, let \(w_h \in V_h\) arbitrary. Then

\[\begin{split} \| u - u_h \|_a^2 &= a \left( u - u_h, u - u_h \right) \\ &= a \left( u - w_h, u - u_h \right) + a \left( w_h - u_h, u - u_h \right) \qquad (\textit{adding and subtracting the same term}) \\ &= a \left( u - w_h, u - u_h \right) \qquad (\textit{Galerkin orthonormality of $u_h$}) \\ &\le \left\| u - w_h \right\|_a \left\| u - u_h \right\|_a \qquad (\textit{Cauchy-Schwarz}) \end{split}\]

i.e. in the energy norm the following inequality holds:

\[ \| u - u_h \|_a \le \| u - w_h \|_a, \]

and using the continuity of \(a\) the following equation holds:

\[ \| u \|_a = \sqrt{a(u , u)} \le \sqrt{M} \| u \|_V. \]

\[\begin{split} \alpha \| u - u_h \|_V^2 &\le a \left( u - u_h, u - u_h \right) \qquad (\textit{coercivity of $a$}) \\ &= \| u - u_h \|_a^2 \qquad (\textit{definition of $\|\cdot\|_a$}) \\ &\le \| u - w_h \|_a^2 \qquad (\textit{Energy norm inequality}) \\ &\le M \| u - w_h \|_V^2 \qquad (\textit{Conversion to $V$-norm}). \end{split}\]

As \(w_h\) is chosen arbitrarily, we can take the inf:

\[ \| u - u_h \|_V \le \sqrt{\frac{M}{\alpha}} \inf_{w_h \in V_h} \| u - w_h \|_V. \]

■