在Arbitrary Lagrangian Eulerian架構下發展一具有守恆形式的有限元素法

(1)

國立臺灣大學工學院工程科學及海洋工程學系博士論文

Department of Engineering Science and Ocean Engineering College of Engineering

National Taiwan University Doctoral Dissertation

在Arbitrary Lagrangian Eulerian 架構下發展一具守恆形式的有限元素法

Development of an Arbitrary Lagrangian Eulerian Finite Element Formulation in Conservative Form

伊菲利 Filip Ivanˇci´c

指導教授: 許文翰博士、馬克沁博士 Advisors: Tony Wen–Hann Sheu, Ph.D., and

Maxim Solovchuk, Ph.D.

中華民國109年8月 August, 2020

(2)

(3)

doi:10.6342/NTU202003676

(4)

Acknowledgements

I would like to express my sincere gratitude to my advisors, Prof. Tony Sheu and Prof.

Maxim Solovchuk, for their continuous support and guidance, patience, motivation, knowl- edge and time invested in my research and academic development. I thank them for always encouraging me to conduct individual research, yet always being available with useful advice and guidelines. I have learned a lot from them.

I would like to thank professors Hong–Yueh Lo, Chien–Chen Chang, An–Bang Wang, Yi–Ju Chou, I–Liang Chern and Shiu–Wu Chau for accepting to read and review my thesis and for being a part my dissertation defense committee. I thank them sincerely for their valuable comments and suggestions.

My research has been made possible thanks to the financial support from National Health Research Institutes (NHRI Taiwan) and Ministry of Science and Technology (MOST Taiwan).

I owe my deepest gratitude to my former master thesis advisor, Prof. Boris Muha, who was the first one to introduce and guide me to the world of mathematical fluid dynamics.

Thanks to Neo Shih–Chao Kao who helped me a lot during my stay in Taiwan, espe- cially in my first year.

I wish to thank Ms. Wei–Hsuan Huang for her precious help with the administrative work during my study at NTU. I would have been lost without her.

(5)

doi:10.6342/NTU202003676 Thanks to Bela and Melanie for helping me get Visa prior to coming to Taiwan.

Finally, I owe my deepest gratitude to my parents, my brothers and my grandparents for their unconditional support and help they selflessly provided me during my whole life.

All of my success I equally consider theirs as well.

(6)

摘要

本論文目的主要為發展一數值方法用以模擬在時變域上的多物理場系統。考慮此類問題的動機大部分來自於以時變域的偏微分方程式觀點所描述的生醫及生物流體力學問題。為此，我們將採用有限元素法來求解此類問題。此外，我們只考慮在演化過程中其場域拓樸不變的問題，此限制允許我們採用以對齊網格來顯式描述場域介面的任意拉格朗日-歐拉架構。因此，整個數值方法屬於顯式界面追蹤類別方法中的移動網格架構。論文的第一部分主要在守恆形式的ALE架構下推導出一個新穎的有限元數值方法，提供一個系統的方法用以消除由於移動網格下而產生的人工沉降及源項。即便此類人工數值沉降及源項已被眾所皆知，此問題仍是一個開放性且具挑戰性的主題。質量及離散空間律的守恆則為另外兩個需要被解決的問題，而所本論文的方法正是在結合此兩個特徵所發展出的。論文的第二部分將採用所提出的數值方法來解決真實的流體問題，將會著重在自由液面流跟流

固耦合問題這兩類主題上，所選取的驗證問題中將會驗證所開發的方法具有良好

靈活性及可信賴性。

關鍵字: 有限元方法任意拉格朗日歐拉方法移動網格人工沉降源自由液面流流固耦合

(7)

doi:10.6342/NTU202003676

(8)

Abstract

The purpose of this thesis is to develop a numerical method for simulations of multiphysical systems on evolving domains. Motivation for the problems considered in this work comes largely from the field of bio–medicine and bio–fluid mechanics. These multiphysical systems are described in terms of systems of partial differential equations (PDEs) posed on time dependent domains. Finite element method (FEM) is employed for numerical approximation of such problems. Furthermore, only a special class of ”domain–evolving” problems is considered – problems in which domain’s topology does not change during its evolution. This restriction allows to work within the so–called arbitrary Lagrangian–Eulerian(ALE) framework in which the interface of domain is described explicitly by the aligned mesh. Thus, the complete numerical method employed falls under a moving mesh category within an explicit, so called interface tracking, approach.

The first part of the thesis deals with derivation of a novelty approach in finite element method within ALE framework focused on conservative formulations. This approach offers a systematic way to eliminate artificial sinks and sources arising from the moving mesh. Although the numerical origins of these artificial sinks and sources are well known, this problematics still remains to be an active and challenging topic. The mass conser- vationproblem and the discrete space conservation law (SCL) are the two major issues

(9)

doi:10.6342/NTU202003676 resolved; actually, the novelty approach is integrated upon these two characteristics.

In the second part of the thesis, the newly proposed approach is applied to (academic) problems arising from the real world situations. The attention is on two particular class of problems: free–surface flows and fluid–structure interaction (FSI) problems. The flex- ibility and credibility of the methodology derived in the first part are demonstrated on selected examples.

Keywords: finite element method, arbitrary Lagrangian–Eulerian, moving mesh, artificial sink/source, free–surface flow, fluid–structure interaction

(10)

I FEM approximation of differential problems within ALE frame-

work 1

1 Parabolic equations in time–dependent domains 3

1.1 The Arbitrary Lagrangian–Eulerian framework . . . 4

1.1.1 Fundamentals of ALE framework . . . 5

1.1.2 The ALE temporal derivative . . . 8

1.1.3 Euler expansion formula . . . 9

1.1.4 Test function spaces in ALE framework . . . 10

1.1.5 Strong forms of conservation laws . . . 12

1.1.6 Weak formulations of conservation laws . . . 14

1.1.7 The transformation of configurations . . . 16

1.1.8 Pullback of weak formulation to reference configuration . . . 19

(11)

doi:10.6342/NTU202003676

1.2 ALE finite element formulation . . . 20

1.2.1 Finite element discretization of the ALE map . . . 23

1.2.2 Finite element formulation . . . 24

1.2.3 Remark on notation . . . 25

1.2.4 Example . . . 28

1.3 Artificial sinks/sources on moving meshes . . . 33

2 Volume preserving moving mesh method 37 2.1 Motivation . . . 37

2.2 Construction of volume preserving deformation . . . 41

2.3 FEM formulation with Lagrange multiplier . . . 48

2.4 Numerical validation . . . 52

2.4.1 Volume gain/loss . . . 53

2.4.2 Accumulated volume oscillation during the simulation . . . 53

2.5 Discussion . . . 56

3 Space Conservation Law 57 3.1 Space conservation law . . . 59

3.1.1 SCL in finite element method . . . 62

3.2 Mesh velocity calculation and vanishing discrete SCL . . . 64

3.2.1 Mesh velocity piecewise constant in time . . . 64

3.2.2 Mesh velocity continuous in time . . . 67

3.3 Discretization schemes . . . 70

3.3.1 Implicit Euler scheme . . . 72

3.3.2 Crank–Nicolson scheme . . . 73

3.3.3 Backward differentiation formula – BDF . . . 73

3.4.1 Stability . . . 77

3.4.2 Convergence . . . 79

3.4.3 Accuracy . . . 82

(12)

4 Stabilization methods for FEM on moving meshes 87 4.1 S–SS decomposition of parabolic equations . . . 89

4.1.1 S-SS decomposition of differential operator on time–dependent domain . . . 90

4.2 Numerical diffusion based stabilizations . . . 92

4.2.1 Selection of the stabilization parameter on time–dependent domain 98 4.3 Temporal discretization of stabilized conservative formulation . . . 99

4.4 Numerical validation for scalar conservation laws . . . 102

4.4.1 Heat equation on an oscillating domain . . . 102

4.4.2 Convergence of stabilized methods on moving meshes . . . 108

4.5 Stabilization of the Navier–Stokes equations . . . 111

4.5.1 S-SS decomposition of Navier–Stokes differential operator . . . . 111

4.5.2 Ladyˇzenskaya–Babuˇska–Brezzi (inf–sup) condition . . . 113

4.5.3 Flow past an oscillating cylinder . . . 113

4.6 Beyond convection stabilization . . . 118

4.6.1 Convection–diffusion equation on domain with a moving cylinder 118 4.7 Discussion . . . 120

5 Curvature evaluation of mesh–fitted interface in FEM 123 5.1 Curvature in weak form: employment of the Laplace–Beltrami operator . 127 5.2 Introduction of spurious velocities due to curvature approximation . . . . 129

5.2.1 Model problem setup . . . 129

5.2.2 Finite element formulation . . . 131

5.3 Detour framework for Laplace–Beltrami operator in finite elements . . . . 136

5.3.1 Finite element formulation for discrete curvature calculation . . . 136

5.3.2 Effect of finite element spaces on numerical curvature . . . 138

5.3.3 Beyond linear meshes . . . 141

5.4 FEM formulation with the numerically corrected curvature . . . 145

(13)

doi:10.6342/NTU202003676

5.4.1 Decoupling the curvature evaluation from the primary problem . . 147

5.4.2 Numerical validation – FEM on polygnal meshes . . . 149

5.4.3 Numerical validation – isoparametric concept . . . 151

II Applications 161

6 Dynamic contact line problem – sliding droplet 163 6.1 Introduction . . . 164

6.2 Moving contact line problem . . . 166

6.3 Non–dimensionalization . . . 168

6.4 Weak and FEM formulation . . . 170

6.5 Numerical results . . . 171

6.5.1 Mesh adaptation . . . 171

6.5.2 Droplet on a horizontal solid surface . . . 172

6.5.3 Droplet on an inclined solid surface . . . 174

7 Chemotaxis 179 7.1 Introduction . . . 180

7.2 Chemotaxis–diffusion–convection (CDC) coupling system with fixed free surface . . . 183

7.2.1 The dimensional CDC system . . . 183

7.2.2 The dimensionless CDC system . . . 186

7.3 Chemotaxis–diffusion–convection (CDC) coupling system with dynamic free surface . . . 188

7.3.1 The generalized Navier boundary conditions . . . 188

7.3.2 CDC system with dynamic free surface . . . 191

7.4 FEM formulation . . . 194

7.4.1 Weak formulation of system (7.14,7.15) . . . 194

(14)

7.4.2 Numerical (FEM) approach . . . 196

7.4.3 Multiscale to singlescale formulation . . . 198

7.5 Numerical simulations . . . 199

7.5.1 Two–dimensional setup . . . 199

7.5.2 Three–dimensional setup . . . 202

7.5.3 Bacterial chemotaxis in bacterial droplets . . . 207

7.6 Free (thermal) convection . . . 209

7.6.1 Mathematical model . . . 209

7.6.2 Non–dimensionalization . . . 211

7.6.3 Weak and FEM formulation . . . 213

7.6.4 Numerical results . . . 214

8 Fluid–Structure Interaction 219 8.1 Introduction . . . 220

8.2 Mathematical models for blood and vessel wall . . . 222

8.2.1 Mathematical model for the blood . . . 222

8.2.2 Mathematical model for the vessel wall . . . 223

8.3 Fluid–structure interaction modeling . . . 226

8.3.1 Fluid–structure coupling . . . 227

8.3.2 Weak formulation and implicit coupling . . . 228

8.3.3 Implicit coupling in FEM formulation . . . 231

8.3.4 Implicit coupling through Lagrange multipliers . . . 232

Conclusions 239

(15)

doi:10.6342/NTU202003676

(16)

List of Figures

1.1 Transformation between two configurations. . . 6 1.2 Example of a C¹–smooth domain Ω ⊂ R² (a) and its discretized coun-

terpart Ω_h (b). In figure (b) a triangulation T_h of Ω_h is shown. Ω_h is polygonal approximation of Ω. . . 21 1.3 A reference simplex eK , its image bK = fM^K₁^b ⊂ bΩ_h in ALE reference

domain, and its image in physical domain Ω, K = bA_t( bK ). . . 23 1.4 The domain volume gain and loss during the simulation due to mesh

movement. . . 33 1.5 Variation of u over time for various FEM formulations. Finite element

space is chosen as V = P1. f (t) = R

Ωu dx denotes the variation of u over Ω at time t. . . 33

2.1 Domain Ω has a fixed rigid bottom and walls at all times. Free surface Σ is moving in time. . . 40 2.2 Target velocity field ϑh (a), corrected velocity field q¹_h (b), and their dif-

ference (c). Domains obtained by A^[1,0]_h constructed from ϑh and q_h are shown in figure (d). . . 54

(17)

doi:10.6342/NTU202003676 2.3 Convergence of the artificial velocity field (q^k_h)_k towards the velocity

which results in volume preserving deformation A^[0,1]_h . Ordinate is shown in log–scale. . . 54 2.4 Artificial volume gain/loss for various constructions of deformation maps

A^[n+1,n]_h and different choices of time step ∆t. . . 55

3.1 Sketch of the ALE transformation K_hⁿ 7→ K_hⁿ⁺¹with ∆t = 0.1. . . 61 3.2 Evolution of discrete configurations on a time interval [t_n−1, t_n+1]. . . 65 3.3 Sketch of evolution of the mesh node denoted by x = x(x, t) on intervalb

[t_n−1, t_n+1]. . . 68 3.4 Sketch of the implicit Euler method on [t_n−1, t_n+1]. Evaluating uⁿ⁺¹_h takes

place in configuration on time interval [t_n, t_n+1]. Test functions involved in SCL ”carry the same weight” uⁿ⁺¹_h on [t_n, t_n+1] (pure implicit method) so there is no violation of SCL. . . 72 3.5 Sketch of the Crank–Nicolson method on [t_n−1, t_n+1]. Evaluating uⁿ⁺¹_h

takes place in configuration on time interval [t_n, t_n+1]. Test functions involved in SCL ”carry the same weight” ¹₂(uⁿ_h+uⁿ⁺¹_h ) on [t_n, t_n+1] so there is no violation of SCL. . . 73 3.6 Sketch of the BDF2 method on [tn−1, t_n+1]. Evaluating uⁿ⁺¹_h takes place

in configuration on time interval [t_n−1, t_n+1]. In order not to violate the SCL, the whole evolution of configuration on interval [t_n−1, t_n+1] has to be taken into account. Test functions involved in SCL ”carry different weights on different time intervals”, namely, for the case of BDF2, ³₂uⁿ⁺¹_h on [tn, tn+1] and −¹₂uⁿ⁺¹_h on [tn−1, tn]. . . 74 3.7 The L²(Ω(t)) norms of the discrete solutions for different time steps and

different methods: implicit Euler method (a) (mIE-dc, m denoting modified), Crank–Nicolson method (b) (mCN–dc), BDF2 (mBDF2–dc) (c) and BDF3 (c) (mBDF3–dc) methods. Grid velocity is piecewise constant in time calculated by the first proposed approach (3.12) (discontinuous in time reconstruction). . . 80

(18)

3.8 The L²(Ω(t)) norms of discrete solutions for different time steps and different methods: implicit Euler method (mIE-c, m denoting modified), Crank–Nicolson method (mCN–c), BDF2 (mBDF2–c) and BDF3(mBDF3–

c) methods. Grid velocity is continuous in time calculated by the second proposed approach (3.20) (continuous in time reconstruction). . . 81 3.9 Rates of convergence for different time–stepping methods in log–log scale:

implicit Euler method (mIE-c, mIE-dc), Crank–Nicolson method (mCN–

c, mCN–dc), BDF2 (mBDF2–c, mBDF2–dc) and BDF3(mBDF3–c, mBDF3–

dc) methods. m denotes modified, c continuous and dc discontinuous in time grid velocity reconstruction. x–axis represents the time step ∆t in discretization scheme (∆t = 0.001, 0.005, 0.01, 0.05), y–axis represents the kuⁿ⁺¹_h − u_h(t_n+1)k_L²_(Ω(t_n+1₎₎with n + 1 such that t_n+1 = 0.3. Dashed black lines denote the slopes. . . 83 3.10 The L²(Ω_t)–norms of errors between the exact and numerical solutions

for different schemes and time step ∆t = 0.05. In the legend ”fixed”

refers to solutions obtained on fixed grids, while ”mov-wSCL” and ”mov- nSCL” for ones obtained on moving grids with the proposed non–violating SCL (wSCL, w denoting with) schemes (”dc” and ”c” standing for the discontinuous and continuous in time reconstruction of grid velocity), and the classical, SCL–violating (nSCL, n denoting no) schemes, respectively. . . 85

4.1 Energy of the solution for problem (4.20) obtained by the implicit Euler method during the time interval [0, 0.4] with two different choices of finite element space: P1 on the left and P2 on the right. In legend, numerical value next to the ”stabilized Galerkin” represents the choice of % in S^%_IE. . 104 4.2 Implicit Euler method. Solution of problem (4.20) at time t = 0.1 in case

of P1 finite element space. Standard Galerkin FEM (left) produces spurious oscillations in the solution, while stabilized FEM produces smooth solution. . . 105

(19)

doi:10.6342/NTU202003676 4.3 Implicit Euler method. Solution of problem (4.20) at time t = 0.1 in case

of P2 finite element space. Choice of P2 finite element space produces a more stable solution in comparison with choice of P1 finite element space (expected). Still, some spurious oscillations can be observed which vanish in all cases of stabilized FEM. . . 105

4.4 L²–energy of the solution of problem (4.20) obtained by the Crank–Nicolson method during the time interval [0, 0.4] for P1 finite element spaces and two different choices of the time step (∆t = 0.01 on the left, and ∆t = 0.001 on the right). Energy of the solution is non–oscillating only for a sufficiently small time step. In legend, numerical value next to the ”stabilized Galerkin” represents the choice of % in S^%_IE. . . 106

4.5 Crank–Nicolson method. Solution of problem (4.20) at time t = 0.15 in case of P1 finite element space. Standard Galerkin FEM (left) produces spurious oscillations in the solution, while stabilized FEM produces smoother solution. . . 107

4.6 Crank–Nicolson method. Solution of problem (4.20) at time t = 0.15 in case of P2 finite element space. Choice of P2 finite element space produces more stable solution in comparison with choice of P1finite element space (expected). Still, some spurious oscillations can be observed which vanish in case of stabilization with SUPG and GLS technique. . . 107

4.7 Evolution of error in L²(Ω_t) norm computed with time step ∆t = 10⁻³ and Crank–Nicolson method for the temporal discretization (problem (4.25)).

P²finite element space is employed with GLS stabilization technique. . . 110

4.8 Convergence rate for GLS stabilization technique on moving mesh in the sense of L²(0, T ; L²(Ω_t)) norm. Crank–Nicolson method is employed for temporal discretization with time step ∆t = 10⁻³ (problem (4.25)). . . 110

(20)

4.9 Computational mesh for the flow past oscillating cylinder problem. Cylin- der oscillates along {y = 0} line and rotates around its center. Initial shape and position of the cylinder are given by parametric description:

(x(t), y(t)) = (¹₂ + 0.08 cos(t), 0.1 sin(t)), t ∈ [0, 2π). . . 116 4.10 Number of Newton’s iterations on two different meshes where finite el-

ement space is chosen as [P^b1]² for velocity and P1 for pressure. ”Stan- dard” denotes the method without stabilization but on the very fine mesh (h = 0.004). . . 116 4.11 Pressure field in area near the cylinder obtained by standard Galerkin

and GLS method on the same mesh with characteristic size h = 0.03 at time t = 0.15. Newton method does not converge for standard Galerkin method for later times. . . 117 4.12 x–component of the velocity field in area near the cylinder obtained by

standard Galerkin and GLS method on the same mesh with characteristic size h = 0.03 at time t = 0.15. Newton method does not converge for standard Galerkin method for later times. . . 117 4.13 Energies of velocity field produced by stabilized methods on two different

meshes. ”Standard” denotes the method without stabilization but on very fine mesh (h = 0.004), i.e. our reference solution. . . 117 4.14 Numerical solution of equation (4.31) obtained by stabilized FEM (left)

and modified discontinuity–capturing stabilized FEM (right) employing P1finite element space and implicit Euler method. Small undershoots are still present when β is linearized using the previous step solution uⁿ_h as described in subsection 4.6.1. If one considers an iterative algorithm and uses the current iteration u^n+1,k_h for linearization of β, undershoots vanish. 120

5.1 Sketches of the domain Ω and its discrete counterparts, Ω^r1_h and Ω^r2_h . Tri- angulation T^r1_h of Ω^r1_h is coarser than triangulation T^r2_h of Ω^r2_h near the interface. . . 130

(21)

doi:10.6342/NTU202003676 5.2 Mesh on which computations resulting in Figures 5.3, 5.4 and 5.5 were

performed. Mesh parameters are: h = 0.943 (the largest triangle diam- eter) and h^e_Σ

h = 0.416 (the longest edge on the interface Σ_h). Formal definitions of mesh parameters are given in equations (5.20) and (5.21). . 133 5.3 Numerical solution (vh, ph) ∈ Xh× Mh in case of Crouzeix–Raviart fi-

nite elements space. Solution is obtained using the FEM formulation (5.10).133 5.4 Numerical solution (vh, p_h) ∈ X_h× M_hin case of the Taylor–Hood finite

elements space. Solution is obtained using the FEM formulation (5.11). . 134 5.5 Numerical solution (vh, p_h) ∈ X_h× M_h in case of the Mini finite ele-

ments space. Solution is obtained using the FEM formulation (5.11). . . . 135 5.6 Three different triangulations of the same discrete domain Ω_h. All three

triangulations approximate the circle with the same precision, i.e. by the same polygon. . . 137 5.7 Mean curvature and mean curvature vector of the unit circle obtained by

finite element method using FEM formulation (5.13) for Wh = P²× P² and Kh = P2. . . 139 5.8 Sketch of DOFs for P1, P2and P3 finite element spaces on one part of Σh

that consists of three edges. Dotted line illustrates how Laplace–Beltrami operator smoothens the discrete curve in various cases. . . 140 5.9 Mean curvature and mean curvature vector in case W_h(T_h,1) = [P1, P1]

and K_h(T_h,1) = P1. This figure should be compared with Figure 5.8 (a). . 142 5.10 Mean curvature and mean curvature vector in case W_h(T_h,1) = [P2, P2],

K_h(T_h,1) = P2 and W_h(T^s2_h,1) = [P1, P1], K_h(T^s2_h,1) = P1. This figure should be compared with Figure 5.8 (b). . . 142 5.11 Mean curvature and mean curvature vector in case Wh(Th,1) = [P³, P³],

K_h(T_h,1) = P3 and Wh(T^s3_h,1) = [P1, P1], K_h(T^s3_h,1) = P1. This figure should be compared with Figure 5.8 (c). . . 143

(22)

5.12 Mean curvature and mean curvature vector for the cases W_h(Σ_h, P2) = [P2, P2, P2], K_h(Σ_h, P2) = P2 (subfigures (a) and (b)) and W_h(Σ_h, P1) = [P1, P1, P1], K_h(Σ_h, P1) = P1 (subfigures (c) and (d)). Σ_h is the discrete counterpart of unit sphere Σ = {x²+ y²+ z² = 1}. . . 144 5.13 Mean curvature vector for all of the combinations of spaces (A^k_h = [Pk, Pk], W^l_h =

[Pl, Pl]), k = 2, 3, l = 1, 2, 3. For the cases where k ≥ l the mean curvature vector is credibly evaluated and in all of these cases 0.9 ≤ | h_h| ≤ 1.1. For the case of k < l (k = 2, l = 3) spurious oscillations appear (c). . 146 5.14 L² and L^∞–error in velocity and pressure fields with respect to the mesh

parameter h^e_Σ

h(denoted h_Σ,ein plot) in log–log scale. Taylor–Hood finite element spaces are employed for the unknown and A¹_h space for the geometry construction (linear mesh). Π_h(v_e) denotes the projection of exact solution to the corresponding finite element space (f = v, p). . . 155 5.15 L² and L^∞–error in velocity and pressure fields with respect to the mesh

parameter h^e_Σ

h (denoted h_Σ,e in plot) in log–log scale. Crouzeix–Raviart finite element spaces are employed for the unknown and A¹_hspace for the geometry construction (linear mesh). Π_h(f_e) denotes the projection of exact solution to the corresponding finite element space (f = v, p). . . 155 5.16 L² and L^∞–error in velocity and pressure fields with respect to the mesh

parameter h^e_Σ

h (denoted h_Σ,e in plot) in log–log scale. Scott–Vogelius finite element spaces are employed for the unknown and A¹_hspace for the geometry construction (linear mesh). Π_h(f_e) denotes the projection of exact solution to the corresponding finite element space (f = v, p). . . 156 5.17 L² and L^∞–error in velocity and pressure fields with respect to the mesh

parameter h^e_Σ

h (denoted h_Σ,e in plot) in log–log scale. Isoparametric Crouzeix–Raviart and isoparametric Taylor–Hood finite element spaces are employed for the unknown and A²_h space for the geometry construction (linear mesh). Πh(f_e) denotes the projection of exact solution to the corresponding finite element space (f = v, p). . . 158

(23)

doi:10.6342/NTU202003676 6.1 Liquid droplet on an inclined plane with the inclination angle α. . . 165

6.2 Sketch of the vectors describing the geometry of droplet interface in 2D (a) and 3D (b) cases. . . 167 6.3 Initial mesh is pre–adapted – it is denser near the free surface and, espe-

cially, near the contact points. . . 172 6.4 Initial configuration (light gray) and configuration near the final equilib-

rium state (dark gray). Volume is preserved within the order of 10⁻¹⁰. . . 173 6.5 Velocity field at a time near the initial time (start of simulations) when

droplet starts to deform (a) and at a time when the steady state will soon be reached (b). . . 174 6.6 Magnitude of velocity field at time near the initial time (start of simula-

tions) when droplets starts to deform (left) and at time when the steady state will soon be reached (right). Volume is preserved within the order of 10⁻¹⁰. . . 175 6.7 Velocity magnitude (lower figures) and pressure field (upper figures) at

different times. . . 177 6.8 The droplet mesh at time t = 4 (a) and enlarged caption of the mesh near

the contact points (b) and (c). . . 178 7.1 The sketch of the time–independent domain occupied by the suspension

of bacteria in water. . . 184 7.2 The sketch of part of the domain near the contact line with the various

unit vectors necessary for formulating the boundary conditions. . . 189 7.3 An example of pre–adapted mesh with finer triangulation near the free

surface and contact points. Figure shows the clip of the mesh near the right contact point. The meniscus position is obtained from the Young–

Laplace equation with contact angle θ = 5π/3. . . 197 7.4 Formation of the bacteria depletion layer near the free surface at early

stages of chemotaxis phenomenon, at times t = 0.005 (a) and 0.05 (b), for the case of Ra = 2000. ”n” denotes bacteria concentration. . . 202

(24)

7.5 State of chemotaxis phenomenon at time t = 0.175 for Ra = 2000. De- veloped physics, bacteria (a) and oxygen (b) concentration and velocity (c) properties, is mirror symmetric with respect to line {x = 0}. Con- centration of bacteria is in logarithmic scale. Black curves represent the induced velocity streamlines. ”n” denotes the bacteria concentration. . . . 203 7.6 Bacterial plumes from Figure 7.5 (a) at a time when they hit the bottom

of the container (t > 0.175). The shape changes into a mushroom–like because of the head of the plume being convected with the fluid. ”n”

denotes bacteria concentration. . . 203 7.7 State of chemotaxis phenomenon at time t = 0.395 for Ra = 2000.

Developed physics, bacteria (a) and oxygen (b) concentration and velocity (c) properties, is mirror symmetric with respect to line {x = 0}.

Concentration of bacteria is in logarithmic scale. Black curves represent the induced velocity streamlines. ”n” denotes the bacteria concentration.

Concentration of oxygen in lower layers of the container has reached its critical value (c = 0.3) at which bacteria become inactive. Comparison with Figure 7.5 (a) reveals that plume merging occured. . . 204 7.8 Layer of the domain near the free surface scaled with respect to the y–

direction in order to emphasize distortion of the free surface. ”n” denotes bacteria concentration. . . 204 7.9 The formation of the bacterial plumes (at time t = 2.6). The domain in

(a) is clipped in order to illustrate the plume formation inside the domain – on the cross section. . . 205 7.10 Figure shows an enlarged cubic clip of the domain in order to get a better

insight of the details. In (a) the velocity streamlines are shown. In (b) velocity streamlines in context of bacterial plumes are shown. . . 205 7.11 Velocity streamlines in 3D context. Streamline plot is clipped at front in

order to illustrate the vortex–like patterns. . . 206

(25)

doi:10.6342/NTU202003676 7.12 Surface streamlines corresponding to state shown at Figure 7.10. Convec-

tion cells can be observed. . . 206 7.13 State of the two dimensional bacterial chemotaxis at a time when biocon-

vection patterns are already developed. . . 208 7.14 State of the three dimensional bacterial chemotaxis at time when biocon-

vection patterns are already developed. . . 208 7.15 Temperature distribution (a) and vertical component of the velocity field

(b) on a clip of the domain. From the vertical component of the velocity field it can be observed that convection from top to bottom is of small magnitude. Hence, the primary mechanism for the heat transport from bottom to top is diffusion. . . 215 7.16 Temperature distribution and surface velocity streamlines on the free sur-

face. Tangential velocities induced by the surface tension forces arising from the surface tension gradient can be observed. . . 215 7.17 Enlarged clip of the surface streamlines on the free surface (a), and scaled

deformation of the free surface (b). . . 216 8.1 Simplified blood vessel geometry. . . 222 8.2 Mesh of the discretized domain (a) and enlarged mesh around the elastic

beam (b). . . 235 8.3 Velocity (on the left) and corresponding pressure (on right) fields at four

different time instants, at the time when beam has already started oscillating. Deformation of the (enlarged) beam can be observed. . . 236

(26)

List of Tables

5.1 Errors in L²–norm in velocity field with respect to the finite element space employed, the technique for numerical evaluation of the curvature and the mesh parameter h^e_Σ_h. Πhdenotes the projection to the corresponding finite element space. ”LB–embedded” refers to formulations (5.10,5.11) where curvature computation is embedded directly into the formulation. . . 151 5.2 Errors in L²–norm in pressure field with respect to the finite element space

employed, the technique for numerical evaluation of the curvature and the mesh parameter h^e_Σ

h. Π_hdenotes the projection to the corresponding finite element space. ”LB–embedded” refers to formulations (5.10,5.11) where curvature computation is embedded directly into the formulation. . . 152 5.3 Errors in L^∞–norm in velocity field with respect to the finite element

space employed, the technique for numerical evaluation of the curvature and the mesh parameter h^e_Σ

h. Π_h denotes the projection to the corresponding finite element space. ”LB–embedded” refers to formulations (5.10,5.11) where curvature computation is embedded directly into the formulation. . . 153

(27)

doi:10.6342/NTU202003676 5.4 Errors in L^∞–norm in pressure field with respect to the finite element

space employed, the technique for numerical evaluation of the curvature and the mesh parameter h^e_Σ

h. Π_h denotes the projection to the corresponding finite element space. ”LB–embedded” refers to formulations (5.10,5.11) where curvature computation is embedded directly into the formulation. . . 154 7.1 Nomenclature description. . . 185 7.2 Definitions of characteristic values and dimensionless numbers governing

the system (7.14)–(7.15). σ₀ denotes the surface tension for pure water (σ₀ ≈ 72 mN/ m at 25^◦C). . . 192

(28)

Introduction

Time–dependent domains appear in many multiphysics models governing various phenomena arising from the physics and engineering, commonly related to fluid flow problems. Deformation of the free interface in free surface flows and deformations of the structure in fluid–structure interaction (FSI) problems correspond to the deformation of the fluid domain.

Fluid dynamics is most naturally described in terms of fluid velocity field in Eulerian coordinates. In this work, a special class of ”domain–evolving” problems is considered – problems in which domain’s topology does not change during its evolution. This class of problems, for example, excludes the breaking waves and splashing problems. Restriction to such class of problems allows to work within the so–called arbitrary Lagrangian–

Eulerian(ALE) framework in which the interface of domain is described explicitly by the aligned mesh in the numerical method. Hence, the numerical method employed falls under a moving mesh category within explicit, so called interface tracking, approach.

Finite element method (FEM) for simulating moving mesh problems is employed in this work (see, e.g., [1, 3, 2]). Within ALE framework, FEM is employed for the spatial discretization of the specific phenomena governing equations (see, e.g., [4, 5, 6]). For the time discretization, typically finite difference method is employed. Fundamentals of ALE framework are recalled in Chapter 1. General idea consists of an interplay between the

(29)

doi:10.6342/NTU202003676 fixed referential domain and the current physical domain occupied by the medium. The

interplay between these two domains is realized through the so called ALE map which maps the referential domain into the physical one:

Ab_t: bΩ → R^d, bΩ 7→ Ω. (0.1)

In the same manner, the discrete referential triangulation bT_h(bΩ_h) is mapped into the physical triangulation T (Ω_h) through the discrete ALE map bA_h,t.

An important technical problem arises in the context of incompressible fluid flows.

Due to the incompressiblity constraint, it is clear that triangulations at two different times, T^t_hⁿ and T^t_hⁿ⁺¹, must have the same volume. Assume that T^t_hⁿ⁺¹ is obtained from T^t_hⁿ by the ALE map

A^[n+1,n]_h = x + uⁿ⁺¹_h,n , for x ∈ Tⁿ_h,

where uⁿ⁺¹_h,n denotes the displacement field. The displacement field has to be somehow constructed from the fluid velocity field, which is divergence free only on Ωⁿ_h. This proved to be a hard and non–trivial task, and is one of the main drawbacks of the ALE approach for describing the moving mesh methods. This issue is addressed in detail in Chapter 2 where a method for the construction of a volume preserving ALE map is derived. The newly proposed method consists of solving a constrained optimization problem for the mesh displacement field. An artificially derived constraint is constructed from the fluid velocity and the discrete time step, and it ensures the volume preservation.

In the context of ALE framework, the mesh velocitywbhis defined as

wb_h = ∂

∂tx(x, t), x ∈ Ωb _h. (0.2)

Let f : QT → R be an Eulerian field, and bf = f ◦ bAtits ALE counterpart, where QT = {(x, t) | x ∈ Ω(t), t ∈ (0, T )}. The time derivative of an Eulerian field f in the ALE framework, i.e. time derivative of f written from the viewpoint of reference configuration,

(30)

is given by the relation

∂

∂t

bx

f (x, t) = ∂

∂tf (x, t) + w(x, t) · ∇f (x, t). (0.3)

Consider for a moment a generic conservation law for the scalar quantity u:

∂

∂tu − div B(u) = f in Q_T, (0.4)

where B is a first order, linear or non–linear, differential operator. In ALE framework, employing the ALE time derivative defined above, equation (0.4) is rewritten in the form

∂

∂t

xb

u − w ·∇u − div B(u) = f in Q_T. (0.5)

Finite element method is based on the weak formulation of the considered partial differential equation(PDE). Weak formulation of problem (0.5) reads

Z

Ω

ψ ∂

∂t

xb

u − ψ w ·∇u + ψ · B(u) − ψf

dx +

Z

∂Ω

ψ B(u) · n dS = 0, (0.6)

where ψ is a smooth test function, such that _∂t^∂

xbψ = 0. Then, the transient term in equation (0.6) can be rewritten in the conservative form:

Z

Ω

ψ ∂

∂t

xb

u − ψ w ·∇u

dx = d dt

Z

Ω

ψu dx − Z

Ω

ψ div(u w) dx . (0.7)

In the context of ALE framework, weak formulations in which the time derivative is for- mally extracted in front of the integral sign are referred to as ”conservative formulations”.

Such formulations enjoy better conservation properties in the discretized form than their non–conservative counterpart, where the temporal derivative is kept under the integral sign. However, it is also very well known that if temporal discretization is not handled carefully, artificial sinks and/or sources may appear in the discrete scheme. The reason behind this quite problematic and unwanted property lies in the most simple form of the

(31)

doi:10.6342/NTU202003676 Reynolds transport theorem which gives the relation between the volume change and the

domain velocity w:

d dt

Z

K

dx = Z

∂K

div w dx , for any control volume K ⊂ Ω. (0.8)

However, on the discrete level, between two time steps tnand tn+1, above identity holds only approximately in general:

Z

Kⁿ⁺¹

dx − Z

Kⁿ

dx ≈ ∆t Z

K

div w dx, (0.9)

where ∆t = tn+1 − tn, and the right–hand–side in the above equations is evaluated at some point t ∈ [tn+1, t_n]. Hence, the change in volume of K is not preserved in the mesh velocity w for an arbitrary discretization scheme. This issue is investigated in detail in Chapter 3 where a systematic way for constructing SCL preserving time–discretization schemes is developed for PDEs on time–dependent domains within the ALE FEM framework. Most of the original work on this topic has been done in the context of the finite difference method ([7, 8, 9, 10, 11, 12]) and finite volume method ([13, 14, 15, 16, 17, 18]), and, lately, extended to the finite element method ([15, 16, 5, 6, 19, 20]). The material presented in Chapter 3 has already been published by the author in [21], coauthored by T. W. H. Sheu and M. Solovchuk.

Conservation laws considered in this work are typically of a parabolic type. When the character of the system of equations to be solved is of elliptic or parabolic type, yet close to the hyperbolic type, a numerical scheme may produce nonphysical oscillations in the numerical solution if computational mesh is not sufficiently dense ([1]). Typical examples are convection–diffusion (CD) equations with dominating convection term. In these situations, continuous problem is well posed and it has a unique solution based on the Lax–Milgram lemma, yet numerical problem obtained by standard FEM is not stable.

Loss of stability is a consequence of too small coercivity constant of the bilinear form in the weak formulation (see e.g. [3]). Greater insights on these issues as well as some popular techniques on handling them can be found in [3] for the case of problems posed in the

(32)

time–independent (stationary) domains. Essentially, the so called stabilization methods are introduced with the aim of stabilizing the numerical scheme for the parabolic PDEs with dominating convection. Stabilized schemes take the form

C(u_h, ψ_h) + S^%(u_h, ψ_h) = 0

where C(u_h, ψ_h) is the weak formulation of the considered PDE and S^%(u_h, ψ_h) is the stabilization term. Desirable property of the stabilization term is that it vanishes when the exact solution is plugged in. In that case, stabilization scheme is called strongly consistent. In Chapter 4, three popular stabilization methods commonly found in the literature, which are strongly consistent methods in stationary domain scenario, are extended to ALE framework: Streamline Upwind/Petrov Galerkin (SUPG) method (introduced in [22] and extended in [23] for conservative formulations in ALE framework), Galerkin Least Squares(GLS) method (introduced in [24]) and Douglas–Wang/Galerkin (DWG) method (introduced in [25, 26, 27] and occasionally referred to as unusual stabilized finite element method). It has to be noted that extension of stabilization methods for the conservative formulations in ALE framework is no trivial work when one wants to preserve the strong consistency of the method. Yet, it is shown in Chapter 4 that the strong consistency can indeed be achieved in the spirit of approach introduced in Chapter 3.

In chapters 3 and 4 a novel SCL–preserving approach for moving mesh problems is introduced. This approach is then applied in Part 2 for simulating some (academic) multiphysics problems.

In Chapter 5, reconstruction of curvature of discrete surface is investigated. Fluid flow problems in which curvature plays an important role typically include multiphase and multifluid flows. Immiscible multifluid flow problems are typical source of inspira- tion for the moving domain problems. In such problems, surface tension on fluid–fluid interface generates surface force which is a function of the interface curvature which de- pends on the geometry of the interface. Two essentially different techniques have been used to describe the interface in the literature: implicit and explicit. In the implicit approach, a fixed computational mesh is used and an additional scalar field is introduced to

(33)

doi:10.6342/NTU202003676 describe the interface. This approach is often referred to as ”interface capturing” and its

main advantage is that it can easily handle topological changes. It is often combined with mesh adaptation techniques in order to ensure credible interface capturing. Level–set ([28, 29]) and volume–of–fluid ([30, 31]) methods, for example, fall within this approach. In the explicit approach, often referred to as ”interface tracking”, the interface is described explicitly with an aligned mesh i.e. ”mesh fits the interface”. In this environment, when the interface moves, the mesh has to be moved accordingly with it. Lagrange and ALE approaches fall into this category. The most common approximation of geometry in FEM is with linear interpolation functions. This means that interface is approximated with piecewise linear edges in 2D and triangles or quadrangles in 3D. In this case, a popular choice for curvature calculation that can be found in the literature is the higher order interpolation of the interface. It allows the use of the curvature formula that involves second derivatives of the boundary parametrization. A spline interpolation of the interface is reconstructed from the linear computational mesh and it is then used solely for the curvature calculation.

This was, for example, studied in [32] where they used cubic splines and in [33] where non–uniform rational B–splines (NURBS) were used. In [34] the authors used a simple finite difference version of Frenet–Serret formula to calculate curvature and surface tension force in two–phase flow and in [35] the authors computed the curvature for interfacial tension by least squares parabola fitting method. A somewhat different but particularly attractive approach for curvature calculation within interface tracking FEM employs the Laplace–Beltramioperator. It is used in both standard FEM with linear meshes and with isoparametric FEM (both of these approaches are studied in [32]). Laplace–Beltrami operator falls into machinery from the discrete differential geometry where it plays an important role in discrete surface modeling (see e.g. [36]). In the context of FEM in fluid dynamics, the Laplace–Beltrami operator technique was already employed for problems with free surfaces. Weak form can be derived from the mathematical expression for the curvature which involves the Laplace–Beltrami operator. Thus, the curvature can be very easily and naturally incorporated into the FEM formulation using this technique. In [32], they noticed the appearance of spurious oscillations in velocity field due to the introduced

(34)

numerical error in the evaluation of surface tension force. Two different approaches for curvature calculation on interface fitted meshes were studied – cubic interpolation of the interface and Laplace–Beltrami operator technique. It has been reported that the cubic interpolation technique performed much better in general. In Chapter 5, the issue why Laplace–Beltrami operator performs poorly on boundary fitted meshes in general cases is examined and resolved. It turns out that when finite element space is not chosen carefully, the Laplace–Beltrami operators ”viewpoint” of the discrete surface (curve) is ”distorted”.

This results in locally nonphysical oscillations of the curvature which, in turn, introduce the local spurious surface forces.

This concludes Part 1 of this work.

In Part 2, methodology derived in Part 1 is applied to solve complex multiphysics problems. In particular, Part 2 consists of three chapters in which dynamic contact line problem, chemotaxis phenomenon and fluid–structure interaction problem are simulated with ALE FEM methods derived in Part 1. Each problem carries specific problematics which is addressed.

In Chapter 6, a sliding droplet problem is considered. A small liquid droplet is placed on an inclined plane where surface tension, gravity and friction forces compete. When gravity force is stronger the friction forces, droplet starts to move downwards along the plane. The most common approach for describing a viscous fluid flow in contact with some solid surface is to prescribe the so called no-slip boundary condition on the fluid–

solid interface. This condition ensures that the fluid velocity is equal to the solid velocity and, in general, describes the physics of such flows credibly. However, it is well known that the contact line (solid–liqudid–gas interface) is able to move in real world examples.

If one employs no–slip boundary condition, physics of the flow in the numerical simulations is ruined, at least near the contact line. Hence, boundary condition with roots from the molecular dynamics approach has been derived for the continuum modeling approach in [37, 38]. The so called generalized Navier boundary condition (GNBC) credibly describes the fluid behavior near the contact line, and the no–slip boundary condition can be

(35)

doi:10.6342/NTU202003676 derived from the GNBC limiting case. Volume preserving method derived in Chapter 2

and proper curvature evaluation described in Chapter 5 are of particular importance for this problem.

In Chapter 7 bacterial chemotaxis phenomenon is considered. Suspension of an oxytactic bacteria (e.g. the species Bacillus subtilis) placed in a container with its upper surface open to the atmosphere results in the formation of complex bioconvection patterns. The bacteria consume the oxygen diluted in the water, thereby causing the decrease of oxygen concentration everywhere except on the free surface. Through the free surface, which is in direct contact with air, oxygen diffuses into the water. Slightly denser than water, the oxytactic bacteria are able to swim towards the higher concentration of oxygen (i.e. upwards) and they concentrate in a thin layer below the free surface. This causes the change of the suspension density and Rayleigh–Taylor type instabilities to oc- cur. The chemotaxis phenomenon has been successfully modeled in the literature within continuum mechanics approach under the assumption that domain is fixed in time. In this chapter, this model is extended to credibly model the phenomenon when free surface is allowed to move, as is the case in the realistic situation. The chemotaxis phenomenon exhibits the similar behavior as the free thermal convection, which is a well studied problem due its significance in engineering and industry. Hence, the governing system of equations is constrained with fewer assumptions and approximations. For example, the dependence of the surface tension of water on the temperature has been estimated. There- fore, one is able to consider the thermal gradients on the free surface accompanied with the (tangential) Marangoni flows. Similar behavior is expected for the bacterial chemotaxis, however, the physics of the surface tension depending on bacteria concentration is still under the research and hence is neglected. SCL preserving method derived in Chap- ter 3 is of particular importance for the chemotaxis phenomenon since the bacteria has to be preserved at all times. Majority of material presented in Chapter 7 has been published recently by the author in [39, 40], coauthored by T. W. H. Sheu and M. Solovchuk.

In Chapter 8, methods derived in Part 1 are illustrated on fluid–structure interaction (FSI) problems arising from the field of bio–medicine. FSI plays a major role in mathe-

(36)

matical modeling and numerical simulations of the blood flow (in large arteries). Failure to take into account the changes in dynamics of the blood vessels due to change in blood pressure may result in a bad estimation of the wall shear stress. This is particularly important in modeling of aneurysms which are prone to rupture when aneurysm wall is weakened by the effects of shear stress. A monolithic approach for solving the FSI FEM formulation is employed which, when the finite element spaces are appropriately chosen, ensures the implicit (strong) coupling of the boundary conditions on the fluid–structure interface. Monolithic approaches, generally, enjoy very good stability properties. In this regard, they are more desirable than the partitioned approaches which require unstable explicit coupling. FSI is a large field in computational fluid dynamics, and details on derivations, numerical approaches and particular models outreach the scope of this work.

The idea of this chapter is only to demonstrate the ability of adapting the methodology developed in Part 1 for this class of problems.

Finally, in the concluding Chapter 8, the content of the thesis is summarized and conclusions are drawn.

Publication list

Parts of this thesis have already been published as a requirement by National Taiwan Universityfor PhD defense. Publication list is given below:

1 S. Chakraborty, F. Ivanˇci´c, M. Solovchuk, T. W.–H. Sheu: Stability and dynamics of a chemotaxis system with deformed free-surface in a shallow chamber, Phys.

Fluids 30 (2018), 071904.

2 F. Ivanˇci´c, T. W.–H. Sheu, M. Solovchuk: Arbitrary Lagrangian Eulerian–type finite element methods formulation for PDEs on time–dependent domains with vanishing discrete space conservation law, SIAM J. Sci. Comput. 41 (2019), No. 3, pp. A1548–A1573.

3 F. Ivanˇci´c, T. W.–H. Sheu, M. Solovchuk: The free surface effect on a chemotaxis–

diffusion–convection coupling system, Comput. Methods Appl. Mech. Engrg. 356

(37)

doi:10.6342/NTU202003676 (2019), pp. 387–406.

4 F. Ivanˇci´c, T. W.–H. Sheu, M. Solovchuk: Bacterial chemotaxis in thin fluid layers with free surface, Phys. Fluids 32 (2020), 061902.

5 F. Ivanˇci´c, T. W.–H. Sheu, M. Solovchuk: Elimination of spurious velocities gen- erated by curvature dependent surface force in finite element flow simulation with mesh–fitted interface, to appear @ Comput. Methods Appl. Mech. Engrg..

(38)

Part I

FEM approximation of differential

problems within ALE framework

(39)

(40)

CHAPTER 1 Parabolic equations in time–dependent domains

This chapter is devoted to describe the ALE approach for a generic conservation law of the form

∂

∂tu + div B(u) = f in Ω(t), t ∈ (0, T ). (1.1) In equation (1.1), u = u(x, t) is an unknown scalar field representing some physical quantity, for example concentration of some substance or heat. Ω(t) ⊂ R^dis the physical domain at time instant t ∈ (0, T ) in which the equation (1.1) is posed. Notice that the domain Ω is itself possibly a function of time, Ω = Ω(t). Term B(u) is a vector field with d–components which are linear or non–linear functions of u. The most simple example is B(u) = −∇u in which case the equation (1.1) is the classic heat equation and can be written in the form

∂

∂tu − ∆u = f. (1.2)

Another classical example is B(u) = −∇u + u v, where v is a vector field typically representing the fluid velocity. In this case, equation (1.1) is the convection diffusion equation governing the transport of the concentration of some quantity in the fluid. It can

(41)

doi:10.6342/NTU202003676 be written in the form

∂

∂tu − ∆u + v ·∇u + u div v = f. (1.3)

The incompressible Navier–Stokes equations can be viewed as a system of non–linearly coupled convection–diffusion equations subjected to the incompressibility constraint. De- note by v : Ω → R^dand p : Ω → R the velocity and pressure field. Then, in their vector form, the Navier–Stokes equations governing the flow of incompressible Newtonian fluid can be written as

∂

∂tv +(v ·∇) v −∆ v +∇p = f in Ω(t), t ∈ (0, T ), div v = 0 in Ω(t), t ∈ (0, T ).

(1.4)

In this chapter a general ALE approach for equations governing the conservation laws on time–dependent domains is revisited. Most of the content in this chapter is already well known and can be found in the literature. The main idea is to introduce the problematics one faces when employing finite element method on moving meshes and establish the notation to be used throughout this work. A detailed description of the ALE framework can be found in [41] and references therein. The model equations on which the methodology is described are taken in their dimensionless (and normalized) forms for simplicity since generalizations to more specific cases are straightforward.

1.1 The Arbitrary Lagrangian–Eulerian framework

The general idea of the ALE framework consists of an interplay between the fixed referential domain and the current physical domain occupied by the medium. ALE framework is most often employed for fluid flow problems with free boundaries, such as free surface flows. Referential domain most often coincides with the initial domain, but does not necessarily have to. The interplay between these two domains is realized through the so called ALE map which maps the referential domain into the physical one. In order to per- form the necessary calculus, a minimal smoothness of the ALE map has to be demanded – e.g. a kind of inverse of the ALE map has to exist in order to ensure the correspondence

(42)

between the referential domain and the current one (in sense that one can be obtained from the other).

1.1.1 Fundamentals of ALE framework

Let bΩ ⊂ R^d, d = 2, 3, be a fixed referential domain and let Ω ≡ Ω(t) ⊂ R^da physical domain occupied by the fluid. It is assumed that boundaries of the domain are sufficiently smooth – this usually refers to the Lipschitz continuous boundary – and that the domain evolution can be followed through a one–parameter family of mappings ( bA_t)_{t∈[0,T ]}, T <

∞,

Ab_t: bΩ → R^d, t ∈ [0, T ], x 7→ x ,b x ∈ bb Ω , x ∈ Ω(t).

(1.5)

For the sake of compact notation, it is denoted

QT = {(x, t) | x ∈ Ω(t), t ∈ (0, T )}, and Qb_T = {(x, t) |b x ∈ bb Ω(t), t ∈ (0, T )}.

(1.6)

ALE map is often defined as a single vector field on bQ_T, bA : bQ_T → R^d, rather than a one–parameter family of mappings introduced a moment ago. However, slightly abusing the notation, these two terminologies are usually identified for convenience,

( bA_t)_{t∈[0,T ]} ≡ ( bA(·, t))_{t∈[0,T ]}.

Ab_t maps the referential into the physical domain (see Figure 1.1), bΩ 7→ Ω ≡ Ω(t) = Ab_t(bΩ). In this context,x ∈ bb Ω is referred to as the ALE coordinate while x = bA_t(x) ∈b Ω(t) is referred to as an Eulerian (or spatial) coordinate. It is often of interest to write the ALE map bA_tin terms of displacementu =b u(b x, t):b

A(b x, t) =b x +b u(b x, t).b (1.7)

(43)

doi:10.6342/NTU202003676

Ωb Ω(t)

Ab_t

Ab^-1_t

xb x

ub

Figure 1.1: Transformation between two configurations.

More precisely,u(b x, t) is the displacement ofb x at time t,b

u(b x, t) = x −b x , x = bb A_t(bx) , x ∈ bb Ω , t ∈ [0, T ]. (1.8)

Figure 1.1 illustrates the two configurations, reference and physical, and the maps between them.

Let f : Q_T → R andbg : bQ_T → R be two scalar fields defined on the physical and the referential configurations, respectively. Their ALE and Eulerian counterparts are defined respectively by

bf : bQ_T → R , bf = f ◦ bA_t, g : Q_T → R , g =bg ◦ bA^-1_t .

(1.9)

Therefore, the ”hat” operator is just an abbreviation for a composition with the ALE map. Dropping the ”hat” operator on the functions defined on bQ_T is then understood as the composition with the inverse of ALE map. Note that the ”hat” operator notation only makes sense in the context of physical/referential configurations interplay realized through the family of ALE maps. To make a clear difference between functions on physical and referential domains, the ”hat” operator is used to identify functions which ”live”

on the referential domain, while it is dropped for the functions on the physical domain.

The same convention will be employed for (differential) operators: if the operator operates with respect to the referential configuration, this is emphasized by using the ”hat”

symbol, and if it operates with respect to the current configuration, the ”hat” symbol is

(44)

dropped; e.g.

∇ =b ∂

∂xb and ∇ = ∂

∂ x.

Gradient and Jacobian of the ALE map play important roles in the interplay between configurations. The are given by

Fb_t = b∇ bA_t , bJ_t= det( b∇ bA_t). (1.10)

In a strictly mathematical sense, the ALE map is just a family of coordinate transforma- tions. Therefore, the information on volume element between the physical and reference domains is kept in the Jacobian. Specially, from this observation it is easily deduced that J (t) > 0 must hold for all t ∈ [0, T ] for a physically reasonable transformation. Dueb to the regularity of bA_t, it can be shown that bA_t(∂ bΩ) = ∂Ω, i.e. Ab_t maps ”boundary to boundary”. As was the case with the ALE map, the notation throughout this work is slightly abused for Jacobian and gradient of the ALE map. The following is identified

J (b x, t) ≡ bb J_t(x) and bb F (x, t) ≡ bb F_t(x).b

Consider again, for a moment, a fieldbg : bQ_T → R defined on reference configuration and its physical configuration counterpart g : QT → R , g = bg ◦ bA^-1_t . It is not obvious whether regularity of bg in some norm on bΩ implies regularity of g on Ω. A result on sufficient condition on ALE map which ”preserves regularity” of configurations is given in the following proposition.

PROPOSITION 1.1.1 (ALE map regularity condition) Let bΩ be a bounded domain with Lipschitz continuous boundary and let bA_t be a C⁰–diffeomorphism¹ and assume that∀t ∈ (0, T )

(i) Ωt= bAt(bΩ) is bounded and ∂Ω is Lipschitz continuous,

(ii) bAt∈ W^1,∞(bΩ; R^d) and bA^-1_t ∈ W^1,∞(Ω; R^d).

1Differentiable map bAt: bΩ → Ω is called a diffeomorphism if it has a differentiable inverse bA^-1_t . bAtis a C^k–diffeomorphism if bAtand bA^-1_t are k times continuously differentiable.

(45)

doi:10.6342/NTU202003676 Then, g ∈ H¹(Ω) if and only if bg ∈ H¹(bΩ). Furthermore, ∀bg ∈ H¹(bΩ) kgk_H¹_(Ω) is

equivalent tokbg k_H1(bΩ).

1.1.2 The ALE temporal derivative

In the ALE framework, the temporal derivative of an Eulerian field can be considered from different viewpoints. Let f : QT → R be an Eulerian field, and bf = f ◦ bA_tits ALE counterpart. The time derivative of an Eulerian field f in the ALE framework, i.e. time derivative of f written from the viewpoint of reference configuration, is defined as

∂

∂t

xb

f : QT → R , ∂

∂t

xb

f (x, t) = ∂

∂tbf (x, t),b x = bb A^-1_t (x). (1.11)

The time derivative of an Eulerian field in the spatial (physical) framework is just the temporal partial derivative in classical sense,

∂

∂t _x

f (x, t) = ∂

∂tf (x, t). (1.12)

From the practical point of view, discretization of the ALE temporal derivative makes more sense than discretization of the Eulerian temporal derivative on moving mesh. Val- ues of a discrete field f are in correspondance with grid nodes which vary in time but have

”fixed numbering” – therefore, it is possible that a point in control volume at time t is not inside the control volume at time t + ∆t. In that case, the discrete temporal derivative from Eulerian viewpoint doesn’t even make sense.

At this point, one is able to define the domain velocity as

w(x, t) = ∂

∂t

xb

x , x = bA(x, t),b (1.13)

i.e.

w(x, t) =w(b x, t) =b ∂

∂tA(b bx, t) , x = bA(x, t).b (1.14) Note that

∂ub

∂t =wb (1.15)