Robust Iterative Solvers for Gao Type Nonlinear Beam Models in Elasticity

: 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.


Introduction
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 Du IV = q and D∆ 2 u = q, (1.1) where D = EI 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.

Nonlinear Gao Beam Models
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: with the following constants: E > 0 is the elastic modulus, I > 0 is the moment of inertia for the beam's crosssection, α = 3h(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 where P is the axial force, assumed to be below the Euler critical load P E cr , and μ = (1 + ν)(1 − ν 2 ).

Some Iterative Methods
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] 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.
The formulation of the third theorem uses the energy * -norm ‖v‖ * := ⟨F (u * ) −1 v, v⟩ 1/2 . 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 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 ).
Remark 2.1. (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.

Estimates in Sobolev spaces
Here let Ω ⊂ R 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.

Proposition 2.1 (Generalized Hölder Inequality). If the numbers p
Then for some constant C p > 0 independent of v.

Weak Formulation and Finite Elements
We rewrite (2.1) by dividing with EI and letting further, we impose clamped boundary conditions. Then our problem becomes The weak formulation uses the Sobolev space H 2 (J); moreover, using the boundary conditions, we work in the space , which has the standard inner product and corresponding norm Then the weak solution of problem (3.1) is defined as The well-posedness of the problem will be discussed in the next section. The weak solution minimizes the energy (3.4) The formulation and properties of the problem will be given for a general FEM subspace V h ⊂ H 2 0 (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
in the space V h (also equipped with the H 2 0 inner product). Then, using standard calculations (see, e.g., [8, Chapter 6.1]), one can derive that F has a Gâteaux derivative satisfying 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.

Properties of the Linearized Operator
The following ellipticity properties hold.

Proposition 3.1. There exists a continuous increasing function
Proof. Owing to (3.7), the quadratic form reads as Since β, k ≥ 0, the first inequality of (3.8) holds.
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. One can also establish local Lipschitz continuity.
Proof. 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, As in the proof of Proposition 3.1 above, we use Proposition 2.1 and the fact that Since ‖z ‖ L 2 = ‖z‖ H 2 0 , this readily entails Repeating the technique used in (3.9) gives Combining (3.10)-(3.11) yields

Finite Elements Using Splines
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 = hk.

The Iterative Solvers: Construction, Convergence and Mesh-Independence
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.
Construction 3.4. The three recurrences can be summarized as follows: 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, Here and henceforth, let us define the right-hand side as a functional ℓ n , that is, by (3.5), Then the z n ∈ V h is the solution of the auxiliary FEM problem • Newton's method: F (u n )z n = F(u n ), σ := 1.
To determine z n , let us again involve test functions, Using (3.7) and (3.13), z n ∈ V h is the solution of the auxiliary FEM problem b ∫ 0 (z n v + 3β(u n ) 2 z n v + kz 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 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 That is, by (3.13), z n ∈ V h is the solution of the auxiliary FEM problem It is informative to summarize also the strong versions of the auxiliary problems, which are formal equations involving fourth derivatives. Namely, let us denote It is readily seen by integration that (for smooth u n ) 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 IV n = r n , • Newton's method: z IV n − 3β((u n ) 2 z n ) + kz n = r n , • quasi-Newton/variable preconditioning: z IV n − w n z n + kz 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, 0 ≤ w n ≤ 3β max(u n ) 2 . (

3.15)
A convenient choice is the arithmetic mean (3.15). Then

Proposition 3.3. For given u n ∈ V h , let the constant w n satisfy
Proof. Firstly, to obtain the upper bound, by (3.7), one can write on the other hand, owing to (3.15), we have 1 ≤ M n := 3β max Ω (u n ) 2 w n for each n ∈ ℕ; hence To obtain the lower bound, by (3.14) and Proposition 2.2, we have writing this back to integral form and adding the term 3β(u n ) 2 (h ) 2 ≥ 0 yields This readily entails ⟨B n h, h⟩ H 2 0 ≤ (1 + w n C 2 2 )⟨F (u n )h, h⟩ H 2 0 . Now we can formulate the convergence results.

Theorem 3.1. The iterative methods, defined in Construction 3.4, provide the following convergence estimates:
• simple iteration/gradient method: • Newton's method: (b) Global convergence for the Newton and quasi-Newton methods can be achieved via damped versions (see, e.g., [5,8] 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.

Other Boundary Conditions
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 This problem can be posed in the Sobolev space H 2 (J) ∩ H 1 0 (J), which is endowed with the same inner product (3.2) as used in H 2 0 (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).

A Modified Equation
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 E cr , 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.

Extension to Plane Problems
For the 2D plate problem in (1.1), the analogue of equation (3.1) is for a thin plate Ω ⊂ R 2 . The weak solution minimizes the energy in the Sobolev space H 2 0 (Ω). It is easy to see that our 1D results can be readily extended to this situation. The problem is posed in H 2 0 (Ω) endowed with the inner product ⟨u, v⟩ H 2 0 := ∫ Ω D 2 u : D 2 v. The main analytic point is that the required Sobolev embedding H 1 0 (Ω) ⊂ 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 1 0 (Ω). Altogether, the calculations can be repeated to obtain the same convergence results as before.

Numerical Experiments
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 [11,14].
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. 8  3  4  5  3  4  5  80  3  4  5  3  4  5  800  3  4  5  3  4  5  8000  4  4  5  3 4 5 Table 1: Number of iterations using the quasi-Newton method.   3  3  4  3  3  4  80  3  3  4  3  3  4  800  3  3  4  3  3  4  8000  3  3  4  3  4 4 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.
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.      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.
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.

Conclusions
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.