A New Finite Element Formulation for Nonlinear Vibration Analysis of the Hard-Coating Cylindrical Shell

Abstract

In this paper, a four-node composite cylindrical shell finite element model based on Love’s first approximation theory is proposed to solve the nonlinear vibration of the hard-coating cylindrical shell efficiently. The developed model may have great significance for vibration reduction of the cylindrical shell structures of the aero engine or aircraft. The influence of the strain dependence of the coating material on the complex stiffness matrix is considered in this model. Nonlinear iterative solution formulas with a unified iterative method are theoretically derived for solving the resonant frequency and response of the composite cylindrical shell. Then, a cylindrical shell coated with a thin layer of NiCoCrAlY + yttria-stabilized zirconia (YSZ) is chosen to demonstrate the proposed formulation, and the rationality is validated by comparing with the finite element iteration method (FEIM). Results show that the developed finite element method is more efficient, and the hard-coating cylindrical shell has the characteristics of soft nonlinearity due to the strain dependence of the coating material.

1. Introduction

The cylindrical shell has been widely used in the aviation and aerospace fields, especially in the aero engine and aircraft, which have an urgent demand to restrain their excessive vibration [1]. It is well known that the hard-coating treatment has found application in many engineering applications, such as the blades of gas turbines, new turbofans, cylindrical shell structures of the aero engine and other structures. The hard-coating plays an important role in improving the performance of the thermal barrier, anti-friction and anti-erosion of composite structures [2,3,4,5]. Recently, both Ivancic and Palazotto [6] and Patsias et al. [7] found that the hard-coating with special microstructural features can remarkably enhance the additional damping for vibration reduction of structures. Results also revealed that the storage modulus and the loss factor of the coating material change regularly with the strain amplitude, called the strain dependence, which makes the nonlinear vibration analysis of the hard-coating structures become more challenging.

To clearly explain the strain dependence of hard-coating, researchers have conducted many experiments. Blackwell et al. [8] experimentally studied the damping effects of hard-coating coated on the titanium plate and found its significant nonlinear damping characterization. Several research works, e.g., [9,10], indicated that both the elastic and damping properties of hard-coatings exhibit obvious nonlinearity, as well. Furthermore, Torvik [11] successfully extracted the mechanical parameters of hard-coating, such as the storage modulus and the loss factor, which are dependent on the equivalent strain amplitude and mainly responsible for the nonlinearity of the hard-coating structures. The above indicates that the hard-coating has nonlinear characteristics due to the strain dependence, which makes the establishment of an analysis model of the hard-coating structures become very difficult and time consuming.

To better apply the hard-coating for vibration reduction and solve the vibration problems of the hard-coating structures, several analysis models have been developed. Sun and Liu [12] developed an analytical model for the resonant frequencies and response analysis of the hard-coating beam with considering the strain dependence of the coating material. Chen et al. [13] proposed an analytical approach to calculate the free vibration characteristics and damping effect of the hard-coating on blades based on the constitutive model of the complex modulus and Rayleigh–Ritz method. Sun et al. [14] built an analytical model to solve the free vibration of the hard-coating cantilever cylindrical shell by Love’s first approximation theory and the Rayleigh–Ritz method, but did not consider the strain dependence of the hard-coating. In addition to the analytical method, the finite element method with practicality and effectiveness has been used to solve the nonlinear vibration of the hard-coating structures. Filippi and Torvik [15] proposed a finite element method to solve the nonlinear vibration response of the hard-coating blade. However, it is difficult to exactly determine the nonlinear characteristics by using only the linear analytic expressions. Then, Li et al. [16] developed a finite element iterative method (FEIM), which successfully solved the resonant frequencies and responses of the hard-coating plate by satisfactorily characterizing the nonlinear storage modulus and loss modulus of hard-coating with the polynomial method. The practicability and reliability of the FEIM have been verified through comparing with the experimental results. The above research works provide an important support for the nonlinear vibration analysis of the hard-coating structures, such as beams, plates, blades and shells. However, it is quite inadequate for the nonlinear vibration analysis of the hard-coating cylindrical shell implemented by the finite element method with considering the strain dependence of the coating material. Moreover, the nonlinear vibration analysis of the hard-coating cylindrical shell may have great significance for vibration reduction of the cylindrical shell structures of the aero engine or aircraft.

Generally, the available FEIM needs to apply the commercial FEM software to obtain the element matrices of stiffness, mass and damping, and accurate results can only be obtained when with very small meshes, which makes the FEIM inefficient. Moreover, two different iteration methods have been used to solve the resonant frequencies and responses, respectively, leading to more inefficiency. The problems become even serious when solving the nonlinear vibration characteristics of the hard-coating cylindrical shell, which have more degrees of freedom. Unfortunately, no applicable elements exist to solve the problems more efficiently. The vast majority of elements are only used to deal with geometric nonlinearity [17,18,19,20]. Although Ferreira et al. [21] and Masti and Sainsbury [22] presented the composite cylindrical shell elements, which can describe the nonlinear behaviors of the hyperelastic coating and the viscoelastic coating, respectively, the elements have no ability to deal with the strain dependence of hard-coating, especially the strain dependence of the loss factor, which determines the value of the global material damping matrix.

Thus, this paper presents a new finite element formulation and a unified iterative method to solve the nonlinear vibration of the hard-coating cylindrical shell with considering the strain dependence of the coating material more efficiently. In Section 2, the geometry and element of the hard-coating cylindrical shell and Love’s first approximation theory are briefly introduced, and then, the finite element formulation is deduced in detail. Meanwhile, the strain-dependent storage modulus and loss factor of hard-coating are characterized by the polynomial method in Section 3, and the solution formulas of the nonlinear vibration of the hard-coating cylindrical shell in Section 4 are theoretically derived. Finally, in Section 5, a cylindrical shell coated with a thin layer of NiCoCrAlY + yttria-stabilized zirconia (YSZ), which is a common hard-coating material [3], is set as the study subject; both the accuracy and convergence of the developed finite element method are compared with the FEIM, which indicates the rationality of the proposed method. Moreover, the nonlinear vibration analysis of the cylindrical shell coated with a thin layer of NiCoCrAlY + YSZ is implemented.

2. Derivation of the Finite Element Formulation for the Hard-Coating Cylindrical Shell

2.1. Geometry and Element of the Hard-Coating Cylindrical Shell

The hard-coating cylindrical shell with thickness h, axial length L and mean radius R is shown in Figure 1a. An orthogonal curvilinear coordinate system is established at the middle surface of the hard-coating cylindrical shell along the x (axial), θ (circumferential) and z (radial) directions, respectively. The finite element of the hard-coating cylindrical shell is shown in Figure 1b, which has four nodes and nine degrees of freedom at each node.

2.2. Love’s First Approximation Theory

Based on Love’s first approximation theory [23], the displacement fields u, v and w of the hard-coating cylindrical shell can be represented by the following relationships:

$$\begin{array}{l}u\left(x,\mathsf{\theta },z\right)=u^0\left(x,\mathsf{\theta }\right)+z{\mathsf{\psi }}_x\left(x,\mathsf{\theta }\right)\\ v\left(x,\mathsf{\theta },z\right)=v^0\left(x,\mathsf{\theta }\right)+z{\mathsf{\psi }}_{\mathsf{\theta }}\left(x,\mathsf{\theta }\right)\\ w\left(x,\mathsf{\theta },z\right)=w^0\left(x,\mathsf{\theta }\right)\end{array}$$

(1)

where u⁰, v⁰ and w⁰ are the orthogonal components of displacement of the mid-surface in the x, θ and z directions, respectively, and ψ_x and ψ_θ are the rotations of the mid-surface about the θ and x axes, respectively.

The magnitude of an arbitrary infinitesimal changing in the mid-surface of the hard-coating cylindrical shell ds is determined by:

$${\left(\mathrm{d}s\right)}^2=g_x\mathrm{d}{x}^2+g_{\mathsf{\theta }}\mathrm{d}{\mathsf{\theta }}^2+g_z\mathrm{d}{z}^2$$

(2)

where g_x, g_θ, g_z, are Lamb constants [24]. For thin cylindrical shell (h/R$\ll 1$):

$$g_x=1,\hspace{1em}{g}_{\mathsf{\theta }}=R,\hspace{1em}{g}_z=1$$

(3)

Thus, ${\mathsf{\epsilon }}_x^0$, ${\mathsf{\epsilon }}_{\mathsf{\theta }}^0$, ${\mathsf{\epsilon }}_{x\mathsf{\theta }}^0$ and κ_x, κ_θ, κ_xθ, the strains and bending strains of the middle surface, can be given as follows [25]:

$$\begin{array}{l}{\epsilon }_x^0=\frac{\partial u^0}{\partial x},\hspace{1em}{\epsilon }_{\mathsf{\theta }}^0=\frac{1}{R}\frac{\partial v^0}{\partial \mathsf{\theta }}+\frac{w^0}{R},\hspace{1em}{\epsilon }_{x\mathsf{\theta }}^0=\frac{\partial v^0}{\partial x}+\frac{1}{R}\frac{\partial u^0}{\partial \mathsf{\theta }}\\ {\kappa }_x=\frac{{\partial }^2{w}^0}{\partial x^2},\hspace{1em}{\kappa }_{\mathsf{\theta }}=\frac{1}{R^2}\left(\frac{{\partial }^2{w}^0}{\partial {\mathsf{\theta }}^2}-\frac{\partial v^0}{\partial \mathsf{\theta }}\right),\hspace{1em}{\kappa }_{x\mathsf{\theta }}=\frac{2}{R}\left(\frac{{\partial }^2{w}^0}{\partial x\partial \mathsf{\theta }}-\frac{\partial v^0}{\partial x}\right)\end{array}$$

(4)

The stress-strain equations of the hard-coating cylindrical shell can be expressed as:

$$\mathsf{\sigma }=\left\{\begin{array}{c}{\mathsf{\sigma }}_x\\ {\mathsf{\sigma }}_{\theta }\\ {\mathsf{\sigma }}_{x\theta }\end{array}\right\}=\left[\begin{array}{ccc}{Q}_{11}^k& Q_{12}^k& 0\\ Q_{21}^k& Q_{22}^k& 0\\ 0& 0& Q_{66}^k\end{array}\right]\left\{\begin{array}{c}{\mathsf{\epsilon }}_x\\ {\mathsf{\epsilon }}_{\mathsf{\theta }}\\ {\mathsf{\epsilon }}_{x\mathsf{\theta }}\end{array}\right\}$$

(5)

where ${\mathrm{Q}}_{ij}^k$ are the reduced stiffness coefficients of the k-th lamina, given by:

$$Q_{11}^k=Q_{22}^k=\frac{E_k}{1-v_k^2},\hspace{1em}{Q}_{12}^k=Q_{21}^k=\frac{E_k{v}_k}{1-v_k^2},\hspace{1em}{Q}_{66}^k=\frac{E_k}{2\left(1+v_k\right)}$$

(6)

The material parameters E_k and v_k are the Young’s modulus and Poisson’s ratio of the k-th lamina, respectively. For a double-layered cylindrical shell, the force and moment resultants are defined by:

$$\left\{\begin{array}{c}{N}_x\\ N_{\mathsf{\theta }}\\ N_{x\mathsf{\theta }}\\ M_x\\ M_{\mathsf{\theta }}\\ M_{x\mathsf{\theta }}\end{array}\right\}={\displaystyle \sum _{k=1}^2{\displaystyle {\int }_{z_{k-1}}^{z_k}\left\{\begin{array}{c}{\mathsf{\sigma }}_x\\ {\mathsf{\sigma }}_{\mathsf{\theta }}\\ {\mathsf{\sigma }}_{x\mathsf{\theta }}\\ z{\mathsf{\sigma }}_x\\ z{\mathsf{\sigma }}_{\mathsf{\theta }}\\ z{\mathsf{\sigma }}_{x\mathsf{\theta }}\end{array}\right\}}}\mathrm{d}z=\mathit{S}\left\{\begin{array}{c}{\mathsf{\epsilon }}_x^0\\ {\mathsf{\epsilon }}_{\mathsf{\theta }}^0\\ {\mathsf{\epsilon }}_{x\mathsf{\theta }}^0\\ {\mathsf{\kappa }}_x\\ {\mathsf{\kappa }}_{\mathsf{\theta }}\\ {\mathsf{\kappa }}_{x\mathsf{\theta }}\end{array}\right\}$$

(7)

where:

$$\mathit{S}=\left[\begin{array}{cccccc}{A}_{11}& A_{12}& 0& B_{11}& B_{12}& 0\\ A_{21}& A_{22}& 0& B_{21}& B_{22}& 0\\ 0& 0& A_{66}& 0& 0& B_{66}\\ B_{11}& B_{12}& 0& D_{11}& D_{12}& 0\\ B_{21}& B_{22}& 0& D_{21}& D_{22}& 0\\ 0& 0& B_{66}& 0& 0& D_{66}\end{array}\right]$$

(8)

in which:

$$\begin{array}{l}{A}_{ij}={\displaystyle \sum _{k=1}^2{Q}_{ij}^k\left(z_k-z_{k-1}\right)}\\ B_{ij}=\frac{1}{2}{\displaystyle \sum _{k=1}^2{Q}_{ij}^k\left(z_k^2-z_{k-1}^2\right)}\\ D_{ij}=\frac{1}{3}{\displaystyle \sum _{k=1}^2{Q}_{ij}^k\left(z_k^3-z_{k-1}^3\right)}\end{array}$$

(9)

The three coordinates z₀, z₁ and z₂ are defined as follows:

$$z_0=-T_0, z_1=T_1-T_0, z_2=T_1+T_2-T_0$$

(10)

where T₁ and T₂ are the thickness of metal substrate and hard-coating, respectively (see Figure 1b). T₀ is the distance from the middle surface to the inner surface, defined by:

$$T_0=\frac{E_1{T}_1^2+2E_2{T}_1{T}_2+E_2{T}_2^2}{2\left(E_1{T}_1+E_2{T}_2\right)}$$

(11)

2.3. Finite Element Formulation

The degrees of freedom at each node are defined as follows (see Figure 1b):

$$u^0,u_x^0,u_{\theta }^0,v^0,v_x^0,v_{\theta }^0,w^0,w_x^0,w_{\theta }^0$$

(12)

where:

$$u_x^0=\frac{\partial u^0}{\partial x},u_{\mathsf{\theta }}^0=\frac{\partial u^0}{\partial \mathsf{\theta }},v_x^0=\frac{\partial v^0}{\partial x},v_{\mathsf{\theta }}^0=\frac{\partial v^0}{\partial \mathsf{\theta }},w_x^0=\frac{\partial w^0}{\partial x},w_{\mathsf{\theta }}^0=\frac{\partial w^0}{\partial \mathsf{\theta }}$$

(13)

The displacement field for u⁰ in the coordinate directions x and θ is assumed to be of the form:

$$\begin{array}{ll}{u}^0=& a_1+a_{2}x+a_3\mathsf{\theta }+a_4{x}^2+a_{5}x\mathsf{\theta }+a_6{\mathsf{\theta }}^2+a_7{x}^3+a_8{x}^2\mathsf{\theta }+a_{9}x{\mathsf{\theta }}^2+a_{10}{\mathsf{\theta }}^3\\ & +a_{11}{x}^3\mathsf{\theta }+a_{12}x{\mathsf{\theta }}^3\end{array}$$

(14)

The same form of expression is used for the displacements v⁰ and w⁰. The equations of the displacement field for u⁰, $u_x^0$ and $u_{\mathsf{\theta }}^0$ can be written in matrix form:

$$U=XA$$

(15)

in which U, X and A can be expressed as follows:

$$\mathit{U}={\left[\begin{array}{cccc}{\overline{\mathit{U}}}_1& {\overline{\mathit{U}}}_2& {\overline{\mathit{U}}}_3& {\overline{\mathit{U}}}_4\end{array}\right]}^{\mathrm{T}}$$

(16)

$$\mathit{X}={\left[\begin{array}{cccc}{\overline{\mathit{X}}}_1& {\overline{\mathit{X}}}_1& {\overline{\mathit{X}}}_1& {\overline{\mathit{X}}}_1\end{array}\right]}^{\mathrm{T}}$$

(17)

$$A={\left(\begin{array}{cccccccccccc}{a}_1& a_2& a_3& a_4& a_5& a_6& a_7& a_8& a_9& a_{10}& a_{11}& a_{12}\end{array}\right)}^{\mathrm{T}}$$

(18)

where:

$${\overline{U}}_i={\left[\begin{array}{ccc}{u}_i^0& u_{xi}^0& u_{\theta i}^0\end{array}\right]}^{\mathrm{T}}$$

(19)

$${\overline{X}}_i=\left[\begin{array}{cccccccccccc}1& x_i& {\mathsf{\theta }}_i& x_i^2& x_i{\mathsf{\theta }}_i& {\mathsf{\theta }}_i^2& x_i^3& x_i^2{\mathsf{\theta }}_i& x_i{\mathsf{\theta }}_i^2& {\mathsf{\theta }}_i^3& x_i^3{\mathsf{\theta }}_i& x_i{\mathsf{\theta }}_i^3\\ 0& 1& 0& 2x_i& {\mathsf{\theta }}_i& 0& 3x_i^2& 2x_i{\mathsf{\theta }}_i& {\mathsf{\theta }}_i^2& 0& 3x_i^2{\mathsf{\theta }}_i& {\mathsf{\theta }}_i^3\\ 0& 0& 1& 0& x_i& 2{\mathsf{\theta }}_i& 0& x_i^2& 2x_i{\mathsf{\theta }}_i& 3{\mathsf{\theta }}_i^2& x_i^3& 3x_i{\mathsf{\theta }}_i^2\end{array}\right]$$

(20)

By applying inverse matrix method, the displacement interpolation functions can be calculated as follows:

$$\left(\begin{array}{cccccccccccc}{N}_1& N_{x1}& N_{\mathsf{\theta }1}& N_2& N_{x2}& N_{\mathsf{\theta }2}& N_3& N_{x3}& N_{\mathsf{\theta }3}& N_4& N_{x4}& N_{\mathsf{\theta }4}\end{array}\right)=X^{-1}U$$

(21)

where $N_i$, $N_{xi}$ and $N_{\theta i}$ are the displacement interpolation functions at node i.

Then, the shape function matrix N can be expressed as:

$$N=\left[\begin{array}{cccc}{\overline{N}}_1& {\overline{N}}_2& {\overline{N}}_3& {\overline{N}}_4\end{array}\right]$$

(22)

where ${\overline{N}}_i$ is the interpolation functions corresponding to node i, given by:

$${\overline{N}}_i=\left[\begin{array}{ccccccccc}{N}_i& N_{xi}& N_{\mathsf{\theta }i}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& N_i& N_{xi}& N_{\mathsf{\theta }i}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& N_i& N_{xi}& N_{\mathsf{\theta }i}\end{array}\right]$$

(23)

The general strain-displacement relations for the hard-coating cylindrical shell can be expressed as:

$$\mathsf{\epsilon }=\left\{\begin{array}{c}{\mathsf{\epsilon }}_x^0\\ {\mathsf{\epsilon }}_{\mathsf{\theta }}^0\\ {\mathsf{\epsilon }}_{x\mathsf{\theta }}^0\\ {\mathsf{\kappa }}_x\\ {\mathsf{\kappa }}_{\mathsf{\theta }}\\ {\mathsf{\kappa }}_{x\mathsf{\theta }}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial u^0}{\partial x}\\ \frac{1}{R}\frac{\partial v^0}{\partial \mathsf{\theta }}+\frac{w^0}{R}\\ \frac{\partial v^0}{\partial x}+\frac{1}{R}\frac{\partial u^0}{\partial \mathsf{\theta }}\\ \frac{{\partial }^2{w}^0}{\partial x^2}\\ \frac{1}{R^2}\left(\frac{{\partial }^2{w}^0}{\partial {\mathsf{\theta }}^2}-\frac{\partial v^0}{\partial \mathsf{\theta }}\right)\\ \frac{2}{R}\left(\frac{{\partial }^2{w}^0}{\partial x\partial \mathsf{\theta }}-\frac{\partial v^0}{\partial x}\right)\end{array}\right\}=\mathbf{\Gamma }\left\{\begin{array}{c}{u}^0\\ v^0\\ w^0\end{array}\right\}=\mathbf{\Gamma }\mathbf{N}\mathit{d}=\mathbf{B}\mathit{x}$$

(24)

where the matrix differential operator Г, the strain matrix B and the displacement vector x are defined, respectively, by:

$$\mathbf{\Gamma }=\left[\begin{array}{ccc}\frac{\partial }{\partial x}& 0& 0\\ 0& \frac{1}{R}\frac{\partial }{\partial \mathsf{\theta }}& \frac{1}{R}\\ \frac{1}{R}\frac{\partial }{\partial \mathsf{\theta }}& \frac{\partial }{\partial x}& 0\\ 0& 0& \frac{{\partial }^2}{\partial x^2}\\ 0& -\frac{1}{R^2}\frac{\partial }{\partial \mathsf{\theta }}& \frac{1}{R^2}\frac{{\partial }^2}{\partial {\mathsf{\theta }}^2}\\ 0& -\frac{2}{R}\frac{\partial }{\partial x}& \frac{2}{R}\frac{{\partial }^2}{\partial x\partial \mathsf{\theta }}\end{array}\right]$$

(25)

$$\mathbf{B}=\left[\begin{array}{cccc}{\overline{\mathbf{B}}}_1& {\overline{\mathbf{B}}}_2& {\overline{\mathbf{B}}}_3& {\overline{\mathbf{B}}}_4\end{array}\right]$$

(26)

$$\mathit{x}={\left(\begin{array}{cccc}{\overline{\mathit{x}}}_1& {\overline{\mathit{x}}}_2& {\overline{\mathit{x}}}_3& {\overline{\mathit{x}}}_4\end{array}\right)}^{\mathrm{T}}$$

(27)

The matrix ${\overline{\mathbf{B}}}_i$ is:

$${\overline{\mathbf{B}}}_i=\left[\begin{array}{ccccccccc}{N_i|}_x& {N_{xi}|}_x& {N_{\theta i}|}_x& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{{N_i|}_{\theta }}{R}& \frac{{N_{xi}|}_{\theta }}{R}& \frac{{N_{\theta i}|}_{\theta }}{R}& \frac{N_i}{R}& \frac{N_{xi}}{R}& \frac{N_{\theta i}}{R}\\ \frac{{N_i|}_{\theta }}{R}& \frac{{N_{xi}|}_{\theta }}{R}& \frac{{N_{\theta i}|}_{\theta }}{R}& {N_i|}_x& {N_{xi}|}_x& {N_{\theta i}|}_x& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& {N_i|}_x^2& {N_{xi}|}_x^2& {N_{\theta i}|}_x^2\\ 0& 0& 0& \frac{-{N_i|}_{\theta }}{R^2}& \frac{-{N_{xi}|}_{\theta }}{R^2}& \frac{-{N_{\theta i}|}_{\theta }}{R^2}& N_i|& \frac{{N_{xi}|}_{\theta }^2}{R^2}& \frac{{N_{\theta i}|}_{\theta }^2}{R^2}\\ 0& 0& 0& \frac{-2{N_i|}_x}{R}& \frac{-2{N_{xi}|}_x}{R}& \frac{-2{N_{\theta i}|}_x}{R}& \frac{2{N_i|}_{x\theta }^2}{R}& \frac{2{N_{xi}|}_{x\theta }^2}{R}& \frac{2{N_{\theta i}|}_{\theta }^2}{R}\end{array}\right]$$

(28)

where:

$${|}_x=\frac{\partial }{\partial x},\hspace{1em}{|}_x^2=\frac{{\partial }^2}{\partial x^2},\hspace{1em}{|}_{\theta }=\frac{\partial }{\partial \theta },\hspace{1em}{|}_{\theta }^2=\frac{{\partial }^2}{\partial {\theta }^2},\hspace{1em}{|}_{x\theta }^2=\frac{{\partial }^2}{\partial x\partial \theta }$$

(29)

The displacement vector ${\overline{\mathit{x}}}_i$ corresponding to node i is defined as:

$${\overline{\mathit{x}}}_i={\left(\begin{array}{ccccccccc}{u}_i& u_{xi}& u_{\theta i}& v_i& v_{xi}& v_{\theta i}& w_i& w_{xi}& w_{\theta i}\end{array}\right)}^{\mathrm{T}}$$

(30)

Thus, the element stiffness matrix K^el, the element mass matrix M^el and the external load vector F^el introduced by the base excitation can be expressed, respectively, by:

$${\mathbf{K}}^{\mathrm{el}}=\frac{1}{2}{\displaystyle \iiint {\mathbf{B}}^{\mathrm{T}}\mathbf{S}\mathbf{B}\left(R+z\right)\mathrm{d}z\mathrm{d}\mathsf{\theta }\mathrm{d}x}$$

(31)

$${\mathbf{M}}^{\mathrm{el}}=\frac{1}{2}{\displaystyle \iiint \mathsf{\rho }{\mathbf{N}}^{\mathrm{T}}\mathbf{N}\left(R+z\right)\mathrm{d}z\mathrm{d}\mathsf{\theta }\mathrm{d}x}$$

(32)

$${\mathbf{F}}^{\mathrm{el}}=\frac{1}{2}{\displaystyle \iiint {\mathbf{N}}^{\mathrm{T}}\mathit{f}\left(R+z\right)\mathrm{d}z\mathrm{d}\mathsf{\theta }\mathrm{d}x}$$

(33)

where the mass density per unit area ρ is defined by [22]:

$$\mathsf{\rho }={\mathsf{\rho }}_1{T}_1+{\mathsf{\rho }}_2{T}_2$$

(34)

where ρ₁ and ρ₂ are the densities of metal substrate and hard-coating, respectively.

The f is a 3 × 1 order vector composed of the x, θ and z components of the base excitation per unit volume, defined by:

$$\mathit{f}=\mathsf{\rho }a{\left\{{\mathsf{\mu }}_x,{\mathsf{\mu }}_{\mathsf{\theta }},{\mathsf{\mu }}_z\right\}}^{\mathrm{T}}$$

(35)

in which a is the absolute acceleration of the base and μ_x, μ_θ and μ_z are the influence coefficients of the base excitation on the x, θ and z directions of the hard-coating cylindrical shell.

The global stiffness matrix K, the global mass matrix M and the global external load vector F of the hard-coating cylindrical shell can be assembled by K^el, M^el and F^el in the corresponding sequence, respectively.

3. Characterizing the Strain Dependence Using the High Order Polynomial

To describe the strain dependence of the hard-coating cylindrical shell, the complex modulus E₂* of the hard-coating is considered, given by:

$$E_2^{\ast }=E_2+{\mathrm{i}\mathsf{\eta }}_2$$

(36)

where * denotes the complex value, $\mathrm{i}=\sqrt{-1}$, and E₂ (namely the E₂ in Equation (6)) and η₂ are the storage modulus and the loss factor of hard-coating, respectively, which can be expressed by the p-order polynomials,

$$E_2=E_{20}+E_{21}{\mathsf{\epsilon }}_{\mathrm{e}}+E_{22}{\mathsf{\epsilon }}_{\mathrm{e}}^2+\cdot \cdot \cdot +E_{2p}{\mathsf{\epsilon }}_{\mathrm{e}}^p$$

(37)

$${\mathsf{\eta }}_2={\mathsf{\eta }}_{20}+{\mathsf{\eta }}_{21}{\mathsf{\epsilon }}_{\mathrm{e}}+{\mathsf{\eta }}_{22}{\mathsf{\epsilon }}_{\mathrm{e}}^2+\cdot \cdot \cdot +{\mathsf{\eta }}_{2p}{\mathsf{\epsilon }}_{\mathrm{e}}^p$$

(38)

in which E_2j and η_2j (j = 0, 1, …, p) are specific j-order coefficients of the storage modulus and the loss factor, respectively. For convenience, one can define:

$$E_2^{\ast }=\left(E_{20}+\mathrm{i}{\mathsf{\eta }}_{20}\right)+E_{2,\mathrm{inc}}^{\ast }$$

(39)

where $E_{2,\mathrm{inc}}^{*}$ is the increment of the complex modulus caused by the equivalent strain ε_e, given by:

$$E_{2,\mathrm{inc}}^{\ast }=E_{21}^{\ast }{\mathsf{\epsilon }}_{\mathrm{e}}+E_{22}^{\ast }{\mathsf{\epsilon }}_{\mathrm{e}}^2+\cdot \cdot \cdot +E_{2p}^{\ast }{\mathsf{\epsilon }}_{\mathrm{e}}^p$$

(40)

$$E_{21}^{\ast }=E_{21}\left(1+{\mathrm{i}\mathsf{\eta }}_{21}\right),E_{22}^{\ast }=E_{22}\left(1+{\mathrm{i}\mathsf{\eta }}_{22}\right),\cdot \cdot \cdot ,E_{2p}^{\ast }=E_{2p}\left(1+{\mathrm{i}\mathsf{\eta }}_{2p}\right)$$

(41)

According to the principle of equal strain energy density [20], the equivalent strain ε_e of an element can be calculated as:

$${\mathsf{\epsilon }}_{\mathrm{e}}=\sqrt{{\mathit{x}}^{\mathrm{T}}\frac{{\mathbf{K}}_{20}^{\mathrm{el}}}{E_{20}{V}_2^{\mathrm{el}}}\mathit{x}}$$

(42)

where ${\mathbf{K}}_{20}^{\mathrm{el}}$ and E₂₀ are the element stiffness matrix and the storage modulus of the hard-coating at zero strain, respectively. $V_2^{\mathrm{el}}$ is the element volume of hard-coating. In addition, the displacement vector x is defined as:

$$\mathit{x}={\mathsf{\phi }}_0\mathit{q}$$

(43)

in which q is the response vector in the normal coordinate and ${\mathsf{\phi }}_0$ is the normal mode shape matrix at zero strain without damping. Submitting Equation (43) to Equation (42), one can obtain:

$${\mathsf{\epsilon }}_{\mathrm{e}}=\sqrt{\frac{{\mathit{q}}^{\mathrm{T}}{\mathsf{\phi }}_0^{\mathrm{T}}{\mathit{K}}_{20}^{\mathrm{el}}{\mathsf{\phi }}_0\mathit{q}}{E_{20}{V}_2^{\mathrm{el}}}}=\sqrt{\frac{{\mathit{q}}^{\mathrm{T}}{\mathbf{\Lambda }}_{\mathrm{el}}\mathit{q}}{V_2^{\mathrm{el}}}}$$

(44)

where Λ is the normalized element stiffness matrix per unit storage modulus of the hard-coating.

Usually, the storage modulus and the loss factor of the hard-coating can be identified by the relevant experiments. For the NiCoCrAlY + YSZ hard-coating considered here, the expressions of the strain-dependent storage modulus E₂ and loss factor η₂ obtained from the experiment [26] are:

$$E_2=54.494-0.06402{\mathsf{\epsilon }}_{\mathrm{e}}+8.09137\times {10}^{-5}{\mathsf{\epsilon }}_{\mathrm{e}}^2-6.1681\times {10}^{-8}{\mathsf{\epsilon }}_{\mathrm{e}}^3$$

(45)

$${\mathsf{\eta }}_2=0.0212-0.43521{\mathsf{\epsilon }}_{\mathrm{e}}-0.37997{\mathsf{\epsilon }}_{\mathrm{e}}^2-0.5441{\mathsf{\epsilon }}_{\mathrm{e}}^3$$

(46)

4. Solution of the Nonlinear Vibration of the Hard-Coating Cylindrical Shell

The dynamic equation of the hard-coating cylindrical shell in the normal coordinate is defined as:

$$\left({\mathit{K}}_{\mathsf{\epsilon }\mathrm{N}}^{\ast }-{\omega }^2{\mathit{M}}_{\mathrm{N}}\right)\mathit{q}={\mathit{F}}_{\mathrm{N}}$$

(47)

in which:

$${\mathbf{K}}_{\mathsf{\epsilon }\mathrm{N}}^{*}={\mathsf{\phi }}_0^{\mathrm{T}}{\mathbf{K}}_{\mathsf{\epsilon }}^{*}{\mathsf{\phi }}_0,\hspace{1em}{\mathbf{M}}_{\mathrm{N}}={\mathsf{\phi }}_0^{\mathrm{T}}\mathbf{M}{\mathsf{\phi }}_0,\hspace{1em}{\mathbf{F}}_{\mathrm{N}}={\mathsf{\phi }}_0^{\mathrm{T}}\mathbf{F}$$

(48)

where ω is the angular frequency of excitation, ${\mathbf{K}}_{\mathsf{\epsilon }}^{*}$ is the global complex stiffness matrix at ε strain and M and F are the global mass matrix and external force vector, respectively. The normalized response vector q can be obtained from Equation (47) by applying the QR-method.

It is assumed that the global material damping matrix D is independent of the strain level. Neglecting the coupling effect between the layers of hard-coating and the metal substrate of the cylindrical shell element, the global complex stiffness matrix at ε strain can be defined by:

$${\mathbf{K}}_{\mathsf{\epsilon }}^{*}\approx {\mathbf{K}}_0^{\ast }+{\displaystyle \sum _{\mathrm{el}}\frac{{\mathbf{K}}_{20}^{\mathrm{el}}}{E_{20}}\left(\mathrm{i}{\mathsf{\eta }}_{20}+E_{2,\mathrm{inc}}^{\ast }\right)}$$

(49)

in which:

$${\mathbf{K}}_0^{\ast }=\mathbf{K}+\mathrm{i}\mathsf{D}$$

(50)

The global material damping matrix D of the hard-coating cylindrical shell can be calculated by the same method as the global stiffness matrix K, which only requires multiplying Equation (6) by the zero-strain loss factor of the metal substrate and hard-coating η₁₀ and η₂₀, respectively

$$Q_{11}^1=Q_{22}^1=\frac{{\eta }_{10}{E}_1}{1-v_1^2},\hspace{1em}{Q}_{12}^1=Q_{21}^1=\frac{{\eta }_{10}{E}_1{v}_1}{1-v_1^2},\hspace{1em}{Q}_{66}^1=\frac{{\eta }_{10}{E}_1}{2\left(1+v_1\right)}$$

(51)

$$Q_{11}^2=Q_{22}^2=\frac{{\eta }_{20}{E}_2}{1-v_2^2},\hspace{1em}{Q}_{12}^2=Q_{21}^2=\frac{{\eta }_{20}{E}_2{v}_2}{1-v_2^2},\hspace{1em}{Q}_{66}^2=\frac{{\eta }_{20}{E}_2}{2\left(1+v_2\right)}$$

(52)

It should be noted that to improve the solution precision, the zero-strain loss factor of metal substrate η₁₀ is considered.

The Newton-Raphson solution method [27,28] is employed to solve the nonlinear Equation (47). Further, the corresponding residual vector r can be expressed as:

$$\mathit{r}=\left({\mathit{K}}_{\mathsf{\epsilon }\mathrm{N}}^{*}-{\mathrm{\omega }}^2{\mathit{M}}_{\mathrm{N}}\right)\mathit{q}-{\mathit{F}}_{\mathrm{N}}$$

(53)

Submitting Equation (49) to Equation (53), one can obtain:

$$\mathit{r}={\mathit{J}}_{\mathrm{N}}\mathit{q}-{\mathit{F}}_{\mathrm{N}}+{\displaystyle \sum _{\mathrm{el}}\frac{{\mathit{K}}_{20}^{\mathrm{el}}}{E_{\mathsf{\epsilon }0}^{\mathrm{c}}}\left(\mathrm{i}{\eta }_{20}+E_{2,\mathrm{inc}}^{*}\right)}\mathit{q}$$

(54)

in which:

$${\mathit{J}}_{\mathrm{N}}={\mathit{K}}_{0\mathrm{N}}^{*}-{\omega }^2{\mathit{M}}_{\mathrm{N}}$$

(55)

Since the normalized response vector q is complex-valued, it is necessary to separate the real and imaginary parts of q in the iterative formula derived by the Newton–Raphson solution method. The separated iterative formula can be expressed as:

$$\left\{\begin{array}{c}R\left({\mathit{q}}_{\mathrm{new}}\right)\\ S\left({\mathit{q}}_{\mathrm{new}}\right)\end{array}\right\}=\left\{\begin{array}{c}R\left({\mathit{q}}_{\mathrm{old}}\right)\\ S\left({\mathit{q}}_{\mathrm{old}}\right)\end{array}\right\}-{\left[\begin{array}{cc}R\left(\partial \mathit{r}/\partial {\mathit{q}}_{\mathrm{R}}\right)& R\left(\partial \mathit{r}/\partial {\mathit{q}}_{\mathrm{S}}\right)\\ S\left(\partial \mathit{r}/\partial {\mathit{q}}_{\mathrm{R}}\right)& S\left(\partial \mathit{r}/\partial {\mathit{q}}_{\mathrm{S}}\right)\end{array}\right]}^{-1}\left\{\begin{array}{c}R\left(\mathit{r}\right)\\ S\left(\mathit{r}\right)\end{array}\right\}$$

(56)

where R( ) and S( ) are the functions of extracting the real and imaginary parts of the bracketed expressions, respectively. The derivatives of the residual vector r with respect to the real and imaginary parts, q_R and q_S, are:

$$\frac{\partial \mathit{r}}{\partial {\mathit{q}}_{\mathrm{R}}}={\mathit{J}}_{\mathrm{N}}+\frac{{\mathit{K}}_{0\mathrm{N}}^{\mathrm{c}}}{E_{20}}\left[\frac{{\mathit{q}}_{\mathrm{R}}^{\mathrm{T}}\left({\mathrm{\Lambda }}_{\mathrm{el}}^{\mathrm{T}}+{\mathrm{\Lambda }}_{\mathrm{el}}\right)\mathit{q}}{V_2^{\mathrm{el}}}\left(\frac{E_{21}^{\ast }}{2{\epsilon }_{\mathrm{e}}}+E_{22}^{\ast }+\frac{3E_{23}^{\ast }{\epsilon }_{\mathrm{e}}}{2}\right)+E_{2,\mathrm{inc}}^{\ast }\right]$$

(57)

$$\frac{\partial \mathit{r}}{\partial {\mathit{q}}_S}=\mathrm{i}{\mathit{J}}_{\mathrm{N}}+\frac{{\mathit{K}}_{0\mathrm{N}}^{\mathrm{c}}}{E_{20}}\left[\frac{{\mathit{q}}_{\mathrm{S}}^{\mathrm{T}}\left({\mathrm{\Lambda }}_{\mathrm{el}}^{\mathrm{T}}+{\mathrm{\Lambda }}_{\mathrm{el}}\right)\mathit{q}}{V_2^{\mathrm{el}}}\left(\frac{E_{21}^{\ast }}{2{\epsilon }_{\mathrm{e}}}+E_{22}^{\ast }+\frac{3E_{23}^{\ast }{\mathsf{\epsilon }}_{\mathrm{e}}}{2}\right)+\mathrm{i}{E}_{2,\mathrm{inc}}^{\ast }\right]$$

(58)

The termination of the Newton–Raphson solution method can be defined by:

$${\Vert \mathit{r}\Vert }_2\le \zeta$$

(59)

where ${\Vert \mathit{r}\Vert }_2$ is the two-norm of the residual vector r and ζ is the solution precision. When the residual vector r of Equation (56) satisfies the condition of the solution precision, the newest global complex stiffness matrix at ε strain ${\mathbf{K}}_{\mathsf{\epsilon },\mathrm{new}}^{*}$ (in Equation (49)) and the normalized response vector q_new can be output. Then, the n-order nonlinear resonant circular frequency ω_n and the nonlinear harmonic response x can be obtained, respectively, by:

$$\left|\Re \left({\mathit{K}}_{\mathsf{\epsilon },\mathrm{new}}^{*}\right)-{\mathrm{\omega }}_n^2\mathit{M}\right|=0$$

(60)

$$\mathit{x}={\mathsf{\phi }}_0^{\mathrm{T}}{\mathit{q}}_{\mathrm{new}}$$

(61)

The whole solution procedure of the nonlinear vibration of the hard-coating cylindrical shell is shown in Figure 2.

5. Case Study

For this study, the NiCoCrAlY + YSZ hard-coating with strain dependence consists of NiCrAlY and yttria-stabilized zirconia (YSZ), which is the most widely used in the hottest turbine sections of the aero engine and aircraft due to its optimal performance in high temperature and vibration damping applications [3]. The relevant geometrical parameters and material parameters of the hard-coating cylindrical shell are listed in

Table 1. Geometrical parameters of the hard-coating cylindrical shell.

Parameters	L (mm)	R (mm)	T1 (mm)	T2 (mm)
Value	95	142	2	0.31

and

Table 2. Material parameters of the hard-coating cylindrical shell.

Lamina	Material	Young’s Modulus (GPa)	Loss Factor	Density (kg/m3)	Poisson’s Ratio
Metal substrate	Ti-6Al-4V	110.32	0.0007	4420	0.3
Hard coating	NiCoCrAlY + YSZ	E2	η2	5600	0.3

. It should be noted that E₂ and η₂ in Table 2 are functions of the equivalent strain ε_e (see Equations (45) and (46)).

To simulate the vibration during the aero engine or aircraft working, the hard-coating cylindrical shell is forced in the horizontal direction by a sinusoidal base excitation. Assume that the absolute acceleration of the base is along the z direction in the Cartesian coordinate system (see Figure 3). Converting it into the orthogonal curvilinear coordinate system, one can obtain the values of the influence coefficients μ_x, μ_θ and μ_z in different quadrants. The μ_x is equal to zero; the values of μ_θ and μ_z are listed in

Table 3. Values of the influence coefficients in different quadrants.

Quadrant	μθ	μz
First	cosθ	sinθ
Second	sinθ	−cosθ
Third	−cosθ	−sinθ
Fourth	cosθ	−sinθ

.

5.1. Validation Analysis

The finite element models with element size 5 mm corresponding to the proposed method and the FEIM are shown in Figure 4, which contains 3496 and 6992 elements in total, respectively. The initial base excitation level is set to 1 g. Besides, all nodes at the bottom of the hard-coating cylindrical shell are constrained in all degrees of freedom to represent the clamped-free boundary condition. It should be noted that the finite element model of the FEIM contains two layers of the element, which are coupled together to describe the connection of the hard-coating and metal substrate [16].

Since the classical shell element (Shell281) applied in the FEIM contains eight nodes and six degrees of freedom at each node and the proposed cylindrical shell element contains four nodes and nine degrees of freedom at each node, the total degrees of freedom of the finite element model of the FEIM is about 2.67-times that of the present model under the same element size, which indicates the advantage of the lower computing cost of the present model briefly.

Based on the new finite element formulation, the whole process of nonlinear vibration analysis of the hard-coating cylindrical shell is implemented by a MATLAB (MathWorks Inc., Natick, MA, USA) self-made finite element program. The numerical results of the hard-coating cylindrical shell will now be compared with the FEIM, which applied the Shell281 element of ANSYS (ANSYS Inc., Canonsburg, PA, USA) to calculate the element matrices of stiffness, mass and damping. Numerical experiences indicate that the appropriate results can only be obtained with fine enough meshing. Thus, different element sizes are used to solve the nonlinear resonant frequency of the hard-coating cylindrical shell by the FEIM and the present method, respectively. The comparison of the convergences of the first order nonlinear resonant frequencies under different element sizes is shown in Figure 5, which indicates that the developed finite element method is less affected by the element size and has lower computing cost. The advantages become more obvious when solving the nonlinear resonant response spectrums of the hard-coating cylindrical shell, because its calculation process can be regarded as a combination of a number of nonlinear resonant frequency calculation processes, which have a basically consistent iterative method.

In addition, the comparisons of the first six order nonlinear resonant frequencies and responses of the hard-coating cylindrical shell with the FEIM under element size of 1 mm are listed in

Table 4. Comparisons of the first 6 order nonlinear resonant frequencies and responses with the FEIM.

Modal Order	Nonlinear Resonant Frequencies			Nonlinear Resonant Responses
	Present Method A (Hz)	FEIM B (Hz)	Difference |A − B|/A (%)	Present Method C (10−2 mm)	FEIM D (10−2 mm)	Difference |C − D|/C (%)
1	1274.561	1281.912	0.577	1.786	1.738	2.732
2	1283.778	1294.985	0.873	3.199	3.048	4.940
3	1443.842	1446.837	0.207	1.753	1.705	2.793
4	1513.250	1526.021	0.844	1.409	1.345	4.753
5	1738.061	1736.589	0.085	1.172	1.164	0.702
6	1990.964	2002.791	0.594	1.239	1.200	3.259

. The difference represents the relative error between present method and the FEIM. As can be seen from Table 4, the nonlinear resonant frequencies and responses calculated by the present method and the FEIM show a good agreement, and the differences are within 0.873% and 4.940%, respectively, which indicate the rationality of the developed finite element method.

5.2. Nonlinear Vibration Analysis and Results Discussion

The finite element model for nonlinear vibration analysis of the hard-coating cylindrical shell is shown in Figure 4a, and the element size is changed to 1 mm. To analyze the nonlinear vibration characteristics of the hard-coating cylindrical shell, the linear resonant frequencies and responses without considering the strain dependence of the hard-coating (the equivalent strain ε_e = 0) need to be calculated, which are taken as the reference.

The first six order linear and nonlinear resonant frequencies and responses of the hard-coating cylindrical shell under a 5 g excitation level are listed in

Table 5. The first 6 order linear and nonlinear resonant frequencies of the hard-coating cylindrical shell under a 5 g excitation level.

Modal Order	Linear (Hz) E	Nonlinear (Hz) F	Descent (Hz) E − F
1	1274.899	1274.561	0.338
2	1283.922	1283.778	0.144
3	1444.074	1443.842	0.232
4	1513.410	1513.25	0.160
5	1738.293	1738.061	0.232
6	1991.071	1990.964	0.107

and

Table 6. The first 6 order linear and nonlinear resonant responses of the hard-coating cylindrical shell under a 5 g excitation level.

Modal Order	Linear (10−2 mm) G	Nonlinear (10−2 mm) H	Descent (10−2 mm) G − H
1	8.926	8.529	0.397
2	16.003	15.472	0.531
3	8.740	8.444	0.296
4	7.067	6.794	0.273
5	5.869	5.714	0.155
6	6.195	6.067	0.128

. Then, the three-order linear and nonlinear resonant frequencies and responses of the hard-coating cylindrical shell under different excitation levels are listed in

Table 7. The 3 order linear and nonlinear resonant frequencies of the hard-coating cylindrical shell under different excitation levels.

Excitation Level (g)	Linear (Hz) E	Nonlinear (Hz) F	Descent (Hz) E − F
1	1444.074	1444.026	0.048
3	1444.074	1443.932	0.142
5	1444.074	1443.842	0.232
7	1444.074	1443.755	0.319
9	1444.074	1443.672	0.402

and

Table 8. The 3 order linear and nonlinear resonant responses of the hard-coating cylindrical shell under different excitation levels.

Excitation Level (g)	Linear (10−2 mm) G	Nonlinear (10−2 mm) H	Descent (10−2 mm) G − H
1	1.752	1.741	0.011
3	5.246	5.136	0.110
5	8.740	8.444	0.296
7	12.234	11.641	0.593
9	15.738	14.728	1.010

. Moreover, the three-order frequency response spectrums of the hard-coating cylindrical shell are shown in Figure 6. It should be noted that the legend L represents the linear vibration response, and the legend Nl represents the nonlinear vibration response.

According to the results listed in Table 5 and Table 6, it can be found that the nonlinear resonant frequencies and responses decrease to a certain degree compared with the linear calculation results, which indicate that the strain dependence of the hard-coating would make the resonant frequencies and responses present a certain degree of decline. As can be seen from Table 7 and Table 8, the descents of the nonlinear resonant frequencies and responses of the hard-coating cylindrical shell increase continually compared with the linear calculation results, while the excitation level increases from 1 to 9 g; that is, the increase of the excitation level would make the strain dependence of the hard-coating more remarkable, which reveals the characteristics of the soft stiffness nonlinearity or “strain softening”. This phenomenon may be more obvious and prominent in the frequency response spectrums (see Figure 6), which may indicate that the strain dependence of the hard-coating is beneficial for vibration reduction.

6. Conclusions

With considering the influence of the strain dependence of the coating material on the complex stiffness matrix, the nonlinear modeling, validation example and analysis are studied, and the following conclusions can be drawn.

• Based on Love’s first approximation theory, a four-node composite cylindrical shell finite element model is proposed. Then, the nonlinear iterative solution formulas with a unified iterative method are theoretically derived for solving the resonant frequency and response of the hard-coating cylindrical shell.
• A cylindrical shell coated with a thin layer of NiCoCrAlY + YSZ is chosen to demonstrate the proposed formulation. The nonlinear resonant frequencies and responses calculated by the present method and the FEIM show a good agreement, which indicates the rationality of the developed finite element method. Moreover, the developed finite element method is less affected by the element size and has lower computing cost.
• Moreover, the nonlinear vibration analysis of the cylindrical shell coated with a thin layer of NiCoCrAlY + YSZ is implemented. Compared with the linear calculation results, the nonlinear resonant frequencies and responses of each order decrease to a certain degree, and the descents increase continually with the increase of the excitation level; that is, the increase of the excitation level would make the strain dependence of the hard-coating more remarkable, which reveals the characteristics of the soft stiffness nonlinearity or “strain softening”.