# Analysis and numerical simulation of the nonlinear beam equation with moving ends

The numerical analysis for the small amplitude motion of an elastic beam with internal damping is investigated in domain with moving ends. An efficient numerical method is constructed to solve this moving boundary problem. The stability and convergence of the method is studied, and the errors of both the semi-discrete and fully-discrete schemes are derived, using Hermite's polynomials as a base function, proving that the method has order of quadratic convergence in space and time. Numerical simulations using the finite element method associated with the finite difference method (Newmark's method) are employed, for one-dimensional and two-dimensional cases. To validate the theoretical results, tables are shown comparing approximate and exact solutions. In addition, numerically the uniform decay rate for energy and the order of convergence of the approximate solution are also shown.

## Authors

• 1 publication
• 1 publication
05/05/2020

### Convergence analysis of the scaled boundary finite element method for the Laplace equation

The scaled boundary finite element method (SBFEM) is a relatively recent...
09/03/2020

### A variational analysis for the moving finite element method for gradient flows

By using the Onsager principle as an approximation tool, we give a novel...
03/11/2021

### Convergence of Computed Dynamic Models with Unbounded Shock

This paper studies the asymptotic convergence of computed dynamic models...
02/04/2022

### Convergence Analysis of Virtual Element Method for Nonlinear Nonlocal Dynamic Plate Equation

In this article, we have considered a nonlinear nonlocal time dependent ...
10/15/2020

### Crack tip fields and fracture resistance parameters based on strain gradient plasticity

The crack tip mechanics of strain gradient plasticity solids is investig...
07/15/2019

### Numerical study of vanishing and spreading dynamics of chemotaxis systems with logistic source and a free boundary

The current paper is to investigate the numerical approximation of logis...
04/15/2022

### Convergence of the Discrete Minimum Energy Path

The minimum energy path (MEP) describes the mechanism of reaction, and t...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The equation of motion of a thin beam with weak-internal damping undergoing cylindrical bending can be written as

 u′′(x,t)−(ζ0+ζ1∫Ωt∣∣∇u(x,t)∣∣2dx)Δu(x,t)+Δ2u(x,t)+νu′(x,t)=0, in Qt, (1)

where is the transverse displacement,

is the vector of spatial coordinates,

is the time. The aerodynamic damping term is denoted by is the nonlinear stiffness of beam, is an in-plane tensile load, and belongs to the noncylindrical domain, defined by

 Qt={(x,t)∈Rn+1;x=K(t)y, y∈Ω, 0

with is an open, limited and regular and is a - real function. All quantities are physically nondimensionalized, are fixed positive and is not necessarily positive.

We will study the equation (1) with zero boundary conditions , as follows,

 u(x,t)=0and∇u≡0, on ∂Ωt×[0,T[, (3)

where .

In the absence of the damping term , many authors have already studied Cauchy’s problems and the mixed problems associated with the equation (1) on bounded and unlimited cylinders. The existence and uniqueness of solution was demonstrated in [9].

The Cauchy problem associated with the (1) equation in abstract terms on a Hilbert space is studied, among other authors, by [4, 18]. In addition, in the last were established results on existence, unicity and asymptotic stability of the solution.

For the one-dimensional case, omitting the term and considering the cylinder , with open and limited, we obtain the Kirchhoff equation with internal damping, namely,

 u′′(x,t)−(ζ0+ζ1∫Ω∣∣ux(t)∣∣2dx)uxx(x,t)+νu′(x,t)=0, in Q. (4)

The equation (4) has been extensively studied by a wide variety of authors on -dimensional cases and generic mathematical models defined in Hilbert spaces. The existence of a local and global solution can be found in several physical-mathematical contexts, such as [14, 15].

In the article [15] a new model associated with the Kirchhoff equation (4) was introduced, called the Medeiros-Kirchhoff equation, given by

 u′′(x,t)−(a(t)+b(t)∫Ω∣∣ux(t)∣∣2dx)uxx(x,t)=0, in Qt.

There are several areas that work with the application of the Beam equation, for example, in civil and naval engineering, as well as in the aerospace area as the reaction of rocket and satellite metal structures to various space situations, as can be observed in the articles [12, 13]. In addition to these areas, it is possible to associate Maxwell’s equations, which govern electromagnetic theory, and the Euler-Bernoulli Beam Equation, in the study presented by article [16]

The equation (1) follows the Euler-Bernoulli Beam model. In addition to this model, there are other models, extensions of the Euler-Bernoulli model, which consider rotational movements besides transverse ones, being possible to create models closer to reality. The article [11] presents a more detailed study on three of these models.

The existence and uniqueness of the solution, the asymptotic decay and the numerical simulations of the beam equation with movinf boundary were studied in [Raquel] for the one-dimensional case.

In this paper, we present a family of numerical method for the one-dimensional and two-dimensional cases, based on the finite element method and the finite difference method. In addition, we are doing numerical analysis for both the semi-discrete problem and the fully discrete problem, showing that the order of convergence is quadratic in space and time. For this we will use the Faedo-Galerkin method, using the Hermite polynomials as base functions. In addition, the nonlinear system associated with the system of ordinary differential equations is being solved by Newton’s method and for temporal discretization is applying the Neumark’s method, where we will show the convergence of the approximate solution to a family of numerical methods depending on the choice of

. Numerical simulations are presented as examples, with different types of boundary and tables showing the order of numerical convergence. In addition, we present the asymptotic decay, and the graphs of the solutions for the one-dimensional and two-dimensional cases.

## 2 Analytical results

Our section will present, without demonstrating, the existence and uniqueness results of the mixed problem solution for the equation (1), considering the boundary conditions (3), given by:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩u′′(x,t)−(ζ0+ζ1∫Ωt∣∣∇u(x,t)∣∣2dx)Δu(x,t)+Δ2u(x,t)+νu′(x,t)=0,%in Qt u(x,0)=u0(x),  u′(x,0)=u1(x), in Ω0, u(x,t)=0 and ∇u≡0, in ∂Ωt×[0,T[, (5)

where is defined in 2 and satisfying the hypotheses

 {Qt⊂Ωt×]0,T[, with Ωt⊂Rn open and limited with 0∈Ωt,K∈C2([0,+∞[), % with 0

where and are positive real constants.

Let us consider the diffeomorphism , defined by , with . Therefore, the change of variable transforms the problem(5) into the equivalent problem with cylindrical boundary

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩v′′(y,t)−b1(t)∣∣∇v∣∣20Δv(y,t)+b2(t)Δ2v(y,t)+ν v′(y,t)+a(2)ij(y,t)∇yi,yjv(y,t)−a(1)i(y,t)Δyiv(y,t)−a(3)i(y,t)∇yiv(y,t)−a(4)i(y,t)∇yiv′(y,t)=0, in% Q v(y,0)=v0(y),  v′(y,0)=v1(y)in  Ω, v(y,t)=0 and ∇v≡0, in ∂Ω×[0,T[, (6)

where, for simplicity, we denote and

 b1(t)=ζ1K4, b2(t)=1K4, a(1)i(y,t)=1K2[ζ0−4(yiK′)2], a(2)ij(y,t)=4yiyj(K′/K)2, (7) a(3)i(y,t)=1K2[2yi(K′)2−yiK(νK′+K′′)], a(4)i(y,t)=−2yi(K′/K).

In the following, we present the results on the existence and uniqueness of the weak solution of problems (5) and (6) and asymptotic decay, whose demonstrations can be found in [8] for the one-dimensional case and [5] for the two-dimensional case.

###### Theorem 1.

(Existence and Uniqueness of solution in nocylindrical domain.)

Let , and the hypothesis (H1). Then there is a unique weak solution for the problem (5) satisfying the following conditions:

 1.u∈L2(0,T;H20(Ωt)), 2.u′∈L2(0,T;L2(Ωt))

and the weak formulation of (5) is verified in the sense of .

###### Theorem 2.

(Existence and Uniqueness of solution in cylindrical domain)

Let , and the hypothesis (H1). Then , Then there is a unique weak solution for the problem (5) satisfying the following conditions:

 1.v∈L2(0,T;H20(Ω)), 2.v′∈L2(0,T;L2(Ω))

and the weak formulation of (6) is verified in the sense of .

###### Theorem 3.

(Asymptotic Decay: Two-dimensional case)

Consider the hypothesis H1, then we have the energy (5), given by

satisfie

 E(t)≤A0 e−A1t, ∀t≥0,

where and are positive constants.

From the important theoretical results, we can now develop numerical methods to determine an approximate numerical solution and verify the consistency of the numerical methods used.

Due to the equivalence of the problems, the methods used will be for the problem (6) that is defined in a cylindrical domain.

## 3 Numerical Method

The purpose of this section is to develop a numerical method, based on the finite element method and the finite difference method to make the numerical simulations in order to obtain a approximate solution of the problem (6).

Denoting the Hilbert space , then weak formulation of the problem (6) can be written in the follows form, when we integrate by parts the second, third, fifth and sixth terms:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(v′′(t),w)+b1(t)∣∣∇v(t)∣∣20(∇v(t),∇w)+b2(t)(Δv(t),Δw)+ν (v′(t),w)+(a(1)i(t)∇yiv(t),∇yiw)+(a(2)ij(t)∇yi,yjv(t),w)+(a(4)i(t)∇yiv′(t),w)+(a(5)i(t)∇yiv(t),w)=0, in Q(v(y,0),w)=(v0(y),w),  (v′(y,0),w)=(v1(y),w), in  Ω,(v,w)=0 and (∇v,w)≡0, in ∂Ω×[0,T[, ∀w∈V, (8)

where .

### 3.1 Approximated Problem

Let be a subspace generated by the first vectors of the basis of Hilbert space . Then we want to determine an approximate solution of the following approximate system, given by:

 (9) +(a(1)i(t)∇yivh(t),∇yiwh)−(a(2)ij(t)∇yivh(t),∇yjwh) +(a(4)i(t)∇yiv′h(t),wh)+(a(5)i(t)∇yivh(t),wh)=0, ∀wh∈Vh,

where we integrate by parts the sixth term.

Since that, , then the function can be represented by

 vh(y,t)=m∑i=1di(t)φi(y),φi(y)∈Vh. (10)

Taking , for some , and substituting (10) in (9), we obtain

 m∑k=1{d′′k(t)(φk(y),φl(y))+b1(t)dk(t)∣∣m∑k=1dk(t)∇φk(y)∣∣20(∇φk(y),∇φl(y)) +b2(t)dk(t)(Δφk(y),Δφl(y))+d′k(t)[ν (φk(y),φl(y))+(a(4)i(t)∇yiφk(y),φl(y))] (11) +dk(t)[(a(1)i(t)∇yiφk,∇yiφl)−(a(2)ij(t)∇yiφk,∇yjφl)+(a(5)i(t)∇yiφk,φl)]}=0.

In order to validate the method (11), we add the term source on the right side of the equation (9) and, taking , we have the following definitions for the matrices and vectors:

 Akl=(φk(y),φl(y)),K1(kl)=(∇φk(y),∇φl(y)),K2(kl)=(Δφk(y),Δφl(y)), (12) B1(ikl)(t)=(a(1)i(t)∇yiφk(y),∇yiφl(y)),B2(ijkl)(t)=(a(2)ij(t)∇yiφk(y),∇yjφl(y)), B3(ikl)(t)=(a(4)i(t)∇yiφk(y),φl(y)),B4(ikl)(t)=(a(5)i(t)∇yiφk(y),φl(y)),

Considering the definitions (12) and varying , we obtain the following nonlinear system of second order ordinary differential equations:

 (13)

with

 G(t,d)=b1(t)∣∣m∑k=1dk(t)φk∣∣20,L1(t)=νA+ˆB3(i)(t), (14) L2(t)=b2(t)K2+(ˆB1(i)(t)+ˆB4(i)(t)−ˆB2(ij)(t)),

where we denote by the transpose of the matrix .

Now, to solve the system of ordinary nonlinear equations (13) in discrete time, we will use the finite difference method.

### 3.2 The Newmark’s method

In order, consider the discretize the fixed time interval , using a uniform time step of size , then ,  for and . Consider the Neumark’s approximations:

 χη+θ=θχη−1+(1−2θ)χη+θχη+1.

Substituting the approximations in the system (13) and taking the equation in the discrete time , we obtain

 (Mη+11+θ (Δt)2Gη+1(d)K1)dη+1+Mη2dη (15) +(Mη−13+θ (Δt)2Gη−1(d)K1)dη−1=(Δt)2Fη+θ,

or equivalently

 (Mη+11+θ (Δt)2Gη+1(d)K1)dη+1=−Mη2dη−(Mη−13+θ (Δt)2Gη−1(d)K1)dη−1+(Δt)2Fη+θ, (16)

where

 Mη+11 =A+Δt2Lη+11+θ(Δt)2 Lη+12,Mη2=(Δt)2(1−2θ)(Gη(d)K1+Lη2)−2A, Mη−13 =A−Δt2Lη−11+θ(Δt)2 Lη−12,Fη+θ=θFη−1+(1−2θ)Fη+θFη+1.

#### Remark.

For the initial step, that is, for , appears in the equation the term , that does not exist. Proceeding formally, we use a quadratic order approximation for the initial data . Therefore, we consider , so we have . So, for in (16), we have

 [M11+M03+θ (Δt)2(G1(d)+G−1(d))K1]d1=−M02d0+2Δt M03d1+(Δt)2F0+θ, (17)

where we assume the approximation and

 G−1(X)=b01∣∣m∑k=1(d1k−2Δt d1(k))∇φk(y)∣∣20,

considering the approximation .

To determine the solution vector , we must now solve the algebraic system (17) and (16), for . As the system is not linear, we will employ the Newton Method, which we briefly describe below, aiming at an optimization of the computational cost in the calculations.

### 3.3 Newton’s method

For each , consider the left side of the equation in (17). In this way, solve (17) is equivalent to finding the root of . By Newton’s method, given the initial approximation and denoting ,, we obtain a new approximate solution of , given by , where is the linear system solution . The key point of the method is the calculation of the Jacobian matrix . This is because, if it has a high computational cost, it can make the method unfeasible.

Define

 Fη(X)=(Mη+11+θ(Δt)2Gη+1(X)K1)X+Γη, (18)

where .

For , follows

 F0(X)=[M11+M03+θ(Δt)2(G1(X)+G−1(X))K1]X+2θ(Δt)3G−1(X)K1d1+Γ0, (19)

where .

For :

 JFη(X)=d FηdX(X) =Mη+11+θ(Δt)2ddX[Gη+1(X)K1X] (20) =Mη+11+θ(Δt)2(K1X dGη+1dX(X)+Gη+1(X)K1),

and for :

 JF0(X) =M11+M03+θ(Δt)2[K1X (dG1dX(X)+dG−1dX(X))+(G1(X)+G−1(X))K1] (21) +2θ(Δt)3 K1d1 dG−1dX(X),

where

 dGηdX(X)=[∂Gη∂X1(X)  ⋯  ∂Gη∂Xm(X)]m×1,∀η∈{−1}∪{1,⋯,N−1}.

Using (14), we obtain

 ∂Gη∂Xk(X)=2 bη1 Xk [∫Ω(∇yiφk(y))2 dx], ∂G−1∂Xk(X)=2 b01 (Xk−2Δt d1k) [∫Ω(∇yiφk(y))2 dx],

where is the -th coordinate of the initial velocity vector .

From there, we apply Newton’s refinement algorithm, for , until one of the stopping criteria is satisfied: or , where is the solution obtained in the -th iteration of Newton. We also define a maximum limit for iterations in order to avoid infinite refinements.

### 3.4 Finite Element Method- Hermite Base

Let be a polygonalization family of , where , satisfying the minimum angle condition (see [7]) and indexed by the parameter , representing the maximum diameter of the elements . Given an integer , we introduce the following finite element space

spanned by the Hermite interpolation polynomials basis functions. More specifically, we let

 Vh={ψh∈V3h∩C1(Ω); ψh∣∣K∈Pl(K), ∀ K∈Th}⊂V,

where represents a finite-dimensional function space generated by piecewise polynomial functions of degree less or equal associated with a discretization by finite elements of size . In particular, in this paper we will use the Hermite polynomials of degree equal .

From the Lemma of Douglas-Dupont (see in [7]), it follows that, given a function , there is an interpolator such that

 ∥w(t)−ˆwh(t)∥m≤C1 hl+1−m∥w(t)∥l+1,

where denotes the seminorm over Hilbert space and .

Knowing the Hermite polynomials, then all the matrices of the system (16) can be calculated explicitly and solving the nonlinear system by Newton’s method gives the approximate solution, given by

 vh(y,t)=m∑i=1di(t)φi(y),φi(y)∈Vh.

We note that since the Hermite polynomials are orthogonal then

## 4 Numerical analysis

In this section we present error estimates for the finite element semi-discrete and fully discrete approximations of the nonlinear beam equation (8) or equivalently for the weak formulation of the nonlinear beam equation with moving ends.

Remarks on regularity: In order to get an error estimate, we require more regularity for the function and its derivatives in time:

 v′,v′′∈L2(0,T;H4(Ω)),  v′′′∈L∞(0,T;H1(Ω))  and  v(iv)∈L∞(0,T;L2(Ω)). (H2)

### 4.1 Error estimates for the semi-discrete problem

The semi-discrete problem consists of finding such that satisfies the problem (8) restricted to , that is, to find for all such that

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(v′′h(t),wh)+b1(t)∣∣∇vh(t)∣∣20(∇vh(t),∇wh)+b2(t)(Δvh(t),Δwh)+ν (v′h(t),wh)+(a(1)i(t)∇yivh(t),∇yiwh)+(a(2)ij(t)∇yi,yjvh(t),wh)+(a(4)i(t)∇yiv′h(t),wh)+(a(5)i(t)∇yivh(t),wh)=0, in Q, (vh(y,0),wh)=(v0h(y),wh),  (v′h(y,0),wh)=(v1h(y),wh)in  Ω, (vh,wh)=0 and (∇vh,wh)≡0, on ∂Ω×[0,T[. (22)

For the estimates, we will need the following hypothesis about the projections of in , denoted by and defined in (27), of the initial data

 |v1h−Rhv1|0 ≤ c1h4,|∇v0h−Rh∇v0|0 ≤ c2h3and|Δv0h−RhΔv0|0 ≤ c3h2, (H3)

where and are the projections of and in , respectively, and are constants independent of .

Let us assume that,

 maxt∈[0,T](K′(t))2<ζ04 (H4)

which requires that the speed of the end points be smaller than the characteristic speed of the equation.

Let us denote by , the error between the exact and approximate solution. Then, we have the following error estimate:

###### Theorem 4.

(Semi-discrete problem - Estimate)

Let the solution of the problem (8). Under hypotheses H1– (H4) and given the initial data , there exists a constant , such that

 ∥e′∥L∞(0,T;L2(Ω))+∥e∥L∞(0,T;H20(Ω)) ≤ Ch2. (23)

where the constant is dependent on and independent of .

###### Proof.

We note that, as , taking in particular in (8) and doing the difference with (22), denoting , and , we obtain:

 (24) +ν (e′,wh)+(a(1)i(t)∇yie ,∇yiwh)+(a(2)ij(t)∇y