MTH4323/5323 Advanced Numerical analysis of partial differential equations Assignment 2 (2023)

To suBM1T: upload on moodle a pdf file (preferably typed, but a clear and legible scanned copy would do) of your solutions AND a .zip archive of your codes . Any numerical result should be accompanied with comments (which may be brief, but meaningful). Each scheme should correspond to a different file. As usual, the marker should be able to reproduce your results, without modifying anything in your code.

LATE suBM1ss1oN: the usual penalty for a late submission is 10% per day.

Total marks: 60 + 5 bonus

1.  A mixed formulation for linear elasticity in terms of displacement-pressure.

Consider a nearly incompressible elastic material (that is, a solid for which A > u > 0) occupying the open bounded domain Ω c R2  with polygonal boundary Γ . Let f e L2 (Ω) a given body force and consider that the solid is clamped everywhere on Γ .  The displacement u satisfies the linear elasticity equation

_2 u div e(u) _ A V(div u) = f    in Ω ,    u = 0   on Γ ,                           (1)

where e(u) :=  (V u + (V u)t ) is the infinitesimal strain tensor.

(a)  [3  marks] Define the additional unknown p =  _Adivu in Ω and show that a variational

formulation for (1) is given by: Find (u, p)  e H0(1)(Ω)  × L0(2)(Ω) (where L0(2)(Ω) = }u e L2 (Ω) : Ω u(z) dz = 0|) such that

2 u      e(u) : e(u) _      p div u  =      f . u Ω                                    Ω                          Ω

u div u +  1        pu = 0

Ω                             A    Ω

V u  e H0(1)(Ω) ,

V u e L0(2)(Ω) .

(2a)

(2b)

(b)  [7 marks] Define suitable bilinear forms and show that this variational problem satisfies the

conditions of Theorem 3.4 in the lecture notes (for perturbed saddle-point problems), therefore establishing the unique solvability and stability of (2).  Comment on the differences you note between the conditions (3.21)– (3.22) of that theorem and the Babuˇska–Brezzi conditions (BB2) and (BB4) in Theorem 3.1.

(c)  [3 marks] Modify the boundary conditions in (1) to be of mixed displacement-traction type (as in the lecture notes, but considering them homogeneous g 1   = g2   =  0) and write the discrete version of (2) using the generic discrete space xh  c HΓ(1)祉 (Ω) for displacement, and yh  c L2 (Ω) for discrete pressure.  Note we are no longer asking yh  c L0(2)(Ω).  In which cases you think this might be not appropriate?

(d)  [7 marks] Implement in FEniCS (or another finite element library of your choice) the method defined in part (c) using the following finite element pairs for discrete displacement and pres- sure: [P2]2 _P1 , [P1]2 _P1 , [P2]2 _P0 (where [P2]2 are piecewise quadratic and overall continuous vector-valued polynomials and P0  are piecewise constants). Use that implementation to solve numerically the linear elasticity problem when Ω = (0 , 1)2 , f = (0, _1)t , and considering that Γu  consists of the bottom and left sides of the square where the domain is clamped, and no- stress conditions are considered on the remainder of the boundary, Γ 扌  (top and right sides). Use the model parameters E = 1 and 页 = 0.499 (Young modulus and Poisson ratio). Employ a coarse mesh (say, no more than 10 elements per side) and plot the solutions (displacement and pressure on the deformed domain) obtained with the three finite element methods.  Compare also with the solution obtained with the [P1]2 finite element method for the primal formulation of linear elasticity seen in Weeks 7-8.

2.  Augmenting a mixed formulation of the Poisson equation.

Let Ω be an open bounded domain in Rd  with smooth boundary Γ, and for a given f e L2 (Ω), consider the problem

_ ∆ o = f   in Ω ,    o = 0   on Γ .

A  weak  formulation  for  this  problem,  seen  in  the  lectures,  writes  as  follows:   Find  (口 , o)   e H(div, Ω) × L2 (Ω) such that

口 . T dz +      o div(T) dz _      u div(口) dz =      f u dz,

Ω                             Ω                                     Ω                                     Ω

for all (T , u) e H(div, Ω) × L2 (Ω).

(a)  [3 marks] Given constants 61 , 62  > 0, prove that the following relations hold

61        (Vo _ 口) . (Vu + T) dz = 0   V (T , u) e H := H(div, Ω) × H0(1)(Ω),

Ω

(3)

(4a)

62         div(口) div(T) dz  =  _ 62         f div(T) dz   V T  e H(div, Ω).           (4b)

Ω                                                              Ω

(b)  [4  marks] Combine  (3),  (4a) and  (4b) to obtain a different variational formulation:  Find

(口 , o) e H such that

A((口 , o), (T , u))  =  F (T , u)       V (T , u) e H ,                                (5)

specifying the form of the bilinear form A : H×H → R and of the linear functional F : H → R.

(c)  [7 marks] Show that, choosing adequately 61  and 62 , problem (5) has a unique solution in H (defining the norm of the product space), depending continuously on the datum f .

(d)  [4 marks] Define a finite element method associated with (5), using (BDM1 , P1 ) elements to construct the discrete space Hh  c H. This differs from the usual mixed finite element method one would consider for (3). Discuss how and why.

(e)  [7 marks] Implement in FEniCS (or other finite element library) the method defined in part (e)

and perform a convergence study plotting (or tabulating) the errors obtained for different mesh refinements of a uniform triangulation of the domain Ω = (_π, π)2 , using the manufactured exact solution o = sin(z) cos(y) and a suitable f .  Comment on the robustness of the error decay with respect to the values of the augmentation constants 61 , 62 .

3.  Error estimate for the primal variable.

Consider the abstract saddle-point problem (3.2) from the lecture notes and its discrete counterpart (3.13).

[5 marks] Assuming that the following C´ea estimate holds for the error associated with 口h

l口 _ 口h l尸  s /1 + 、/ 1 + 、T h(i)f尸h  l口 _ Th l尸 +  亿 h(i)，(f)h  lo _ uh l， ,

show that the error associated with oh  satisfies

lo _oh l，  s  /1+ 、/ 1+ 、T h(i)f尸h  l口 _Th l尸 + /1+  + 、亿 h(i)，(f)h  lo _uh l， .

4.  Mixed formulation for Darcy equations.

Consider the Darcy equations for flow in porous media, solving for the vector field u and the scalar field p (they represent filtration flux on the domain, and fluid pressure, respectively)

u + KVp = g    in Ω ,       div u = f   in Ω,       p = 0   on aΩ,                      (6)

where g e L2 (Ω), f e L2 (Ω) are given data, u > 0 is the constant viscosity of the interstitial fluid, and K is a matrix of permeability (it can depend on the position z e Ω).

(a)  [5 marks] Write (6) in the form: Find (u, p) e x × y such that

a(u, u) + b(u, p) = l(u)

b(u, u) = m(u)

Au e x,

Au e y,

(7a) (7b)

and state a condition on K to guarantee that the problem (7) has a unique solution.

(b)  [5 marks] Define a mixed method for (7) using RTk _ Pk  elements, implement the method in

FEniCS (or other finite element library), and verify numerically its convergence for k = 0, 1 when we take K = 0.001, Ω = (0, 1)2 , and defining the exact pressure p = sin(πzy).

(c)  [5 marks – BONUS (on top of the 60 points)] Let Πh be the Raviart-Thomas interpolator

satisfying

div(u _ Πh u)uh = 0

Ω

and assume that

Auh  e P0 , u e H(div, Ω) n H1 (T), T e Th ,

div (RT0 ) = P0 .