On Newton-Raphson formulation and algorithm for displacement based structural dynamics problem with quadratic damping nonlinearity

Quadratic damping nonlinearity is challenging for displacement based structural dynamics problem as the problem is nonlinear in time derivative of the primitive variable. For such nonlinearity, the formulation of tangent stiffness matrix is not lucid in the literature. Consequently, ambiguity related to kinematics update arises when implementing the time integration-iterative algorithm. In present work, an Euler-Bernoulli beam vibration problem with quadratic damping nonlinearity is addressed as the main source of quadratic damping nonlinearity arises from drag force estimation, which is generally valid only for slender structures. Employing Newton-Raphson formulation, tangent stiffness components associated with quadratic damping nonlinearity requires velocity input for evaluation purpose. For this reason, two mathematically equivalent algorithm structures with different kinematics arrangement are tested. Both algorithm structures result in the same accuracy and convergence characteristic of solution.


Introduction
In finite element analysis of time dependent problems, time integration is of primary importance. Other than satisfying accuracy and numerical stability, the selected time integration scheme has to be practical in term of computational effort. In general, there are implicit and explicit time integration schemes. From a computational perspective, the selection of a time integration scheme depends on the type of problem and analysis interest. When the desired overall time span is extensive, it is natural to adopt a time step that is as large as possible to save computing time but small enough to ensure an accurate solution. Such situation is the case for low speed dynamics such as structural dynamics where the structural response can be characterized by a few low frequency modes with sufficient accuracy [1]. Therefore, implicit time integration is normally preferred over explicit time integration in structural dynamics since a relatively larger time step can be used.
When applying implicit time integration onto nonlinear problem, an iterative procedure is required within each time step [1]. Since convergence is a must for each time step, proper implementation of the nonlinear algorithm is crucial. Misunderstanding of the nonlinear formulation can lead to a theoretically incompatible algorithm. Under certain circumstances, theoretically incompatible algorithm may ironically yield a converged solution. However, in most cases such algorithm will yield an inaccurate solution or give rise to numerical instability. Discussion on the implementation of nonlinear algorithm is vast in the literature of finite element analysis [1,2,3]. The numerical behavior of various nonlinear solvers is not within present work's scope. Rather, the full Newton-Raphson method which is widely practised in finite element analysis, is focused. On the basis of a full Newton-Raphson method, a specific source of nonlinearity arised from quadratic damping, also known as quadratic velocity or hydrodynamic damping, is deemed challenging for a few reasons as explained in the forthcoming paragraph.
For quadratic damping nonlinearity, the nonlinear source arises from the presence of quadratic structural velocity. In other words, the problem is nonlinear in term of the time derivative of structural displacement. In a standard displacement based finite element formulation, the primitive variable is structural displacement. Thus, such nonlinearity appears to be more complicated than the standard case which is nonlinear in term of structural displacement. A representative example of nonlinear quadratic damping exists in the estimation of fluid drag force. In particular, drag component of the Morison equation [4,5] widely used in offshore applications contains the quadratic structural velocity. Many researchers in the past had presented the solution for structural analysis involving Morison equation, to name a few [6,7,8,9]. Owing to the diverse focuses and objectives, details on the contribution of nonlinear quadratic damping in the Morison equation towards the formulation of tangent stiffness matrix are often omitted in these works. In ref. [10], the iterative procedure was neglected within each time step when dealing with quadratic damping nonlinearity. To address the lack of lucidity in the Newton-Raphson formulation, the main objective in present work is to show how the quadratic damping nonlinearity contributes to the tangent stiffness matrix. As quadratic damping nonlinearity involves time derivative of displacement, the time integration-iterative algorithm appears to be challenging. In particular, a question arises such that should acceleration and/ or velocity be iterated simultaneously with displacement? To the best knowledge of the authors, such confusion involving quadratic damping nonlinearity has not been addressed before. Therefore, in present work the stated question about nonlinear algorithm is addressed.
The finite element discretization involved in present work is based on a semi-discrete formulation. For structural dynamics, following spatial discretization the semi-discrete form is of a set of hyperbolic ordinary differential equations [11]. A variety of implicit time integration schemes are available for temporal discretization and their implementations had been investigated rigorously in the past, to name a few [12,13,14,15]. Similarly, it is not within the scope of the present work to compare and contrast the various time integration schemes. Rather, the trapezoidal rule or Newmark's constant average acceleration method [1] which is unconditionally stable in linear analysis is adopted. However, the nonlinear formulation can be analogously applied for other implicit time integration schemes.
As mentioned in the third paragraph of section 1, the lack of lucidity in the Newton-Raphson formulation may be due to the diverse focuses and objectives of works in the past. Thus, the authors learn that the best approach to achieve clarity is to present a specific structural dynamics problem where quadratic damping nonlinearity is involved. A suitable choice is the Euler-Bernoulli beam vibration problem because in general, slender structures are prone to bending thus can be modelled as Euler-Bernoulli beam. Meanwhile, drag force estimation which exhibits quadratic damping nonlinearity is normally valid only for slender structures [5]. Therefore, in present work finite element formulation for the Euler-Bernoulli beam vibration problem with quadratic damping nonlinearity is illustrated.
The present work is organized as such. In section 2, a specific Euler-Bernoulli beam vibration problem involving quadratic damping is described. In section 3-4, spatial discretization using Galerkin weighted residual method is performed followed by temporal discretization using trapezoidal rule. In section 5, the Newton-Raphson formulation is given with emphasis on the contribution of quadratic damping nonlinearity towards formulation of the tangent stiffness matrix. Following that, in section 6, discussion of the nonlinear algorithm is given such that the stated question about iterating acceleration and/ or velocity is addressed. Lastly, in section 7, the finite element solution is verified with an established semianalytical solution, while the convergence characteristic of the nonlinear algorithm is discussed. The presented example is taken from Patel and Park [16].

Euler-Bernoulli beam vibration with quadratic damping nonlinearity
In their work, a marine tether subjected to a harmonic support excitation (forced excitation) and a harmonic tension fluctuation (parametric excitation) is idealized as a 1D Euler-Bernoulli beam with constant axial effect as shown in Fig. 1. The quadratic damping nonlinearity arises from the hydrodynamic damping. Patel and Park [16] had solved the problem using a semi-analytical approach. Herein, the same problem is solved using finite element method. The chosen problem is deemed representative because the only source of nonlinearity comes from the hydrodynamic damping. The partial differential equation is given as where ‫ݕ‬ is the displacement variable, a function of 1D space ‫ݔ‬ and time ‫,ݐ‬ ݉ is the mass per unit length of beam, ‫ܫܧ‬ is the flexural rigidity of beam, ܶ is the constant axial tension, ܵ and ߱ are the amplitude and angular frequency of tension fluctuation respectively, ‫ܥ‬ is the drag coefficient, ߩ ௪ is the external fluid density, ‫ܦ‬ is the outer diameter of an annular cross-section and ‫ܮ‬ ் is the total beam length. Similar to the original work [16], only a pinnedpinned boundary condition is studied. Thus, the Dirichlet boundary conditions and Neumann boundary conditions are given respectively as where ‫ݕ‬ and ߱ are the amplitude and angular frequency of support excitation respectively. Due to the nonhomogeneous Dirichlet boundary conditions, the initial displacement and initial velocity of the beam are not necessarily zero. For example, in this particular problem the initial displacement is zero but the initial velocity is not. For the finite element solution, the initial displacement and initial velocity are calculated based on the semi-analytical solution [16] at ‫ݐ‬ = 0.
domain, eq. (1) is weighted using the Bubnov-Galerkin method. Following integration by parts and substitution of cubic Hermitian interpolation function for the displacement variable, the elemental semi-discrete form is given in index notation as where ‫ܮ‬ is the element length, ܰ is the shape function for a Hermitian beam element and ‫ݑ‬ is the nodal degree of freedom. The Neumann boundary conditions at the elemental level can be identified more clearly by rewriting eq. (4) as eq. (5). where Equivalently, eq. (5) can be written in matrix form as where ‫]ܯ[‬ is the mass matrix, ‫]ܭ[‬ is the stiffness matrix, ‫})̇ݑ(ܨ{‬ is the hydrodynamic damping force vector which is nonlinear in velocity ‫,̇ݑ‬ {ℎ} is the Neumann boundary load vector, ‫}ݑ{‬ is the nodal displacement vector and ‫}̈ݑ{‬ is the nodal acceleration vector. Components for the matrices and vectors are given as In this particular problem, the Neumann boundary load vector {ℎ} in eq. (7) vanishes as the forcing effect is due to the nonhomogeneous Dirichlet boundary condition as stated in eq. (2) and the harmonic tension fluctuation. Hence {ℎ} is omitted hereinafter. Note that although the stiffness matrix ‫]ܭ[‬ is time-dependent, the same temporal discretization and Newton-Raphson formulation hold as in the case of a zero tension fluctuation i.e. when ܵ = 0.

4
Temporal discretization using trapezoidal rule To proceed, the equation of motion in eq. (7) is expressed at time ‫ݐ‬ + ‫ݐ∆‬ which is to be solved as Eq. (9) signifies that equilibrium condition should hold at any time and is mathematically equivalent to the notion that the external load vector is equal to the internal load vector at any time [1]. For temporal discretization using trapezoidal rule, the acceleration and velocity at time ‫ݐ‬ + ‫ݐ∆‬ are given implicitly with respect to time ‫ݐ‬ + ‫ݐ∆‬ and previous time ‫ݐ‬ as where ‫ݐ∆‬ is the time step used. Substituting eq. (10) and eq. (11) to eq. (9) and rearranging, one obtains Apparently, eq. (12) is nonlinear in velocity hence an iterative procedure is required within each time step.

Newton-Raphson formulation
All terms in eq. (9) or equivalently eq. (12) can be collected on one side such that the residual function desired at present iteration ‫ݎ‬ is defined as ‫})̇ݑ(ܴ{‬ ௧ା∆௧ = 0 (13) In Newton-Raphson formulation, the residual function at present iteration ‫ݎ‬ is expanded using Taylor series with respect to previous iteration ‫ݎ‬ − 1 . Terms involving second order function and above are truncated.
Thus, the tangent stiffness equations are obtained as where [ܶ] ௧ା∆௧ ିଵ is the tangent stiffness matrix and ିଵ . Both tangent stiffness matrix and residual function are expressed at the known iteration ‫ݎ‬ − 1 hence eq. (15) can be readily solved in an iterative manner.
In particular, components of the tangent stiffness matrix can be evaluated as Utilizing relations given in eq. (10) and eq. (11), eq. (16) becomes The third term on the right hand side of eq. (17) is elaborated as By differentiating using product rule in eq. (18), it can be shown that Accordingly, components of the tangent stiffness matrix can be finalized as The spatial discretization, temporal discretization and Newton-Raphson formulation shown above are operated on the elemental level. However, this is immaterial for the global level operation as global matrices and vectors can be constructed using direct stiffness method.

Implementation of nonlinear algorithm
In solving nonlinear time dependent problems using implicit time integration, there are essentially two loops in the algorithm structure. The outer loop is for time integration whereas the inner loop is for Newton-Raphson iteration. Since input of velocity is required for evaluation of the nonlinear hydrodynamic damping force vector, an important question herein is that should the acceleration and/ or velocity be updated in the outer loop or the inner loop? The confusion may be further demonstrated through literature for the standard case where the problem is nonlinear in displacement. For instance, update in the outer loop is seen in section 16.3.3 of ref. [17] whereas update in the inner loop is seen in section 6.3.7 of ref. [18] and ref. [19].
In present work, both algorithm structures had been tested for quadratic damping nonlinearity. To clarify, this section shows that both algorithm structures are mathematically equivalent thus either algorithm structure can be implemented. However, when dealing with quadratic damping nonlinearity the velocity should be iterated in the inner loop for either algorithm structure to ensure convergence. Otherwise, the nonlinear hydrodynamic damping vector may be inappropriately evaluated thus producing divergence of solution.
In eq. (15), both the tangent stiffness matrix as given in eq. (9) or eq. (12) and residual function as given in eq. (22) are expressed at the known iteration ‫ݎ‬ − 1 . If the next iteration is required, the matrices and vectors are naturally expressed at the succeeding iteration. For the residual function, the general picture is clearer by rewriting eq. (9) and eq. (12) respectively with the iteration superscript included as Accordingly, for the standard case where the problem is nonlinear in displacement, the following statements shall apply. If the acceleration and velocity were updated in the outer loop, this indicates that the perspective is set such that eq. (22) and eq. (24) are meant. On the other hand, the perspective for the case where acceleration and velocity are updated in the inner loop is set such that eq. (22) and eq. (23) are meant. The latter can be implemented as following.
By substituting ‫}ݑ{‬ ௧ା∆௧ = ‫}ݑ{‬ ௧ା∆௧ ିଵ + ‫}ݑ∆{‬ ௧ା∆௧ into eq. (10) and eq. (11) at the desired iteration ‫ݎ‬ [1], the acceleration and velocity to be iterated in the inner loop are given respectively as As shown in eq. (25) and eq. (26), once ‫}ݑ∆{‬ ௧ା∆௧ is solved, acceleration and velocity can be updated in the inner loop. Note that in this particular problem however, the velocity should be iterated in the inner loop regardless of the algorithm structure used. For this particular problem, both algorithm structures are given in Table 1 and Table 2.

IN-3
Evaluate tangent stiffness matrix and residual function using eq. (22) and eq. (24) respectively and assemble into global matrices

IN-4
Impose stationary condition on ‫}ݑ∆{‬ ௧ା∆௧ * for all known DOFs associated with Dirichlet boundary conditions

IN-6
Check if tolerance is acceptable while iterating: x If pass: go to OUT-4 x If fail: set current iteration as ‫ݎ‬ − 1 and perform

OUT-4
Calculate acceleration at time ‫ݐ‬ + ‫ݐ∆‬

OUT-6
Set current time as ‫ݐ‬ and go to OUT-1 * For any iteration in the same time instant ‫ݐ‬ + ‫,ݐ∆‬ the known DOFs associated with the nonhomogeneous Dirichlet boundary conditions do not change.

OUT-3
Initialization for kinematics at time ‫ݐ‬ + ‫ݐ∆‬ to be iterated in the inner loop

IN-2
Impose stationary condition on ‫}ݑ∆{‬ ௧ା∆௧ * for all known DOFs associated with Dirichlet boundary conditions

IN-4
Check if tolerance is acceptable while iterating: x If pass: go to OUT-3 x If fail: set current iteration as ‫ݎ‬ − 1 and perform Set current time as ‫ݐ‬ and go to OUT-1 * For any iteration in the same time instant ‫ݐ‬ + ‫,ݐ∆‬ the known DOFs associated with the nonhomogeneous Dirichlet boundary conditions do not change.

Numerical example
To verify the nonlinear formulation, both finite element solution given in present work and semi-analytical solution given by Patel and Park [16] are computed for comparison. In the semi-analytical solution, the solution is expressed in the form of Fourier sine series. To obtain the unknown Fourier coefficients, eq. (1) is reduced to a system of 2 nd order nonlinear ordinary differential equations using Galerkin mode decomposition method, which is subsequently solved using fourth order Runge-Kutta (RK4) method. With respect to Fig. 1, the problem parameters are given in Table 3. Similar to the original work [16], three cases are studied herein to illustrate the effect of different forcing excitations. Likewise, the solution parameters are given in Table 4. In the finite element solution, a rather strict tolerance criterion has been set as to test the feasibility of both algorithm structures in handling quadratic damping nonlinearity.  [16] at ‫ݐ‬ = 0 B. Semi-analytical solution (for notation see [16]) Initial conditions for the system of 2 nd order ODEs The displacement time histories of the beam at midlength for all three cases are given in Fig. 2-4. As shown in the figures, both algorithm structures agree well with the semi-analytical solution. To further ascertain the feasibility of both algorithm structures, the number of iterations required and tolerance upon converged are analysed. As shown in Fig. 5-6, both algorithm structures used are effective in the numerical example with only three to four iterations required in each time step.  Fig. 7-9. The similarity of convergence properties in both algorithm structures is due to mathematical manipulations as explained earlier.

Conclusion
In present work, nonlinear formulation for an Euler-Bernoulli beam vibration problem with quadratic damping nonlinearity is demonstrated. The finite element solution shows a good agreement with the established semi-analytical solution by Patel and Park [16]. Through Newton-Raphson formulation, it is lucid that the tangent stiffness components associated with the quadratic damping nonlinearity require input of velocity for evaluation purpose in the iterative procedure. As a result, velocity has to be iterated in the inner (iterative) loop when implementing the time integration-iterative algorithm. The nonlinear formulation has been tested with two mathematically equivalent algorithm structures with different kinematics arrangement. The first algorithm structure which is termed "update in the outer loop", involves iteration of displacement and velocity, whereas the second algorithm which is termed "update in the inner loop", involves iteration of displacement, velocity and acceleration. As expected, both mathematically equivalent algorithm structures exhibit the same accuracy and convergence characteristic of solution.