A mixed finite element method for solving coupled wave equation of Kirchhoff type with nonlinear boundary damping and memory term

This paper is concerned with the numerical approximation of the solution of the coupled wave equation of Kirchhoff type with nonlinear boundary damping and memory term using a mixed finite element method. The Raviart‐Thomas mixed finite element method is one of the most prominent techniques to discretize the second‐order wave equations; therefore, we apply this scheme for space discretization. Furthermore, an L2‐in‐space error estimate is presented for this mixed finite element approximation. Finally, the efficiency of the method is verified by a numerical example.


INTRODUCTION
Mixed finite element methods are among the most efficient techniques to solve second-order differential equations; therefore, this method has become an interesting research area during recent years. In particular, in Pani and Yuan, 1 a mixed finite element Galerkin method is used for a strongly damped wave equation. Liu et al. 2 presented an H 1 -Galerkin mixed finite element method for a class of second-order Schrödinger equation. Liu et al. 3 considered a new numerical scheme based on the H 1 -Galerkin mixed finite element method for a class of second-order pseudo-hyperbolic equations. Liu et al. 4 presented two splitting mixed finite element schemes for the pseudo-hyperbolic equation where a mixed finite element method applied for approximating the solution of nearly incompressible elasticity and Stokes equations.
Wang and Ye 5 introduced a novel mixed finite element for the second-order elliptic equation formulated as a system of two first-order linear equations. In Monk et al., 6 a hybridized mixed method was presented to discretize the Helmholtz problem using Raviart-Thomas finite elements. Discretization of the shallow-water equations using the Raviart-Thomas finite element method was presented in Rostand and Le Roux. 7 Ainsworth 8 obtained an a posteriori error estimator for the lowest order Raviart-Thomas mixed finite element which provided the computable upper bounds on the error in the flux variable regardless of jumps in the material coefficients across interfaces. The numerical approximation of the displacement form of the acoustic wave equation using mixed finite elements was explained in Jenkins. 9 In Glowinski and Lapin, 10 the space-time discretization using a combination of mixed finite element method (Raviart-Thomas) and a finite difference scheme was applied for the time discretization of the wave equation. A splitting positive definite mixed finite element method was used for the second-order viscous elasticity wave equation in Liu et al. 11 In Parvizi and Eslahchi, 12 an extended Raviart-Thomas mixed finite element method was applied to approximate the solution of damped Boussinesq equation. In Gao et al., 13 the numerical solution of a viscoelasticity wave equation using mixed finite element approximations is studied. Jenkins et al. 14 presented a priori error estimation for the mixed finite element displacement formulations of the acoustic wave equation. In He et al., 15 mixed finite element methods, explicit and implicit in time, for a fourth-order wave equation were studied. Liu and Li 16 applied Galerkin mixed finite element methods to approximate the solution of a class of second-order pseudo-hyperbolic equations. In Bécache et al., 17 a new family of quadrangular (two-dimensional) or cubic (three-dimensional) mixed finite elements for the approximation of elastic wave equation was used. A new approach for stabilization of low-order velocity-pressure pairs for the incompressible Stokes equations using mixed finite element was presented in Bochev et al. 18 For several model problems, the general inf-sup condition for the mixed finite element was reviewed in Bathe. 19 In Shi et al., 20 a new stabilized mixed finite element method for the Poisson equation was presented. For more details, we refer to the following paper and corresponding equations, Cahn-Hilliard equation, 21,22 phase-field fracture, 23 drift-diffusion-Poisson system, 24 and damped Boussinesq equation. 25 A crucial step in the analysis of the mixed finite element schemes is to study the error analysis of the discretized solutions obtained from the mixed finite element method. But first, we need to verify conditions such as consistency, ellipticity, and the so-called inf-sup (or LBB) condition. 26 A main advantage of using the mixed finite element method is the freedom to choose suitable approximation spaces for different types of equations. Robey 27 shows that in the mixed framework, it is not always necessary to meet the LBB condition, and instead, the idea of a weak LBB condition is introduced in this paper.
This paper is concerned with the Raviart-Thomas mixed finite element method for the following coupled Kirchhoff-type wave equation with nonlinear boundary damping and memory term 28,29 : ( ( where U and V represent the transverse displacements, Ω is a polygon domain of R 2 with a boundary Γ ∶= Ω such that Γ = Γ 0 ∪ Γ 1 and Γ 0 and Γ 1 have positive measures, and Here, (1a) has its origin in the mathematical description of small amplitude vibrations of an elastic string. 29 We define U ∶= ∇U · where is the unit outer normal vector pointing toward Ω and for some m 0 , m 1 , m 2 > 0, m 1 > 2( + 2), and 1 − ∫ ∞ 0 g(r)dr. The function f satisfies the mentioned conditions of Ferreira et al., 30 which guaranties its meaningfulness. Finally, the existence and the uniqueness of the system Equation 1 can be found in Ferreira et al. 30 A fully discretization analysis is done by Parvizi et al. 31 Let us define W ∶= {u ∈ H 1 (Ω)|u = 0onΓ 1 } and assume U 0 , U 1 , V 0 , and V 1 belong to H 3/2 (Ω) ∩ W should satisfy the following assumptions: Now, we review the history and the background of (1). The first Kirchhoff equation was of the form where this equation extends the classical d'Alembert's wave equation by considering the effects of the changes in the length of the strings during the vibrations. Equation 3 is called the wave equation of the Kirchhoff type because Kirchhoff was the first one who introduced this equation in the study of oscillation of stretched strings and plates. 32 The equation is a mathematical Kirchhoff model of nonlinear transverse vibration, neglecting the displacements along the string's axis and averaging tension N over its length L. In this equation, is the density of the string material, E is the Young modulus of the latter, C is the viscous damping parameter, 0 is the initial string tension value, and q(x, t) is the transverse load intensity. 33 Here, we consider the following equation 34 : In the special case of (4), the dynamics of the moving string in Figure 1 can be described by where U = U(x, t) is the lateral displacement at the space coordinate x and the time t. Additionally, E, , h, 0 , a, and f are the Young modulus of the latter, the mass density, the cross-section of the area, the length, the initial axial tension, and the For a freely vibrating fixed-fixed string with two polarized displacements U 1 and U 2 , the Kirchhoff-Carrier equations are described by 35,36 where N is the axial tension created by the large amplitude motions and the coupling with the transverse motion and defined by In these equations, is density, E is Young's modulus, h is the crosse section, 0 is the tension, and L is the length. Equation 1 is a nonlinear PDE with a gradient term. In order to obtain the convergence rate for the approximation technique, the efficient strategy is converting the equation into a system of two nonlinear first-order equations (using an auxiliary variable). Then, the Raviart-Thomas mixed finite element can be employed to solve the system. The paper is organized as follows. In Section 2, the necessary definitions and lemma that will be used in the next sections are introduced. In Section 3, the semi-discrete Raviart-Thomas mixed finite element method 37 for solving (1) is presented in detail, and the main theorem of this paper is mentioned. In Section 4, we provide a proof for the main result mentioned in Section 3. In Section 5, we demonstrate the accuracy of theoretical results using a numerical example.

PRELIMINARY AND NOTATIONS
This section is devoted to the necessary definitions and some preliminaries. We also mention a lemma regarding the regularity of the solutions of Equation 1 and introduce two required lemmas that will be used later in our analysis.

Lemma 1. Under the mentioned assumptions for Equation 1, these equations have unique solutions
then, satisfies the following estimate

Lemma 2.3. (Houston et al 39 and Liu & Barrett 40 ). Let ∈ C([0, T] ×Ω) satisfies the following estimate:
where m f and M f are positive constants. Then, there exist the constants

MIXED FINITE ELEMENT METHOD (SPATIAL DISCRETIZATION)
In this section, we present a mixed variational formulation for Equation 1 and describe the Raviart-Thomas mixed finite element scheme associated with this variational formulation. Furthermore, we mention the main theorem of the paper which is related to an a priori estimate of the error of the semi-discrete scheme.

The weak formulation
This subsection is devoted to introduce the variational formulation of Equation 1. We denote X ∶= L 2 (Ω) and M ∶= H(div, Ω); then, the mixed formulation is obtained by splitting each of (1a) and (1b) into two equations where W ∶= ∇U ∈ M and Z ∶= ∇V ∈ M. Thus, Equation 1 can be written in the following form: Then, the weak formulation associated with (11) gives rise to the following variational formulation of (1): find

Mixed finite element discretization
. Also, the mesh  h assume to be regular in the Ciarlet sense. Furthermore, we assume the elements are -shape regular in the sense that Then, we define the following spaces: Finally, based on the variational formulation (12) and the mentioned finite element spaces, we define the Raviart-Thomas mixed finite element method for (1): Considering suitable regularity assumptions on the exact solutions, the following theorem provides the theoretical rate of convergence for the semi-discretized Raviart-Thomas method (13).
be the unique solutions of the continuous formulation (12) and u h , v h ∈ X h and w h , z h ∈ M h be the unique solutions of the discrete formulation (13). Then, the following error estimation holds: where c is a constant depending only on Ω and  h .
The mentioned theorem will be proved in Section 4. Before we start to prove the theorem, some necessary lemmas should be presented and proved.

PROOF OF THE MAIN RESULTS
The main purpose of this section is to prove Theorem 3.1, which is mentioned in Section 3. Hence, for the sake of simplicity and preventing the complexity in the proof, we separate this section into three subsections and prove the required lemma and theorems. Then, we employ these results to prove the main theorem.
where c is a constant depending only on Ω and  h .
Let us define the error terms of U, V, W, and Z as Subtracting (12a-12d) from (13a-13d) gives us ( where From now on, we consider the following regularity assumptions.

Assumptions 4.1. We assume U(t), V(t) ∈ H l (Ω) and W(t), Z(t)
∈ (H l (Ω)) 2 , 1 ≤ l ≤ k + 1, are the exact solutions of (11) and u h (t), v h (t) ∈ X h and w h , z h ∈ M h are the discretized solutions of (11) obtained from (13). Additionally, we assume U 0 , In order to prove Theorem 3.1, the challenging part is to deal with L 1 and L 2 . Therefore, the following subsections are devoted to obtain upper bounds for L 1 and L 2 .

An upper bound for
Proof. Taking derivative with respect to t from (21) and (22), and putting = e 3,h and = e 4,h into these equations as test functions, respectively, yield d dt Considering the orthogonality estimate from Equation 14, we have Then, the combination of (26)- (29) gives us the desired result.
then, we have where c i , i = 1, ·, 4 are constants depending only on Ω and  h .
Proof. First, we should note that L 1, 0 can be written in the following form: Then, rearranging the terms of (32) gives rise to Employing Lemma 4.2 and using the following estimates result in the following inequality: In order to simplify (34), we rewrite it as where Finally, applying the Cauchy-Schwarz and Young inequalities and approximate properties of Π h and R h from Lemma 4.1 yields which gives us the desired results.
then, under Assumption 4.1, L 1, 2 can be written in the following form: where L 1, 2, 0 and L 1, 2, 1 satisfy the following properties: and Proof. Using the definition of L 1, 2 as well as adding and subtracting some terms, one arrives at Substituting (26) and (27) from Lemma 4.2 into the above formula gives us where and Employing the Cauchy-Schwarz and Young inequalities and approximate property of R h from Lemma 4.1, we conclude and The combination of (41), (44), and (45) completes the proof.

Theorem 4.1. Let W(t), Z(t)
∈ (H l (Ω)) 2 , 1 ≤ l ≤ k + 1; then, L 1 satisfies the following inequality: and L 1, 2, 0 and L 1, 2, 1 are defined in (42) and (43), respectively, and M 0 can be found in (36). Furthermore, L 1, 1 and L 1, 3 satisfy the following estimates: and Proof. Rearranging L 1 gives us where L 1, 0 and L 1, 2 are defined in Lemmas 4.3 and 4.4, respectively, and we denote L 1, 4 by Considering (48) and Lemmas 4.2-4.4, L 1 satisfies the following estimate: Adding and subtracting some terms allow us to rewrite L 1, 1 as Applying the Cauchy-Schwarz and Young inequalities and Equation 18, we conclude the following estimate for L 1, 1 : Adding and subtracting some terms allow us to rewrite L 1, 3 in the following form: Using the Cauchy-Schwarz and Young inequalities and the approximate property of R h from Lemma 4.1, the following inequality can be deduced: Finally, the combination of (49), (50), and (51) completes the proof.

An upper bound for
where From (55) and (56), we get the following estimate for L 2 : Finally, the Trace theorem and Cauchy-Schwarz inequality and Equation 15 imply that (58)

Proof of Theorem 3.1
Proof. Let us denote L 3, 1 and L 3, 2 by and Using the Cauchy-Schwarz and Young inequalities, Trace theorem, and Lemma 2.3, we get and Now, from Lemma 2, we get  .

=1
and  = { } N 2 =1 are basis for X h and M h , respectively; therefore, from (13), we have ( Using the above equations and applying (15)-(18), we find that In this section, we validate the theoretical results obtained in Theorem 3.1 by a test problem. The computational geometry (Ω = [0, 1] × [0, 1]) including the boundaries is shown in Figure 2. In this section, we use the following spaces: The initial conditions for the model problem (1) are given as ) .
In this example, the boundary Γ consists of Γ 0 and Γ 1 (see Figure 1) such that the homogenous Neumann and Dirichlet boundary conditions are defined on Γ 1 . For the boundary conditions (1d) and (1e) defined on Γ 0 , we set g(t) ∶= 0.1 exp(−t), p ∶= 3.5, ∶= 2.5. Finally, the exact solutions of Equation 1 considering the mentioned boundary and initial conditions are given by ) . In order to discretize in the space, we use the mixed variational formulation given in (13). The time discretization is done by an implicit finite difference (Crank-Nicolson method) scheme. Furthermore, to treat the nonlinearity, a Newton scheme is used. The numerical results for two different time step sizes (Δt = 0.005 and Δt = 0.003) are presented in Figure 3. Moreover, for the final time step, the discretization errors (U(T) − u h (T) and V(T) − v h ) with respect to the exact solutions for two specific mesh sizes are illustrated in Figure 4.