5.1 Errors in L2–norm in velocity field with respect to the finite element space employed, the technique for numerical evaluation of the curvature and the mesh parameter heΣ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 L2–norm in pressure field with respect to the finite element space
employed, the technique for numerical evaluation of the curvature and the mesh parameter heΣ
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 heΣ
h. Πh denotes the projection to the corre-sponding finite element space. ”LB–embedded” refers to formulations (5.10,5.11) where curvature computation is embedded directly into the formulation. . . 153
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 heΣ
h. Πh denotes the projection to the corre-sponding 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). σ0 denotes the surface tension for pure water (σ0 ≈ 72 mN/ m at 25◦C). . . 192
Introduction
Time–dependent domains appear in many multiphysics models governing various phe-nomena arising from the physics and engineering, commonly related to fluid flow prob-lems. 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
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:
Abt: bΩ → Rd, bΩ 7→ Ω. (0.1)
In the same manner, the discrete referential triangulation bTh(bΩh) is mapped into the phys-ical triangulation T (Ωh) through the discrete ALE map bAh,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, Tthn and Tthn+1, must have the same volume. Assume that Tthn+1 is obtained from Tthn by the ALE map
A[n+1,n]h = x + un+1h,n , for x ∈ Tnh,
where un+1h,n denotes the displacement field. The displacement field has to be somehow constructed from the fluid velocity field, which is divergence free only on Ωnh. 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
wbh = ∂
∂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,
is given by the relation
Consider for a moment a generic conservation law for the scalar quantity u:
∂
∂tu − div B(u) = f in QT, (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
∂
Finite element method is based on the weak formulation of the considered partial differ-ential equation(PDE). Weak formulation of problem (0.5) reads
Z equation (0.6) can be rewritten in the conservative form:
Z
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
doi:10.6342/NTU202003676 Reynolds transport theorem which gives the relation between the volume change and the
domain velocity w:
However, on the discrete level, between two time steps tnand tn+1, above identity holds only approximately in general:
where ∆t = tn+1 − tn, and the right–hand–side in the above equations is evaluated at some point t ∈ [tn+1, tn]. 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 frame-work. 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 mate-rial 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 pop-ular techniques on handling them can be found in [3] for the case of problems posed in the
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(uh, ψh) + S%(uh, ψh) = 0
where C(uh, ψh) is the weak formulation of the considered PDE and S%(uh, ψ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 con-sistent. In Chapter 4, three popular stabilization methods commonly found in the litera-ture, 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 conser-vative 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 ap-proach, a fixed computational mesh is used and an additional scalar field is introduced to
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 ap-proaches fall into this category. The most common approximation of geometry in FEM is with linear interpolation functions. This means that interface is approximated with piece-wise 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 ten-sion 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
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 lems. In particular, Part 2 consists of three chapters in which dynamic contact line prob-lem, 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 simula-tions 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 de-scribes the fluid behavior near the contact line, and the no–slip boundary condition can be
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 oxy-tactic 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 pat-terns. 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 sur-face, 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 prob-lem 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 chemo-taxis, 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 bacChap-teria 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-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 im-portant 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.