**Three-Dimensional Simulation of Fluid–Structure Interaction Problems Using Monolithic Semi-Implicit Algorithm**

#### **Cornel Marius Murea**

Département de Mathématiques, IRIMAS, Université de Haute Alsace, 6, rue des Frères Lumière, 68093 Mulhouse CEDEX, France; cornel.murea@uha.fr

Received: 17 March 2019; Accepted: 13 May 2019; Published: 22 May 2019

**Abstract:** A monolithic semi-implicit method is presented for three-dimensional simulation of fluid–structure interaction problems. The updated Lagrangian framework is used for the structure modeled by linear elasticity equation and, for the fluid governed by the Navier–Stokes equations, we employ the Arbitrary Lagrangian Eulerian method. We use a global mesh for the fluid–structure domain where the fluid–structure interface is an interior boundary. The continuity of velocity at the interface is automatically satisfied by using globally continuous finite element for the velocity in the fluid–structure mesh. The method is fast because we solve only a linear system at each time step. Three-dimensional numerical tests are presented.

**Keywords:** fluid–structure interaction; monolithic method; Updated Lagrangian; Arbitrary Lagrangian Eulerian

#### **1. Introduction**

There exists a rich literature on solving numerically fluid–structure interaction problems. Some methods are based on partitioned procedures, the fluid and structure sub-problems are solved separately using iterative process: fixed point iterations [1–3], Newton-like methods [4–6] or optimization techniques [7–9]. Monolithic methods solve the fluid–structure interaction problem as a single system of equations, [10–13], or more recently [14–17].

In some methods such as the Arbitrary Lagrangian Eulerian (ALE) framework, the fluid equations are written over a moving mesh which follows the structure displacement (see [18,19]). Other methods use a fixed mesh for fluid domain: immersed boundary method [20], distributed Lagrange multiplier [21,22], penalization [23,24], extended finite element method (XFEM) [25,26], and Nitsche-XFEM [27]. Distributed Lagrange multiplier strategy with remeshing is used in [28] and a monolithic fictitious domain without Lagrange multiplier is employed in [29,30].

Most of these methods are implicit. For a long time, the explicit methods were considered not suitable because of the lack of stability, but these methods are successfully applied in [31,32]. A third class of methods are so-called semi-implicit methods, where the domain is computed explicitly while the fluid velocity and pressure as well as the structure displacement are computed implicitly, [33,34]. In [35], it is proved that a schema of this kind is unconditionally stable.

In this paper, a monolithic semi-implicit method is employed for three-dimensional simulation. For the structure modeled by the linear elasticity equations, we use the updated Lagrangian framework and, for the fluid governed by the Navier–Stokes equations, we employ the ALE method. A similar strategy is used in [36] for a bi-dimensional compressible neo-Hookean model or in [37] for a bi-dimensional linear elasticity model for the structure. As in [38], we employ a global mesh for the fluid–structure domain where the fluid–structure interface is an interior boundary. Using globally continuous finite element for the velocity in the fluid–structure mesh, the continuity of velocity at the

interface is automatically satisfied. The method is fast because we solve only a linear system at each time step. Three-dimensional numerical tests are presented.

In [14–17], a global mesh obtained from the deformed structure mesh and a fluid mesh generated at each time step, compatible at the interface with the structure mesh are used. Remeshing the fluid domain improves the quality of the mesh in the case of large deformation. The non-linear structure equation written in the Eulerian coordinates is obtained by using Cayley–Hamilton theorem. The fluid equations are solved by the characteristics method. The weak formulation of the fluid–structure interaction problem is written in the Eulerian domain, which is unknown, and a fixed-point algorithm solves the global non-linear problem at each time step.

In [30], it is assumed that the structure is viscoelastic with the same viscosity as the fluid. Based on fictitious domain without Lagrange multiplier, the fluid is solved in a fixed mesh of the fluid–structure domain. The weak formulation contains integrals over the unknown Eulerian domain of the structure. At each time step, a fixed-point algorithm is employed.

In [36], by using the Updated Lagrangian framework for a compressible Neo-Hookean structure, the weak formulation is written in the known configuration obtained at the precedent time step. By linearization around this configuration, at each time step, only a linear system is solved for the fluid–structure coupled equations and consequently the computing time is reduced. In the present paper, we follow this approach for three-dimensional simulations using linear elastic model for the structure.

If at each time step of the monolithic implicit methods, the fixed-point algorithm does not converge quickly or the computational time by fixed-point iteration is very expensive, thus the monolithic semi-implicit methods, which have similar stability properties and a reduced computational time, could be a good alternative.

#### **2. Problem Statement**

The initial fluid domain Ω*<sup>F</sup>* <sup>0</sup> is a right circular cylinder of bases Σ1, Σ<sup>3</sup> and lateral surface Γ0. We denote by Ω*<sup>S</sup>* <sup>0</sup> the initial structure domain and we assume that it is a right annular cylinder of bases Σ5, Σ7, interior lateral surface Γ<sup>0</sup> and exterior lateral surface Γ*<sup>N</sup>* <sup>0</sup> (see Figure 1). We suppose that the initial structure domain is undeformed (stress-free). The boundary Γ<sup>0</sup> is common of both domains and it represents the initial position of the fluid–structure interface.

**Figure 1.** Initial geometrical configuration.

At the time instant *t*, the fluid occupies the domain Ω*<sup>F</sup> <sup>t</sup>* bounded by the moving interface Γ*<sup>t</sup>* and by the fixed boundaries Σ1, Σ3, while the structure occupies the domain Ω*<sup>S</sup> <sup>t</sup>* bounded by the moving lateral surfaces Γ*t*, Γ*<sup>N</sup> <sup>t</sup>* and by the fixed boundaries Σ5, Σ7.

We denote by **U***<sup>S</sup>* : Ω*<sup>S</sup>* <sup>0</sup> <sup>×</sup> [0, *<sup>T</sup>*] <sup>→</sup> <sup>R</sup><sup>3</sup> the displacement of the structure. A particle of the structure whose initial position was the point **X** will occupy the position **x** = **X** + **U***<sup>S</sup>* (**X**, *t*) in the deformed domain Ω*<sup>S</sup> <sup>t</sup>* . At the time instant *<sup>t</sup>*, the interface <sup>Γ</sup>*<sup>t</sup>* is the image of <sup>Γ</sup><sup>0</sup> via the map **<sup>X</sup>** <sup>→</sup> **<sup>X</sup>** <sup>+</sup> **<sup>U</sup>***<sup>S</sup>* (**X**, *<sup>t</sup>*). The same relationship exists between Γ*<sup>N</sup> <sup>t</sup>* and Γ*<sup>N</sup>* <sup>0</sup> . On <sup>Γ</sup>*<sup>D</sup>* <sup>0</sup> = Σ<sup>5</sup> ∪ Σ7, we impose zero displacements.

We set <sup>∇</sup>**XU***<sup>S</sup>* the gradient of the displacement **<sup>U</sup>***<sup>S</sup>* <sup>=</sup> - *U<sup>S</sup>* <sup>1</sup> , *<sup>U</sup><sup>S</sup>* <sup>2</sup> , *<sup>U</sup><sup>S</sup>* 3 *<sup>T</sup>* with respect to the Lagrangian coordinates **X** = (*X*1, *X*2, *X*3) *<sup>T</sup>*. We denote by **<sup>F</sup>** (**X**, *<sup>t</sup>*) <sup>=</sup> **<sup>I</sup>** <sup>+</sup> <sup>∇</sup>**XU***<sup>S</sup>* (**X**, *<sup>t</sup>*) the gradient of the deformation, where **I** is the unity matrix and we set *J* (**X**, *t*) = det **F** (**X**, *t*). The second Piola–Kirchhoff stress tensor is denoted by **Σ**.

We assume that the fluid is governed by the Navier–Stokes equations. For each time instant *<sup>t</sup>* <sup>∈</sup> [0, *<sup>T</sup>*], we denote the fluid velocity by **<sup>v</sup>***F*(*t*) = - *vF* <sup>1</sup> (*t*), *<sup>v</sup><sup>F</sup>* <sup>2</sup> (*t*), *<sup>v</sup><sup>F</sup>* <sup>3</sup> (*t*) *<sup>T</sup>* : Ω*<sup>F</sup> <sup>t</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> and the fluid pressure by *pF*(*t*) : Ω*<sup>F</sup> <sup>t</sup>* <sup>→</sup> <sup>R</sup>. Let *-* - **v***F* = <sup>1</sup> 2 <sup>∇</sup>**v***<sup>F</sup>* <sup>+</sup> - ∇**v***FT* be the fluid rate of strain tensor and let *<sup>σ</sup><sup>F</sup>* <sup>=</sup> <sup>−</sup>*pF***<sup>I</sup>** <sup>+</sup> <sup>2</sup>*μS-* - **v***F* be the fluid stress tensor. To simplify the notation, we write <sup>∇</sup>**v***<sup>F</sup>* in place of <sup>∇</sup>**xv***F*, when the gradients are computed with respect to the Eulerian coordinates **<sup>x</sup>**.

The problem is to find the structure displacement **U***S*, the fluid velocity **v***<sup>F</sup>* and the fluid pressure *p<sup>F</sup>* such that:

$$
\rho\_0^S(\mathbf{X}) \frac{\partial^2 \mathbf{U}^S}{\partial t^2}(\mathbf{X}, t) - \nabla\_\mathbf{X} \cdot (\mathbf{F} \mathbf{\Sigma}) \left(\mathbf{X}, t\right)
\quad = \quad \rho\_0^S(\mathbf{X}) \mathbf{g}\_r \quad \text{in } \Omega\_0^S \times (0, T) \tag{1}
$$

$$\mathbf{U}^S(\mathfrak{X}, t) \quad = \quad 0, \quad \text{on } \Gamma^D\_0 \times (0, T) \tag{2}$$

$$\left(\mathbf{F}\mathbf{\dot{z}}\right)\left(\mathbf{X},t\right)\mathbf{N}^S\left(\mathbf{X}\right) \quad = \quad 0, \quad \text{on } \Gamma\_0^N \times \left(0,T\right) \tag{3}$$

$$\rho^F \left( \frac{\partial \mathbf{v}^F}{\partial t} + (\mathbf{v}^F \cdot \nabla) \mathbf{v}^F \right) - 2\mu^F \nabla \cdot \boldsymbol{\epsilon} \left( \mathbf{v}^F \right) + \nabla \left. p^F \right|\_{} = \rho^F \mathbf{g}\_{\prime} \,\forall t \in (0, T), \forall \mathbf{x} \in \Omega^F\_t \tag{4}$$

$$\nabla \cdot \mathbf{v}^F \quad = \quad 0, \; \forall \mathbf{t} \in (0, T), \forall \mathbf{x} \in \Omega^F\_t \tag{5}$$

$$
\sigma^F \mathbf{n}^F \quad = \quad \mathbf{h}\_{in\prime} \quad \text{on } \Sigma\_1 \times (0, T) \tag{6}
$$

$$
\sigma^F \mathbf{n}^F \quad = \quad \mathbf{h}\_{\text{out}} \quad \text{on } \Sigma\_{\mathcal{B}} \times (0, T) \tag{7}
$$

$$\mathbf{v}^{\mathcal{F}}\left(\mathbf{X} + \mathbf{U}^{\mathcal{S}}\left(\mathbf{X}, t\right), t\right) \\ \quad = \left.\frac{\partial \mathbf{U}^{\mathcal{S}}}{\partial t}\left(\mathbf{X}, t\right), \text{ on } \Gamma\_0 \times (0, T) \tag{8}$$

$$\left(\sigma^{F}\mathbf{n}^{F}\right)\_{\left(\mathbf{X}+\mathbf{U}^{S}(\mathbf{X},t),t\right)} = -\left(\mathbf{F}\mathbf{\bar{\omega}}\right)\left(\mathbf{X},t\right)\mathbf{N}^{S}\left(\mathbf{X}\right),\mathbf{on}\,\Gamma\_{0}\times\left(0,T\right) \tag{9}$$

$$\mathbf{U}^{\mathcal{S}}(\mathfrak{X},0) \;=\;\mathsf{U}^{\mathcal{S},0}(\mathfrak{X})\;,\;\text{in }\Omega\_0^{\mathcal{S}}\tag{10}$$

$$\frac{\partial \mathbf{U}^S}{\partial t}(\mathbf{X}, 0) \quad = \quad \mathbf{V}^{S,0}\left(\mathbf{X}\right) \text{, in } \Omega\_0^S \tag{11}$$

$$\mathbf{v}^{\mathrm{F}}\left(\mathfrak{X},0\right) \;=\;\mathrm{v}^{\mathrm{F},0}\left(\mathfrak{X}\right)\;\mathrm{in}\,\Omega\_{0}^{\mathrm{F}}\tag{12}$$

where *ρ<sup>S</sup>* <sup>0</sup> : <sup>Ω</sup>*<sup>S</sup>* <sup>0</sup> <sup>→</sup> <sup>R</sup> is the initial mass density of the structure, **<sup>g</sup>** is the acceleration of gravity vector and it is assumed to be constant, **N***<sup>S</sup>* is the unit outer normal vector along the boundary *∂*Ω*<sup>S</sup>* <sup>0</sup> , **<sup>U</sup>***S*,0 and **V***S*,0 are the initial displacement and velocity of the structure, *ρ<sup>F</sup>* > 0 and *μ<sup>F</sup>* > 0 are constants representing the mass density and the viscosity of the fluid, **h***in* and **h***out* are prescribed boundary stress, **n***<sup>F</sup>* is the unit outer normal vector along the boundary *∂*Ω*<sup>F</sup> <sup>t</sup>* , and **v***F*,0 is the initial velocity of the fluid.

#### **3. Updated Lagrangian Framework for the Structure Approximation**

Introducing **V***S*, the velocity of the structure in the Lagrangian coordinates, Equation (1) can be rewritten as

$$\rho\_0^S(\mathbf{X}) \frac{\partial \mathbf{V}^S}{\partial t}(\mathbf{X}, t) - \nabla\_\mathbf{X} \cdot (\mathbf{F} \mathbf{\Sigma})(\mathbf{X}, t) \quad = \quad \rho\_0^S(\mathbf{X}) \mathbf{g}\_\prime \quad \text{in } \Omega\_0^S \times (0, T) \tag{13}$$

$$\frac{\partial \mathbf{U}^S}{\partial t} \left( \mathbf{X}, t \right) \quad = \quad \mathbf{V}^S \left( \mathbf{X}, t \right), \quad \text{in } \Omega\_0^S \times (0, T). \tag{14}$$

Let *<sup>N</sup>* <sup>∈</sup> <sup>N</sup><sup>∗</sup> be the number of time steps and <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>T</sup>*/*<sup>N</sup>* the time step. We set *tn* <sup>=</sup> *<sup>n</sup>*Δ*<sup>t</sup>* for *n* = 0, 1, ... , *N*. Let **V***S*,*<sup>n</sup>* (**X**) and **U***S*,*<sup>n</sup>* (**X**) be approximations of **V***<sup>S</sup>* (**X**, *tn*) and **U***<sup>S</sup>* (**X**, *tn*). In the sequel,

$$\mathbf{F}^{\mathrm{n}} = \mathbf{I} + \nabla\_{\mathbf{X}} \mathbf{U}^{\mathrm{S,n}}, \quad \mathbf{\Sigma}^{\mathrm{n}} = \mathbf{\Sigma}(\mathbf{F}^{\mathrm{n}}), \ n \ge 0.$$

Using the implicit Euler scheme, we approach the system in Equations (13) and (14) by

$$\rho\_0^\mathbb{S}\left(\mathbf{X}\right)\frac{\mathbf{V}^{\mathbb{S},n+1}\left(\mathbf{X}\right)-\mathbf{V}^{\mathbb{S},n}\left(\mathbf{X}\right)}{\Delta t}-\nabla\_\mathbf{X}\cdot\left(\mathbf{F}^{n+1}\Sigma^{n+1}\right)\left(\mathbf{X}\right) \\ = \quad \rho\_0^\mathbb{S}\left(\mathbf{X}\right)\mathbf{g}\_\nu\text{ in }\Omega\_0^\mathbb{S}\tag{15}$$

$$\frac{\mathbf{U}^{\mathbb{S},n+1}\left(\mathbf{X}\right)-\mathbf{U}^{\mathbb{S},n}\left(\mathbf{X}\right)}{\Delta t} \\ = \quad \mathbf{V}^{\mathbb{S},n+1}\left(\mathbf{X}\right), \text{ in }\Omega\_0^\mathbb{S}\tag{16}$$

Using Equation (16), we get **<sup>F</sup>***n*+<sup>1</sup> <sup>=</sup> **<sup>F</sup>***<sup>n</sup>* <sup>+</sup> <sup>Δ</sup>*t*∇**XV***<sup>S</sup>*,*n*+<sup>1</sup> and, consequently, **<sup>F</sup>***n*+<sup>1</sup> and **<sup>Σ</sup>***n*+<sup>1</sup> depend on the velocity **V***S*,*n*<sup>+</sup>1. We have eliminated the unknown displacement and we have now an equation of unknown **V***S*,*n*<sup>+</sup>1.

We can put Equation (15) in a weak form: find **V***S*,*n*+<sup>1</sup> : Ω*<sup>S</sup>* <sup>0</sup> <sup>→</sup> <sup>R</sup>3, **<sup>V</sup>***S*,*n*+<sup>1</sup> <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* <sup>0</sup> , such that

$$\begin{split} \int\_{\Omega\_{0}^{S}} \rho\_{0}^{S} \frac{\mathbf{V}^{S,n+1} - \mathbf{V}^{S,n}}{\Delta t} \cdot \mathbf{W}^{S} \, d\mathbf{X} + \int\_{\Omega\_{0}^{S}} \mathbf{F}^{n+1} \boldsymbol{\Sigma}^{n+1} : \nabla\_{\mathbf{X}} \mathbf{W}^{S} \, d\mathbf{X} \\ = \int\_{\Omega\_{0}^{S}} \rho\_{0}^{S} \mathbf{g} \cdot \mathbf{W}^{S} \, d\mathbf{X} + \int\_{\Gamma\_{0}} \mathbf{F}^{n+1} \boldsymbol{\Sigma}^{n+1} \mathbf{N}^{S} \cdot \mathbf{W}^{S} \, d\mathbf{S} \end{split} \tag{17}$$

for all **W***<sup>S</sup>* : Ω*<sup>S</sup>* <sup>0</sup> <sup>→</sup> <sup>R</sup>3, **<sup>W</sup>***<sup>S</sup>* <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* <sup>0</sup> . For instance, we have assumed that the forces **<sup>F</sup>***n*+1**Σ***n*+1**N***<sup>S</sup>* on the interface Γ<sup>0</sup> are known.

The above formulation is also called the total Lagrangian framework, since the equations are written in the undeformed domain Ω*<sup>S</sup>* <sup>0</sup> . Now, we present the updated Lagrangian framework, where the equations are written in the domain Ω*<sup>S</sup> <sup>n</sup>*. We follow a similar method as in [36], where the structure is a bi-dimensional compressible neo-Hookean material, or in [37], where the bi-dimensional linear elasticity model is used.

We set Ω*<sup>S</sup> <sup>n</sup>* the image of Ω*<sup>S</sup>* <sup>0</sup> via the map **<sup>X</sup>** <sup>→</sup> **<sup>X</sup>** <sup>+</sup> **<sup>U</sup>***S*,*<sup>n</sup>* (**X**) and we denote by <sup>Ω</sup> *<sup>S</sup>* <sup>=</sup> <sup>Ω</sup>*<sup>S</sup> <sup>n</sup>* the computational domain for the structure. The map from Ω*<sup>S</sup>* <sup>0</sup> to <sup>Ω</sup>*<sup>S</sup> <sup>n</sup>*+<sup>1</sup> defined by **<sup>X</sup>** <sup>→</sup> **<sup>x</sup>** <sup>=</sup> **<sup>X</sup>** <sup>+</sup> **<sup>U</sup>***S*,*n*+<sup>1</sup> (**X**) is the composition of the map from Ω*<sup>S</sup>* <sup>0</sup> to <sup>Ω</sup> *<sup>S</sup>* given by **<sup>X</sup>** <sup>→</sup> **<sup>x</sup>** <sup>=</sup> **<sup>X</sup>** <sup>+</sup> **<sup>U</sup>***S*,*<sup>n</sup>* (**X**) with the map from <sup>Ω</sup> *<sup>S</sup>* to Ω*<sup>S</sup> <sup>n</sup>*+<sup>1</sup> given by

$$
\widehat{\mathbf{x}} \to \mathbf{x} = \widehat{\mathbf{x}} + \mathbf{U}^{\mathrm{S}, \mathbb{W}+1} \left( \mathbf{X} \right) - \mathbf{U}^{\mathrm{S}, \mathbb{W}} \left( \mathbf{X} \right) = \widehat{\mathbf{x}} + \widehat{\mathbf{u}} \left( \widehat{\mathbf{x}} \right) \dots
$$

Using the notations **<sup>F</sup>** <sup>=</sup> **<sup>I</sup>** <sup>+</sup> <sup>∇</sup> **xu** and *<sup>J</sup>* <sup>=</sup> det **<sup>F</sup>** , *<sup>J</sup><sup>n</sup>* <sup>=</sup> det **<sup>F</sup>***n*, we get

$$\mathbf{F}^{\mathbb{n}+1}\left(\mathbf{X}\right) = \widehat{\mathbf{F}}\left(\widehat{\mathbf{x}}\right)\mathbf{F}^{\mathbb{n}}\left(\mathbf{X}\right), \quad \mathsf{J}^{\mathbb{n}+1}\left(\mathbf{X}\right) = \widehat{\mathsf{J}}\left(\widehat{\mathbf{x}}\right)\mathbf{J}^{\mathbb{n}}\left(\mathbf{X}\right). \tag{18}$$

We have the relation between the Cauchy stress tensor of the structure *σ<sup>S</sup>* and the second Piola–Kirchhoff stress tensor **Σ**,

$$\boldsymbol{\sigma}^{\mathcal{S}}\left(\mathbf{x},t\right) = \left(\frac{1}{\mathcal{J}} \mathbf{F} \mathbf{E} \mathbf{F}^{\mathrm{T}}\right)\left(\mathbf{X},t\right), \quad \mathbf{x} = \mathbf{X} + \mathbf{U}^{\mathcal{S}}\left(\mathbf{X},t\right).$$

The mass conservation assumption gives *<sup>ρ</sup><sup>S</sup>* (**x**, *<sup>t</sup>*) <sup>=</sup> *<sup>ρ</sup><sup>S</sup>* <sup>0</sup> (**X**) *<sup>J</sup>*(**X**,*t*), where *<sup>ρ</sup><sup>S</sup>* (**x**, *<sup>t</sup>*) is the mass density of the structure in the Eulerian framework.

For the semi-discrete scheme, we use the notations

$$\sigma^{\mathbf{S}, \mathbf{n}+1}\left(\mathbf{x}\right) = \left(\frac{1}{J^{\mathbf{n}+1}} \mathbf{F}^{\mathbf{n}+1} \mathbf{L}^{\mathbf{n}+1} \left(\mathbf{F}^{\mathbf{n}+1}\right)^{T}\right)(\mathbf{X}) , \quad \mathbf{x} = \mathbf{X} + \mathbf{U}^{\mathbf{S}, \mathbf{n}+1}\left(\mathbf{X}\right)$$

and *<sup>ρ</sup>S*,*<sup>n</sup>* ( **<sup>x</sup>**) <sup>=</sup> *<sup>ρ</sup><sup>S</sup>* <sup>0</sup> (**X**) *<sup>J</sup>n*(**X**) , **<sup>x</sup>** <sup>=</sup> **<sup>X</sup>** <sup>+</sup> **<sup>U</sup>***S*,*<sup>n</sup>* (**X**).

Let us introduce **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> and **<sup>v</sup>***S*,*<sup>n</sup>* : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> defined by **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> ( **<sup>x</sup>**) <sup>=</sup> **<sup>V</sup>***S*,*n*+<sup>1</sup> (**X**) and **<sup>v</sup>***S*,*<sup>n</sup>* ( **<sup>x</sup>**) <sup>=</sup> **<sup>V</sup>***S*,*<sup>n</sup>* (**X**). In addition, for **<sup>W</sup>***<sup>S</sup>* : <sup>Ω</sup>*<sup>S</sup>* <sup>0</sup> <sup>→</sup> <sup>R</sup>3, we define **<sup>w</sup>** *<sup>S</sup>* : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> and **<sup>w</sup>***<sup>S</sup>* : <sup>Ω</sup>*<sup>S</sup> <sup>n</sup>*+<sup>1</sup> <sup>→</sup> <sup>R</sup><sup>3</sup> by **<sup>w</sup>** *<sup>S</sup>* ( **<sup>x</sup>**) <sup>=</sup> **<sup>w</sup>***<sup>S</sup>* (**x**) <sup>=</sup> **<sup>W</sup>***<sup>S</sup>* (**X**).

Now, we rewrite Equation (17) over the domain <sup>Ω</sup> *<sup>S</sup>*. For the first term of Equation (17), we get

$$\int\_{\Omega\_0^S} \rho\_0^S \frac{\mathbf{V}^{\rm S, n+1} - \mathbf{V}^{\rm S, n}}{\Delta t} \cdot \mathbf{W}^{\rm S} \, d\mathbf{X} = \int\_{\hat{\Omega}^S} \rho^{\rm S, n} \frac{\hat{\mathbf{v}}^{\rm S, n+1} - \mathbf{v}^{\rm S, n}}{\Delta t} \cdot \hat{\mathbf{w}}^{\rm S} \, d\hat{\mathbf{x}}$$

and similarly

$$\int\_{\Omega\_0^S} \rho\_0^S \mathbf{g} \cdot \mathbf{W}^S \, d\mathbf{X} = \int\_{\widehat{\Omega}^S} \rho^{S,\mathfrak{n}} \mathbf{g} \cdot \widehat{\mathbf{w}}^S \, d\widehat{\mathbf{x}}.$$

Using the identity - <sup>∇</sup>**w***<sup>S</sup>* (**x**) **<sup>F</sup>***n*+<sup>1</sup> (**X**) <sup>=</sup> <sup>∇</sup>**XW***<sup>S</sup>* (**X**) and the definition of *<sup>σ</sup>S*,*n*<sup>+</sup>1, we get

$$\int\_{\Omega\_0^S} \mathbf{F}^{n+1} \mathbf{L}^{n+1} : \nabla\_{\mathbf{X}} \mathbf{W}^S \, d\mathbf{X} = \int\_{\Omega\_{n+1}^S} \sigma^{S,n+1} : \nabla \mathbf{w}^S \, d\mathbf{x}.$$

Details about this kind of transformation can be found in [39], Chapter 1.2. To write the above integral over the domain <sup>Ω</sup> *<sup>S</sup>*, let us introduce the tensor

$$
\hat{\mathbf{E}}\left(\hat{\mathbf{x}}\right) = \hat{f}\left(\hat{\mathbf{x}}\right)\hat{\mathbf{F}}^{-1}\left(\hat{\mathbf{x}}\right)\sigma^{S,n+1}\left(\mathbf{x}\right)\hat{\mathbf{F}}^{-T}\left(\hat{\mathbf{x}}\right).
\tag{19}
$$

Since - <sup>∇</sup>**w***<sup>S</sup>* (**x**) **<sup>F</sup>** ( **<sup>x</sup>**) <sup>=</sup> <sup>∇</sup> **xw** *<sup>S</sup>* ( **<sup>x</sup>**) (see [39], Chapter 1.2) and taking into account Equation (19), we get

$$\int\_{\Omega\_{n+1}^{\mathbb{S}}} \sigma^{\mathbb{S},n+1} : \nabla \mathbf{w}^{\mathbb{S}} \, d\mathbf{x} = \int\_{\widehat{\Omega}^{\mathbb{S}}} \widehat{\mathbf{F}} \widehat{\mathbf{z}} : \nabla\_{\widehat{\mathbf{x}}} \widehat{\mathbf{w}}^{\mathbb{S}} \, d\widehat{\mathbf{x}}.$$

Now, it is possible to present the updated Lagrangian version of Equation (17). Knowing **U***S*,*<sup>n</sup>* : Ω*<sup>S</sup>* <sup>0</sup> <sup>→</sup> <sup>R</sup>3, <sup>Ω</sup> *<sup>S</sup>* <sup>=</sup> <sup>Ω</sup>*<sup>S</sup> <sup>n</sup>* and **<sup>v</sup>***S*,*<sup>n</sup>* : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup>3, we try to find **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup>3, **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* <sup>0</sup> such that

$$\begin{split} & \int\_{\widehat{\Omega}^{S}} \rho^{\widehat{\mathbf{S}},\boldsymbol{\mu}} \frac{\widehat{\mathbf{v}}^{\widehat{\mathbf{S}},\boldsymbol{\mu}+1} - \mathbf{v}^{S,\boldsymbol{\mu}}}{\Delta t} \cdot \widehat{\mathbf{w}}^{S} \, d\widehat{\mathbf{x}} + \int\_{\widehat{\Omega}^{S}} \widehat{\mathbf{F}} \widehat{\mathbf{E}} : \nabla\_{\widehat{\mathbf{x}}} \widehat{\mathbf{w}}^{S} \, d\widehat{\mathbf{x}} \\ & \qquad = \int\_{\widehat{\Omega}^{S}} \rho^{\widehat{\mathbf{S}},\boldsymbol{\mu}} \mathbf{g} \cdot \widehat{\mathbf{w}}^{S} \, d\widehat{\mathbf{x}} + \int\_{\Gamma\_{0}} \mathbf{F}^{n+1} \boldsymbol{\Sigma}^{n+1} \mathbf{N}^{S} \cdot \mathbf{W}^{S} \, d\mathbf{S} \end{split} \tag{20}$$

for all **<sup>w</sup>** *<sup>S</sup>* : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup>3, **<sup>w</sup>** *<sup>S</sup>* <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* <sup>0</sup> . We recall that the forces **<sup>F</sup>***n*+1**Σ***n*+1**N***<sup>S</sup>* on the interface <sup>Γ</sup><sup>0</sup> are assumed known.

Using the identity **<sup>u</sup>** ( **<sup>x</sup>**) <sup>=</sup> **<sup>U</sup>***S*,*n*+<sup>1</sup> (**X**) <sup>−</sup> **<sup>U</sup>***S*,*<sup>n</sup>* (**X**) <sup>=</sup> <sup>Δ</sup>*<sup>t</sup>* **<sup>V</sup>***S*,*n*+<sup>1</sup> (**X**) <sup>=</sup> <sup>Δ</sup>*<sup>t</sup>* **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> ( **<sup>x</sup>**), we obtain

$$
\widehat{\mathbf{F}} = \mathbf{I} + \Delta t \nabla\_{\widehat{\mathbf{x}}} \widehat{\mathbf{v}}^{S, n+1}.\tag{21}
$$

Using Equations (18) and (19), it follows that

$$
\hat{\mathbf{E}}^{\top} = \hat{\mathbf{J}} \hat{\mathbf{F}}^{-1} \sigma^{\mathbf{S}, n+1} \hat{\mathbf{F}}^{-T} = \hat{\mathbf{J}} \hat{\mathbf{F}}^{-1} \frac{1}{l^{n+1}} \mathbf{F}^{n+1} \mathbf{E}^{n+1} \left(\mathbf{F}^{n+1}\right)^{T} \hat{\mathbf{F}}^{-T}
$$

$$
= \frac{1}{l^{n}} \mathbf{F}^{n} \boldsymbol{\Sigma}^{n+1} \left(\mathbf{F}^{n}\right)^{T}.\tag{22}
$$

For the linear elastic material, we have

$$\boldsymbol{\Sigma}(\mathbf{U}) = \boldsymbol{\lambda}^S(\nabla\_{\mathbf{X}} \cdot \mathbf{U}) + \boldsymbol{\mu}^S \left(\nabla\_{\mathbf{X}} \mathbf{U} + (\nabla\_{\mathbf{X}} \mathbf{U})^T\right),$$

where *λ<sup>S</sup>* and *μ<sup>S</sup>* are the Lamé coefficients. We have

$$
\Delta^{n+1} = \Sigma(\mathbf{U}^{\mathbb{S}, n+1}) = \Sigma(\mathbf{U}^{\mathbb{S}, n}) + (\Delta t)\Sigma(\mathbf{V}^{\mathbb{S}, n+1}) = \Sigma^n + (\Delta t)\Sigma(\mathbf{V}^{\mathbb{S}, n+1}).
$$

From Equations (21) and (22) and the above equality, we get

**<sup>F</sup> <sup>Σ</sup>** <sup>=</sup> <sup>1</sup> *<sup>J</sup><sup>n</sup>* **<sup>F</sup>***n***Σ***<sup>n</sup>* (**F***n*) *<sup>T</sup>* <sup>+</sup> <sup>Δ</sup>*t*<sup>∇</sup> **xv** *<sup>S</sup>*,*n*+<sup>1</sup> <sup>1</sup> *<sup>J</sup><sup>n</sup>* **<sup>F</sup>***n***Σ***<sup>n</sup>* (**F***n*) *T* + Δ*t <sup>J</sup><sup>n</sup>* **<sup>F</sup>***n***Σ**(**V***S*,*n*+1)(**F***n*) *<sup>T</sup>* <sup>+</sup> (Δ*t*)<sup>2</sup> *<sup>J</sup><sup>n</sup>* <sup>∇</sup> **xv** *<sup>S</sup>*,*n*+1**F***n***Σ**(**V***S*,*n*+1)(**F***n*) *T* <sup>=</sup> *<sup>σ</sup>S*,*<sup>n</sup>* <sup>+</sup> <sup>Δ</sup>*t*<sup>∇</sup> **xv** *<sup>S</sup>*,*n*+1*σS*,*<sup>n</sup>* + Δ*t <sup>J</sup><sup>n</sup>* **<sup>F</sup>***n***Σ**(**V***S*,*n*+1)(**F***n*) *<sup>T</sup>* <sup>+</sup> (Δ*t*)<sup>2</sup> *<sup>J</sup><sup>n</sup>* <sup>∇</sup> **xv** *<sup>S</sup>*,*n*+1**F***n***Σ**(**V***S*,*n*+1)(**F***n*) *T* .

We introduce **<sup>Σ</sup> <sup>x</sup>**(**<sup>u</sup>** ) = *<sup>λ</sup>S*(<sup>∇</sup> **<sup>x</sup>** · **<sup>u</sup>** ) + *<sup>μ</sup><sup>S</sup>* - <sup>∇</sup> **xu** + (<sup>∇</sup> **xu** )*T* and **<sup>u</sup>** *<sup>S</sup>*,*n*( **<sup>x</sup>**) = **<sup>U</sup>***S*,*n*(**X**). For small deformations, we have **<sup>F</sup>***<sup>n</sup>* <sup>≈</sup> **<sup>I</sup>**, *<sup>J</sup><sup>n</sup>* <sup>≈</sup> 1, then **<sup>Σ</sup>**(**V***S*,*n*+1) could be approached by **<sup>Σ</sup> <sup>x</sup>**(**<sup>v</sup>** *<sup>S</sup>*,*n*+1) and *<sup>σ</sup>S*,*<sup>n</sup>* by **<sup>Σ</sup> <sup>x</sup>**(**<sup>u</sup>** *<sup>S</sup>*,*n*).

Finally, we can approach the map **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> <sup>→</sup> **<sup>F</sup> <sup>Σ</sup>** , by the simplified linear application

$$
\hat{\mathbf{L}}\left(\hat{\mathbf{v}}^{S,n+1}\right) \quad = \quad \Sigma\_{\hat{\mathbf{x}}}(\hat{\mathbf{u}}^{S,n}) + (\boldsymbol{\Lambda}t)\Sigma\_{\hat{\mathbf{x}}}(\hat{\mathbf{v}}^{S,n+1}).\tag{23}
$$

The linearized updated Lagrangian weak formulation of the structure is: knowing **U***S*,*<sup>n</sup>* : Ω*<sup>S</sup>* <sup>0</sup> <sup>→</sup> <sup>R</sup>3, <sup>Ω</sup> *<sup>S</sup>* = <sup>Ω</sup>*<sup>S</sup> <sup>n</sup>* and **<sup>v</sup>***S*,*<sup>n</sup>* : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup>3, find **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup>3, **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* <sup>0</sup> such that

$$\begin{split} \int\_{\hat{\Omega}^{S}} \rho^{\hat{\mathbf{y}},\mu} \frac{\hat{\mathbf{v}}^{S,n+1} - \mathbf{v}^{S,n}}{\Delta t} \cdot \hat{\mathbf{w}}^{S} \, d\hat{\mathbf{x}} + \int\_{\hat{\Omega}^{S}} \hat{\mathbf{L}} \left( \hat{\mathbf{v}}^{S,n+1} \right) : \nabla\_{\hat{\mathbf{x}}} \hat{\mathbf{w}}^{S} \, d\hat{\mathbf{x}} \\ = \int\_{\hat{\Omega}^{S}} \rho^{\hat{\mathbf{y}},\mu} \mathbf{g} \cdot \hat{\mathbf{w}}^{S} \, d\hat{\mathbf{x}} + \int\_{\Gamma\_{0}} \mathbf{F}^{n+1} \boldsymbol{\Sigma}^{n+1} \mathbf{N}^{S} \cdot \mathbf{W}^{S} \, d\boldsymbol{S} \end{split} \tag{24}$$

for all **<sup>w</sup>** *<sup>S</sup>* : <sup>Ω</sup> *<sup>S</sup>* <sup>→</sup> <sup>R</sup>3, **<sup>w</sup>** *<sup>S</sup>* <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* 0 .

**Remark 1.** *We can use a non-linear model for the structure such as St. Venant Kirchhoff, neo-Hookean, etc. By linearization of* **F Σ** *(around the deformed state at the precedent time instant), we obtain in place of Equation (23)*

$$
\hat{\mathbf{L}}\left(\hat{\mathbf{v}}^{S,n+1}\right) = n\ell(\hat{\mathbf{u}}^{S,n}) + \left(\Delta t\right)\ell(\hat{\mathbf{v}}^{S,n+1})
$$

*where n is a non-linear operator and a linear one. Since n*-(**<sup>u</sup>** *<sup>S</sup>*,*n*) *is a known term, we can transfer it to the right-hand side, then the problem to solve is linear, similar to Equation (24).*

#### **4. Arbitrary Lagrangian Eulerian (ALE) Framework for Approximation of Fluid Equations**

The Arbitrary Eulerian Lagrangian (ALE) framework is a successful method to solve fluid equations in moving domain (see [19]). Let <sup>Ω</sup> *<sup>F</sup>* be a reference fluid domain and let <sup>A</sup>*t*, *<sup>t</sup>* <sup>∈</sup> [0, *<sup>T</sup>*] be a family of transformations such that: <sup>A</sup>*t*( **<sup>x</sup>**) = **<sup>x</sup>** for all **<sup>x</sup>** <sup>∈</sup> <sup>Σ</sup><sup>1</sup> <sup>∪</sup> <sup>Σ</sup><sup>3</sup> and <sup>A</sup>*t*(<sup>Ω</sup> *<sup>F</sup>*) = <sup>Ω</sup>*<sup>F</sup> t* , where **<sup>x</sup>** = (*<sup>x</sup>* 1, *<sup>x</sup>* 2, *<sup>x</sup>* <sup>3</sup>)*<sup>T</sup>* <sup>∈</sup> <sup>Ω</sup> *<sup>F</sup>* are the ALE coordinates and **<sup>x</sup>** = (*x*1, *<sup>x</sup>*2, *<sup>x</sup>*3)*<sup>T</sup>* <sup>=</sup> <sup>A</sup>*t*( **<sup>x</sup>**) the Eulerian coordinates.

Let **<sup>v</sup>***<sup>F</sup>* be the fluid velocity in the Eulerian coordinates. We denote by **<sup>v</sup>** *<sup>F</sup>* : <sup>Ω</sup> *<sup>F</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> the corresponding function in the ALE coordinates, which is defined by **<sup>v</sup>** *<sup>F</sup>*( **<sup>x</sup>**, *<sup>t</sup>*) = **<sup>v</sup>***F*(A*t*( **<sup>x</sup>**), *<sup>t</sup>*) = **<sup>v</sup>***F*(**x**, *<sup>t</sup>*).

We denote the mesh velocity by *ϑ*(**x**, *t*) = *<sup>∂</sup>*A*<sup>t</sup> <sup>∂</sup><sup>t</sup>* ( **<sup>x</sup>**) = *<sup>∂</sup>*A*<sup>t</sup> <sup>∂</sup><sup>t</sup>* (A−<sup>1</sup> *<sup>t</sup>* (**x**)) and the ALE time derivative of the fluid velocity by *<sup>∂</sup>***v***<sup>F</sup> ∂t* (**x**, *t*) = *<sup>∂</sup>***<sup>v</sup>** *<sup>F</sup> <sup>∂</sup><sup>t</sup>* ( **<sup>x</sup>**, *<sup>t</sup>*).

 **x** The Navier–Stokes equations in the ALE framework give:

$$\begin{aligned} \rho^{\boldsymbol{F}} \left( \frac{\partial \mathbf{v}^{\boldsymbol{F}}}{\partial t} \bigg|\_{\tilde{\mathbf{x}}} + \left( \left( \mathbf{v}^{\boldsymbol{F}} - \boldsymbol{\theta} \right) \cdot \nabla \right) \mathbf{v}^{\boldsymbol{F}} \right) - 2\mu^{\boldsymbol{F}} \nabla \cdot \boldsymbol{\epsilon} \left( \mathbf{v}^{\boldsymbol{F}} \right) + \nabla p^{\boldsymbol{F}} &= \quad \rho^{\boldsymbol{F}} \mathbf{g}, \text{ in } \Omega^{\mathrm{F}}\_{t} \times (0, T),\\ \nabla \cdot \mathbf{v}^{\boldsymbol{F}} &= \quad 0, \text{ in } \Omega^{\mathrm{F}}\_{t} \times (0, T). \end{aligned}$$

We denote by **<sup>v</sup>***F*,*n*, *<sup>p</sup>F*,*n*, and *<sup>ϑ</sup><sup>n</sup>* the approximations of **<sup>v</sup>***F*(·, *tn*), *<sup>p</sup>F*(·, *tn*), and *<sup>ϑ</sup>*(·, *tn*), respectively, all defined in Ω*<sup>F</sup> <sup>n</sup>*. Here, we set <sup>Ω</sup> *<sup>F</sup>* = <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>*. The time advancing scheme for fluid equations is: find **<sup>v</sup>** *<sup>F</sup>*,*n*+<sup>1</sup> : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> and *<sup>p</sup> <sup>F</sup>*,*n*+<sup>1</sup> : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup> such that

$$\begin{split} \rho^{\mathbb{F}} \left( \frac{\widehat{\mathbf{v}}^{\mathbb{F},n+1} - \mathbf{v}^{\mathbb{F},n}}{\Delta t} + \left( \left( \mathbf{v}^{\mathbb{F},n} - \boldsymbol{\theta}^{n} \right) \cdot \nabla\_{\widehat{\mathbf{x}}} \right) \widehat{\mathbf{v}}^{\mathbb{F},n+1} \right) \\ -2\mu^{\mathbb{F}} \nabla\_{\widehat{\mathbf{x}}} \cdot \boldsymbol{\epsilon} \left( \widehat{\mathbf{v}}^{\mathbb{F},n+1} \right) + \nabla\_{\widehat{\mathbf{x}}} \widehat{\boldsymbol{p}}^{\mathbb{F},n+1} \end{split} \tag{25}$$

$$\begin{array}{rcl} \nabla\_{\hat{\mathbf{x}}} \cdot \hat{\mathbf{v}}^{F, \mathbf{n}+1} &=& 0, \text{ in } \Omega^{F}\_{\mathbf{n}} \\ \cdot \quad \cdot \quad \cdot \quad \cdot \end{array} \tag{26}$$

$$\sigma^{F}(\hat{\mathbf{v}}^{F,n+1}, \hat{p}^{F,n+1})\mathbf{n}^{F} = \mathbf{h}^{n+1}\_{\text{in}} \text{ on } \Sigma\_1 \tag{27}$$

$$\sigma^F(\hat{\mathbf{v}}^{F,n+1}, \hat{p}^{F,n+1}) \mathbf{n}^F = \mathbf{h}\_{\text{out}}^{n+1}, \text{ on } \Sigma\_3 \tag{28}$$

The above time discretization scheme is based on the backward Euler scheme and a linearization of the convective term.

We multiply Equation (25) by a test function **<sup>w</sup>** *<sup>F</sup>* : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> and Equation (26) by a test function *<sup>q</sup>* : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup>. After integrating them over the domain <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* and using the Green's formula and the corresponding boundary conditions, we get the below discrete weak form. Find **<sup>v</sup>** *<sup>F</sup>*,*n*+<sup>1</sup> : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> and *<sup>p</sup> <sup>F</sup>*,*n*+<sup>1</sup> : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup> such that:

$$\begin{split} & \int\_{\Omega\_{n}^{F}} \rho^{F} \frac{\widehat{\mathbf{v}}^{F,n+1}}{\Delta t} \cdot \widehat{\mathbf{w}}^{F} d\widehat{\mathbf{x}} + \int\_{\Omega\_{n}^{F}} \rho^{F} \left( \left( \left( \mathbf{v}^{F,n} - \boldsymbol{\mathfrak{o}}^{\boldsymbol{n}} \right) \cdot \nabla\_{\widehat{\mathbf{x}}} \right) \widehat{\mathbf{v}}^{F,n+1} \right) \cdot \widehat{\mathbf{w}}^{F} d\widehat{\mathbf{x}} \\ & - \int\_{\Omega\_{n}^{F}} \left( \nabla\_{\widehat{\mathbf{x}}} \cdot \widehat{\mathbf{w}}^{F} \right) \widehat{p}^{F,n+1} d\widehat{\mathbf{x}} + \int\_{\Omega\_{n}^{F}} 2\mu^{F} \boldsymbol{\varepsilon} \left( \widehat{\mathbf{v}}^{F,n+1} \right) : \boldsymbol{\varepsilon} \left( \widehat{\mathbf{w}}^{F} \right) d\widehat{\mathbf{x}} \\ & = \mathcal{L}\_{F} (\widehat{\mathbf{w}}^{F}) + \int\_{\Gamma\_{n}} \left( \boldsymbol{\sigma}^{F} \left( \widehat{\mathbf{v}}^{F,n+1}, \widehat{p}^{F,n+1} \right) \mathbf{n}^{F} \right) \cdot \widehat{\mathbf{w}}^{F} d\mathbf{s}, \\ & \int\_{\Gamma\_{n}} (\nabla\_{\widehat{\mathbf{x}}} \cdot \widehat{\mathbf{v}}^{F,n+1}) \widehat{\mathbf{x}} d\widehat{\mathbf{x}} = 0. \end{split} \tag{29}$$

$$\int\_{\Omega\_n^{\mathbb{F}}} (\nabla\_{\hat{\mathbf{x}}} \cdot \hat{\mathbf{v}}^{\mathbb{F}, n+1}) \hat{q} \, d\hat{\mathbf{x}} = 0,\tag{30}$$

for all **<sup>w</sup>** *<sup>F</sup>* : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> and for all *<sup>q</sup>* : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup>, where

$$\mathcal{L}\_F(\widehat{\mathbf{w}}^F) = \int\_{\Omega\_n^F} \boldsymbol{\rho}^F \frac{\widehat{\mathbf{v}}^{F,n}}{\Delta t} \cdot \widehat{\mathbf{w}}^F d\widehat{\mathbf{x}} + \int\_{\Omega\_n^F} \boldsymbol{\rho}^F \mathbf{g} \cdot \widehat{\mathbf{w}}^F + \int\_{\Sigma\_1} \mathbf{h}\_{in}^{n+1} \cdot \widehat{\mathbf{w}}^F + \int\_{\Sigma\_3} \mathbf{h}\_{out}^{n+1} \cdot \widehat{\mathbf{w}}^F.$$

We have assumed, for instance, that the forces *<sup>σ</sup>F*(**<sup>v</sup>** *<sup>F</sup>*,*n*<sup>+</sup>1, *<sup>p</sup> <sup>F</sup>*,*n*+1)**n***<sup>F</sup>* on the interface <sup>Γ</sup>*<sup>n</sup>* are known. The mesh velocity *ϑ <sup>n</sup>*+<sup>1</sup> : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> can be computed from

$$\begin{cases} \begin{array}{rcl} \Delta\_{\hat{\mathbf{x}}} \hat{\boldsymbol{\theta}}^{n+1} &=& 0, & \text{in } \Omega\_{n}^{F} \\ \hat{\boldsymbol{\theta}}^{n+1} &=& 0, & \text{on } \Sigma\_{1} \cup \Sigma\_{3} \\ \hat{\boldsymbol{\theta}}^{n+1} &=& \hat{\mathbf{v}}^{F,n+1}, & \text{on } \Gamma\_{n}. \end{array} \end{cases}$$

For all *<sup>n</sup>* <sup>=</sup> 0, ··· , *<sup>N</sup>* <sup>−</sup> 1, we denote by <sup>A</sup>*tn*+<sup>1</sup> the map from <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* to <sup>R</sup><sup>3</sup> defined by <sup>A</sup>*tn*+<sup>1</sup> (*<sup>x</sup>* 1, *<sup>x</sup>* 2, *<sup>x</sup>* <sup>3</sup>) = (*<sup>x</sup>* <sup>1</sup> <sup>+</sup> <sup>Δ</sup>*tϑn*+<sup>1</sup> <sup>1</sup> , *<sup>x</sup>* <sup>2</sup> <sup>+</sup> <sup>Δ</sup>*tϑn*+<sup>1</sup> <sup>2</sup> , *<sup>x</sup>* <sup>3</sup> <sup>+</sup> <sup>Δ</sup>*tϑn*+<sup>1</sup> <sup>3</sup> ). We set <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>*+<sup>1</sup> <sup>=</sup> <sup>A</sup>*tn*+<sup>1</sup> (Ω*<sup>F</sup> <sup>n</sup>*), Γ*n*+<sup>1</sup> = A*tn*+<sup>1</sup> (Γ*n*) and we remark that **<sup>x</sup>** <sup>=</sup> <sup>A</sup>*tn*+<sup>1</sup> ( **<sup>x</sup>**), for all **<sup>x</sup>** <sup>∈</sup> <sup>Σ</sup><sup>1</sup> <sup>∪</sup> <sup>Σ</sup>3.

We define the fluid velocity **v***F*,*n*+<sup>1</sup> : Ω*<sup>F</sup> <sup>n</sup>*+<sup>1</sup> <sup>→</sup> <sup>R</sup>3, the fluid pressure *<sup>p</sup>F*,*n*+<sup>1</sup> : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>*+<sup>1</sup> <sup>→</sup> <sup>R</sup> and the mesh velocity *ϑn*+<sup>1</sup> : Ω*<sup>F</sup> <sup>n</sup>*+<sup>1</sup> <sup>→</sup> <sup>R</sup><sup>3</sup> by:

$$\mathbf{v}^{\mathsf{F},\mathsf{u}+1}(\mathbf{x}) = \widehat{\mathbf{v}}^{\mathsf{F},\mathsf{u}+1}(\widehat{\mathbf{x}}), \quad \mathbf{p}^{\mathsf{F},\mathsf{u}+1}(\mathbf{x}) = \widehat{\mathbf{p}}^{\mathsf{F},\mathsf{u}+1}(\widehat{\mathbf{x}}), \quad \mathbf{\theta}^{\mathsf{u}+1}(\mathbf{x}) = \widehat{\mathbf{\theta}}^{\mathsf{u}+1}(\widehat{\mathbf{x}}),$$

for all **<sup>x</sup>** <sup>∈</sup> <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>*, **<sup>x</sup>** <sup>=</sup> <sup>A</sup>*tn*+<sup>1</sup> ( **<sup>x</sup>**) <sup>∈</sup> <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>*+1.

#### **5. Monolithic Formulation for the Fluid–Structure Equations**

Let Ω*<sup>n</sup>* = Ω*<sup>F</sup> <sup>n</sup>* <sup>∪</sup> <sup>Γ</sup>*<sup>n</sup>* <sup>∪</sup> <sup>Ω</sup>*<sup>S</sup> <sup>n</sup>* be the global fluid–structure domain at time instant *n* and let us introduce the global velocity and test function

$$
\hat{\mathbf{v}}^{n+1} : \Omega\_{\boldsymbol{n}} \to \mathbb{R}^{\hat{3}}, \quad \hat{\mathbf{w}} : \Omega\_{\boldsymbol{n}} \to \mathbb{R}^{\hat{3}}
$$

$$
\hat{\mathbf{v}}^{n+1} = \begin{cases}
\hat{\mathbf{v}}^{F,n+1} \text{ in } \Omega\_{\boldsymbol{n}}^{F} \\
\hat{\mathbf{v}}^{S,n+1} \text{ in } \Omega\_{\boldsymbol{n}}^{S}
\end{cases} \qquad \hat{\mathbf{w}} = \begin{cases}
\quad \hat{\mathbf{w}}^{F} \text{ in } \Omega\_{\boldsymbol{n}}^{F} \\
\quad \hat{\mathbf{w}}^{S} \text{ in } \Omega\_{\boldsymbol{n}}^{S}.
\end{cases}
$$

At each time step, we solve the linear coupled problem: find **<sup>v</sup>** *<sup>n</sup>*+<sup>1</sup> <sup>∈</sup> - *H*<sup>1</sup> (Ω*n*) 2 , **<sup>v</sup>** *<sup>n</sup>*+<sup>1</sup> <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* 0 and *<sup>p</sup> <sup>F</sup>*,*n*+<sup>1</sup> <sup>∈</sup> *<sup>L</sup>*<sup>2</sup> - Ω*<sup>F</sup> n* , such that:

$$\begin{split} &\int\_{\Omega\_{n}^{F}}\rho^{F}\frac{\widehat{\mathbf{v}}^{n+1}}{\Delta t}\cdot\widehat{\mathbf{w}}d\widehat{\mathbf{x}}+\int\_{\Omega\_{n}^{F}}\rho^{F}\left(\left(\left(\mathbf{v}^{n}-\boldsymbol{\uptheta}^{n}\right)\cdot\nabla\_{\widehat{\mathbf{x}}}\right)\widehat{\mathbf{v}}^{n+1}\right)\cdot\widehat{\mathbf{w}}d\widehat{\mathbf{x}} \\ &= -\int\_{\Omega\_{n}^{F}}\left(\nabla\_{\widehat{\mathbf{x}}}\cdot\widehat{\mathbf{w}}\right)\widehat{\mathbf{y}}^{F,n+1}d\widehat{\mathbf{x}}+\int\_{\Omega\_{n}^{F}}2\mu^{F}\boldsymbol{\varepsilon}\left(\widehat{\mathbf{v}}^{n+1}\right):\boldsymbol{\upepsilon}\left(\widehat{\mathbf{w}}\right)d\widehat{\mathbf{x}} \\ &+\int\_{\Omega\_{n}^{S}}\rho^{S,n}\frac{\widehat{\mathbf{v}}^{n+1}}{\Delta t}\cdot\widehat{\mathbf{w}}d\widehat{\mathbf{x}}+\int\_{\Omega\_{n}^{S}}\widehat{\mathbf{L}}\cdot\left(\widehat{\mathbf{v}}^{n+1}\right):\nabla\_{\widehat{\mathbf{x}}}\widehat{\mathbf{w}}d\widehat{\mathbf{x}} \\ &= \mathcal{L}\_{F}(\widehat{\mathbf{w}})+\int\_{\Omega\_{n}^{S}}\rho^{S,n}\frac{\mathbf{v}^{n}}{\Delta t}\cdot\widehat{\mathbf{w}}d\widehat{\mathbf{x}}+\int\_{\Omega\_{n}^{S}}\rho^{S,n}\mathbf{g}\cdot\widehat{\mathbf{w}}d\widehat{\mathbf{x}}, \\ &\int\_{\Omega\_{n}^{F}}\widehat{\mathbf{w}}\cdot\widehat{\mathbf{w}}d\widehat{\mathbf{x}}d\widehat{\mathbf{x}}. \end{split} \tag{31}$$

$$\int\_{\Omega\_{\rm u}^{\mathbb{F}}} (\nabla\_{\widetilde{\mathbf{x}}} \cdot \widehat{\mathbf{v}}^{n+1}) \widehat{q} \, d\widetilde{\mathbf{x}} = 0,\tag{32}$$

for all **<sup>w</sup>** <sup>∈</sup> - *H*<sup>1</sup> (Ω*n*) 3 , **<sup>w</sup>** <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* <sup>0</sup> and for all *<sup>q</sup>* <sup>∈</sup> *<sup>L</sup>*<sup>2</sup> - Ω*<sup>F</sup> n* .

From the regularity **<sup>v</sup>** *<sup>n</sup>*+<sup>1</sup> <sup>∈</sup> - *H*<sup>1</sup> (Ω*n*) 2 , the traces of **<sup>v</sup>** *<sup>F</sup>*,*n*+<sup>1</sup> and **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> on <sup>Γ</sup>*<sup>n</sup>* are well defined and **<sup>v</sup>** *<sup>F</sup>*,*n*+<sup>1</sup> <sup>|</sup>Γ*<sup>n</sup>* <sup>=</sup> **<sup>v</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> <sup>|</sup>Γ*<sup>n</sup>* which is a discrete form of the continuity of the velocity at the interface (8). Similarly, we get **<sup>w</sup>** *<sup>F</sup>*,*n*+<sup>1</sup> <sup>|</sup>Γ*<sup>n</sup>* <sup>=</sup> **<sup>w</sup>** *<sup>S</sup>*,*n*+<sup>1</sup> <sup>|</sup>Γ*<sup>n</sup>* .

Equation (31) is obtained by adding Equations (24) and (29). By enforcement of the hypothesis of continuity of stress (Equation (9)) at the discrete level, the expression

$$\int\_{\Gamma\_0} \mathbf{F}^{n+1} \boldsymbol{\Sigma}^{n+1} \mathbf{N}^S \cdot \mathbf{W}^S \, dS + \int\_{\Gamma\_n} \left( \sigma^F(\hat{\mathbf{v}}^{F, n+1}, \hat{p}^{F, n+1}) \mathbf{n}^F \right) \cdot \hat{\mathbf{w}}^F ds = 0$$

does not appear anymore in Equation (31).

#### **Algorithm for fluid–structure interaction Time advancing scheme from** *n* **to** *n* + 1

We assume that we know Ω*<sup>n</sup>* = Ω*<sup>F</sup> <sup>n</sup>* <sup>∪</sup> <sup>Γ</sup>*<sup>n</sup>* <sup>∪</sup> <sup>Ω</sup>*<sup>S</sup> <sup>n</sup>*, **<sup>v</sup>***<sup>n</sup>* : <sup>Ω</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup>3, *<sup>ϑ</sup><sup>n</sup>* : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup>3.

**Step 1**: Solve the **linear** system in Equations (31) and (32) and get the velocity **<sup>v</sup>** *<sup>n</sup>*+<sup>1</sup> and the pressure *<sup>p</sup> <sup>F</sup>*,*n*<sup>+</sup>1.

**Step 2**: Compute the fluid mesh velocity *ϑ <sup>n</sup>*+<sup>1</sup> : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup>

$$\begin{cases} \begin{array}{rcl} \Delta\_{\overline{\mathsf{X}}} \widehat{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathcal{\mathsf{\mathsf{\mathsf{\mathcal{\mathsf{\mathcal{\mathsf{\mathcal{\mathsf{\mathcal{\mathsf{\mathcal{\mathsf{\mathcal{\mathcal{\mathsf{\mathcal{\mathsf{\mathcal{\mathcal{\mathcal{\theta}}}}}}}}}}}}}}}}}}}}}} \mathbf}}} \mathbf}} \mathbf}} \mathbf} \mathbf}} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf} \mathbf$$

To improve the quality of the mesh, we can replace in Equation (33) the Laplacian by a linear elasticity operator.

**Step 3**: Define the map <sup>T</sup>*<sup>n</sup>* : <sup>Ω</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup><sup>3</sup> by:

$$\mathbb{T}\_n(\widehat{\mathbf{x}}) = \widehat{\mathbf{x}} + (\Delta t) \widehat{\mathfrak{b}}^{n+1}(\widehat{\mathbf{x}}) \chi\_{\Omega\_n^F}(\widehat{\mathbf{x}}) + (\Delta t) \widehat{\mathbf{v}}^{n+1}(\widehat{\mathbf{x}}) \chi\_{\Omega\_n^S}(\widehat{\mathbf{x}}) $$

where *χ*Ω*<sup>F</sup> <sup>n</sup>* and *χ*Ω*<sup>S</sup> <sup>n</sup>* are the characteristic functions of fluid and structure domains.

**Step 4**: Set Ω*<sup>F</sup> <sup>n</sup>*+<sup>1</sup> <sup>=</sup> <sup>T</sup>*n*(Ω*<sup>F</sup> <sup>n</sup>*), Ω*<sup>S</sup> <sup>n</sup>*+<sup>1</sup> <sup>=</sup> <sup>T</sup>*n*(Ω*<sup>S</sup> <sup>n</sup>*), and <sup>Γ</sup>*n*+<sup>1</sup> = T*n*(Γ*n*); consequently, <sup>Ω</sup>*n*+<sup>1</sup> = T*n*(Ω*n*). Define **<sup>v</sup>***n*+<sup>1</sup> : <sup>Ω</sup>*n*+<sup>1</sup> <sup>→</sup> <sup>R</sup><sup>3</sup> by

$$\mathbf{v}^{n+1}(\mathbf{x}) = \widehat{\mathbf{v}}^{n+1}(\widehat{\mathbf{x}}), \forall \widehat{\mathbf{x}} \in \Omega\_{n\prime} \ \mathbf{x} = \mathbb{T}\_n(\widehat{\mathbf{x}})$$

and *pF*,*n*+<sup>1</sup> : Ω*<sup>F</sup> <sup>n</sup>*+<sup>1</sup> <sup>→</sup> <sup>R</sup>, *<sup>ϑ</sup>n*+<sup>1</sup> : <sup>Ω</sup>*<sup>F</sup> <sup>n</sup>*+<sup>1</sup> <sup>→</sup> <sup>R</sup><sup>3</sup> by:

$$p^{\mathsf{F}, \mathsf{u} + 1}(\mathbf{x}) = \widehat{p}^{\mathsf{F}, \mathsf{u} + 1}(\widehat{\mathbf{x}}), \ \mathfrak{P}^{\mathsf{u} + 1}(\mathbf{x}) = \widehat{\mathfrak{P}}^{\mathsf{u} + 1}(\widehat{\mathbf{x}}), \ \forall \widehat{\mathbf{x}} \in \Omega^{\mathrm{F}}\_{\mathsf{u} \mathsf{v}}.\ \mathbf{x} = \mathbb{T}\_{\mathsf{u}}(\widehat{\mathbf{x}}).\ \mathsf{u}$$

**Remark 2.** *(i) The domain is computed explicitly while the velocity and the pressure are computed implicitly. This kind of schema is also called semi-implicit. The monolithic system in Equations (31) and (32) is linear.*

*(ii) Using globally continuous finite element for the velocity* **<sup>v</sup>** *<sup>n</sup>*+<sup>1</sup> <sup>∈</sup> - *H*<sup>1</sup> (Ω*n*) <sup>2</sup> *defined all over the fluid–structure global mesh, the continuity of the velocity at the interface holds, automatically.*

*(iii) The vertices in the structure mesh are moved using the structure velocity, thus the structure mesh is of updated Lagrangian type.*

**Remark 3.** *In [35], for linear elastic model for the structure and Navier–Stokes equations for the fluid, a semi-implicit monolithic algorithm is introduced. The unconditional stability in time is established. In a forthcoming paper, a proof of the unconditional stability of an algorithm will be presented for a non-linear model of the structure. This stable algorithm is similar to the one presented in this paper, only a stabilization term has been added, as in [35].*

**Remark 4.** *In this paper, the derivative of fluid velocity as well as the derivative of structure velocity are approached by the implicit Euler scheme. It is possible to use different time discretization schemes, for example Newmark for the structure and implicit Euler for the fluid. We have to pay attention to the time advancing algorithm for the interface. One solution is to solve the fluid–structure coupled equations written in the domain at the precedent time instant to find the fluid–structure velocity and the fluid pressure. Then, the structure including the interface is advanced by the Newmark scheme, and finally the fluid mesh velocity is computed using the new position and velocity of the interface.*

#### **6. Numerical Experiments**

The numerical tests were produced using *FreeFem++* (see [40]).

#### *6.1. Straight Cylinder*

We tested the benchmark studied in [3,4] concerning blood flow in artery. The geometrical configuration is represented in Figure 1. The fluid occupies initially the straight cylinder of length

*L* = 5 cm and radius *R* = 0.5 cm. The disk Σ<sup>1</sup> is in the plane *x*1*Ox*<sup>2</sup> and the axis of the cylinder is *Ox*3. The viscosity of the fluid is *μ<sup>F</sup>* = 0.03 <sup>g</sup> cm·<sup>s</sup> and its density is *<sup>ρ</sup><sup>F</sup>* <sup>=</sup> <sup>1</sup> <sup>g</sup> cm<sup>3</sup> .

The fluid is surrounded by a structure of constant thickness *h<sup>S</sup>* = 0.1 cm. The others physical parameters of the structure are: the Young's modulus is *<sup>E</sup>* <sup>=</sup> <sup>3</sup> <sup>×</sup> 106 <sup>g</sup> cm·s<sup>2</sup> , the Poisson ratio is *<sup>ν</sup><sup>S</sup>* <sup>=</sup> 0.3, and the density is *ρ<sup>S</sup>* <sup>0</sup> <sup>=</sup> 1.2 <sup>g</sup> cm<sup>3</sup> . The Lamé parameters are computed by the formulas *<sup>λ</sup><sup>S</sup>* <sup>=</sup> *<sup>ν</sup>SE* (1−2*νS*)(1+*νS*) and *μ<sup>S</sup>* = *<sup>E</sup>* 2(1+*νS*) .

For the volume force in fluid and structure, we put **g** = (0, 0, 0)*T*. The prescribed boundary stress at the inlet <sup>Σ</sup><sup>1</sup> is **<sup>h</sup>***in* = - 0, 0, 1.3332 <sup>×</sup> <sup>10</sup><sup>4</sup> g cm·s2 for *<sup>t</sup>* <sup>≤</sup> 0.005 <sup>s</sup> and **<sup>h</sup>***out* <sup>=</sup> (0, 0, 0) <sup>g</sup> cm·s2 at the outlet Σ3. The structure is clamped at both ends, Σ<sup>5</sup> and Σ7. The remaining boundary conditions are Equations (3), (8) and (9). Initially, the fluid and the structure are at rest.

Using *FreeFem++*, it is possible to construct a global fluid–structure mesh with an "interior boundary" that is the fluid–structure interface. For the finite element approximation of the fluid–structure velocity, we used the finite element P<sup>1</sup> + *bubble* and we employed for the pressure the finite element P1. The parameters of meshes used for the numerical tests are presented in Table 1.

**Table 1.** The number of vertices, tetrahedra and degrees of freedom (DOF) of fluid–structure linear system for each mesh.


The time step was set to Δ*t* = 0.0005 s and the number of time steps to *N* = 40. The radial displacement of the interface was measured at three points *A*(*R*, 0, *L*/4), *B*(*R*, 0, *L*/2), and *C*(*R*, 0, 3*L*/4). We observed that the displacements were small, less than 0.012 cm (see Figure 2, left). Similar results are observed in [41] using non-conforming meshes. In addition, we measured the radial displacement along the line (*R*, 0, *x*3), *x*<sup>3</sup> ∈ [0, *L*] on the interface (see Figure 2, right.)

**Figure 2.** Constant stress at the inlet: the time history of the displacement of three points on the interface (**left**); and the radial displacement at time instants *t* = 0.0060, *t* = 0.0080, *t* = 0.0100 (**right**).

To obtain visible displacement, we also tested with the sin-like stress at the inlet:

$$\mathbf{h}\_{in}(\mathbf{x},t) = \begin{cases} \begin{pmatrix} 0, 0, 5 \times 1.3332 \times 10^4 \times \frac{(1 - \cos(2\pi t/0.001))}{2} \end{pmatrix}, & \mathbf{x} \in \Sigma\_1, 0 \le t \le 0.001\\ \begin{pmatrix} 0, 0, 0 \end{pmatrix}, & \mathbf{x} \in \Sigma\_1, 0.001 \le t \le T \end{cases}$$

where *T* = 0.02 s. The other parameters were the same. The displacements at three points are plotted in Figure 3 and the pressure at three time instants is plotted in Figure 4.

**Figure 3.** The time history of the displacement of three points on the interface using sin-like stress at the inlet.

At each time step, we had to solve a sparse non-symmetric linear system for the fluid–structure velocity and pressure and a sparse symmetric positive definite linear system for the mesh displacement. The linear systems were solved using MUMPS (MUltifrontal Massively Parallel sparse direct Solver), implemented in *FreeFem++*. This method was very efficient; the total CPU time for the first three meshes were: 44 s (1.1 s/iteration) for Mesh 1, 173 s (4.3 s/iteration) for Mesh 2 and 1083 s (27 s/iteration) for Mesh 3, on a computer with a processor of 4 × 3.30 GHz frequency and 16 Go RAM. For Mesh 4, the CPU time was about 10 min/iteration on a noed Intel Sandy-Bridge 16 × 3.30 GHz and 64 Go RAM.

**Figure 4.** The pressure obtained with Mesh 4 using sin-like stress at the inlet.

In [4], the mean number of iterations per time step is: 33.9 for the fixed-point with Aitken acceleration and 8.9 iterations for a quasi-Newton algorithm. At each iteration of quasi-Newton algorithm, a linear system is solved. In [41,42], the average number of Newton iterations per time step is 3.

#### *6.2. Cane Cylinder*

We considered the fluid–structure interaction in a cane-like geometry inspired from [4,43]. The parameters of the fluid domain were: *R* = 0.5 cm, *R*<sup>1</sup> = 1 cm, *L*<sup>1</sup> = 1 cm, *L*<sup>2</sup> = 5 cm (see Figure 5, left). The fluid was surrounded by a structure of constant thickness *h<sup>S</sup>* = 0.1 cm.

**Figure 5.** The parameters of the fluid domain (**left**); and a global mesh for fluid–structure domain (**right**).

The time step was Δ*t* = 0.0005 s, but the final time was *T* = 0.04 s. The other parameters were the same as in the case of the straight cylinder. The details of meshes used for the numerical tests are presented in Table 2.

**Table 2.** The number of vertices, tetrahedra and degrees of freedom (DOF) of fluid–structure linear system for each mesh.


We measured the displacement of the interface ar three points: *U<sup>S</sup>* <sup>1</sup> at *<sup>A</sup>*(*R*<sup>1</sup> + <sup>2</sup>*R*, 0, 0), *<sup>U</sup><sup>S</sup>* <sup>3</sup> at *B*(0, 0, *R*<sup>1</sup> + 2*R*) and *U<sup>S</sup>* <sup>1</sup> at *C*(−*R*<sup>1</sup> − 2*R*, 0, 0) (see Figure 6). The pressure at three time instants is plotted in Figure 7. We observed that, for the sin-like stress at the inlet, the structure displacements were greater than 0.23 (cm) and a non-linear model for the structure could be more appropriate. We chose a sin-like stress at the inlet with maximal value five times greater than the constant case to obtain visible deformations. Even though we used the same sin-like stress at the inlet, the structure displacements were larger than in the straight cylinder case because of the shape as well as because the cane was longer. Recall that the structure was fixed at both ends.

The total CPU time for the three meshes were: 255 s (3.18 s/iteration) for Mesh 1, 502 s (6.2 s/iteration) for Mesh 2 and 915 s (11.2 s/iteration) for Mesh 3, on a computer with a processor of 4 × 3.30 GHz frequency and 16 Go RAM.

In [4], the mean number of iterations per time step is 8.9 for quasi-Newton algorithm. The fixed-point algorithm with Aitken acceleration failed.

**Figure 6.** The time history of the displacement of three points on the interface using constant (**left**) and sin-like (**right**) stress at the inlet.

**Figure 7.** The computed pressure obtained with sin-like stress at the inlet (Mesh 2).

#### *6.3. Artery Stenosis*

Now, we considered a fluid–structure interaction in artery stenosis inspired from the paper [42]. The inlet and outlet surfaces Σ1, Σ<sup>3</sup> were disks of radius *R*, normal to the axis *Ox*<sup>1</sup> of centers (*xmin*, 0, 0) and (*xmax*, 0, 0), respectively. The lateral surface Γ<sup>0</sup> of the initial fluid domain was composed by the bottom half straight cylinder surface

$$\left\{ (\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) \in \mathbb{R}^3 \colon \mathbf{x}\_{\text{min}} < \mathbf{x}\_1 < \mathbf{x}\_{\text{max}}, \sqrt{\mathbf{x}\_2^2 + \mathbf{x}\_3^2} = R, \ x\_3 < 0 \right\}$$

and the stenosis surface (see Figure 8), obtained from the top half straight cylinder surface

$$\left\{ (\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) \in \mathbb{R}^3 ; \ x\_{\text{min}} < x\_1 < x\_{\text{max}}, \sqrt{\mathbf{x}\_2^2 + \mathbf{x}\_3^2} = R, \ x\_3 \ge 0 \right\}$$

via the map

$$g(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) \to \left(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3 - \left(h\_{\max}^S - h\_{\min}^S\right) \frac{\left(\ell^2 - \mathbf{x}\_1^2\right)}{\ell^2} \frac{\left(R^2 - \mathbf{x}\_2^2\right)}{R^2}\right), \ |\mathbf{x}\_1| \le \ell, |\mathbf{x}\_2| \le R.$$

The initial boundary of the structure domain was composed by: the interior lateral surface Γ0, which is the fluid–structure interface; the exterior lateral surface

$$\Gamma\_0^N = \left\{ (\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) \in \mathbb{R}^3 \colon \mathbf{x}\_{\min} < \mathbf{x}\_1 < \mathbf{x}\_{\max}, \sqrt{\mathbf{x}\_2^2 + \mathbf{x}\_3^2} = R + h\_{\min}^S \right\};$$

and two annular surfaces Σ5, Σ7, of radii *R* and *R* + *h<sup>S</sup> min*, normal to the axis *Ox*1, of centers (*xmin*, 0, 0) and (*xmax*, 0, 0), respectively.

**Figure 8.** The stenosis interface.

The numerical values were: *R* = 0.5 cm, *xmin* = −3, *xmax* = 3, - = 1.5 cm, *h<sup>S</sup> min* = 0.1 cm, and *h<sup>S</sup> max* = 0.5 cm. The details of meshes used for the numerical tests are presented in Table 3.

**Table 3.** The number of vertices, tetrahedra and degrees of freedom (DOF) of fluid–structure linear system for each mesh.


We used the same sin-like stress at the inlet similar to the precedent experiments

$$\mathbf{h}\_{in}(\mathbf{x},t) = \begin{cases} \begin{pmatrix} 0, \ 0, \ m \text{ag} \times 1.3332 \times 10^4 \times \frac{(1 - \cos(2\pi t/0.001))}{2} \end{pmatrix}, & \mathbf{x} \in \Sigma\_1, 0 \le t \le 0.001\\ \begin{pmatrix} 0, \ 0, \ 0 \end{pmatrix}, & \mathbf{x} \in \Sigma\_1, 0.001 \le t \le T \end{cases}$$

where *T* = 0.02 s and *mag* > 0 is a parameter.

For *Mesh 1*, the time step was set to <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> s, the number of time steps to *<sup>N</sup>* <sup>=</sup> <sup>40</sup> and *mag* = 5. The total CPU time was 323 s (8.07 s/iteration), on a computer with a processor of 4 × 3.30 GHz frequency and 16 Go RAM. The radial displacement of the interface was measured at three points: *A*(−-, 0, *<sup>R</sup>*), *<sup>B</sup>*(0, 0, *<sup>R</sup>* <sup>−</sup> (*h<sup>S</sup> max* <sup>−</sup> *<sup>h</sup><sup>S</sup> min*)), and *C*(-, 0, *R*) (see Figure 9). We observed that the maximal displacement of point A was more than 0.08 cm, which was more important than in the case of straight cylinder (see Figure 3). The artery deformation was more important in the uphill zone of the stenosis than in the case of healthy artery.

**Figure 9.** The time history of the displacement of three points on the interface using sin-like stress at the inlet (Mesh 1).

For *Mesh 2*, the time step was set to Δ*t* = 10−<sup>4</sup> s, the number of time steps to *N* = 200 and *mag* = 2. The total CPU time was 15,123 s (75.6 s/iteration), on a noed Intel Sandy-Bridge 16 × 3.30 GHz and 64 Go RAM.

For *Mesh 3*, the time step was set to Δ*t* = 10−<sup>4</sup> s, the number of time steps to *N* = 200 and *mag* = 1. The average CPU time was 4 min 55 s by iteration, on a noed Intel Sandy-Bridge 16 × 3.30 GHz and 64 Go RAM. The pressure at three time instants is plotted in Figure 10.

**Figure 10.** The computed pressure obtained with *Mesh 3*.

For *Mesh 4*, the time step was set to <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> s, the number of time steps to *<sup>N</sup>* <sup>=</sup> 400 and *mag* = 1. The average CPU time was 18 min 8 s by iteration, on a noed Intel Sandy-Bridge 16 × 3.30 GHz and 64 Go RAM.

For finer mesh, we were forced to use smaller time steps. This phenomenon was observed for neither the three-dimensional tests presented above nor for the two-dimensional benchmark flow around a flexible thin structure attached to a fixed cylinder (see [37]). The semi-implicit algorithms had

good stability properties (see Remark 3). We suppose that the source of the problem is this particular surface mesh of the stenosis zone obtained from the mesh of the top half straight cylinder surface by vertical projection on the stenosis surface.

#### **7. Conclusions**

We have presented a monolithic semi-implicit method for three-dimensional fluid–structure interaction problems. At each time step, we solve only a linear system to find the fluid–structure velocity and the fluid pressure, thus the method is fast. We use a global mesh for the fluid–structure domain where the fluid–structure interface is an interior boundary. Using globally continuous finite element for the velocity in the fluid–structure mesh, the continuity of velocity at the interface is automatically satisfied.

**Funding:** This research received no external funding.

**Acknowledgments:** I gratefully thank Computing Center, University of Strasbourg for giving me access to their computing resources.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **On the Kutta Condition in Compressible Flow over Isolated Airfoils**

#### **Farzad Mohebbi 1,\*, Ben Evans <sup>1</sup> and Mathieu Sellier <sup>2</sup>**


Received: 10 April 2019; Accepted: 27 May 2019; Published: 1 June 2019

**Abstract:** This paper presents a novel and accurate method to implement the Kutta condition in solving subsonic (subcritical) inviscid isentropic compressible flow over isolated airfoils using the stream function equation. The proposed method relies on body-fitted grid generation and solving the stream function equation for compressible flows in computational domain using finite-difference method. An expression is derived for implementing the Kutta condition for the airfoils with both finite angles and cusped trailing edges. A comparison of the results obtained from the proposed numerical method and the results from experimental and other numerical methods reveals that they are in excellent agreement, which confirms the accuracy and correctness of the proposed method.

**Keywords:** computational aerodynamics; Kutta condition; compressible flow; stream function

#### **1. Introduction**

Nowadays, computational fluid dynamics complements its experimental and theoretical counterpart. Benefiting from high-speed digital computers, the use of sophisticated numerical methods has made possible the numerical solution of fluid flow problems which were heretofore intractable. Compressible flows over airfoils and wings play a vital role in computational aerodynamics, and demand advanced computational techniques. Since introducing the source and vortex panel methods in the late 1960s [1], they have become the standard tools to numerically solve low-speed flows over bodies of arbitrary shape, and have had extensive applications in flow modeling [2–6]. The panel methods used for the simulation of flows over an airfoil are concerned with the vortex panel strength and circulation quantities; the evaluation of such quantities allows one to calculate the velocity distribution over the airfoil surface, and hence, to determine the pressure coefficients. The Kutta condition (a viscous boundary condition based on physical observation used with inviscid theoretical model) states that the flow leaves the sharp trailing edge of an airfoil smoothly [7]. Various methods have been proposed to impose the Kutta condition [8–10]. In panel methods, the Kutta condition is incorporated in the numerical formulation by requiring that the strength of vortex sheet is zero at the airfoil trailing edge. Moreover, the use of panel methods along with the compressibility corrections such as the Prandtl-Glauert method [11] allow one to consider the compressible flows over bodies. The panel methods have been extensively investigated in the aerodynamics literature, so these will not be discussed further here. However, these methods often have trouble with accuracy at the trailing edge of airfoils with zero angle cusped trailing edges [12]. Moreover, the compressibility corrections do not give accurate results for the compressible flows over airfoils of any shape at any angle of attack. For example, the Prandtl-Glauert method is based on the linearized perturbation velocity potential equation, and hence, it is limited to thin airfoils at small angles of attack [4].

In this paper, we propose a novel method to numerically solve the steady irrotational compressible flow over an airfoil which is exempt from considerations of quantities such as the vortex panel strength and circulation. This method takes advantage of an O-grid, generated by the elliptic grid generation technique, over the flow field and approximates the flow field quantities such as stream function, density, velocity, pressure, speed of sound, and Mach number at the nodal points. An accurate Kutta condition scheme is proposed and implemented into the computational loop by an exact derived expression for the stream function. The exact expression is general and encompasses both the finite-angle and cusped trailing edges. Finally, the results obtained from the proposed numerical method, and the results from experimental and other numerical methods, are compared to reveal the accuracy of the proposed method (The material given in this article are implemented in the code FOILcom which is freely available at: DOI:10.13140/RG.2.2.36459.64801/1).

It is worth emphasizing that the stream function equation considered in this study is completely equivalent to the full potential equation. Here, the stream function equation is solved in non-conservative form and can be employed to obtain accurate results for subsonic (subcritical) inviscid isentropic compressible flow over isolated airfoils. The conservative stream function equation, along with upwinding schemes such as artificial compressibility [13], can be used to solve transonic flows over airfoils which is not considered in this study.

#### **2. Governing Equations**

The stream function equation for two-dimensional, irrotational, steady, and isentropic flow of a compressible fluid in non-conservative form is as follows [14–16]

$$(c^2 - u^2)\psi\_{xx} + (c^2 - v^2)\psi\_{yy} - 2uv\psi\_{xy} = 0\tag{1}$$

where ψ is the stream function, *u* and *v* are the components of the velocity vector **V**, i.e.,**V** = *u***i** + *v***j** (**i** and **j** are the unit vectors in *x* and *y* directions, respectively). For a two-dimensional compressible flow,

$$
\mu = \frac{\rho\_0}{\rho} \psi\_y \tag{2}
$$

$$w = -\frac{\rho\_0}{\rho} \psi\_x \tag{3}$$

*c* is the local sound speed [4]

$$c^2 = c\_0^2 - \frac{\gamma - 1}{2}V^2 = c\_0^2 - \frac{\gamma - 1}{2}(u^2 + v^2) \tag{4}$$

*c*<sup>0</sup> is the stagnation speed of sound, ρ is the density, ρ<sup>0</sup> is the stagnation density, and γ = *cp*/*cv* (*cv* and *cp* are the specific heats at constant volume and constant pressure, respectively) is the ratio of specific heats of gas (for air at standard conditions, γ = 1.4). Dividing both sides of Equation (1) by *c*<sup>2</sup> and then substituting the expressions in Equations (2) and (3) into Equation (1) gives the stream function equation as

$$[1 - \frac{1}{c^2}(\frac{\rho\_0}{\rho})^2 \psi\_y^2] \psi\_{\text{xx}} + [1 - \frac{1}{c^2}(\frac{\rho\_0}{\rho})^2 \psi\_x^2] \psi\_{yy} + \frac{2}{c^2}(\frac{\rho\_0}{\rho})^2 \psi\_x \psi\_y \psi\_{xy} = 0 \tag{5}$$

by substituting Equation (4) into Equation (5), we get

$$\begin{cases} 1 - \frac{1}{c\_0^2 - \frac{\gamma - 1}{2} \left(\frac{\rho\_0}{\rho}\right)^2 (\psi\_x^2 + \psi\_y^2)} \left(\frac{\rho\_0}{\rho}\right)^2 \psi\_y^2 \vert \psi\_{xx} + \left[1 - \frac{1}{c\_0^2 - \frac{\gamma - 1}{2} \left(\frac{\gamma \mu}{\rho}\right)^2 (\psi\_x^2 + \psi\_y^2)} (\frac{\rho\_0}{\rho})^2 \psi\_x^2\right] \psi\_{xy} + \frac{\gamma - 1}{c\_0^2 - \frac{\gamma - 1}{2} \left(\frac{\gamma \mu}{\rho}\right)^2 (\psi\_x^2 + \psi\_y^2)} (\frac{\rho\_0}{\rho})^2 \psi\_x \psi\_y \psi\_{xy} = 0 \end{cases} \tag{6}$$

and from the local isentropic stagnation properties of an ideal gas, we have

$$\frac{\rho\_0}{\rho} = \left(1 + \frac{\mathcal{V} - 1}{2} M^2\right)^{\frac{1}{\mathcal{V} - 1}} = \left(1 + \frac{\mathcal{V} - 1}{2} \frac{u^2 + v^2}{c^2}\right)^{\frac{1}{\mathcal{V} - 1}}\tag{7}$$

*M* = *<sup>V</sup> <sup>c</sup>* is the local Mach number. Equations (6) and (7) should be solved simultaneously to obtain the local values of ψ and ρ in the flow field.

#### *2.1. Transformation*

The solution of the governing PDE is based on transformation of the physical domain (*x*, *y*) and the governing equations into the regular computational domain (ξ, η) (see Figure 1). Therefore, the derivatives such as ψ*x*, ψ*y*, ψ*xx*, ψ*yy*, and ψ*xy* in the stream function equation, Equation (1), should be transformed from the physical domain (*x*, *y*) to the computational domain (ξ, η) [17–19]. This transformation can be stated as

$$
\xi \equiv \xi(x, y) \tag{8}
$$

$$
\eta \equiv \eta(\mathbf{x}, \mathbf{y}) \tag{9}
$$

and the inverse transformation is given as below.

$$\mathbf{x} \equiv \mathbf{x}(\boldsymbol{\xi}, \boldsymbol{\eta}) \tag{10}$$

$$y \equiv y(\xi, \eta) \tag{11}$$

Since the stream function equation involves first and second derivatives, relationships are needed to transform such derivatives from the (*x*, *y*) system to the (ξ, η) one. In order to do this, the Jacobian of the transformation is needed, which is given below

$$\text{2D}: J = J(\frac{\mathbf{x}\_{\prime}y}{\xi\_{\prime}\eta}) = \begin{vmatrix} \mathbf{x}\_{\xi} & y\_{\xi} \\ \mathbf{x}\_{\eta} & y\_{\eta} \end{vmatrix} = \mathbf{x}\_{\xi}y\_{\eta} - \mathbf{x}\_{\eta}y\_{\xi} \neq \mathbf{0} \tag{12}$$

As will be shown, the transformation relations involve the Jacobian in denominator. Hence, it cannot be zero. Since we deal with the stream function equation, it is necessary to find relationships for the transformation of the first and second derivatives of the variable ψ with respect to the variables *x* and *y*. By using the chain rule, it can be concluded that

$$\frac{\partial\psi}{\partial\mathbf{x}} = \frac{\partial\psi}{\partial\xi}\frac{\partial\xi}{\partial\mathbf{x}} + \frac{\partial\psi}{\partial\eta}\frac{\partial\eta}{\partial\mathbf{x}} = \frac{\partial\psi}{\partial\xi}\xi\_{\mathbf{x}} + \frac{\partial\psi}{\partial\eta}\eta\_{\mathbf{x}}\tag{13}$$

$$\frac{\partial\psi}{\partial y} = \frac{\partial\psi}{\partial\xi}\frac{\partial\xi}{\partial y} + \frac{\partial\psi}{\partial\eta}\frac{\partial\eta}{\partial y} = \frac{\partial\psi}{\partial\xi}\xi\_y + \frac{\partial\psi}{\partial\eta}\eta\_y \tag{14}$$

By interchanging *x* and ξ, and *y* and η, the following relationships can also be derived

$$\frac{\partial\psi}{\partial\xi} = \frac{\partial\psi}{\partial x}\frac{\partial x}{\partial\xi} + \frac{\partial\psi}{\partial y}\frac{\partial y}{\partial\xi} = \frac{\partial\psi}{\partial x}\mathbf{x}\_{\xi} + \frac{\partial\psi}{\partial y}\mathbf{y}\_{\xi} \tag{15}$$

$$\frac{\partial \psi}{\partial \eta} = \frac{\partial \psi}{\partial \mathbf{x}} \frac{\partial \mathbf{x}}{\partial \eta} + \frac{\partial \psi}{\partial y} \frac{\partial y}{\partial \eta} = \frac{\partial \psi}{\partial \mathbf{x}} \mathbf{x}\_{\eta} + \frac{\partial \psi}{\partial y} y\_{\eta} \tag{16}$$

By solving Equations (15) and (16) for ∂ψ <sup>∂</sup>*<sup>x</sup>* and ∂ψ <sup>∂</sup>*<sup>y</sup>* , we finally obtain

$$\frac{\partial \psi}{\partial \mathbf{x}} = \frac{1}{l} (y\_{\eta} \frac{\partial \psi}{\partial \xi} - y\_{\xi} \frac{\partial \psi}{\partial \eta}) \tag{17}$$

$$\frac{\partial \psi}{\partial y} = \frac{1}{f} (-\mathbf{x}\_{\eta} \frac{\partial \psi}{\partial \xi} + \mathbf{x}\_{\xi} \frac{\partial \psi}{\partial \eta}) \tag{18}$$

where *J* = *x*ξ*y*<sup>η</sup> − *x*η*y*<sup>ξ</sup> is Jacobian of the transformation. By comparing Equations (13) and (17), and (14) and (18), it can be shown that

$$\mathbf{x}\_x = \frac{1}{\mathbf{J}} y\_\eta, \mathbf{x}\_y = -\frac{1}{\mathbf{J}} \mathbf{x}\_\eta \tag{19}$$

$$
\eta\_{\rm x} = -\frac{1}{J} y\_{\xi}, \eta\_{\rm y} = \frac{1}{J} x\_{\xi} \tag{20}
$$

To transform terms in the stream function equation (Equation (6)), the second order derivatives are needed. Therefore,

$$\psi\_{\rm xx} = \left(\psi\_{\rm x}\right)\_{\rm x} = \left(\frac{1}{l}(y\_{\eta}\psi\_{\xi} - y\_{\xi}\psi\_{\eta})\right)\_{\rm x} \tag{21}$$

$$\psi\_{\mathcal{Y}\mathcal{Y}} = \left(\psi\_{\mathcal{Y}}\right)\_{\mathcal{Y}} = \left(\frac{1}{I}(-\mathbf{x}\_{\mathcal{U}}\psi\_{\mathcal{E}} + \mathbf{x}\_{\mathcal{E}}\psi\_{\mathcal{V}})\right)\_{\mathcal{Y}} \tag{22}$$

$$\psi\_{xy} = (\psi\_x)\_y = \left(\frac{1}{I}(y\_\eta \psi\_\xi - y\_\xi \psi\_\eta)\right)\_y \tag{23}$$

which result in the following expressions for the second order derivatives

$$\psi\_{\mathbf{xx}} = \frac{(y\_{\eta}^2 \psi\_{\xi\xi} - 2y\_{\xi}y\_{\eta}\psi\_{\xi\eta} + y\_{\xi}^2 \psi\_{\eta\eta})}{l^2} + \frac{(y\_{\eta}^2 y\_{\xi\xi} - 2y\_{\xi}y\_{\eta}y\_{\xi\eta} + y\_{\xi}^2 y\_{\eta\eta})(x\_{\eta}\psi\_{\xi\xi} - x\_{\xi}\psi\_{\eta})}{l^3} + \dotsb}{\frac{(y\_{\eta}^2 x\_{\xi\xi} - 2y\_{\xi}y\_{\eta}x\_{\xi\eta} + y\_{\xi}^2 x\_{\eta\eta})(y\_{\xi}\psi\_{\eta} - y\_{\eta}\psi\_{\xi})}{l^3}} + \dotsb} \tag{24}$$

$$\psi\_{yy} = \frac{(\mathbf{x}\_{\eta}^{2}\psi\_{\xi\xi} - 2\mathbf{x}\_{\xi}\mathbf{x}\_{\eta}\psi\_{\xi\eta} + \mathbf{x}\_{\xi}^{2}\psi\_{\eta\eta})}{l^{2}} + \frac{(\mathbf{x}\_{\eta}^{2}y\_{\xi\xi} - 2\mathbf{x}\_{\xi}\mathbf{x}\_{\eta}y\_{\xi\eta} + \mathbf{x}\_{\xi}^{2}y\_{\eta\eta})(\mathbf{x}\_{\eta}\psi\_{\xi} - \mathbf{x}\_{\xi}\psi\_{\eta})}{l^{3}} + \mathbf{x}\_{\xi}^{3} \\ \tag{25}$$

$$\begin{split} \psi\_{xy} &= \frac{1}{f} \Big[ -\chi\_{\eta} \frac{(y\_{\xi\eta}\psi\_{\xi} + y\_{\eta}\psi\_{\xi\xi} - y\_{\xi\xi}\psi\_{\eta} - y\_{\xi}\psi\_{\xi\eta})}{f} - \frac{(y\_{\eta}\psi\_{\xi} - y\_{\xi}\psi\_{\eta})(x\_{\xi\xi}y\_{\eta} + x\_{\xi}y\_{\xi\eta} - x\_{\xi\eta}y\_{\xi} - x\_{\eta}y\_{\xi\xi})}{f^{2}} \Big] + \\ &\propto \Big( \frac{y\_{\eta\eta}\psi\_{\xi} + y\_{\eta}\psi\_{\xi\eta} - y\_{\xi\eta}\psi\_{\eta} - y\_{\xi}\psi\_{\eta\eta}}{f} - \frac{(y\_{\eta}\psi\_{\xi} - y\_{\xi}\psi\_{\eta})(x\_{\xi\eta}y\_{\eta} + x\_{\xi}y\_{\eta\eta} - x\_{\eta\eta}y\_{\xi} - x\_{\eta}y\_{\xi\eta})}{f^{2}} \Big) \end{split} \tag{26}$$

The finite-difference method can be used to discretize the above expressions in the regular computational domain (ξ, η). The actual values of ξ and η in the computational domain are immaterial, because they do not appear in the final expressions. Thus, without a loss of generality, we can select the coordinates of the node *A* in the computational domain as ξ = η = 1 and the mesh size as Δξ = Δη = 1 [19]. Therefore, we have [17]

$$(f\_{\vec{\xi}})\_{i,j} = \frac{1}{2}(f\_{i+1,j} - f\_{i-1,j})\tag{27}$$

$$(f\_{\eta})\_{i,j} = \frac{1}{2}(f\_{i,j+1} - f\_{i,j-1})\tag{28}$$

$$(f\_{\mathbb{E}\xi})\_{i,j} = (f\_{i+1,j} - 2f\_{i,j} + f\_{i-1,j})\tag{29}$$

$$(f\_{\eta\eta})\_{i,j} = (f\_{i,j+1} - 2f\_{i,j} + f\_{i,j-1})\tag{30}$$

$$(f\_{\mathbb{E}\eta})\_{i,j} = \frac{1}{4}(f\_{i+1,j+1} - f\_{i-1,j+1} - f\_{i+1,j-1} + f\_{i-1,j-1})\tag{31}$$

where *f* ≡ *x*, *y*,ψ.

**Figure 1.** The physical and the computational domains. (**a**) Physical domain. (**b**) Computational domain.

#### *2.2. Boundary Conditions*

*Conditions at infinity*: Far away from the airfoil surface (toward infinity), in all directions, the flow approaches the free stream conditions. The known free stream conditions are the velocity *V*∞, the pressure *p*∞, the density ρ∞, and the temperature *T*∞. Thus, the free stream Mach number is

$$M\_{\infty} = \frac{V\_{\infty}}{\mathcal{L}\_{\infty}} \prec M\_{\text{critical}} \tag{32}$$

and the speed of sound *c*<sup>∞</sup> is

$$x\_{\circ \circ} = \sqrt{\gamma RT\_{\circ \circ}} \tag{33}$$

where *R* is the specific gas constant, which is a different value for different gases. For air at standard conditions, *<sup>R</sup>*air <sup>=</sup> 287J/(kg.K), <sup>γ</sup>air <sup>=</sup> 1.4, <sup>ρ</sup><sup>∞</sup> <sup>=</sup> 1.23kg/m3, and *<sup>p</sup>*<sup>∞</sup> <sup>=</sup> 1.01 <sup>×</sup> <sup>10</sup><sup>5</sup> N/m2 . The critical Mach number *M*critical is that free stream Mach number at which sonic flow is first achieved on the airfoil surface.

*Condition on the airfoil surface*: The relevant boundary condition at the airfoil surface for the inviscid flow is the no-penetration boundary condition. Thus, the velocity vector must be tangential to the surface. This wall boundary condition can be expressed by

$$\frac{\partial \psi}{\partial s} = 0 \text{or} \psi = \text{constant} \tag{34}$$

where *s* is tangent to the surface.

#### *2.3. Grid Generation*

Here, an O-grid is initially generated around the airfoil using elliptic grid generation method [17] and then the stream function equation is solved in the computational domain. The size of mesh is *M* × *N* where *M* is the number of nodes on the branch cut (a horizontal line connecting the trailing edge to outer boundary) and *N* is the number of nodes on the airfoil surface (and hence, the outer boundary), as shown in Figure 1. The physical domain before and after meshing using an O-grid of size 110 × 101 is shown in Figure 2 (using a NACA 0012 airfoil). Moreover, the magnified view of different parts of domain is shown in Figure 3 to highlight the orthogonality of the gridlines at airfoil surface and outer boundary.

**Figure 2.** The physical domain before and after meshing. (**a**) Physical domain before meshing. (**b**) Physical domain after meshing.

**Figure 3.** Magnified view of different parts of domain.

#### *2.4. Kutta Condition*

In inviscid flows, because the flow cannot penetrate the surface, the velocity vector must be tangential to the surface. In other words, the component of velocity normal to the surface must be zero and only the tangential velocity component must be considered. The unit tangent vector on the airfoil surface can be expressed as

$$
\pi^{(\xi)} = \mathbf{n}^{(\xi)} \times \mathbf{k} \tag{35}
$$

$$\mathbf{r}^{(\eta)} = \mathbf{n}^{(\eta)} \times \mathbf{k} \tag{36}$$

where **n**(ξ) and **n**(η) are the outward-pointing unit normal vector to a airfoil surface in ξ and η directions, respectively, and **k** is the unit vector in *z* direction.

At the airfoil surface, corresponding to surface *S*<sup>3</sup> in Figure 4, we have

$$\mathbf{n}^{(\xi)} = -\frac{\nabla \xi}{|\nabla \xi|} = -\frac{\xi\_x \mathbf{i} + \xi\_y \mathbf{j}}{\sqrt{\xi\_x^2 + \xi\_y^2}}\tag{37}$$

**Figure 4.** The outward—pointing unit normal vectors to ξ = constant and η = constant lines.

From the transformation relationships (Equation (19)), we have

$$\mathbf{n}^{(\xi)} = -\frac{\frac{1}{f}y\_{\eta}\mathbf{i} + \left(-\frac{1}{f}x\_{\eta}\mathbf{j}\right)}{\sqrt{\left(\frac{y\_{\eta}}{f}\right)^2 + \left(-\frac{x\_{\eta}}{f}\right)^2}} = \frac{-y\_{\eta}\mathbf{i} + x\_{\eta}\mathbf{j}}{\sqrt{\alpha}}\tag{38}$$

where α = *x*<sup>η</sup> <sup>2</sup> + *y*<sup>η</sup> 2. Using Equation (35), we get

$$\mathbf{r}\_{\mathbf{S}\_3}^{(\ell)} = \mathbf{n}\_{\mathbf{S}\_3}^{(\ell)} \times \mathbf{k} = \frac{1}{\sqrt{\alpha}} (-y\_\uparrow \mathbf{i} + \mathbf{x}\_\uparrow \mathbf{j}) \times \mathbf{k} = \frac{1}{\sqrt{\alpha}} (-y\_\uparrow (-\mathbf{j}) + \mathbf{x}\_\uparrow \mathbf{i}) = \frac{1}{\sqrt{\alpha}} (\mathbf{x}\_\uparrow \mathbf{i} + y\_\uparrow \mathbf{j}) \tag{39}$$

The velocity component tangential to the airfoil surface (*S*3) is

$$V\_t = \mathbf{V}.\tau\_{\mathbb{S}\_3}^{(\xi)} = (u\mathbf{i} + v\mathbf{j}). \frac{1}{\sqrt{\alpha}}(x\_\eta \mathbf{i} + y\_\eta \mathbf{j}) \tag{40}$$

For compressible flows, the velocity components of *u* and *v* can be expressed in terms of stream function ψ as

$$
\mu = \frac{\rho\_0}{\rho} \psi\_y = \frac{\rho\_0}{\rho} \frac{1}{l} (-\mathbf{x}\_{\eta} \psi\_{\xi} + \mathbf{x}\_{\xi} \psi\_{\eta}) \tag{41}
$$

and

$$
\psi = -\frac{\rho\_0}{\rho} \psi\_x = -\frac{\rho\_0}{\rho} \frac{1}{J} (\mathcal{y}\_\eta \psi\_\xi - \mathcal{y}\_\xi \psi\_\eta) \tag{42}
$$

where ρ<sup>0</sup> is a constant. On the airfoil surface, ψ*<sup>s</sup>* = ψη = 0 where *s* is the distance measured along the airfoil surface because the airfoil contour is a streamline of the flow. Thus Equations (41) and (42) become

$$
\mu = \frac{\rho\_0}{\rho} \frac{1}{J} (-\mathbf{x}\_{\overline{\eta}} \psi\_{\xi}) \tag{43}
$$

and

$$w = -\frac{\rho\_0}{\rho} \frac{1}{J} (y\_\eta \psi\_\xi) \tag{44}$$

By substituting Equations (43) and (44) into Equation (40), we get

$$V\_t = ( -\frac{\rho\_0}{\rho} \frac{1}{\rho} \mathbf{x}\_\eta \boldsymbol{\psi}\_\xi \mathbf{i} - \frac{\rho\_0}{\rho} \frac{1}{\rho} \mathbf{y}\_\eta \boldsymbol{\psi}\_\xi \mathbf{j}) \cdot \frac{1}{\sqrt{a}} (\mathbf{x}\_\eta \mathbf{i} + \mathbf{y}\_\eta \mathbf{j}) = -\frac{\rho\_0}{\rho} \frac{1}{\sqrt{a}} \mathbf{x}\_\eta^2 \boldsymbol{\psi}\_\xi - \frac{\rho\_0}{\rho} \frac{1}{\rho} \frac{1}{\sqrt{a}} \mathbf{y}\_\eta^2 \boldsymbol{\psi}\_\xi = -\frac{\rho\_0}{\rho} \frac{1}{\rho} \frac{1}{\sqrt{a}} (a \mathbf{x})$$

Based on the Kutta condition, the velocity at the trailing edge of the airfoil must be the same when approached from the upstream direction along the upper and lower airfoil surfaces [20] (see Figure 5). Thus,

$$\left. V\_{l} \right|\_{l,1} = \left. V\_{l} \right|\_{l,\mathbb{N}} = \mathbb{C} \begin{cases} \mathbb{C} = 0 \text{ for finite angle trailing edge} \\ \mathbb{C} \neq 0 \text{ for cusped trailing edge} \end{cases} \tag{46}$$

For the finite angle trailing edge, we have

$$\left.V\_t\right|\_{1,1} = 0$$

$$-\frac{\rho\_0}{\rho} \frac{\sqrt{\alpha}}{J} \psi\_\xi \Big|\_{1,1} = 0$$

hence

$$
\left.\psi\_{\xi}\right|\_{1,1} = 0\\\psi\_{2,1} - \psi\_{1,1} = 0
$$

$$
\psi\_{2,1} = \psi\_{1,1} \tag{47}
$$

(The above expression is based on forward finite-difference coefficient with first-order accuracy. The forward finite-difference coefficient with second-order accuracy can be used if more accurate Kutta condition implementation is needed.)

Therefore, the wall boundary condition for the finite angle trailing edge becomes

$$
\psi\_{1,j} = \psi\_{2,1}(j = 1, \dots, N) \tag{48}
$$

because the stream function on the airfoil surface is constant.

For the cusped trailing edge, we have

$$|V\_{l}|\_{1,1} = |V\_{l}|\_{1,\text{N}} - \frac{\rho\_{0}}{\rho} \frac{\sqrt{\alpha}}{J} \psi\_{\xi} \Bigg|\_{1,1} = -\frac{\rho\_{0}}{\rho} \frac{\sqrt{\alpha}}{J} \psi\_{\xi} \Bigg|\_{1,\text{N}} \tag{49}$$

The points (1, 1) and (1, N) as well as the points (2, 1) and (2, N) are the same points (see Figure 5). Therefore, ψ1,1 = ψ1,N and ψ2,1 = ψ2,N

$$\left. \psi\_{\xi} \right|\_{1,1} = \psi\_{2,1} - \psi\_{1,1} = \psi\_{2,\mathsf{N}} - \psi\_{1,\mathsf{N}} = \left. \psi\_{\xi} \right|\_{1,\mathsf{N}}$$

On the other hand, we know that

$$-\frac{\rho\_0}{\rho} \left. \frac{\sqrt{\alpha}}{J} \right|\_{1,1} \neq \left. -\frac{\rho\_0}{\rho} \frac{\sqrt{\alpha}}{J} \right|\_{1,\mathcal{N}} \tag{50}$$

Thus, the only possibility to satisfy Equation (49) is that

$$\left.\psi\_{\xi}\right|\_{1,1} = \left.\psi\_{\xi}\right|\_{1,\mathcal{N}} = 0$$

which this relation again results in the same expression as for the finite angle trailing edge case. So, the general expression to implement the Kutta condition based on the stream function for the 2D compressible flow is

$$
\psi\_{1,1} = \psi\_{2,1} \tag{51}
$$

which is the same as one for the incompressible flow [21]

**Figure 5.** Grid notation and tangential velocity components at trailing edge. (**a**) Finite angle trailing edge. (**b**) Cusped trailing edge.

#### *2.5. Computation Procedure*

According to the mapping scheme adopted in Figure 1, there are four sections where the nodal value of the flow variables *fi*,*<sup>j</sup>* (*f* ≡ ψ, ρ, *p*, *u*, *v*, *V*, *c*, *M*, ...) should be calculated.


The known free-stream variables are *M*∞, *p*∞, ρ∞, *T*∞, and the angle of attack α. Thus, we can write

$$M\_{\rm M,j} = M\_{\rm co}(j = 1, \ldots, N) \tag{52}$$

$$T\_{M,j} = T\_{\infty}(j = 1, \dots, N) \tag{53}$$

$$
\rho\_{\rm M,j} = \rho\_{\rm \circledast}(j = 1, \dots, N) \tag{54}
$$

$$p\_{M,j} = p\_{\circ \circ}(j = 1, \dots, N) \tag{55}$$

$$\mathcal{L}\_{M,j} = \mathcal{c}\_{\circ \circ} = \sqrt{\gamma\_{\text{air}} R\_{\text{air}} T\_{\text{co}}} = \sqrt{(1.4)(287)} T\_{M,j} (j = 1, \dots, N) \tag{56}$$

$$V\_{M,j} = V\_{\infty} = M\_{\infty}c\_{\infty} = M\_{M,j}c\_{M,j}(j = 1, \dots, N) \tag{57}$$

from the local isentropic stagnation properties of an ideal gas, we have

$$\frac{\rho\_{0\_{M,j}}}{\rho\_{M,j}} = \left(1 + \frac{\chi\_{\text{air}} - 1}{2} M\_{M,j}^2\right)^{\frac{1}{\chi\_{\text{air}} - 1}} (j = 1, \dots, N) \tag{58}$$

so the local total density ρ0*M*,*<sup>j</sup>* can be computed as

$$\rho\_{0\mathbf{M},j} = \rho\_{M,j} (1 + \frac{\mathcal{V}\_{\text{air}} - 1}{2} \mathcal{M}\_{M,j}^2)^{\frac{1}{\mathcal{V}\_{\text{air}} - 1}} (j = 1, \dots, N) \tag{59}$$

In a similar fashion, we can write

$$p\_{0\mathbf{M}\_{\mathbf{j}}} = p\_{\mathbf{M},j}(1 + \frac{\gamma\_{\text{air}} - 1}{2} \mathbf{M}\_{\mathbf{M},j}^2)^{\frac{\gamma\_{\text{air}}}{\gamma\_{\text{air}} - 1}} (j = 1, \dots, N) \tag{60}$$

$$T\_{0\_{M,j}} = T\_{M,j}(1 + \frac{\mathcal{V}\_{\text{air}} - 1}{2} M\_{M,j}^2)(j = 1, \dots, N) \tag{61}$$

$$\mathcal{L}\_{\text{0M},j} = \sqrt{\gamma\_{\text{air}} R\_{\text{air}} T\_{\text{0M},j}} = \sqrt{(1.4)(287)} T\_{\text{0M},j} \text{ ( $j = 1, \dots, N$ )}\tag{62}$$

We know that if the general flow field is isentropic throughout, then the local total density ρ0, the local total pressure *p*0, and the local total temperature *T*<sup>0</sup> are constant values throughout the flow. Thus

$$\mathfrak{c}\_{0\_{\vec{i}\vec{j}}} = \mathfrak{c}\_{0\mathcal{M}\_{\vec{i}}}(\vec{i} = 1, \dots, \mathcal{M}, \vec{j} = 1, \dots, \mathcal{N}) \tag{63}$$

$$\rho\_{0\_{i\ne j}} = \rho\_{0M\_{\ne i}}(i = 1, \dots, M, j = 1, \dots, N) \tag{64}$$

$$p\_{0\_{ij}} = p\_{0\_{M,j}}(i = 1, \dots, M, j = 1, \dots, N) \tag{65}$$

$$T\_{0\_{\vec{i},\vec{j}}} = T\_{0\_{M,\vec{j}}}(\mathbf{i} = 1, \dots, M, \mathbf{j} = 1, \dots, N) \tag{66}$$

Also, on the outer boundary (far-field) we have

$$u\_{M,j} = V\_{x\_{\rm av}} = V\_{\rm os} \cos a(j = 1, \dots, N) \tag{67}$$

$$\upsilon\_{\mathcal{M},j} = V\_{\mathcal{Y}\_{\mathcal{O}}} = V\_{\mathcal{O}} \sin a(j = 1, \dots, N) \tag{68}$$

If the number of nodes on the outer boundary (the side *BG* in Figure 1) is *N*, then, as shown in Figure 6, the node number of two points at top (point *E*) and bottom (point *F*) of the circular outer boundary would be (*M*, *<sup>N</sup>*+<sup>3</sup> <sup>4</sup> ) and (*M*, <sup>3</sup>*N*+<sup>1</sup> <sup>4</sup> ), respectively.

**Figure 6.** Representation of node number of top and bottom points on the outer boundary (far-field).

Then we can express the magnitude of the stream function on the outer boundary ψ*M*,*<sup>j</sup>* in terms of far-field velocity *VM*,*<sup>j</sup>* as (by considering the equal magnitudes but opposite in sign for stream functions at top and bottom points of *E* and *F*, that is, ψ*<sup>F</sup>* = −*c* and ψ*<sup>E</sup>* = *c*, *c* is a constant)

$$\begin{split} \mu\_{M,\frac{3\mathbb{N}+1}{4}} &= \frac{\rho\_{0,\frac{3\mathbb{N}+1}{4}}}{\rho\_{M,\frac{3\mathbb{N}+1}{4}}} \frac{\partial \psi}{\partial y}\Big|\_{M,\frac{3\mathbb{N}+1}{4}} = \frac{\rho\_{0,\frac{3\mathbb{N}+1}{4}}}{\rho\_{M,\frac{3\mathbb{N}+1}{4}}} (\frac{\psi\_{M,\frac{3\mathbb{N}+1}{4}} - \psi\_{M,\frac{3\mathbb{N}+1}{4}}}{\rho\_{M,\frac{3\mathbb{N}+1}{4}} - \psi\_{M,\frac{3\mathbb{N}+1}{4}}}) = \\ \frac{\rho\_{M,\frac{3\mathbb{N}+1}{4}}}{\rho\_{M,\frac{3\mathbb{N}+1}{4}}} (\frac{-\psi\_{M,\frac{3\mathbb{N}+1}{4}} - \psi\_{M,\frac{3\mathbb{N}+1}{4}}}{\rho\_{M,\frac{3\mathbb{N}+1}{4}} - \psi\_{M,\frac{3\mathbb{N}+1}{4}}}) = -\frac{\rho\_{M,\frac{3\mathbb{N}+1}{4}}}{\rho\_{M,\frac{3\mathbb{N}+1}{4}}} (\frac{2\psi\_{M,\frac{3\mathbb{N}+1}{4}}}{\psi\_{M,\frac{3\mathbb{N}+1}{4}} - \psi\_{M,\frac{3\mathbb{N}+1}{4}}}) \end{split} \tag{69}$$

Therefore,

$$\psi\_{M\_r \frac{3N+1}{4}} = -\frac{\rho\_{M\_r \frac{3N+1}{4}}}{\rho\_{M\_r \frac{3N+1}{4}}} u\_{M\_r \frac{3N+1}{4}} \frac{y\_{M\_r \frac{3N+3}{4}} - y\_{M\_r \frac{3N+1}{4}}}{2} \tag{70}$$

Now the magnitude of the stream function on the outer boundary can be calculated as a (based on the relation <sup>ψ</sup> <sup>=</sup> <sup>−</sup> <sup>ρ</sup> ρ0 *Vyx* <sup>+</sup> <sup>ρ</sup> ρ0 *Vx y*)

$$\psi\_{M,j} = \psi\_{M,\frac{3N+1}{4}} + \frac{\rho\_{M,j}}{\rho\_{0\_{M,j}}} u\_{M,j} (y\_{M,j} - y\_{M,\frac{3N+1}{4}}) - \frac{\rho\_{M,j}}{\rho\_{0\_{M,j}}} v\_{M,j} (x\_{M,j} - x\_{M,\frac{3N+1}{4}}) \tag{71}$$

In Equation (69), we should note that the two points *E* and *F* have the same *x*-coordinates and hence <sup>ψ</sup> <sup>=</sup> <sup>−</sup> <sup>ρ</sup> ρ0 *Vyx* <sup>+</sup> <sup>ρ</sup> ρ0 *Vx <sup>y</sup>* <sup>=</sup> <sup>0</sup> <sup>+</sup> <sup>ρ</sup> ρ0 *Vx <sup>y</sup>* <sup>=</sup> <sup>ρ</sup> ρ0 *Vx y*.

The iterative process to obtain the value of <sup>ψ</sup>*i*,*<sup>j</sup>* and <sup>ρ</sup>*i*,*<sup>j</sup>* may be initiated by assuming <sup>ρ</sup>0*i*,*<sup>j</sup>* <sup>ρ</sup>*i*,*<sup>j</sup>* = 1 or ρ*i*,*<sup>j</sup>* = ρ0*i*,*<sup>j</sup>* . This initial assumption implies that the magnitude of the speed of sound is infinite ( *<sup>c</sup>* = ∞ ⇒ <sup>1</sup> *<sup>c</sup>*<sup>2</sup> = 0 ) and Equation (5) reduces to Laplace's equation

$$
\psi\_{xx} + \psi\_{yy} = 0 \tag{72}
$$

The first iteration is concerned with the solution of Laplace's equation which is described thoroughly in [21] and we do not elaborate further on it. After solving the Laplace's equation to find the value of the stream function at each node, ψ*i*,*j*, in the flow field, the values of the velocity components of *ui*,*<sup>j</sup>* and *vi*,*<sup>j</sup>* can be computed from Equations (2) and (3). Then, the value of the speed of sound at each node, *ci*,*j*,can be calculated from Equation (4). Finally, the value of the density at each node, ρ*i*,*j*, can be determined from Equation (7). If the flow is compressible, ρ*i*,*<sup>j</sup>* ρ0*i*,*<sup>j</sup>* and from the second iteration onwards, instead of Laplace's equation, the stream function given in Equation (6) must be solved by having the density ρ*i*,*<sup>j</sup>* and the stream function ψ*i*,*<sup>j</sup>* from the last iteration (i.e., iteration 1) to get new values of the stream function ψ*i*,*j*. The iterative process (solving Equations (6) and (7) simultaneously) repeated until successive iterations produce a sufficiently small change in density ρ*i*,*<sup>j</sup>* and the stream function ψ*i*,*j*. Equation (6) is solved by initially discretizing the partial differential terms in Equations (21)–(26) using relations given in Equations (27)–(31) and then substituting the terms into Equation (6), and finally, solving the obtained algebraic equation using an algebraic solver (such as Maple) to get an explicit algebraic expression in terms of ψ*i*,*j*. As stated before, the Kutta condition is implemented as

$$
\psi\_{1,1} = \psi\_{2,1} \tag{73}
$$

$$
\psi\_{1,N} = \psi\_{1,1}(\text{thesamerode})\tag{74}
$$

and the wall boundary condition is implemented as

$$
\psi\_{1,j} = \psi\_{1,j-1}(j=2,\ldots,N-1)\tag{75}
$$

and since the branch cut (*AB* or *HG* in Figure 1) is inside the flow field, the same procedure employed to obtain an algebraic expression for ψ*i*,*<sup>j</sup>* can also be used to obtain an expression for the stream function on the branch cut, ψ*i*,1. However, some changes are needed in the terms discretized by the finite-difference method (Equations (27)–(31)) as follows (see Figure 7)

$$(f\_{\xi})\_{i,1} = \frac{1}{2}(f\_{i+1,1} - f\_{i-1,1})\tag{76}$$

$$(f\_{\eta})\_{i,1} = \frac{1}{2}(f\_{i,2} - f\_{i, \mathcal{N}-1})\tag{77}$$

$$(f\_{\xi\xi})\_{i,1} = (f\_{i+1,1} - 2f\_{i,1} + f\_{i-1,1})\tag{78}$$

$$(f\_{\eta\eta})\_{i,1} = (f\_{i,2} - 2f\_{i,1} + f\_{i,N-1})\tag{79}$$

$$\left(f\_{\mathbb{E}\eta}\right)\_{i,1} = \frac{1}{4} \left(f\_{i+1,2} - f\_{i-1,2} - f\_{i+1,N-1} + f\_{i-1,N-1}\right) \tag{80}$$

where *f* ≡ *x*, *y*,ψ.

$$\begin{array}{c} \stackrel{i,2}{\rightsquigarrow} \\\\ \begin{array}{c} \stackrel{i,2}{\rightsquigarrow} \\\\ \stackrel{i-1,1}{\rightsquigarrow} \\\\ \stackrel{i,N}{\rightsquigarrow} \end{array} \end{array} $$

**Figure 7.** Node notation on the branch cut.

and on the branch cut

$$
\psi\_{i,N} = \psi\_{i,1}(i = 2, \dots, M - 1) \tag{81}
$$

In order to solve the elliptic grid generation equations (for *x* and *y*) and the stream function equation, the iterative method *Successive Over Relaxation* (SOR) is used due to its high convergence rate

$$f\_{i,j}^{(k)} = \omega f\_{i,j}^{(k-1)} + (1 - \omega) f\_{i,j}^{(k-2)} \tag{82}$$

where *f* ≡ *x*, *y*,ψ. In these relation for SOR, *k* is iteration number, and ω is *relaxation factor* which has a value of 1 <ω< 2. A value of ω = 1.8 is chosen in this study and the code FOILcom.

We should define the *stopping criteria* for convergence of solution of the stream function equation (Equation (6) and density equation (Equation (7)). These two equations constitute the system of equations to be solved simultaneously. The stopping criteria are defined as follows

$$\lambda\_{\psi} = \sum\_{i=2}^{M-1} \sum\_{j=2}^{N-1} (\psi\_{i,j}^{(k+1)} - \psi\_{i,j}^{(k)})^2 \tag{83}$$

$$\lambda\_{\mathcal{P}} = \sum\_{i=2}^{M-1} \sum\_{j=2}^{N-1} (\rho\_{i,j}^{(k+1)} - \rho\_{i,j}^{(k)})^2 \tag{84}$$

where *k* is iteration number. A value of λψ = λρ = 10−<sup>4</sup> can be considered to get sufficiently accurate results. The other parameters of interest can also be computed after obtaining the density ρ*i*,*<sup>j</sup>* and the stream function ψ*i*,*j*. The components of the velocity at each node, *ui*,*<sup>j</sup>* and *vi*,*j*, can be computed from

$$\mathbf{u}\_{i,j} = \frac{\rho\_{0,j}}{\rho\_{i,j}} \boldsymbol{\psi}\_{\mathcal{Y}}\Big|\_{\mathbf{i},j} = \frac{\rho\_{0,j}}{\rho\_{i,j}} \frac{\mathbf{1}}{\mathbf{J}} (-\mathbf{x}\_{\eta}\boldsymbol{\psi}\_{\mathcal{E}} + \mathbf{x}\_{\mathcal{E}}\boldsymbol{\psi}\_{\eta})\Big|\_{\mathbf{i},j} \tag{85}$$

$$\left. \psi\_{i,j} = -\frac{\rho\_{0\_{i,j}}}{\rho\_{i,j}} \psi\_{\mathbf{x}} \right|\_{i,j} = -\frac{\rho\_{0\_{i,j}}}{\rho\_{i,j}} \frac{1}{\int} (y\_{\eta} \psi\_{\xi} - y\_{\xi} \psi\_{\eta}) \Big|\_{i,j} \tag{86}$$

Then, the velocity at each node, *Vi*,*j*, is given by

$$\mathcal{V}\_{i,j} = \sqrt{\mu\_{i,j}^2 + \upsilon\_{i,j}^2} \tag{87}$$

and the local (nodal) Mach number can be obtained by

$$M\_{i,j} = \frac{V\_{i,j}}{c\_{i,j}}\tag{88}$$

where the local speed of sound,*ci*,*j*, is calculated from Equation (4). And finally, the local pressure, *pi*,*j*, is calculated from

$$\frac{p\_{0\_{i,j}}}{p\_{i,j}} = \left(1 + \frac{\mathcal{Y}\_{\text{air}} - 1}{2} M\_{i,j}^2\right)^{\frac{\mathcal{Y}\_{\text{air}}}{\mathcal{Y}\_{\text{air}} - 1}} \tag{89}$$

Now the values of drag and lift forces on the airfoil can be computed as follows

$$D = A\cos\alpha + N\sin\alpha\tag{90}$$

$$L = -A\sin\alpha + N\cos\alpha\tag{91}$$

where *A* and *N* are axial and normal forces, respectively (Figure 8) [4]. We can write the drag and lift forces as

$$D = \sum\_{j=1}^{N-1} \left( p\_{1,j} (y\_{1,j} - y\_{1,j+1}) \cos \alpha + p\_{1,j} (\mathbf{x}\_{1,j+1} - \mathbf{x}\_{1,j}) \sin \alpha \right) \tag{92}$$

$$L = \sum\_{j=1}^{N-1} \left( -p\_{1,j}(y\_{1,j} - y\_{1,j+1})\sin\alpha + p\_{1,j}(\mathbf{x}\_{1,j+1} - \mathbf{x}\_{1,j})\cos\alpha \right) \tag{93}$$

and the drag and the lift coefficients may be expressed as

$$c\_d = \frac{D}{\frac{1}{2}\rho\_{\infty}V\_{\infty}^2} \tag{94}$$

$$c\_l = \frac{L}{\frac{1}{2}\rho\_{\infty}V\_{\infty}^2} \tag{95}$$

Furthermore, the flowchart of the computational procedure is presented in Figure 9.

The theoretical and computational details presented here for compressible flows are implemented in the freely available code FOILcom. Interested readers can refer to the code and use it by changing the airfoil shape, the free stream subcritical Mach number, and the angle of attack. The freely available code FOILincom(DOI: 10.13140/RG.2.2.21727.15524) also takes advantage of the same computational procedure but only for incompressible flow.

**Figure 8.** Aerodynamic force *R* and the components into which it splits [4].

**Figure 9.** Flowchart of computational procedure.

#### **3. Results**

A few test cases are given to reveal the accuracy and robustness of the numerical scheme. Here, the results from FOILcom are compared with experimental and numerical results (alternative numerical schemes) to show its accuracy and efficiency.

*Test case 1*: Critical Mach number for the NACA 0012 airfoil at zero angle of attack.

The airfoil surface pressure coefficient distribution using FOILcom is shown in Figure 10. The critical Mach number for the NACA 0012 airfoil at zero angle of attack is obtained as 0.722 using FOILcom with a grid of size 110 × 101 and *xM*,1 = 10 m. The critical Much number in references [4,22] is obtained as 0.725 (the experimental data in [22] is investigated in [4]). A comparison of results is shown in Figure 11 which reveals an excellent agreement.

**Figure 10.** The pressure coefficient distribution for NACA0012 using FOILcom (*M*<sup>∞</sup> = *M*critical = 0.722, α = 0).

**Figure 11.** The pressure coefficient distribution. Comparison of FOILcom result and the one from References [4,22].

*Test case 2*: The airfoil surface pressure coefficient distribution for the NACA 0012 using a grid of size 110 × 101, *M*<sup>∞</sup> = 0.5, and α = 3◦ (Figures 12 and 13). The results from FOILcom (Figure 12) and both finite volume and meshless methods (the results from finite volume and meshless methods in [23] are presented in one plot due to an excellent agreement between the results) are compared and depicted in Figure 13.

**Figure 12.** The pressure coefficient distribution for NACA0012 using FOILcom (*M*<sup>∞</sup> = 0.5, α = 3◦).

**Figure 13.** The pressure coefficient distribution. Comparison of FOILcom result and the one from Reference [23].

*Test case 3*: The airfoil surface pressure coefficient distribution for the NACA 0012 using a grid of size 110 × 101, *M*<sup>∞</sup> = 0.63, and α = 2◦ (Figures 14 and 15). In this test case, the result from FOILcom is compared to the one from the code FLO42. There is excellent agreement between the results.

**Figure 14.** The pressure coefficient distribution for NACA0012 using FOILcom (*M*<sup>∞</sup> = 0.63, α = 2◦).

**Figure 15.** The pressure coefficient distribution. Comparison of FOILcom result and the one from FLO42 [24].

*Test case 4*: The airfoil surface pressure coefficient distribution for the NACA 2414. *M*<sup>∞</sup> = 0.435 and α = 2◦ (Figures 16 and 17). In this test case, 19,740 nodes (mesh size of 140 × 141) and 51,792 nodes are used in FOILcom and ANSYS Fluent, respectively. The drag and lift coefficients are calculated as follows:

$$\text{FOILcom}: \mathbf{c}\_d = -4.362 \times 10^{-3}, \mathbf{c}\_l = 0.573$$

$$\text{ANYSYS Fluent}: \mathbf{c}\_d = 9.985 \times 10^{-4}, \mathbf{c}\_l = 0.574$$

which shows perfect agreement. The results are obtained by a FORTRAN compiler and computations are run on a PC with Intel Core i5 and 6G RAM. A tolerance of 10−<sup>7</sup> is used in iterative loops to increase the accuracy of results. The computation time is about 5 min which is high due to the iterative solution of the elliptic grid generation and the stream function equations.

**Figure 16.** The pressure coefficient distribution for NACA2414 using FOILcom (*M*<sup>∞</sup> = 0.435, α = 2◦).

**Figure 17.** The pressure coefficient distribution. Comparison of FOILcom result and the one from ANSYS Fluent (V19.2).

*Test case 5*: The airfoil surface pressure coefficient distribution using different grid sizes for NACA 2214 airfoil and *M*<sup>∞</sup> = 0.55, α = 2◦.

Grid sizes are:


The results are shown in Figure 18. In this test case, four different grid sizes are used to obtain the airfoil surface pressure coefficient distribution. As shown in Figure 18, the distribution can be obtained with reasonable accuracy using even a course grid (30 × 41). The effect of grid size on the drag and lift coefficients is depicted in Figure 19. Moreover, the drag and lift coefficients using the four different grid sizes are compared to the ones from XFOIL (Figure 20) and are given in Table 1. As can be seen, the results from FOILcom are in excellent agreement with the result from XFOIL.

**Figure 18.** The airfoil surface pressure coefficient distribution using different grid sizes for NACA 2214 airfoil (*M*<sup>∞</sup> = 0.55, α = 2◦). The code FOILcom is used.

**Figure 19.** Effect of grid size on drag coefficient (**a**) and lift coefficient (**b**).

**Figure 20.** The airfoil surface pressure coefficient distribution for NACA 2214 airfoil (*M*<sup>∞</sup> = 0.55, α = 2◦) using XFOIL.

**Table 1.** Drag and lift coefficients using different grid sizes. Codes FOILcom and XFOIL are used.


*Test case 6*: The airfoil surface pressure coefficient distribution for NACA 64-012airfoil and *M*<sup>∞</sup> = 0.3, α = 6◦. The mesh size used in FOILcom is 50 × 141. The computation time is 46 s. The results from two codes FOILcom and FOILincom along with the Karman-Tsien compressibility correction are compared and depicted in Figure 21. Moreover, the convergence history of λψ and λρ are shown in Figure 22. As can be seen, in this test case five iterations are needed to satisfy the stopping criteria (λψ = λρ = 10<sup>−</sup>4).

**Figure 21.** Comparison of surface pressure coefficient distributions for NACA 64-012 airfoil (*M*<sup>∞</sup> = 0.3, α = 6◦) using the codes FOILcom and FOILincom along with the Karman-Tsien compressibility correction.

**Figure 22.** Convergence history of λψ and λρ.

*Test case 7*: The airfoil surface pressure coefficient distribution for NACA 2240 airfoil and *M*<sup>∞</sup> = 0.3 at two different angles of attack α = 3◦ and α = 6◦. The mesh size used in FOILcom is 120 × 141. In this test case, a thick airfoil, NACA 2240, is used to reveal the accuracy of the proposed numerical method in dealing with thick airfoils at high angles of attack. The results from FOILcom, FOILincom along with the Karman-Tsien compressibility correction, and XFOIL are compared and depicted in Figure 23 for the angle of attack α = 3◦ and Figure 24 for the angle of attack α = 6◦. Furthermore, the results from XFOIL are given in Figure 25 for the angle of attack α = 3◦ and Figure 26 for the angle of attack α = 6◦. The drag and lift coefficients using both codes for two different angles of attack are given in Table 2.

**Figure 23.** Comparison of surface pressure coefficient distributions for NACA 2240 airfoil (*M*<sup>∞</sup> = 0.3, α = 3◦) using the codes FOILcom, FOILincom along with the Karman-Tsien compressibility correction, and XFOIL.

**Figure 24.** Comparison of surface pressure coefficient distributions for NACA 2240 airfoil (*M*<sup>∞</sup> = 0.3, α = 6◦) using the codes FOILcom, FOILincom along with the Karman-Tsien compressibility correction, and XFOIL.

**Figure 25.** The airfoil surface pressure coefficient distribution for NACA 2240 airfoil (*M*<sup>∞</sup> = 0.3, α = 3◦) using XFOIL.

**Figure 26.** The airfoil surface pressure coefficient distribution for NACA 2240 airfoil (*M*<sup>∞</sup> = 0.3, α = 6◦) using XFOIL.

**Table 2.** Drag and lift coefficients using FOILcom and XFOIL.


#### **4. Conclusions**

The solution of 2D steady, irrotational, subsonic (subcritical) compressible flow over isolated airfoils using the stream function equation and a novel method to implement the Kutta condition have been presented. The numerical scheme takes advantage of transformation of the flow solver and the boundary conditions from the physical domain to the computational domain. The physical domain was meshed by an O-grid elliptic grid generation method and the transformed flow solver is discretized by finite-difference method, a method chosen for its simplicity and ease of implementation. The numerical scheme is exempt from considering the panels and the quantities such as the vortex panel strength and circulation used in the panel method. An accurate Kutta condition scheme is proposed and implemented into the computational loop by an exact derived expression for the stream function at the airfoil trailing edge. The exact expression is general, and encompasses both the finite-angle and cusped trailing edges. Through several test cases, the proposed numerical scheme was validated by results from the experimental and the other numerical methods. The obtained results revealed that the proposed algorithm is very accurate and robust.

**Author Contributions:** Conceptualization, F.M.; Formal analysis, F.M.; Funding acquisition, F.M.; Investigation, F.M., B.E. and M.S.; Methodology, F.M.; Software, F.M.; Validation, F.M.; Visualization, F.M.; Writing—original draft, F.M.; Writing—review and editing, B.E. and M.S.

**Funding:** This research was supported by funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement number 663830.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **Greek symbols**



#### **Subscripts**


#### **Superscript**

k iteration number

#### **Abbreviations**


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Soliton Solution of Schrödinger Equation Using Cubic B-Spline Galerkin Method**

#### **Azhar Iqbal 1,2,\*, Nur Nadiah Abd Hamid <sup>1</sup> and Ahmad Izani Md. Ismail <sup>1</sup>**


Received: 17 April 2019; Accepted: 7 June 2019; Published: 12 June 2019

**Abstract:** The non-linear Schrödinger (NLS) equation has often been used as a model equation in the study of quantum states of physical systems. Numerical solution of NLS equation is obtained using cubic B-spline Galerkin method. We have applied the Crank–Nicolson scheme for time discretization and the cubic B-spline basis function for space discretization. Three numerical problems, including single soliton, interaction of two solitons and birth of standing soliton, are demonstrated to evaluate to the performance and accuracy of the method. The error norms and conservation laws are determined and found to be in good agreement with the published results. The obtained results show that the approach is feasible and accurate. The proposed method has almost second order convergence. The linear stability of the method is performed using the Von Neumann method.

**Keywords:** non-linear Schrödinger equation; cubic B-spline basis functions; Galerkin method

#### **1. Introduction**

The non-linear Schrödinger (NLS) equation describes how the behavior of quantum states of a physical system changes in time and space. The NLS equation can be used to describe the propagation of optical pulses and waves in water and plasmas, among other things. Due to the presence of nonlinearity and the complex nature of the problem, it is still a challenge for researchers to determine the most suitable method. Many theoretical and numerical studies have been carried out to overcome this difficulty. The B-spline Galerkin method is quite advanced and has been used by researchers to solve other complex problems.

Dag [1] presented the quadratic and cubic B-spline Galerkin finite element method for solving the Burger's equation. This method was found to give satisfactory results for the Burger's equation, particularly when continuity of the solutions was essential. Dag et al. [2] also proposed a cubic B-splines collocation method for solving the one-dimensional Burger's equation. The proposed scheme was easy to implement and did not require any inner iteration to deal with the nonlinear term of the Burger's equation. Gorgulu et al. [3] presented exponential B-splines Galerkin finite element method for the numerical solution of the advection-diffusion equation. In this study, a new algorithm was developed for the numerical solution of differential equations. This algorithm was obtained by utilizing exponential B-spline functions for the Galerkin finite element method. It was reported that the proposed method gave satisfactory results. The exponential B-splines Galerkin method for the numerical solution of the Burger's equation was also applied by Gorgulu et al. [4].

Saka et al. [5] proposed a quartic B-splines Galerkin finite element method for solving the regularized long wave equation. The performance and accuracy of the method were observed. The method has a weakness because of its large number of matrix operations. The quartic B-spline method was found to be very useful for finding the numerical solutions of the differential equation when higher continuity of solutions occurs. Aksan [6] presented the quadratic B-spline finite element

method for finding the approximate solution of the nonlinear Burger's equation. The high accuracy of the proposed method was thoroughly examined. Gardner et al. [7] also applied a cubic B-spline finite element method for the numerical solution of the Burger's equation. They noticed that the proposed finite element method gave more accurate results than other approaches.

Sheng et al. [8] solved the NLS equation using a finite difference method based on quartic spline approximation and semidiscretization. The NLS equation was approximated using a finite element method by Zlotnik and Zlotnik [9]. Ersoy et al. [10] studied the numerical solution of the NLS equation using an exponential B-spline with collocation method. In that paper, they used the Crank Nicolson formulation for time integration and exponential cubic B-spline functions for space integration. Mokhtari et al. [11] solved the NLS equation using delta-shaped basis functions. Saka [12] solved the NLS using the collocation method with quintic B-spline functions. The generating function of the Hermite polynomial and the orthogonal Hermite function were introduced with ordinary Hermite polynomials by Cesarano et al. [13]. The nature of Hermite polynomials is described in the Schrödinger Equation. Analytical solutions of the NLS equation for certain initial and boundary conditions were presented by Zakharov and Shabat [14]. The numerical solution of the NLS equation was obtained by Lin [15] using a septic spline function method with the collocation method. The Galerkin finite element method for the NLS equation was applied by Dag, with quadratic B-spline functions as the weight and trial functions over finite elements [16]. The numerical solution of the NLS equation was also obtained by using the collocation method based on cubic B-spline by Gardner et al. [17]. Hu [18] presented a conservative two-grid finite element method for the numerical solution of the NLS equation and applied one Newton iteration for the linearization of the problem by using the coarse-grid solution as the initial guess. The Galerkin B-spline method is not often used for the numerical solution of the NLS equation. A review of some numerical methods for solving the NLS equation indicates that most methods use the collocation finite element method.

In this paper, we study the use of the Galerkin method with cubic B-spline function as the weight and trial functions over finite elements to solve the NLS equation. This paper is organized as follows. In Section 2, we discussed the governing equation and fundamentals of Cubic B-spline Galerkin method. In Section 3, the stability analysis of the numerical scheme is analyzed using the Von Neumann method. The numerical results and test problems are discussed in Section 4. In the last section, the conclusions are given.

#### **2. Governing Equation and Cubic B-Spline Galerkin Method**

The NLS equation is

$$
\mu \left| \mathbf{u} \mu + \mathbf{u}\_{\text{xx}} + q|\mathbf{u}|^2 \mathbf{u} = \mathbf{0}, \tag{1}
$$

The equation in (1) is called a self-focusing NLS equation when *q* > 0 and a defocusing NLS equation when *q* < 0. The subscripts *x* and *t* denote the spatial variable and time variable respectively; *<sup>q</sup>* is a positive real parameter, and *<sup>i</sup>* <sup>=</sup> <sup>√</sup> −1. *u*(*x*, *t*) is a complex-valued function and defined over the whole real line. The initial condition is

$$u(\mathbf{x},0) = f(\mathbf{x}), \qquad \mathbf{x} \in [a,b],\tag{2}$$

and the boundary conditions are as follows:

$$u(a,t) = u(b,t) = 0.\tag{3}$$

The solution *u*(*x*, *t*) is decomposed into its real and imaginary parts as

$$u(\mathbf{x},t) = r(\mathbf{x},t) + is(\mathbf{x},t),\tag{4}$$

Substituting Equation (4) into Equation (1), the following coupled partial differential equations are obtained:

$$\begin{cases} s\_t - r\_{xx} - q(r^2 + s^2)r = 0, \\ r\_t + s\_{xx} + q(r^2 + s^2)s = 0. \end{cases} \tag{5}$$

The solution domain [*a*, *b*] is divided equally into N finite elements by the nodes *xk* such that *<sup>a</sup>* = *<sup>x</sup>*<sup>0</sup> <sup>&</sup>lt; *<sup>x</sup>*<sup>1</sup> <sup>&</sup>lt; ··· *xN* = *<sup>b</sup>* and *<sup>h</sup>* = *<sup>b</sup>*−*<sup>a</sup> <sup>N</sup>* = *xk*<sup>+</sup><sup>1</sup> − *xk*. In this paper, we choose cubic B-spline function as the weight and trial function. The cubic B-splines <sup>∅</sup>*k*(*x*),(*<sup>k</sup>* <sup>=</sup> <sup>−</sup>1, ... , *<sup>N</sup>* <sup>+</sup> <sup>1</sup>) are defined over the interval [*a*, *b*] as follows [1]:

$$\mathbf{x} = \frac{1}{h^3} \begin{cases} \left(\mathbf{x} - \mathbf{x}\_{k-2}\right)^3 & \mathbf{x} \in \left[\mathbf{x}\_{k-2}, \mathbf{x}\_{k-1}\right], \\\ h^3 + 3h^2(\mathbf{x} - \mathbf{x}\_{k-1}) + 3h(\mathbf{x} - \mathbf{x}\_{k-1})^2 - 3(\mathbf{x} - \mathbf{x}\_{k-1})^3 & \mathbf{x} \in \left[\mathbf{x}\_{k-1}, \mathbf{x}\_k\right], \\\ h^3 + 3h^2(\mathbf{x}\_{k+1} - \mathbf{x}) + 3h(\mathbf{x}\_{k+1} - \mathbf{x})^2 - 3(\mathbf{x}\_{k+1} - \mathbf{x})^3 & \mathbf{x} \in \left[\mathbf{x}\_k, \mathbf{x}\_{k+1}\right], \\\ \left(\mathbf{x}\_{k+2} - \mathbf{x}\right)^3 & \mathbf{x} \in \left[\mathbf{x}\_{k+1}, \mathbf{x}\_{k+2}\right], \\\ 0 & \text{otherwise}. \end{cases} \tag{6}$$

The approximate solution *uN*(*x*, *t*) to the analytical solution *u*(*x*, *t*) for Equation (1) is written in terms of the cubic B-spline function as

$$u\_N(\mathbf{x}, t) = s\_N(\mathbf{x}, t) + r\_N(\mathbf{x}, t) = \sum\_{j=-1}^{N+1} \delta\_j(t) \mathcal{Q}\_j(\mathbf{x}),\tag{7}$$

where δ*j*(*t*) = ρ*j*(*t*) + ψ*j*(*t*),

$$\begin{aligned} s\_N(\mathbf{x}, t) &= \sum\_{j=-1}^{N+1} \rho\_j(t) \oslash\_j(\mathbf{x}), \\ r\_N(\mathbf{x}, t) &= \sum\_{j=-1}^{N+1} \psi\_j(t) \oslash\_j(\mathbf{x}). \end{aligned} \tag{8}$$

The coefficients ρ*j*(*t*) and ψ*j*(*t*) are unknown time-dependent parameters that will be determined from the boundary and weighted residual conditions. Since each cubic B-spline covers 4 elements, each element [*xk*, *xk*<sup>+</sup>1] is covered by 4 splines. Using the local coordinate's transformation η = *x* − *xk* for 0 ≤ η ≤ *h*, the cubic B-spline function over the element [*xk*, *xk*<sup>+</sup>1] can be written as

$$\begin{aligned} \mathcal{D}\_{k-1} &= \frac{1}{h^3} \begin{cases} (h-\eta)^3 \\ 4h^3 - 3h^2\eta + 3h(h-\eta)^2 - 3(h-\eta)^3 \\ h^3 + 3h^2\eta + 3h\eta^2 - 3\eta^3 \\ \eta^3 \end{cases}, \quad 0 \le \eta \le h. \end{aligned} \tag{9}$$

Here, all other cubic B-spline functions are zero over the finite element [*xk*, *xk*<sup>+</sup>1]. The approximation function (Equation (8)) over the finite element can be defined in terms of the basis functions (Equation (9)) as

$$\begin{aligned} s\_N(\mathbf{x}, t) &= \sum\_{j=k-1}^{k+2} \rho\_j(t) \mathcal{Q}\_j(\mathbf{x}). \\ r\_N(\mathbf{x}, t) &= \sum\_{j=k-1}^{k+2} \psi\_j(t) \mathcal{Q}\_j(\mathbf{x}). \end{aligned} \tag{10}$$

Using the B-splines basis defined in Equation (6) and the trial function stated in Equation (8), the nodal values of *sk*, *rk*,*s <sup>k</sup>* and *r <sup>k</sup>* at the knots *xk* are determined in terms of time-dependent element parameters ρ*<sup>j</sup>* and ψ*<sup>j</sup>* as:

$$\begin{aligned} s\_k &= s(\mathbf{x}\_k) = \rho\_{k-1} + 4\rho\_k + \rho\_{k+1}, \\ r\_k &= r(\mathbf{x}\_k) = \psi\_{k-1} + 4\psi\_k + \psi\_{k+1}, \\ h\mathbf{s}'\_k &= h\mathbf{s}'(\mathbf{x}\_k) == 3(\rho\_{k+1} - \rho\_{k-1}), \\ h\mathbf{r}'\_k &= h\mathbf{r}'(\mathbf{x}\_k) == 3(\psi\_{k+1} - \psi\_{k-1}), \\ h^2\mathbf{s}''\_k &= h^2\mathbf{s}''(\mathbf{x}\_k) = \dots = 6(\rho\_{k-1} - 2\rho\_k + \rho\_{k+1}), \\ h^2\mathbf{r}''\_k &= h^2\mathbf{r}''(\mathbf{x}\_k) = 6(\psi\_{k-1} - 2\psi\_k + \psi\_{k+1}), \end{aligned} \tag{11}$$

where the prime denotes differentiation with respect to *x*.

Applying the Galerkin approach to Equation (5) with the weight function *W*(*x*), the weak form of Equation (5) over the finite element [*xk*, *xk*<sup>+</sup>1] is

$$\begin{cases} \int\_{x\_k}^{x\_{k+1}} [\mathcal{W}s\_t + \mathcal{W}\_{\mathcal{X}}r\_{\mathcal{X}} - z\_k \mathcal{W}r] dx = 0, \\\\ \int\_{x\_{k+1}}^{x\_{k+1}} [\mathcal{W}r\_t - \mathcal{W}\_{\mathcal{X}}s\_{\mathcal{X}} + z\_k \mathcal{W}s] dx = 0, \end{cases} \tag{12}$$

where

$$z\_k = q(r^2 + s^2). \tag{13}$$

Using the weight function *W* as the cubic B-spline and substituting Equations (10) into (12) with some manipulations, the following differential equations are obtained:

$$\begin{cases} \displaystyle \sum\_{j=k-1}^{k+2} \Big[ \left( \int\_{0}^{h} \mathcal{Q}\_{m} \mathcal{Q}\_{j} d\eta \right) \dot{\boldsymbol{\rho}}\_{j} + \left( \int\_{0}^{h} \mathcal{Q}\_{m} \mathcal{Q}\_{j} \mathcal{Q}\_{j} d\eta \right) \dot{\boldsymbol{\psi}}\_{j} - \left( \boldsymbol{z}\_{k} \int\_{0}^{h} \mathcal{Q}\_{m} \mathcal{Q}\_{j} d\eta \right) \dot{\boldsymbol{\psi}}\_{j} \Big] = \boldsymbol{0},\\ \displaystyle \sum\_{j=k-1}^{k+2} \Big[ \left( \int\_{0}^{h} \mathcal{Q}\_{m} \mathcal{Q}\_{j} d\eta \right) \dot{\boldsymbol{\psi}}\_{j} + \left( \int\_{0}^{h} \mathcal{Q}\_{m}^{\prime} \mathcal{Q}\_{j}^{\prime} d\eta \right) \boldsymbol{\rho}\_{j} - \left( \boldsymbol{z}\_{k} \int\_{0}^{h} \mathcal{Q}\_{m} \mathcal{Q}\_{j} d\eta \right) \boldsymbol{\rho}\_{j} \Big] = \boldsymbol{0},\end{cases} \tag{14}$$

Equation (14) can be written in matrix form as

$$\begin{cases} A^{\epsilon} \dot{\rho}^{\varepsilon} + B^{\epsilon} \psi^{\varepsilon} - z\_k C^{\epsilon} \psi^{\varepsilon} = 0, \\\ A^{\epsilon} \dot{\psi}^{\varepsilon} - B^{\epsilon} \rho^{\varepsilon} + z\_k C^{\epsilon} \rho^{\varepsilon} = 0, \end{cases} \tag{15}$$

where <sup>ρ</sup>*<sup>e</sup>* <sup>=</sup> (ρ*k*−1, <sup>ρ</sup>*k*, <sup>ρ</sup>*k*<sup>+</sup>1, <sup>ρ</sup>*k*+2) *<sup>T</sup>* and <sup>ψ</sup>*<sup>e</sup>* <sup>=</sup> (ψ*k*−1,ψ*k*,ψ*k*<sup>+</sup>1,ψ*k*+2) *<sup>T</sup>* are the element parameters and dot (·) represents differentiation with respect to *t*. The element matrix *Ae mj*, *Be mj* and *Ce mj* are given by

$$A\_{mj}^{\mathfrak{c}} = \mathbb{C}\_{mj}^{\mathfrak{c}} = \int\_0^h \mathbb{} \mathbb{}\_{\mathfrak{m}} \mathbb{1} \mathbb{}\_{\mathfrak{h}} d\eta, \quad B\_{mj}^{\mathfrak{c}} = \int\_0^h \mathbb{1} \mathbb{1}'\_{\mathfrak{m}} \mathbb{1} \mathbb{1}'\_{\mathfrak{h}} d\eta. \tag{16}$$

where *m*, *j* = *k* − 1, *k*, *k* + 1, *k* + 2. The element matrix (16) is calculated as

$$\begin{aligned} A^{\varepsilon}\_{\mathrm{mj}} &= \mathcal{C}^{\varepsilon}\_{\mathrm{mj}} = \int\_0^h \mathcal{Q}\_{\mathrm{m}} \mathcal{Q}\_{\mathrm{j}} d\eta = \frac{h}{140} \begin{bmatrix} 20 & 129 & 60 & 1.0 \\ 129 & 1188 & 933 & 60 \\ 60 & 933 & 1188 & 129 \\ 1.0 & 60 & 129 & 20 \end{bmatrix} \end{aligned}$$

$$B^{\varepsilon}\_{\mathrm{mj}} = \int\_0^h \mathcal{Q}^{\prime}\_{\mathrm{m}} \mathcal{Q}^{\prime}\_{\mathrm{j}} d\eta = \frac{1}{10h} \begin{bmatrix} 18 & 21 & -36 & -3.0 \\ 21 & 102 & -87 & -36 \\ -36 & -87 & 102 & 21 \\ -3.0 & -36 & 21 & 18 \end{bmatrix} .$$

*zk* values are calculated by writing *<sup>s</sup>* <sup>=</sup> *sm*+*sm*+<sup>1</sup> <sup>2</sup> and *<sup>r</sup>* <sup>=</sup> *rm*+*rm*+<sup>1</sup> <sup>2</sup> in Equation (13). Using the values of *sN* and *rN* at the points *xk*, we obtain

$$z\_k = q \left[ \frac{(\rho\_{k-1} + 5\rho\_k + 5\rho\_{k+1} + \rho\_{k+2})^2}{4} + \frac{(\psi\_{k-1} + 5\psi\_k + 5\psi\_{k+1} + \psi\_{k+2})^2}{4} \right]. \tag{17}$$

Assembling the contributions from all of the elements, Equation (15) becomes

$$\begin{cases} A\dot{\rho} + B\psi - \mathcal{C}(z\_k)\psi = 0, \\ A\dot{\psi} - B\rho + \mathcal{C}(z\_k)\rho = 0, \end{cases} \tag{18}$$

where ρ = (ρ−1, ρ0, ... , ρ*k*, ρ*k*+1) *<sup>T</sup>* and <sup>ψ</sup> <sup>=</sup> (ψ−1,ψ0, ... , <sup>ψ</sup>*k*,ψ*k*+1) *<sup>T</sup>* are global element parameters and the matrices *A*, *B and C*(*zk*), which are global matrices with a generalized kth row, have the following form:

$$\begin{aligned} A &= \frac{\hbar}{140} (1, 120, 1191, 2416, 1191, 120, 1), \\ B &= \frac{1}{10\hbar} (-3, -72, -45, 240, -45, -72, -3), \\ C(z\_k) &= \frac{\hbar}{140} (z\_{1k}, 60z\_{1k} + 60z\_{2k}, 129z\_{1k} + 933z\_{2k} + 129z\_{3k}, 20z\_{1k} + 1188z\_{2k} \\ &\quad + 1188z\_{3k} + 20z\_{4k}, 129z\_{2k} + 933z\_{3k} + 129z\_{4k}, 60z\_{3k} + 60z\_{4k}, z\_{4k}). \end{aligned} \tag{19}$$

where

$$\begin{aligned} z\_{1k} &= \frac{q}{4} \Big[ \left( \rho\_{k-2} + 5\rho\_{k-1} + 5\rho\_{k} + \rho\_{k+1} \right)^2 + \left( \psi\_{k-2} + 5\psi\_{k-1} + 5\psi\_{k} + \psi\_{k+1} \right)^2 \Big] \\\ z\_{2k} &= \frac{q}{4} \Big[ \left( \rho\_{k-1} + 5\rho\_{k} + 5\rho\_{k+1} + \rho\_{k+2} \right)^2 + \left( \psi\_{k-1} + 5\psi\_{k} + 5\psi\_{k+1} + \psi\_{k+2} \right)^2 \Big] \\\ z\_{3k} &= \frac{q}{4} \Big[ \left( \rho\_{k} + 5\rho\_{k+1} + 5\rho\_{k+2} + \rho\_{k+3} \right)^2 + \left( \psi\_{k} + 5\psi\_{k+1} + 5\psi\_{k+2} + \psi\_{k+3} \right)^2 \Big] \\\ z\_{4k} &= \frac{q}{4} \Big[ \left( \rho\_{k+1} + 5\rho\_{k+2} + 5\rho\_{k+3} + \rho\_{k+4} \right)^2 + \left( \psi\_{k+1} + 5\psi\_{k+2} + 5\psi\_{k+3} + \psi\_{k+4} \right)^2 \Big]. \end{aligned} \tag{20}$$

Discretization by the Crank-Nicolson method gives us <sup>ρ</sup> = ρ*<sup>n</sup>* + ρ(*n*+1) /2 and <sup>ψ</sup> = ψ*<sup>n</sup>* + ψ*n*+<sup>1</sup> /2. Similarly, the time discretization by finite difference scheme gives us . <sup>ρ</sup> = <sup>ρ</sup>*n*+<sup>1</sup> <sup>−</sup> <sup>ρ</sup>*<sup>n</sup>* /Δ*<sup>t</sup>* and . <sup>ψ</sup> = (ψ*n*+<sup>1</sup> <sup>−</sup> <sup>ψ</sup>*n*)/Δ*t*. Substituting these discretizations into Equation (18), we obtain the following system:

$$\begin{cases} A\rho^{n+1} + 0.5\Delta t (B - \mathbb{C}(z\_k))\psi^{n+1} = A\rho^n - 0.5\Delta t (B - \mathbb{C}(z\_k))\psi^n, \\ A\psi^{n+1} - 0.5\Delta t (B - \mathbb{C}(z\_k))\rho^{n+1} = A\psi^n + 0.5\Delta t (B - \mathbb{C}(z\_k))\rho^n. \end{cases} \tag{21}$$

Using the boundary conditions in Equation (21) yields a septadiagonal matrix. This system is solved by an appropriate method in MATLAB.

First, the initial vectors of parameters <sup>ρ</sup><sup>0</sup> = ρ0 0, ... , ρ<sup>0</sup> *N* and <sup>ψ</sup><sup>0</sup> = ψ0 0, ... ,ψ<sup>0</sup> *N* are calculated to solve System (21) with the use of the initial condition using the relations

$$s\_N(\mathbf{x}, 0) = \sum\_{j=-1}^{N+1} \mathcal{Q}\_j(\mathbf{x}) \rho\_j^o \text{ and } r\_N(\mathbf{x}, 0) = \sum\_{j=-1}^{N+1} \mathcal{Q}\_j(\mathbf{x}) \psi\_{j'}^o \tag{22}$$

where all parameters ρ<sup>0</sup> and ψ<sup>0</sup> are determined. *sN* and *rN* require the following relations to be satisfied at points *xk*:

$$\begin{aligned} \mathbf{s}\_N(\mathbf{x}\_k, \mathbf{0}) &= \mathbf{s}(\mathbf{x}\_k, \mathbf{0}), \\ r\_N(\mathbf{x}\_k, \mathbf{0}) &= r(\mathbf{x}\_k, \mathbf{0}), k = \mathbf{0}, 1, \dots, N. \end{aligned} \tag{23}$$

The initial vectors ρ<sup>0</sup> and ψ<sup>0</sup> can be calculated using the initial and boundary conditions from the following matrix equations:


and

This system is solved using a suitable method in MATLAB. The approximate numerical solution of *sN*(*x*, *t*) and *rN*(*x*, *t*) is obtained from the ρ*<sup>n</sup>* and ψ*<sup>n</sup>* using Equation (21).

#### **3. Stability Analysis**

Stability analysis of the numerical scheme is carried out using the von Neumann method. The coupled system of Equation (5) is written as

$$\frac{\partial \mathcal{U}}{\partial t} + M \frac{\partial^2 \mathcal{U}}{\partial x^2} + G(\mathcal{U}) \mathcal{U} = 0,\tag{24}$$

where *U* = *s r* ! , *M* = <sup>0</sup> <sup>−</sup><sup>1</sup> 1 0 ! and *G*(*U*) = <sup>0</sup> <sup>−</sup>μ<sup>1</sup> μ<sup>2</sup> 0 ! .

To perform linear stability analysis, we further linearize the nonlinear term in Equation (24) by taking *G*(*U*) = μ*M*. Here μ = max[μ1, μ2]. Therefore, the linearized form of Equation (24) is as follows:

$$
\frac{
\partial \mathcal{U}
}{
\partial t
} + M \frac{
\partial^2 \mathcal{U}
}{
\partial \mathbf{x}^2
} + \mu M \mathcal{U} = 0.
\tag{25}
$$

The discretization of linear Equation (24) by the proposed scheme is given by

$$A\mathcal{U}l\_{m}^{n+1} - M\mathcal{D}l\mathcal{U}\_{m}^{n+1} = A\mathcal{U}l\_{m}^{n} + M\mathcal{D}l\mathcal{U}\_{m}^{n+1} \,. \tag{26}$$

where *D* = *B* − μ*I*.

We use the following Fourier modes analysis for Scheme (25)

$$
\mathcal{U}^{n}\_{m} = P \mathcal{G}^{n} e^{im\upsilon}.\tag{27}
$$

Here *<sup>i</sup>* <sup>=</sup> <sup>√</sup> <sup>−</sup>1, *<sup>P</sup>* <sup>∈</sup> *<sup>R</sup>*<sup>2</sup> and *<sup>G</sup>* <sup>∈</sup> *<sup>R</sup>*2×<sup>2</sup> is the amplification matrix. Equation (26) is substituted into Equation (25) to obtain the value of the amplification *G*. After tedious algebra, we obtain the value of *G* as follows

$$G = \left[\frac{A + MD}{A - MD}\right] \tag{28}$$

Matrix *M* is a skew symmetric matrix, and therefore both the matrices *A* + *MD* and *A* − *MD* have the same eigen values. Thus, the maximum value of |*G*| = 1. Hence, the linearized Scheme (25) is unconditionally stable.

#### **4. Numerical Results and Test Problems**

Three test problems, including the single soliton, interaction of two solitons and birth of standing soliton, are presented to evaluate the effectiveness and performance of the proposed method. The accuracy of the proposed method is examined using the error norms *L*2and*L*<sup>∞</sup> and conservation laws defined as

$$\begin{aligned} L\_2 &= \|\mu^{\text{exact}} - \mu\_N\|\_2 = \sqrt{h \sum\_{k=0}^N \left| \mu\_k^{\text{exact}} - (\mu\_N)\_k \right|^2}, \\ L\_{\text{eq}} &= \|\mu^{\text{exact}} - \mu\_N\|\_{\infty} = \max\_{0 \le k \le N} \left| \mu\_k^{\text{exact}} - (\mu\_N)\_k \right|. \end{aligned} \tag{29}$$

Moreover, Equation (1) must satisfy the two conservation laws

$$\begin{aligned} I\_1 &= \underset{a}{\overset{b}{\underset{a}{\left\|}}} |u|^2 dx}, \\ I\_2 &= \underset{a}{\overset{b}{\left\|}} \left| |u\_x|^2 - \frac{1}{2} q |u|^4 \right| \right. \end{aligned} \tag{30}$$

#### *4.1. Single Soliton Solution to the NLS Equation*

The analytical single soliton solution to the NLS equation is given as [2]:

$$u(\mathbf{x},t) = \beta \left(\sqrt{\frac{2}{q}}\right) \text{exp}\{\frac{1}{2}\mathbf{S}\mathbf{x} - \frac{1}{4}(\mathbf{S}^2 - \beta^2)t\} \text{sech}(\beta\mathbf{x} - \beta\mathbf{S}t),\tag{31}$$

where *S* is the speed of the soliton solution whose magnitude depends on the parameter β.

The numerical solution is computed using the following parameters: *q* = 2, *S* = 4, β = 1, *xL* = −20, and *xR* = 20. The conservative quantities *I*<sup>1</sup> and *I*<sup>2</sup> and error norms *L*<sup>∞</sup> and *L*<sup>2</sup> are calculated. Numerical simulations were carried out at different space steps and time steps for comparison with the published results of the previous methods. The results are documented in Tables 1 and 2. The numerical simulations and the absolute numerical error are shown at different times in Figure 1.

Table 2 displays a comparison between the results obtained by the proposed method and published results. The results of the proposed method are in agreement with analytical solutions within very satisfactory limits, and the proposed method exhibits the same results as the quintic B-spline collocation method and quadratic B-spline method, as can be seen from the comparison in Table 2. We found good results even for large step sizes. It is noted that the error norms of the Galerkin cubic B-spline are lower than Explicit method [3], Implicit/explicit [3], and Split-step Fourier method [3] when we compare the results of the proposed method with those referenced in [1,4,19].


**Table 1.** Error norms and conservation laws for single soliton: with *h* = 0.04, *q* = 2, *S* = 4, β = 1.

**Figure 1.** Single soliton solution and errors |Numerical – Analytical|.


**Table 2.** Comparisons of the present results with those of Taha et al. [20]., amplitude = 1 at time = 1.

The rate of convergence for spatial and temporal directions is calculated using the formula [21]:

$$\text{rate of convergence} \approx \log\_2 \frac{error(2h, 2\Delta t)}{error(h, \Delta t)},$$

where *error*(2*h*, 2Δ*t*) is the error norms *L*<sup>∞</sup> and *L*<sup>2</sup> in spatial and temporal directions. The error norms *L*∞, *L*<sup>2</sup> and order of convergence rate at time *t* = 1 are shown in Table 3.

**Table 3.** Rate of convergence in spatial and temporal directions at *t* = 1 with *h* = Δ*t*.


In Table 4, the results are documented for when the soliton amplitude is equal to 2. It was found that the results of the error norms at different times increase very slightly, as shown in Table 4. One soliton solution with amplitude = 2 is calculated, and simulations are shown in Figure 2.

**Table 4.** Error norms and conservation laws for single soliton with *h* = 0.05, Δ*t* = 0.0025 *h* = 0.05, *q* = 2,*S* = 4, β = 2.


Table 5 displays the approximate results obtained using the proposed method compared with the published methods listed in the reference for single soliton when the amplitude is 2. These results are in good agreement with the analytical solution.

**Figure 2.** Single soliton solution with amplitude 2.

**Table 5.** Comparisons of single soliton with those of Taha et al. [20], amplitude = 2 at time = 1.


#### *4.2. The Interaction of Two Solitons for the NLS Equation*

In the second problem, we discuss the interaction of two solitons moving in opposite directions, have the same amplitude of magnitude 1 with initial conditions as follows [16,17,20]

$$\begin{aligned} u(\mathbf{x},0) &= \sum\_{k=1}^{2} u\_k(\mathbf{x},0), \\ u\_k(\mathbf{x},0) &= \beta\_k \Big(\sqrt{\frac{2}{q}}\Big) \exp[\frac{1}{2}S(\mathbf{x}-\mathbf{x}\_k)] \text{sech}(\beta\_k(\mathbf{x}-\mathbf{x}\_k)). \end{aligned} \tag{32}$$

where β*k*, *q* and *xk* are constants.

 ʹͲ The numerical solution is computed by using the following parameters *x*<sup>1</sup> = 10, *x*<sup>2</sup> = −10, *q* = 2, β<sup>1</sup> = β<sup>2</sup> = 1, *S*<sup>1</sup> = −4 and S2 = 4. The first soliton is placed at *x* = 10 and is moving to the left with speed 4, and the second soliton is placed on the other side at *x* = −10, traveling to the right with

speed 4. Both solitons are moving in opposite directions, and they collide and separate. The interactions of these two solitons are visualized in Figure 3. The error norms *L*<sup>∞</sup> and *L*<sup>2</sup> are computed at different times and at different space steps and time steps. The calculated results are tabulated in Table 6.

**Figure 3.** Interaction of two solitons, amplitude = 1 at different space and time steps.

The error norms *L*<sup>∞</sup> and *L*<sup>2</sup> are very small compared with those of Taha et al. [20], as shown in Table 7. It can be clearly seen that the error norms *L*<sup>∞</sup> and *L*<sup>2</sup> obtained by the present method are smaller than those of previous methods [20], as shown in Table 7.


**Table 6.** Two-solitons, amplitude = 1 at different space and time steps.

**Table 7.** Comparisons of two solitons with Taha et al. [20], amplitude = 1.


*4.3. Birth of Standing Soliton with the Maxwellian Initial Condition*

If

$$\int\_{-\infty}^{\infty} u(x,0) \ge \pi. \tag{33}$$

a soliton should appear over time with initial values greater than π**,** otherwise the soliton will decay away [22].

Consider the birth of soliton with the Maxwellian initial condition given by [16]:

$$
\mu(\mathbf{x},0) = A \exp(-\mathbf{x}^2). \tag{34}
$$

The values of all parameters are chosen to be *h* = 0.08, Δ*t* = 0.004 and *q* = 2 for the domain [−45, 45] to exhibit the birth of a soliton. The numerical simulations are shown at different times for the values of *A* = 1.78 over the domain [−45,45] in Figure 4. With *A* = 1, the soliton decays away as expected.

The conserved quantities *I*<sup>1</sup> and *I*<sup>2</sup> are computed using the Maxwellian initial condition (34). Analytical conserved quantities can be computed as:

$$\begin{aligned} I\_1 &= A^2 \sqrt{\frac{\pi}{2}} = 3.97100051267043, \\ I\_2 &= \frac{1}{4} A^2 \Big(2\sqrt{2} - qA^2\Big) \sqrt{\pi} = -4.92561762132093, \end{aligned}$$

The conserved quantities *I*<sup>1</sup> and *I*<sup>2</sup> are tabulated in Table 8. The numerical results obtained are compared with the published results of Mokhtari et al. [11] and the exact solution, as shown in Table 8. The proposed method conserves *I*<sup>1</sup> to 7 decimal places and *I*<sup>2</sup> to 3 decimal places.

**Figure 4.** Birth of a standing soliton for *A* = 1.78.

**Table 8.** Comparison of conserved quantities of the birth of a standing soliton with Mokhtari et al. [11], *h* = 0.08, Δ*t* = 0.004 and *q* = 2.


#### **5. Conclusions**

The approximate solution of the NLS equation was investigated using the Galerkin finite element method with a cubic B-spline shape function. Three numerical problems, including single soliton, interaction of two solitons, and birth of a standing soliton with the Maxwellian initial condition, were demonstrated to evaluate to the performance and accuracy of the method. Furthermore, we simulated the numerical solution by choosing different parameters for motion of the single soliton, the interaction of two solitons and the birth of the standing soliton. The error norms *L*<sup>∞</sup> and *L*<sup>2</sup> and conservation laws *I*<sup>1</sup> and *I*<sup>2</sup> were determined and compared with published results [11,16,17,20]. It was seen that all of our results for the problems were computed to be reasonably low, and were found to be in good agreement with the analytical solution. The present method was shown to be unconditionally stable. The method has almost a second-order convergence. In conclusion, the proposed scheme with cubic B-spline presents an acceptable result for the NLS equation.

**Author Contributions:** Formal analysis, A.I. and N.N.A.H.; Investigation, A.I., N.N.A.H. and A.I.M.I.; Methodology, A.I.; Supervision, N.N.A.H. and Ahmad A.I.M.I.; Writing – original draft, A.I.; Writing – review & editing, N.N.A.H. and A.I.M.I.

**Funding:** This work was supported by USM Short Term Grant of account number 304.PMATHS.6315077.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


22. Delfour, M.; Fortin, M.; Payne, G. Finite-difference solutions of non-linear Schrödinger equation. *J. Comput. Phys.* **1981**, *44*, 277–288. [CrossRef]

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
