Robust Iterative Solvers for Gao Type Nonlinear Beam Models in Elasticity

Robust Iterative Solvers for Gao Type Nonlinear Beam Models in Elasticity

Abstract

1 Introduction

2 Preliminaries

3 Numerical Solution of the Beam Problem

4 Conclusions

References

Abstract

1 Introduction

2 Preliminaries

3 Numerical Solution of the Beam Problem

4 Conclusions

References

Journal and Issue

Articles in the same Issue

2.1 Nonlinear Gao Beam Models

2.2 Some Iterative Methods

Assumptions \ref{subsiter}.

Theorem 2.1 (Simple Iteration/Gradient Method).

Proof.

Theorem 2.2 (Newton’s Method).

Proof.

Theorem 2.3 (Quasi-Newton/Variable Preconditioning Method).

Proof.

Remark 2.1.

2.3 Estimates in Sobolev spaces

Proposition 2.1 (Generalized Hölder Inequality).

Proposition 2.2 (Sobolev Embedding).

3.1 Weak Formulation and Finite Elements

Corollary 3.1.

3.2 Properties of the Linearized Operator

Proposition 3.1.

Proof.

Corollary 3.2.

Proposition 3.2.

Proof.

3.3 Finite Elements Using Splines

3.4 The Iterative Solvers: Construction, Convergence and Mesh-Independence

Construction \ref{subsmesh}.

Proposition 3.3.

Proof.

Theorem 3.1.

Proof.

Remark 3.1.

3.5 Generalizations

3.6 Numerical Experiments

2.1 Nonlinear Gao Beam Models

2.2 Some Iterative Methods

Assumptions \ref{subsiter}.

Theorem 2.1 (Simple Iteration/Gradient Method).

Proof.

Theorem 2.2 (Newton’s Method).

Proof.

Theorem 2.3 (Quasi-Newton/Variable Preconditioning Method).

Proof.

Remark 2.1.

2.3 Estimates in Sobolev spaces

Proposition 2.1 (Generalized Hölder Inequality).

Proposition 2.2 (Sobolev Embedding).

3.1 Weak Formulation and Finite Elements

Corollary 3.1.

3.2 Properties of the Linearized Operator

Proposition 3.1.

Proof.

Corollary 3.2.

Proposition 3.2.

Proof.

3.3 Finite Elements Using Splines

3.4 The Iterative Solvers: Construction, Convergence and Mesh-Independence

Construction \ref{subsmesh}.

Proposition 3.3.

Proof.

Theorem 3.1.

Proof.

Remark 3.1.

3.5 Generalizations

3.6 Numerical Experiments

3.5.1 Other Boundary Conditions

3.5.2 A Modified Equation

3.5.3 Extension to Plane Problems

3.5.1 Other Boundary Conditions

3.5.2 A Modified Equation

3.5.3 Extension to Plane Problems

Open Access Published by De Gruyter June 9, 2021

Benjámin Borsos and János Karátson

From the journal Computational Methods in Applied Mathematics

https://doi.org/10.1515/cmam-2020-0133

The goal of this paper is to present various types of iterative solvers: gradient iteration, Newton’s method and a quasi-Newton method, for the finite element solution of elliptic problems arising in Gao type beam models (a geometrical type of nonlinearity, with respect to the Euler–Bernoulli hypothesis). Robust behaviour, i.e., convergence independently of the mesh parameters, is proved for these methods, and they are also tested with numerical experiments.

Keywords: Gao Beam; Nonlinear Elasticity; Iterative Solvers

MSC 2010: 65N30; 35J60

The numerical study of the deformation of thin beams and plates is a widespread problem in elasticity theory and engineering practice since such elastic structures regularly appear in several real applications; see, e.g., [2, 12, 16, 17, 18]. These models generally lead to fourth-order equations. The linear models, which are basic in engineering applications, can be used only for small deformations. The most popular linear beam model is the Euler–Bernoulli beam: if the elastic modulus E and the moment of inertia I are constant, then one obtains a fourth-order ODE with constant coefficient. The analogous two-dimensional model for a thin plate involves the biharmonic operator. The obtained simple models read as

(1.1) D ⁢ u I ⁢ V = q and D ⁢ Δ 2 ⁢ u = q ,

where D = E ⁢ I is the flexural rigidity and q is the distributed load. However, such linear models are no more valid unless the deformation can be considered as infinitesimal. In more realistic models of engineering structures involving larger deformations, nonlinear behaviour has to be taken into account. Some of these models have been introduced in [9, 10], leading to so-called Gao beam models; see [11, 14] for more recent applications. Gao beam theory respects the Euler–Bernoulli hypothesis, but involves a geometrical type of nonlinearity [15, 19]. The involved nonlinearity requires the application of efficient solvers for the nonlinear system arising for the finite element approximation (FEM).

The goal of this paper is to present various types of such iterative solvers in the setting of the finite element method (FEM) and, in particular, to show the robust behaviour of these methods, i.e., convergence independently of the mesh parameters. The presentation of these methods, based on a Hilbert space framework, includes the proper forms of the Sobolev gradient iteration and of Newton’s method adapted to the given beam problem. Further, based on the authors’ recent results [4, 5] (see also [13]), a quasi-Newton/variable preconditioning method is presented as an intermediate version between the above two. Here the balance between speed and cost is achieved with auxiliary operators chosen via spectral equivalence.

The paper first summarizes the corresponding theory in Section 2. Section 3 contains the main results. After the weak formulation of the problem, the necessary properties of the corresponding operators are proved in the used finite element subspaces, and robust convergence is derived. Then the results of numerical experiments are presented. A brief conclusion is given in the closing section.

For the description of the bending of a beam resting on a foundation, to treat deformations beyond the infinitesimal region, models have been developed in [9, 11] derived from the presence of lateral stress. The following model considers a beam with classical Winkler foundation, which is a widespread concept in civil engineering, also with a profound effect on the field of adhesion mechanics; see [7]. Here the deflection u of the beam is described by the following equation:

(2.1) E I u I ⁢ V - E α ( u ′ ) 2 u ′′ + k F u = f in J : = [ 0 , b ] ,

with the following constants: E > 0 is the elastic modulus, I > 0 is the moment of inertia for the beam’s cross-section, α = 3 ⁢ h ⁢ ( 1 - ν 2 ) , where h is thickness measured from the axis and ν > 0 denotes the Poisson ratio, and k F > 0 is the foundation stiffness coefficient. Further, the transverse distributed load function is denoted by q, and f = ( 1 - ν 2 ) ⁢ q .

A slightly modified version of the above equation, involved, e.g., in contact problems, is

(2.2) E I u I ⁢ V - E α ( u ′ ) 2 u ′′ + P μ u ′′ = f in J : = [ 0 , b ] ,

where P is the axial force, assumed to be below the Euler critical load P cr E , and μ = ( 1 + ν ) ⁢ ( 1 - ν 2 ) .

Here we summarize some iterative methods in a Hilbert space framework. Since we will apply them to our given beam equation, we formulate these theorems under a common set of conditions that can be verified for the beam problem. For more details, we refer to the monograph [8] and the authors’ recent papers [4, 5] (see also [13]). In this section, let H be a real Hilbert space and F : H → H a nonlinear operator which has a bihemicontinuous Gâteaux derivative.

We formulate the following properties for F ′ .

(Symmetry.) For any u ∈ H , the operator F ′ ⁢ ( u ) is self-adjoint.
(Regularity.) There exists a constant λ > 0 such that

(2.3) λ ∥ h ∥ 2 ≤ 〈 F ′ ( u ) h , h 〉 ( for all u , h ∈ H ) .
(Upper growth.) There exists a continuous increasing function Λ : ℝ + → ℝ + such that

(2.4) 〈 F ′ ( u ) h , h 〉 ≤ Λ ( ∥ u ∥ ) ∥ h ∥ 2 ( for all u , h ∈ H ) .
(Local Lipschitz.) There exists a continuous increasing function L : ℝ + → ℝ + such that

∥ F ′ ( u ) - F ′ ( v ) ∥ ≤ L ( max { ∥ u ∥ , ∥ v ∥ } ) ∥ u - v ∥ ( for all u , h ∈ H ) .

We wish to solve the operator equation F ⁢ ( u ) = 0 . Conditions (i)–(ii) imply that there is a unique solution u ∗ ∈ H ; see, e.g., [8]. We study three types of iterative methods.

Let Assumptions 2.2 (i)–(iii) hold. Let u 0 ∈ H be arbitrary and Λ 0 : = Λ ( ∥ u 0 ∥ + 1 λ ∥ F ( u 0 ) ∥ ) . Then the sequence defined by

u n + 1 : = u n - 2 Λ 0 + λ F ( u n ) ( for all n ∈ ℕ )

converges to u ∗ , namely,

∥ u n - u ∗ ∥ ≤ 1 λ ∥ F ( u 0 ) ∥ ( Λ 0 - λ Λ 0 + λ ) n ( for all n ∈ ℕ ) .

Moreover, in fact, ∥ u n - u ∗ ∥ ≤ 1 λ ⁢ ∥ F ⁢ ( u n ) ∥ and

∥ F ⁢ ( u n + 1 ) ∥ ∥ F ⁢ ( u n ) ∥ ≤ Λ 0 - λ Λ 0 + λ .

See [8, Theorem 5.5]. ∎

Let Assumptions 2.2 (i), (ii) and (iv) hold. Let u 0 be in a properly small neighbourhood of u ∗ and L 0 : = L ( ∥ u 0 ∥ + 2 λ ∥ F ( u 0 ) ∥ ) . Then the sequence defined by

u n + 1 : = u n - F ′ ( u n ) - 1 F ( u n ) ( for all n ∈ ℕ )

converges to u ∗ , that is, ∥ u n - u ∗ ∥ ≤ 1 λ ⁢ ∥ F ⁢ ( u n ) ∥ → 0 and

∥ F ( u n + 1 ) ∥ ≤ L 0 2 ⁢ λ 2 ∥ F ( u n ) ∥ 2 ( n ∈ ℕ ) .

This follows from [8, Theorem 5.9 and Remark 5.17]. ∎

The formulation of the third theorem uses the energy * -norm ∥ v ∥ * : = 〈 F ′ ( u * ) - 1 v , v 〉 1 / 2 .

Let Assumptions 2.2 (i)–(iv) hold. Let u 0 be in a sufficiently small neighbourhood of u ∗ . Let the sequence ( u n ) be defined by

u n + 1 : = u n - 2 M n + m n B n - 1 F ( u n ) ( n ∈ ℕ ) ,

where 0 < m n ≤ M n and the symmetric linear operators B n : H → H satisfy

(2.5) m n 〈 B n h , h 〉 ≤ 〈 F ′ ( u n ) h , h 〉 ≤ M n 〈 B n h , h 〉 ( n ∈ ℕ , h ∈ H ) .

We require ( m n ) to be positively bounded from below and ( M n ) bounded from above. Then ( u n ) converges to u ∗ , that is, ∥ u n - u ∗ ∥ ≤ 1 λ ⁢ ∥ F ⁢ ( u n ) ∥ → 0 and

(2.6) lim sup ⁡ ∥ F ⁢ ( u n + 1 ) ∥ * ∥ F ⁢ ( u n ) ∥ * ≤ lim sup ⁡ M n - m n M n + m n < 1 .

This follows from [4, Theorem 2.1] and [5, Theorem 2.5], which extend [8, Theorem 5.16] to the upper non-uniform case (2.4). ∎

We note that a reasonable requirement is to outperform the simple iteration, for which one should choose m n and M n closer than the original spectral bounds of F ′ ⁢ ( u n ) .

(i) (Efficiency.) Obviously, as is well known, Newton’s method is the fastest and the simple iteration is the slowest of the above methods (quadratic speed vs linear speed), but Newton’s method is more costly. The quasi-Newton method is an intermediate version, where the variable preconditioners B n may be cheaper than the full linearized operators; still, its convergence can be superlinear when the lim sup in (2.6) equals 0.

(ii) (Damped versions, global convergence.) The local convergence of the above versions of the Newton and quasi-Newton methods can be made global by damping with proper parameters τ n ≤ 1 . These are not formulated here for simplicity. The convergence speed in the limit is the same as in the above versions.

Here let Ω ⊂ 𝐑 d be a bounded domain. In most of the paper, we will let d = 1 and Ω = J (the studied interval). The following statements, which are well known (see, e.g., [1]), will be useful later.

If the numbers p i ≥ 1 satisfy ∑ i = 1 s 1 p i = 1 , then, for any functions u i ∈ L p i ⁢ ( Ω ) , we have ∫ Ω ∏ i = 1 s | u i | ≤ ∏ i = 1 s ∥ u i ∥ L p i .

Let 1 ≤ p (if d ≤ 2 ) or 1 ≤ p ≤ 2 ⁢ d d - 2 (if d > 2 ). Then

(2.7) H 0 1 ( Ω ) ⊂ L p ( Ω ) , ∥ v ∥ L p ≤ C p ∥ ∇ v ∥ L 2 ( for all v ∈ H 0 1 ( Ω ) )

for some constant C p > 0 independent of v.

We rewrite (2.1) by dividing with EI and letting

β : = α 3 ⁢ I , k : = k F E ⁢ I , g : = f E ⁢ I ;

further, we impose clamped boundary conditions. Then our problem becomes

(3.1) u I ⁢ V - β ⁢ ( ( u ′ ) 3 ) ′ + k ⁢ u = g , u ⁢ ( 0 ) = u ′ ⁢ ( 0 ) = u ⁢ ( b ) = u ′ ⁢ ( b ) = 0 .

The weak formulation uses the Sobolev space H 2 ⁢ ( J ) ; moreover, using the boundary conditions, we work in the space

H 0 2 ⁢ ( J ) = { u ∈ H 2 ⁢ ( J ) : u ⁢ ( 0 ) = u ′ ⁢ ( 0 ) = u ⁢ ( b ) = u ′ ⁢ ( b ) = 0 } ,

which has the standard inner product and corresponding norm

(3.2) 〈 u , v 〉 H 0 2 : = ∫ 0 b u ′′ v ′′ , ∥ u ∥ H 0 2 2 = ∫ 0 b ( u ′′ ) 2 .

Then the weak solution of problem (3.1) is defined as u ∗ ∈ H 0 2 ⁢ ( J ) satisfying

(3.3) ∫ 0 b ( u ∗ ′′ v ′′ + β ( u ∗ ′ ) 3 v ′ + k u ∗ v ) = ∫ 0 b g v ( for all v ∈ H 0 2 ( J ) ) .

The well-posedness of the problem will be discussed in the next section. The weak solution minimizes the energy

E ( u ) : = ∫ 0 b ( 1 2 ( u ′′ ) 2 + β 4 ( u ′ ) 4 + k 2 u 2 - g u )

in H 0 2 ⁢ ( J ) .

The finite element method (FEM) can be used to solve (3.3) numerically. Let V h ⊂ H 0 2 ⁢ ( J ) be a given FEM subspace. Then we seek for u h ∈ V h satisfying

(3.4) ∫ 0 b ( u h ′′ v ′′ + β ( u h ′ ) 3 v ′ + k u h v ) = ∫ 0 b g v ( for all v ∈ V h ) .

The formulation and properties of the problem will be given for a general FEM subspace V h ⊂ H 0 2 ⁢ ( J ) . Later, for the realization, we will specify V h to consist of suitable spline functions.

Rearranging (3.3) and using Riesz representation, one can define an operator F : V h → V h satisfying

(3.5) 〈 F ( u ) , v 〉 H 0 2 = ∫ 0 b ( u ′′ v ′′ + β ( u ′ ) 3 v ′ + k u v - g v ) ( for all v ∈ V h )

so that (3.3) simply becomes an equation

(3.6) F ⁢ ( u ) = 0

in the space V h (also equipped with the H 0 2 inner product). Then, using standard calculations (see, e.g., [8, Chapter 6.1]), one can derive that F has a Gâteaux derivative satisfying

(3.7) 〈 F ′ ( u ) h , v 〉 H 0 2 = ∫ 0 b ( h ′′ v ′′ + 3 β ( u ′ ) 2 h ′ v ′ + k h v ) ( for all u , h , v ∈ V h ) ,

further, that F ′ is bihemicontinuous. In addition, the role of h and v is interchangeable in formula (3.7), which readily implies the following corollary.

For any u ∈ V h , the operator F ′ ⁢ ( u ) is self-adjoint.

The following ellipticity properties hold.

There exists a continuous increasing function Λ : R + → R + , independently of h, such that

(3.8) ∥ h ∥ H 0 2 2 ≤ 〈 F ′ ( u ) h , h 〉 H 0 2 ≤ Λ ( ∥ u ∥ H 0 2 ) ∥ h ∥ H 0 2 2 ( for all u , h ∈ V h ) .

Namely, Λ ⁢ ( t ) = 1 + k ⁢ C 2 4 + 3 ⁢ β ⁢ C 4 4 ⁢ t 2 .

Owing to (3.7), the quadratic form reads as

〈 F ′ ( u ) h , h 〉 H 0 2 = ∫ 0 b ( ( h ′′ ) 2 + 3 β ( u ′ ) 2 ( h ′ ) 2 + k h 2 ) ( for all u , h ∈ V h ) .

Since β , k ≥ 0 , the first inequality of (3.8) holds.

Furthermore, (2.7) implies that ∥ z ′ ∥ L 4 ≤ C 4 ⁢ ∥ z ′′ ∥ L 2 = C 4 ⁢ ∥ z ∥ H 0 2 . Thus we have

∫ 0 b 3 β ( u ′ ) 2 ( h ′ ) 2 ≤ 3 β ∥ ( u ′ ) 2 ∥ L 2 ∥ ( h ′ ) 2 ∥ L 2 = 3 β ∥ u ′ ∥ L 4 2 ∥ h ′ ∥ L 4 2 ≤ 3 β C 4 4 ∥ u ∥ H 0 2 2 ∥ h ∥ H 0 2 2 ( for all u , h ∈ V h ) ,

and similarly,

∫ 0 b k h 2 = k ∥ h ∥ L 2 2 ≤ k C 2 2 ∥ h ′ ∥ L 2 2 ≤ k C 2 4 ∥ h ′′ ∥ L 2 2 = k C 2 4 ∥ h ∥ H 0 2 2 ( for all u , h ∈ V h ) ;

hence the result follows. ∎

The above shows that Assumptions 2.2 (i)–(ii) hold for problem (3.6). As seen in Section 2.2, this implies well-posedness.

Problem (3.3) has a unique solution u ∗ ∈ H 0 2 ⁢ ( J ) ; further, for any FEM subspace V h ⊂ H 0 2 ⁢ ( J ) , problem (3.4) has a unique solution u h ∈ V h .

One can also establish local Lipschitz continuity.

There exists a continuous increasing function L : R + → R + , independently of h, such that

∥ F ′ ( u ) - F ′ ( v ) ∥ ≤ L ( max { ∥ u ∥ H 0 2 , ∥ v ∥ H 0 2 } ) ∥ u - v ∥ H 0 2 ( for all u , h ∈ V h ) .

Namely, L ⁢ ( t ) = 6 ⁢ C 4 4 ⁢ β ⁢ t .

Corollary 3.1 implies the symmetry of F ′ ⁢ ( u ) - F ′ ⁢ ( v ) for each u , v ∈ V h ; therefore, it is possible to obtain its norm using its quadratic form,

∥ F ′ ⁢ ( u ) - F ′ ⁢ ( v ) ∥ = sup ∥ h ∥ H 0 2 = 1 ⁡ | 〈 ( F ′ ⁢ ( u ) - F ′ ⁢ ( v ) ) ⁢ h , h 〉 | = sup ∥ h ∥ H 0 2 = 1 ⁡ | ∫ 0 b 3 ⁢ β ⁢ ( ( u ′ ) 2 - ( v ′ ) 2 ) ⁢ ( h ′ ) 2 | .

As in the proof of Proposition 3.1 above, we use Proposition 2.1 and the fact that ∥ z 2 ∥ L 2 = ∥ z ∥ L 4 2 . With Proposition 2.2, these yield

(3.9) ∥ F ′ ⁢ ( u ) - F ′ ⁢ ( v ) ∥ ≤ sup ∥ h ∥ H 0 2 = 1 ⁡ 3 ⁢ β ⁢ ∥ ( u ′ ) 2 - ( v ′ ) 2 ∥ L 2 ⁢ ∥ ( h ′ ) 2 ∥ L 2 ≤ sup ∥ h ∥ H 0 2 = 1 ⁡ 3 ⁢ β ⁢ ∥ ( u ′ ) 2 - ( v ′ ) 2 ∥ L 2 ⁢ C 4 2 ⁢ ∥ h ′′ ∥ L 2 2 .

Since ∥ z ′′ ∥ L 2 = ∥ z ∥ H 0 2 , this readily entails

(3.10) ∥ F ′ ⁢ ( u ) - F ′ ⁢ ( v ) ∥ ≤ 3 ⁢ C 4 2 ⁢ β ⁢ ∥ ( u ′ ) 2 - ( v ′ ) 2 ∥ L 2 = 3 ⁢ C 4 2 ⁢ β ⁢ ∥ ( u ′ - v ′ ) 2 ⁢ ( u ′ + v ′ ) 2 ∥ L 1 .

Repeating the technique used in (3.9) gives

(3.11) ∥ ( u ′ - v ′ ) 2 ⁢ ( u ′ + v ′ ) 2 ∥ L 1 ≤ ∥ ( u ′ - v ′ ) 2 ∥ L 2 ⁢ ∥ ( u ′ + v ′ ) 2 ∥ L 2 = ∥ u ′ - v ′ ∥ L 4 ⁢ ∥ u ′ + v ′ ∥ L 4 ≤ C 4 2 ⁢ ∥ u ′′ - v ′′ ∥ L 2 ⁢ ∥ u ′′ + v ′′ ∥ L 2

Combining (3.10)–(3.11) yields

∥ F ′ ⁢ ( u ) - F ′ ⁢ ( v ) ∥ ≤ 3 ⁢ C 4 4 ⁢ β ⁢ ∥ u - v ∥ H 0 2 ⁢ ∥ u + v ∥ H 0 2

Finally, since ∥ u + v ∥ H 0 2 ≤ ∥ u ∥ H 0 2 + ∥ v ∥ H 0 2 ≤ 2 ⁢ max ⁡ ( ∥ u ∥ H 0 2 , ∥ v ∥ H 0 2 ) , this yields

∥ F ′ ⁢ ( u ) - F ′ ⁢ ( v ) ∥ ≤ 6 ⁢ C 4 4 ⁢ β ⁢ max ⁡ ( ∥ u ∥ H 0 2 , ∥ v ∥ H 0 2 ) ⁢ ∥ u - v ∥ H 0 2 . ∎

The straightforward way to implement finite elements for fourth-order beam problems is to use piecewise cubic splines, that is, Hermitian elements with four degrees of freedom, which satisfy C 1 -continuity at the node points. These functions are in H 2 ⁢ ( J ) , and accordingly, the integrals in (3.2) are finite. Alternatively, one could use quadratic B-splines and still achieve C 1 -continuity.

We apply a uniform mesh, where the piecewise cubic basis functions are constructed below for mesh parameter h. Two such functions are obtained for each interior node. These nodes are denoted n k , where k ∈ K : = { 1 , 2 , … , b / h - 1 } , and n k is at location x k = h ⁢ k .

Let us consider piecewise cubic functions f 1 , f 2 : [ - 1 , 1 ] → ℝ , which are defined as

f i ( x ) = { f i * ⁢ ( x ) , x ∈ [ 0 , 1 ] , ( - 1 ) ( i - 1 ) ⁢ f i * ⁢ ( - x ) , x ∈ [ - 1 , 0 ) , ( i = 1 , 2 ) ,

where f 1 * ⁢ ( x ) = 2 ⁢ x 3 - 3 ⁢ x 2 + 1 , f 2 * ⁢ ( x ) = x 3 - 2 ⁢ x 2 + x ; see Figure 1.

Figure 1

Piecewise cubic functions on [ - 1 , 1 ] .

The basis functions are obtained via affine transformations L k ( k ∈ K ) such that the domain of L k ⁢ ( f i ) is [ h ⁢ ( k - 1 ) , h ⁢ ( k + 1 ) ] for i = 1 , 2 .

Now we come to the main point of the numerical method, which is the iterative solution of the FEM problem (3.5)–(3.6). We consider the iterative methods of Section 2.2. For notational simplicity, we do not indicate in which subspace V h our sequences ( u n ) n ∈ ℕ run in. In fact, we think of V h as being fixed.

The three recurrences can be summarized as follows:

(3.12) u n + 1 : = u n - σ z n ,

where the function z n ∈ V h and the constant σ > 0 are defined in the way described below. We note that, for the Newton and quasi-Newton methods, one formally has to solve an operator equation; further, the quasi-Newton method will be uniquely defined if we choose proper operators B n .

Simple iteration/gradient method: z n : = F ( u n ) , σ : = 2 Λ 0 + λ .

To determine z n , let us write the equality with test functions,

〈 z n , v 〉 H 0 2 = 〈 F ( u n ) , v 〉 H 0 2 ( for all v ∈ V h ) .

Here and henceforth, let us define the right-hand side as a functional ℓ n , that is, by (3.5),

(3.13) ℓ n v : = 〈 F ( u n ) , v 〉 H 0 2 ≡ ∫ 0 b ( u n ′′ v ′′ + β ( u n ′ ) 3 v ′ + k u n v - g v ) ( for all v ∈ V h ) .

Then the z n ∈ V h is the solution of the auxiliary FEM problem

∫ 0 b z n ′′ v ′′ = ℓ n v ( for all v ∈ V h ) .
Newton’s method: F ′ ⁢ ( u n ) ⁢ z n = F ⁢ ( u n ) , σ : = 1 .

To determine z n , let us again involve test functions,

〈 F ′ ( u n ) z n , v 〉 H 0 2 = 〈 F ( u n ) , v 〉 H 0 2 ( for all v ∈ V h ) .

Using (3.7) and (3.13), z n ∈ V h is the solution of the auxiliary FEM problem

∫ 0 b ( z n ′′ v ′′ + 3 β ( u n ′ ) 2 z n ′ v ′ + k z n v ) = ℓ n v ( for all v ∈ V h ) .
Quasi-Newton/variable preconditioning: B n ⁢ z n = F ⁢ ( u n ) , σ : = 2 M n + m n .

Now the equation for z n with test functions is

〈 B n z n , v 〉 H 0 2 = 〈 F ( u n ) , v 〉 H 0 2 ( for all v ∈ V h ) .

Here we propose to choose the operator B n as an approximation of F ′ ⁢ ( u n ) such that the term 3 ⁢ β ⁢ ( u n ′ ) 2 is replaced by proper constants w n > 0 : let

(3.14) 〈 B n h , v 〉 H 0 2 : = ∫ 0 b ( h ′′ v ′′ + w n h ′ v ′ + k h v ) ( for all v ∈ V h ) .

That is, by (3.13), z n ∈ V h is the solution of the auxiliary FEM problem

∫ 0 b ( z n ′′ v ′′ + w n z n ′ v ′ + k z n v ) = ℓ n v ( for all v ∈ V h ) .

It is informative to summarize also the strong versions of the auxiliary problems, which are formal equations involving fourth derivatives. Namely, let us denote

r n : = u n I ⁢ V - β ( ( u n ′ ) 3 ) ′ + k u n - g .

It is readily seen by integration that (for smooth u n )

ℓ n v = ∫ 0 b r n v ( for all v ∈ V h ) .

Using similar formal integration for the left-hand side, the auxiliary equations become the following, where, in each problem, the boundary conditions are z n ⁢ ( 0 ) = z n ′ ⁢ ( 0 ) = z n ⁢ ( b ) = z n ′ ⁢ ( b ) = 0 :

simple iteration/gradient method: z n I ⁢ V = r n ,
Newton’s method: z n I ⁢ V - 3 ⁢ β ⁢ ( ( u n ′ ) 2 ⁢ z n ′ ) ′ + k ⁢ z n = r n ,
quasi-Newton/variable preconditioning: z n I ⁢ V - w n ⁢ z n ′′ + k ⁢ z n = r n .

That is, the auxiliary problems correspond to the solution of proper linear fourth-order ODEs. In reality, of course, we consider the FEM solution of the weak versions of these problems, where u n is only an H 2 function.

Now we derive the convergence of the iterations, based on the properties in Section 3.2. In addition, we need the spectral equivalence property (2.5) for the above defined operators B n and F ′ ⁢ ( u n ) . A reasonable choice of w n is in the range of the function 3 β max Ω ¯ ( u n ′ ) 2 that it approximates,

(3.15) 0 ≤ w n ≤ 3 β max ( u n ′ ) 2 .

A convenient choice is the arithmetic mean

(3.16) w n : = 3 ⁢ β 2 max ( u n ′ ) 2 .

For given u n ∈ V h , let the constant w n satisfy (3.15). Then

m n 〈 B n h , h 〉 H 0 2 ≤ 〈 F ′ ( u n ) h , h 〉 H 0 2 ≤ M n 〈 B n h , h 〉 H 0 2 ( for all h ∈ V h ) ,

where

(3.17) m n = 1 1 + w n ⁢ C 2 2 , M n = 3 β max ( u n ′ ) 2 w n .

Firstly, to obtain the upper bound, by (3.7), one can write

〈 F ′ ⁢ ( u n ) ⁢ h , h 〉 H 0 2 = ∫ 0 b ( ( h ′′ ) 2 + 3 ⁢ β ⁢ ( u n ′ ) 2 ⁢ ( h ′ ) 2 + k ⁢ h 2 ) ≤ ∫ 0 b ( ( h ′′ ) 2 + 3 β max Ω ¯ { ( u n ′ ) 2 } ( h ′ ) 2 + k h 2 ) ( for all h ∈ V h ) ;

on the other hand, owing to (3.15), we have 1 ≤ M n : = 3 β max Ω ¯ ( u n ′ ) 2 w n for each n ∈ ℕ ; hence

〈 F ′ ( u n ) h , h 〉 H 0 2 ≤ M n ∫ 0 b ( ( h ′′ ) 2 + w n ( h ′ ) 2 + k h 2 ) = M n 〈 B n h , h 〉 H 0 2 ( for all h ∈ V h ) .

To obtain the lower bound, by (3.14) and Proposition 2.2, we have

〈 B n h , h 〉 H 0 2 = ∥ h ′′ ∥ L 2 2 + w n ∥ h ′ ∥ L 2 2 + k ∥ h ∥ L 2 2 ≤ ( 1 + w n C 2 2 ) ∥ h ′′ ∥ L 2 2 + k ∥ h ∥ L 2 2 ( for all h ∈ V h ) ;

writing this back to integral form and adding the term 3 ⁢ β ⁢ ( u n ′ ) 2 ⁢ ( h ′ ) 2 ≥ 0 yields

〈 B n h , h 〉 H 0 2 ≤ ∫ 0 b ( ( 1 + w n C 2 2 ) ( h ′′ ) 2 + 3 β ( u n ′ ) 2 ( h ′ ) 2 + k h 2 ) ( for all h ∈ V h ) .

This readily entails 〈 B n ⁢ h , h 〉 H 0 2 ≤ ( 1 + w n ⁢ C 2 2 ) ⁢ 〈 F ′ ⁢ ( u n ) ⁢ h , h 〉 H 0 2 . ∎

Now we can formulate the convergence results.

The iterative methods, defined in Construction 3.4, provide the following convergence estimates:

simple iteration/gradient method:

∥ F ⁢ ( u n + 1 ) ∥ H 0 2 ∥ F ⁢ ( u n ) ∥ H 0 2 ≤ Λ 0 - λ Λ 0 + λ , 𝑤ℎ𝑒𝑟𝑒 ⁢ Λ 0 = 1 + k ⁢ C 2 4 + 3 ⁢ β ⁢ C 4 4 ⁢ ( ∥ u 0 ∥ + 1 λ ⁢ ∥ F ⁢ ( u 0 ) ∥ ) 2 ,
Newton’s method:

∥ F ⁢ ( u n + 1 ) ∥ H 0 2 ∥ F ⁢ ( u n ) ∥ H 0 2 2 ≤ L 0 2 ⁢ λ 2 , 𝑤ℎ𝑒𝑟𝑒 ⁢ L 0 = 6 ⁢ C 4 4 ⁢ β ⁢ ( ∥ u 0 ∥ + 2 λ ⁢ ∥ F ⁢ ( u 0 ) ∥ ) ,
quasi-Newton/variable preconditioning: if ( 3.15 ) holds, then

lim sup ⁡ ∥ F ⁢ ( u n + 1 ) ∥ * ∥ F ⁢ ( u n ) ∥ * ≤ lim sup ⁡ M n - m n M n + m n

with the constants in ( 3.17 ).

These hold globally for the simple iteration and locally for the Newton and quasi-Newton methods.

This follows from Theorems 2.1–2.3, Propositions 3.1–3.2 and Proposition 3.3. ∎

(a) The estimates in Theorem 3.1 are uniform, i.e., the constant on the right-hand side of the inequalities are mesh-independent.

(b) Global convergence for the Newton and quasi-Newton methods can be achieved via damped versions (see, e.g., [8, 5] for the abstract theorems), which are not detailed here. In this case, the above estimates are ultimate, i.e., they hold in lim sup sense also for Newton’s method.

(c) In the above, we considered the iterations on a fixed mesh. The methods might be generalized to a multilevel setting to increase the efficiency of preconditioning (see related work in [3, 6]), but such extensions are beyond the scope of the present paper.

Instead of the rigidly clamped beam in (3.1), one can consider a freely supported beam. Then the first derivatives at the endpoints are replaced by second derivatives, i.e., the problem becomes

(3.18) u I ⁢ V - β ⁢ ( ( u ′ ) 3 ) ′ + k ⁢ u = g , u ⁢ ( 0 ) = u ′′ ⁢ ( 0 ) = u ⁢ ( b ) = u ′′ ⁢ ( b ) = 0 .

This problem can be posed in the Sobolev space H 2 ⁢ ( J ) ∩ H 0 1 ⁢ ( J ) , which is endowed with the same inner product (3.2) as used in H 0 2 ⁢ ( J ) . Hence the calculations can be repeated, and the same convergence results hold as in Section 3.4. The auxiliary problems are obviously solved with the freely supported boundary conditions as in (3.18).

The above results can also be reproduced for the modified version (2.2). Owing to the fact that the axial force P is below the Euler critical load P cr E , it follows that the energy functional is uniformly convex (see, e.g., [14]), which implies that the uniform monotonicity (2.3) holds. The other conditions remain unchanged; hence the calculations can be repeated to obtain the same convergence results.

For the 2D plate problem in (1.1), the analogue of equation (3.1) is

Δ 2 ⁢ u - β ⁢ div ⁡ ( | ∇ ⁡ u | 2 ⁢ ∇ ⁡ u ) + k ⁢ u = g in ⁢ Ω , u | ∂ Ω = ∂ ⁡ u ∂ ⁡ ν | ∂ ⁡ Ω = 0

for a thin plate Ω ⊂ 𝐑 2 . The weak solution minimizes the energy

E ( u ) : = ∫ Ω ( 1 2 | D 2 u | 2 + β 4 | ∇ u | 4 + k 2 u 2 - g u )

in the Sobolev space H 0 2 ⁢ ( Ω ) . It is easy to see that our 1D results can be readily extended to this situation. The problem is posed in H 0 2 ⁢ ( Ω ) endowed with the inner product 〈 u , v 〉 H 0 2 : = ∫ Ω D 2 u : D 2 v . The main analytic point is that the required Sobolev embedding H 0 1 ⁢ ( Ω ) ⊂ L 4 ⁢ ( Ω ) also holds in 2D owing to (2.7). The other used techniques are independent of the dimension of the domain. In the case of the freely supported plate, the boundary conditions become u | ∂ Ω = ∂ 2 ⁡ u ∂ ⁡ ν 2 | ∂ ⁡ Ω = 0 and the problem is posed in the Sobolev space H 2 ⁢ ( Ω ) ∩ H 0 1 ⁢ ( Ω ) . Altogether, the calculations can be repeated to obtain the same convergence results as before.

The model described by (3.1) was used for simulation. The results are presented with the original physical parameters used in (2.1) for the sake of convenience. The physical and mesh parameters where chosen with the help of [14, 11].

The investigated problems included steel and concrete beams, with length L = 2 m, and we applied contact stiffness k = 3 ⋅ 10 8 ⁢ N m 2 . The steel and concrete beams have modulus of elasticity E 1 = 2.1 ⋅ 10 11 Pa and E 2 = 3 ⋅ 10 10 Pa, respectively, and Poisson’s ratio ν 1 = 0.3 and ν 2 = 0.2 , respectively. As second moment of area, I = 2 / 3 ⋅ 10 - 3 ⁢ m 4 was used. This results from a rectangular cross-section, namely, the beam height is h = 0.1 m, and the beam width is considered as a unit. The total vertical loadings are

F 1 = - 1.5 ⋅ 10 8 N, F 2 = - 3 ⋅ 10 8 N, F 3 = - 5 ⋅ 10 8 N for the steel beam and
F 4 = - 1 ⋅ 10 7 N, F 5 = - 4 ⋅ 10 7 N, F 6 = - 8 ⋅ 10 7 N for the concrete beam.

These loads are distributed uniformly along the beam, and the distributed force is q = F L (that is, q i = F i L in the distinct experiments). The number of finite elements is denoted by NE.

The simulations were carried out in Matlab. For all simulations, the initial guess was constant 0, and the stopping criterion was the relative residual decreased below 10 - 4 . A direct linear solver was used in case of the obtained linear problems in all iteration steps. For the quasi-Newton method, the suggested choice (3.16) was used for the preconditioner. It has been found that all the three methods converge and are robust.

The constant σ in (3.12) was replaced by different values in an attempt to achieve faster convergence. The investigation showed that σ = 1 is suitable for all methods; hence all three methods were considered with this constant. Damping was not necessary for these methods.

Tables 1, 2 and 3 show the number of iterations with the quasi-Newton method, the Sobolev gradient method and the full Newton method, respectively, to illustrate the robustness of the methods.

Table 1

Number of iterations using the quasi-Newton method.

Table 2

Number of iterations using the Sobolev gradient method.

Table 3

Number of iterations using the full Newton method.

The total runtimes (i.e., the runtimes of the whole simulation, not only an individual iteration step) for the quasi-Newton method, the Sobolev gradient method and the Newton method are denoted by t qN , t g and t N , respectively. These values were obtained by averaging the total runtimes of multiple simulations for each parameter combination. Namely, for the cases NE = 8 , 80 , 800 , 8000 , the number of simulations used to measure the total runtimes were 50000 , 5000 , 500 , 50 , respectively, and the averages of the measured runtimes were used. We have compared Newton’s method (considered as a standard nonlinear solver) with the two other ones in Tables 4 and 5.

In Table 4, all of the values t qN / t N are smaller than 1 for the investigated problems, i.e., quasi-Newton method is more efficient than full Newton method for the investigated problems with respect to computational cost. Furthermore, increasing mesh density further improves the relative performance of the quasi-Newton method, owing to the increasing benefit from the simplification of the stiffness matrix.

In Table 5, one can observe that (especially in the case of coarse meshes, e.g., 8 and 80 elements) the Sobolev gradient method may underperform the full Newton method in runtimes, e.g., for the concrete beam, though it is apparently faster for the steel beam. The relative computational cost of Sobolev gradient method also improves with increasing mesh density. Altogether, the Sobolev gradient method often performs well due to no assembling required in the steps.

In 10 of the cases under investigation, the quasi-Newton method is the most effective with respect to computational cost, while the Sobolev gradient method is favoured in the remaining 14 cases.

Table 4

Comparison of quasi-Newton and full Newton runtimes: the values of t qN / t N .

	E = E 1 , ν = ν 1			E = E 2 , ν = ν 2
NE	q = q 1	q = q 2	q = q 3	q = q 4	q = q 5	q = q 6
8	0.656	0.814	0.739	0.733	0.831	0.757
80	0.570	0.679	0.603	0.654	0.722	0.648
800	0.508	0.586	0.500	0.571	0.611	0.527
8000	0.458	0.454	0.370	0.451	0.362	0.372

Table 5

Comparison of Sobolev gradient and full Newton runtimes: the values of t g / t N .

	E = E 1 , ν = ν 1			E = E 2 , ν = ν 2
NE	q = q 1	q = q 2	q = q 3	q = q 4	q = q 5	q = q 6
8	0.567	0.662	0.735	1.377	1.541	1.945
80	0.466	0.541	0.599	1.165	1.328	1.664
800	0.357	0.417	0.440	0.890	0.921	1.114
8000	0.197	0.210	0.196	0.317	0.267	0.456

For one particular case, namely, E = E 2 , ν = ν 2 , q = q 6 , NE = 8000 , for the quasi-Newton method, the apparent fulfilment of Theorem 3.1 is illustrated for each iteration step n in Table 6.

Table 6

One quasi-Newton iteration.

n	1	2	3	4
∥ F ⁢ ( u n + 1 ) ∥ * / ∥ F ⁢ ( u n ) ∥ *	0.055	0.052	0.056	0.056
( M n - m n ) / ( M n + m n )	0.333	0.526	0.511	0.512

Additional experiments have been carried out to determine whether a change in Poisson’s ratio ν affects these results. For this task, Poisson’s ratio of the concrete beam was changed to ν = 0.3 and ν = 0.1 , while other parameters were left unchanged. For all methods, it has been found that the robustness result holds for the new parameter combinations as well; moreover, the corresponding iteration numbers remain almost unchanged. The relative total runtimes also qualitatively coincided with the original results.

Other experiments included replacing the contact stiffness k of the concrete beam with lower values, k = 2 ⋅ 10 7 ⁢ N m 2 and k = 0 ⁢ N m 2 , in the latter case, the model effectively lacking contact stiffness. For the large deformation of these soft models, Gao beam model (2.1) might not hold due to yielding of the material; however, the simulations can be carried out. The robustness result holds for these simulations. The supremacy of the quasi-Newton method over the full Newton method is also sustained, though the Sobolev gradient method exhibits relative improvement.

One can conclude that all three examined methods are robust, and quasi-Newton method can replace full Newton method for this nonlinear model. The Sobolev gradient method is very efficient for thousands of elements and more; otherwise, one should use quasi-Newton method. These coarse meshes appear to be of significance, as [14] states that 32 elements already suffice for accurate computations.

The present paper provides a detailed description of three iterative methods for nonlinear Gao beam models using finite elements. The description includes the Sobolev gradient method, Newton’s method and a quasi-Newton method with a recently developed framework for elliptic problems with non-uniformly monotone bounds. The results of both theoretical and practical work are shown. It has been found that all three methods are robust for the problems, and a comparison has been made between their behaviour under varying parameters.

Funding source: Nemzeti Kutatási Fejlesztési és Innovációs Hivatal

Award Identifier / Grant number: BME NC TKP2020

Award Identifier / Grant number: SNN125119

Funding source: Eötvös Loránd Tudományegyetem

Award Identifier / Grant number: 1783-3/2018/FEKUTSRAT

Funding statement: This research has been supported by the BME NC TKP2020 grant of NKFIH Hungary and also carried out in the ELTE Institutional Excellence Program (1783-3/2018/FEKUTSRAT) supported by the Hungarian Ministry of Human Capacities, and further, it was supported by the Hungarian Scientific Research Fund OTKA SNN125119.

[1] R. A. Adams and J. J. F. Fournier, Sobolev Spaces, 2nd ed., Pure Appl. Math. (Amsterdam) 140, Elsevier/Academic, Amsterdam, 2003. Search in Google Scholar

[2] K. T. Andrews, Y. Dumont, M. F. M’Bengue, J. Purcell and M. Shillor, Analysis and simulations of a nonlinear elastic dynamic beam, Z. Angew. Math. Phys. 63 (2012), no. 6, 1005–1019. 10.1007/s00033-012-0233-9Search in Google Scholar

[3] O. Axelsson and S. Margenov, On multilevel preconditioners which are optimal with respect to both problem and discretization parameters, Comput. Methods Appl. Math. 3 (2003), 6–22. 10.2478/cmam-2003-0002Search in Google Scholar

[4] B. Borsos and J. Karátson, Variable preconditioning for strongly nonlinear elliptic problems, J. Comput. Appl. Math. 350 (2019), 155–164. 10.1016/j.cam.2018.10.004Search in Google Scholar

[5] B. Borsos and J. Karátson, Quasi-Newton variable preconditioning for non-uniformly monotone elliptic problems posed in Banach spaces, IMA J. Numer. Anal. (2021), 10.1093/imanum/drab024. 10.1093/imanum/drab024Search in Google Scholar

[6] P. Deuflhard and M. Weiser, Global inexact Newton multilevel FEM for nonlinear elliptic problems, Multigrid Methods V, Lect. Notes Comput. Sci. Eng. 3, Springer, Berlin (1998), 71–89. 10.1007/978-3-642-58734-4_4Search in Google Scholar

[7] D. A. Dillard, B. Mukherjee, P. Karnal, R. C. Batra and J. Frechette, A review of Winkler’s foundation and its profound influence on adhesion and soft matter applications, Soft Matter 14 (2018), 3669–3683. 10.1039/C7SM02062GSearch in Google Scholar

[8] I. Faragó and J. Karátson, Numerical Solution of Nonlinear Elliptic Problems via Preconditioning Operators: Theory and Applications, Adv. Comput. 11, Nova Science, Hauppauge, 2002. Search in Google Scholar

[9] D. Y. Gao, Nonlinear elastic beam theory with application in contact problems and variational approaches, Mech. Res. Comm. 23 (1996), no. 1, 11–17. 10.1016/0093-6413(95)00071-2Search in Google Scholar

[10] D. Y. Gao, Finite deformation beam models and triality theory in dynamical post-buckling analysis, Int. J. Non-Linear Mech. 35 (2000), no. 1, 103–131. 10.1016/S0020-7462(98)00091-2Search in Google Scholar

[11] D. Y. Gao, J. Machalová and H. Netuka, Mixed finite element solutions to contact problems of nonlinear Gao beam on elastic foundation, Nonlinear Anal. Real World Appl. 22 (2015), 537–550. 10.1016/j.nonrwa.2014.09.012Search in Google Scholar

[12] J. Haslinger, I. Hlaváček and J. Nečas, Numerical methods for unilateral problems in solid mechanics, Handbook of Numerical Analysis. Vol. IV, North-Holland, Amsterdam (1996), 313–485. 10.1016/S1570-8659(96)80005-6Search in Google Scholar

[13] J. Karátson and I. Faragó, Variable preconditioning via quasi-Newton methods for nonlinear problems in Hilbert space, SIAM J. Numer. Anal. 41 (2003), no. 4, 1242–1262. 10.1137/S0036142901384277Search in Google Scholar

[14] J. Machalová and H. Netuka, Solution of contact problems for Gao beam and elastic foundation, Math. Mech. Solids 23 (2018), no. 3, 473–488. 10.1177/1081286517732382Search in Google Scholar

[15] A. H. Nayfeh and P. F. Pai, Linear and Nonlinear Structural Mechanics, John Wiley & Sons, Hoboken, 2004. 10.1002/9783527617562Search in Google Scholar

[16] J. N. Reddy, An Introduction to Nonlinear Finite Element Analysis, Oxford University, Oxford, 2004. 10.1093/acprof:oso/9780198525295.001.0001Search in Google Scholar

[17] S. Stoykov, C. Hofreither and S. Margenov, Isogeometric analysis for nonlinear dynamics of Timoshenko beams, Numerical Methods and Applications, Lecture Notes in Comput. Sci. 8962, Springer, Cham (2015), 138–146. 10.1007/978-3-319-15585-2_16Search in Google Scholar

[18] G. Strang, Computational Science and Engineering, Wellesley-Cambridge, Wellesley, 2007. Search in Google Scholar

[19] C. M. Wang, J. N. Reddy and K. H. Lee, Shear Deformable Beams and Plates, Elsevier Science, Oxford, 2000. Search in Google Scholar

Received: 2020-08-25

Revised: 2021-03-18

Accepted: 2021-05-26

Published Online: 2021-06-09

Published in Print: 2022-01-01

This work is licensed under the Creative Commons Attribution 4.0 International License.

E = E 1 , ν = ν 1

E = E 2 , ν = ν 2