A hybrid meshless–statistical energy analysis method for complex structure vibration analysis

A new hybrid deterministic–statistical energy analysis (SEA) formulation is presented by introducing a meshless method for modeling the deterministic components. Moving least square Ritz (MLSR) meshless method is applied, in which MLS is used to build the discrete model and the Ritz method allows to obtain variational formulation of the deterministic components of the governing equations. Such governing equations can be formulated via boundary conditions by penalty method and Lagrange multipliers. The hybrid model by penalty method keeps a similar formulation with the framework of the finite element SEA, while the model by the latter increases the size of the dynamic stiffness matrix and the expanded components are determined by the constraints. For validation purpose, three case studies are provided, including beam–coupled plates and plate–coupled plates built-up structure. The results by the hybrid MLSR-SEA model are compared with those by FE-SEA and Monte Carlo simulation. Good agreements of responses between the methods demonstrate the reliability of the MLSR-SEA formulation.


Introduction
Manufacturing and assembly imperfections exist in all engineering structural systems. The uncertainties caused by the imperfections lead to difficulties in predicting the vibration responses of complex structures. Moreover, responses of such structures in high frequency range are even more sensitive to the uncertainties. Traditional deterministic methods, e.g., finite element (FE) methods, cannot achieve accurate predictions, unless considering either large number of basis functions or very fine meshes, and both options are computationally expensive. In order to overcome this shortcoming, statistical energy analysis (SEA) was proposed for predicting the high frequency range dynamic response of structure ensembles in the presence of uncertainties [1,2]. A few decades of investigation has demonstrated that the SEA is a reliable method for this class of problems [3,4]. The SEA was presented initially by Lyon to describe the energy transmission and dissipation in the structural system with uncertainties [2,3]. The energy input into the statistical component of the structure (subsystem) is considered to be equal to the sum of the energy dissipated and transmitted to the other components. One key factor when applying the SEA analysis is to calculate the coupling loss factor, and several publications contribute to this respect [5][6][7][8][9]. The traditional SEA has proven successful in predicting the ensemble average energy in high frequency range where the uncertainties result in modal overlapping and lead to statistical modes. However, the low-and mid-frequency modes are less influenced by uncertainties and present less modal overlapping, which causes inapplicability of SEA in these frequency ranges. To improve the robustness in wider frequency range, hybrid models mixing SEA with other deterministic methods (e.g., FE method) were developed. Under the mode-based approach category, Langley presented a hybrid SEA model in which the global and the local modes are considered to capture the characteristics of deterministic and statistical components, and it can realize all-frequency-range analysis [10]. Targeting the ensemble average response in mid-frequency range, Langley proposed a hybrid FE-SEA model featuring a wave-based approach [11,12]. The deterministic components of the structure are modeled by FE method. The statistical field in the subsystem can be seen as the superposition of direct and reverberant field, and the reciprocity relationship plays an important role in providing connection between energy level and reverberant field [13]. Moreover, the transient SEA and its hybrid formulation were developed further to solve the time-domain problem with the impact or time-varying loads [14][15][16]. Regarding nonlinear structures including SEA subsystems, the equivalent stiffness method were applied to linearize the nonlinear joints [17][18][19][20][21]. The variance of the ensembles was also investigated under the framework of SEA, transient SEA and its hybrid formulation [22][23][24][25][26].
The hybrid FE-SEA has been widely explored and proved to be reliable as long as the FE method is applicable to deterministic components. In fact, some other numerical methods for deterministic problem can be more effective than FE method in certain areas. Hence, developing other hybrid deterministic-SEA formulation is a natural extension and improvement in the FE-SEA research field. Vergote applied a wavebased (WB) method as a substitution for the FE in the hybrid model, and it can achieve faster convergence than the FE method based one [27]. Considering the advantage of the boundary element (BE) method in the area of acoustics, Gao developed a hybrid BE-SEA formulation which reduces the number of degrees of freedom (DoFs) compared to the FE-SEA [28]. Wu introduced an embedding of the gradient smoothing technique to the FE-SEA framework, which leads to much more accurate results in numerical cases [29]. Clot proposed a hybrid FE-SEA-experimental model in which some of the structural parameters are determined by experimental data [30,31]. Besides the deterministic methods mentioned above, another numerical method, i.e., the meshless method (MM), has rapidly increased popularity in the last 30 years and developed in a large class of numerical methods. The investigation into MMs is mainly based on two aspects: one related to the construction of the shape functions (discretization) and the other to the formulation procedures. Comprehensive discussion about these two topics was given by Liu [32]. Monaghan [33] proposed the smoothed particle hydrodynamics (SPH) which was initially applied to astrophysical motion. Belytschko developed Element-free Galerkin (EFG) methods based on the moving least square (MLS) which allows to construct meshless shape functions by nodes in least square sense [34]. Liu employed polynomial and radial basis point interpolation combined with Galerkin method or boundary integral equation to realize meshless analysis [35,36]. Collocation methods as a class of procedure constructing the formulation are also developed in meshless sense [37,38]. The corrected collocation method was then derived to enforce the essential boundary conditions, and it has been successfully applied to varieties of problems [39]. A MLS-Ritz (MLSR) method and its improved formulation were proposed to analyze the vibration of beams and plates [40,41]. With the development of meshless formulation, the method is also applied to many research fields, e.g., nanomaterials, thermoelastics [42][43][44][45][46]. Moreover, some other numerical methods for structural vibration and their applications are recently discussed [45,[47][48][49]. MMs has several advantages over the FE method, for instance, MMs do not require a mesh-which is a major time-saving in modeling structures with complex geometry. Moreover, the mesh-free nature means that it can provide more accurate analysis for cases where structural dynamic deformation causes low-quality mesh. In the field of large-deformation structures, traditional FE method often fails to provide smooth interpolation functions when discretizing the model, which may decrease the response accuracy. On the contrary, MMs can provide highly smooth ones which may even be multi-order differentiable, which guarantees more accurate analysis.
In order to further extend the current FE-SEA framework, this research work features a typical MM to describe the motion of deterministic components instead of FE method. Among all the MMs, the MLS-Ritz (MLSR) method is broadly employed and proved robust and effective in engineering structural analysis. Specifically, this work targets the development and validation of the hybrid meshless-SEA model, focusing on deriving the MLSR-SEA formulation. The deterministic components are modeled by MLSR and the statistic components refer to SEA. The couple between them is based on the wave-based approach from Langley's work [11,13]. It is noted that the shape functions by MLS do not include information from the displacement boundary conditions, and this is different from the FE method in which the boundary conditions are considered when applying the polynomial interpolation. Hence, classic methods targeting constraint conditions such as the penalty method and Lagrange multiplier method are firstly introduced in the process of developing the hybrid MLSR-SEA formulation. The formulation regarding Lagrange multiplier method can be seen as one similar form to the work by Ref. [50], but the difference is that the derivation and the numerical results are under the meshless condition. The penalty method and Lagrange multiplier method possess different features and lead to different forms of governing equations. The present investigation derives the hybrid formulation regarding both of them.
Referring to the validation of deterministic-SEA model, previous investigation provided a reliable Lagrange-Rayleigh-Ritz method (LRRM) in which the randomly distributed lumped masses are employed to break the system symmetry of modes, and Monte Carlo simulation (MCS) is applied to obtain the response ensembles [18,19]. LRRM-MCS has seen as a standardized method for verification of SEA and hybrid SEA models. In terms of validating the MLSR-SEA formulation, it is natural to apply the LRRM-MCS as a benchmark model for verification. For deterministic modeling, MLS is used to build the shape functions, while the traditional polynomial is used for statistical components. The Ritz method as equivalent to the Lagrange equation method is applied as a variational formulation to build the governing equations. This paper regards this benchmark model as the Ritz-MCS, and it should be noted that different types of shape functions are used in subsystems and deterministic components, and the couplings between them are realized also by penalty method and Lagrange multiplier method in the Ritz procedure. With respect to every ensemble in the MCS, lumped masses, which allow to randomize the systems, are set to randomly locations on the subsystem. The range of proper magnitude and number of lumped masses has been discussed in our previous work [18]. In the section of numerical results, three cases are presented for validation and the results obtained with the Ritz-MCS method. To further explore the verification of proposed formulation, the FE-SEA and the pure FE by Abaqus are introduced to calculate the results as comparisons.

Hybrid MLSR-SEA formulation
This section derives the hybrid MLSR-SEA formulations regarding both penalty method and Lagrange multiplier method. To achieve the formulation of hybrid model, the first step is obtaining the governing equation by MLSR method. This procedure includes the construction of shape functions by MLS and the Ritz procedure modeling. The formulation with the enforcement of boundary conditions by penalty method and Lagrange multiplier method has been achieved in Ref. [32,34]. The procedure is summarized in the Appendix 4 for integration of the formulation and for readers' easier understanding of the origin of the additional terms in the governing equations. The second step is to build the hybrid MLSR-SEA formulation by coupling the deterministic fields with statistical fields in subsystems.
The MLSR governing equation by penalty method can be obtained as follows: where M is referred as the global mass matrix; K is the global stiffness matrix; F ext is the external force vector; K α is the global penalty stiffness matrix; and F α is the penalty force vector. Enforcing the boundary condition by Lagrange multiplier method, the governing equation by where M, K and F ext have the same meaning and form as those by penalty method; λ is the Lagrange multiplier vector, and the dimension of the vector equals the number of constraints to the deterministic components; the matrices G and Q are related to the parameters of the constraints. For the penalty method formulation Eq. (1), the number of variables, i.e., the generalized coordinates q, is equal to the number of DoFs. Its mathematical form is similar to the one obtained with conventional FE methods and the shape functions by both MLSR and FE are localized. But the governing equation Eq. (2), obtained via Lagrange multiplier method, has more variables than the number of DoFs due to the multiplier vector λ. Hence, it is expected that the final formulation is different from the previous one.

Hybrid MLSR-SEA formulation using the penalty method
In the framework of hybrid deterministic-SEA model by Shorter and Langley, the structural system is divided into the deterministic component, which is modeled by deterministic methods, and the subsystem, which is described by statistical fields [11]. The statistical displacement field in the subsystem can be considered as the superposition of the direct field and the reverberant field from the wave approach perspective. The current MLSR-SEA formulation applies a different numerical method (MLSR), comparing with standard FE-SEA, but does not change the assumption for the statistical fields (subsystems). Hence, its statistical field assumption is proper to be applied in the present new formulation.
Transforming Eq. (1) from time domain to the frequency one and defining F = F ext + F α , the governing equation of deterministic components becomes where D d is the dynamic stiffness matrix. By including the force induced by the direct and the reverberant field to the deterministic components, the equation of motion in a hybrid formulation can be written as where in which f k rev is the reverberant force and D k dir is the direct field matrix. The governing equation Eq. (4) has the same mathematical form and assumption for subsystems as the in FE-SEA [11]. Following an analogous process to derive the SEA equation, the energy equilibrium equation for each subsystem j and the cross-spectral matrix S qq of deterministic displacement field are given as where S F F is the cross-spectral of F; η j is the loss factor of jth subsystem; η d, j corresponds to the power dissipation in jth deterministic component; η jk is the coupling loss factor; n j is the modal density; E j is the ensemble average energy of jth subsystem; P ext in, j represent the power input from the loading to deterministic components, respectively.

Hybrid MLSR-SEA formulation using the Lagrange multiplier method
This subsection presents the hybrid MLSR-SEA formulation in which the deterministic components are governed by Eq. (2). The main difference regarding the formulation between this equation and the one by FE-SEA method is that Eq. (2) contains Lagrange multiplier vector λ which refers to the boundary constraints. Equation (2) can be rewritten as where Defining N n as the number of DoFs of deterministic components by MLSR method and N λ as the length of Lagrange multiplier λ, the dimension of D det,λ is (N n + N λ ) × (N n + N λ ) and that of p and F is (N n + N λ ) × 1.
Considering the coupling between the direct and reverberant fields (subsystems) to the deterministic components, the governing equation above can be written as where D k dir,λ is the direct field matrix regarding kth subsystem and f k rev,λ is the force loaded by reverberant field of kth subsystem. It is noted that . It is assumed that the essential boundary conditions of the deterministic components are independent to the coupling between the deterministic components and subsystems because of the local properties of the junction [11]. Hence, D k dir,λ is an expanded matrix from D k dir and the expanded items are all zero. It is expressed as Using Eq. (12) the cross-spectra matrix can be obtained as where According to the reciprocity relation, the cross-spectra of force by reverberant field can be written as [13] Due to the expanding component in S dir,λ being zeros, Eq. (16) can be rewritten as Considering the ensemble and time average power flow in reverberant fields of the subsystems, the power balance requires that the power into a subsystem should equal the sum of the power dissipated within the subsystem and the power flowing out to others. The power balance in the jth subsystem can then be written as where P in, j is the power injected into the jth subsystem; P rev out, j is the power flowing out of jth reverberant field; P diss, j is the power dissipated in jth subsystem. P in, j is given as The transformation from Eqs. (19) to (20) is based on the matrix block expansion in Eq. (13). Due to the zero blocks, the terms regarding λ are also zero, which means introducing Lagrange multipliers do not change the form of the power input into the subsystems. This is a proper inference because both boundary constraints (referring to Lagrange multipliers) and junctions (referring to power transferred to subsystem) have the property of localization. Combining Eqs. (14), (17) and (20), the input power can be rewritten as where The power flowing out of the reverberant field from jth subsystem can be written as Substituting Eqs. (15) and (17) into the above equation gives where η d, j is power dissipated in deterministic components. P diss, j is the power loss within subsystem and can be given as where η j is the damping loss factor. It should be noted that the external input power P ext, j , the coupling loss factor η jk and the energy transmitted to other subsystem P rev out, j , share a similar form with those in original FE-SEA formulation in Ref. [11]. The differences are that all the matrix expressing the energy power follow in present derivation is expanded due to the existence of Lagrange multiplier λ. Substituting Eqs. (21), (26) and (28) into power balance equation Eq. (18) gives Substitute the reciprocity equation Eq. (17) into Eq. (14), then the response cross-spectra of deterministic component can be written as It can be seen from the governing equations that the formulation regarding Lagrange multiplier has mathematically similar derivation form with those by penalty method. However, the Lagrange multiplier technique enlarge the dimensions of dynamic stiffness matrix. And if the boundary constraints are complex and the number of constraint equations is large, the computational time may increase obviously comparing to the formulation by penalty method.

Numerical results
This section is devoted to the validation of the proposed hybrid MLSR-SEA formulation by three case studies. The case studies feature the meshless MLSR formulation for deterministic components, which are seen as the joints, and the SEA model for statistical components (subsystems). The penalty method and Lagrange multipliers representing two different techniques to take account boundary constraints lead to different formulations but in fact tend to yield nearly the same analysis results for the following cases. Hence, the analysis results in this section by MLSR-SEA model only refer to the formulation regarding penalty method. In addition to numerical results obtained with MLSR-SEA, those realized by several other methods (such as FE-SEA) are computed and compared and are used to validate the new formulation. The Ritz-MCS method which has been mentioned in the Introduction and is used to estimate the ensemble average responses in all three cases.
Considering the importance of the beam and plate in engineering structure analysis, they are introduced as structural components of current case studies. The subsystems in all three cases are the out-of-plane motion of the plates. The modes regarding this motion are more likely to mix and veer over each other within the midand high frequency range. Hence, the out-of-plane motion is proper to be seen as the statistical component. For deterministic component, both one-dimensional (beam) and two-dimensional (plate, in-plane motion) structures are modeled by MLSR method. In the following cases, Timoshenko beams are modeled as the deterministic components (joint); Kirchhoff thin plate is applied and the out-of-plane motion is considered as statistical, while the in-plane motion is robust to uncertainties, within the investigated frequency range, and is considered as deterministic. All the materials in the structures are assumed to be homogeneous, isotropic and linear elastic with Young's modulus E = 70 GPa, Poisson's ratio ν = 0.3 and density ρ = 2700 kg/m 3 . The magnitudes of the beams and plates are referred in Table 1.  The first case study (Case 1) is composed by two plates (statistical subsystems) joined by a beam as shown in Fig. 1. The system is excited by a harmonic force P loaded on the plate 1. The three edges of both plates are simply supported, except for the one connected to the beam. The dimensions of plate 1 and plate 2 are given in Table 1, the 1.2 m left edge of plate 1 is connected to the right edge of plate 2 via the beam component. For every realization of the MCS sample, the random masses on the plates are distributed to obtain the energy response. Hundred Monte Carlo samples are computed to calculate the ensemble average energy. In the formulation of the MLSR-SEA, the plates are considered as subsystems while the beam is the deterministic component modeled by MLS-Ritz method. Figures 2, 3 and 4 show the energy response level of the plates 1, 2 and beam by the methods. It can be seen generally for plate 2 and beam (as the receiving structures in the system) from the figures that the results by the MLSR-SEA formulation are consistent with the one by reference methods. The energy peak at around 2000 rad/s is well predicted by MLSR-SEA model for both plate 2 and the beam. Even in low-frequency range below 1500 rad/s, the dramatic fluctuations of ensemble average energy exist by Ritz-MCS analysis, for the reason that the uncertainties are not large enough to achieve full overlaps between the modes within this frequency range. However, the ensemble average energy by MLSR-SEA can still predict the trend for plate 1, plate 2 and beam. Moreover, obvious differences of the There are two reasons causing this energy difference. First is that some of the energy output from the plate 1 is stored in the beam. Second is that the boundary between the beam and the plate causes reflection of the wave (mostly from reverberant fields) plate 1, and just a part of energy in the wave fields is transferred to the beam. Figure 5 includes curves of ensemble average energy of plate 1 and plate 2 by MLSR-SEA, and corresponding 99% confidence intervals by Ritz-MCS samples. It is clear that energy responses within most frequency range are included in the confidence intervals in the frequency range over 1000 rad/s. Figures 6,  7 8, and 9 show the response convergence and accuracy in selected frequencies of plate 2 and the beam. The red dotted lines represent the energy response convergence when increasing the number of MLS meshless nodes. When the number of nodes reaches 15, the energy level tends to stabilize, and the response is within the 3dB error in comparison with those by the benchmark models. These figures show the good convergence and accuracy of the new MLSR-SEA formulation.
For the second test case (Case 2), the size of the cross section of the beam is set as 0.04 m by 0.06 m, which is larger than the value of the previous one, and the receiving plate is set with an inclination angle β = π/4. This case, by introducing inclination angle, represents a more general study than case 1 when referring to beam- coupled-plate structures. The second test case structure is depicted in Fig. 10. After applying the MLSR-SEA formulation and the validation models, the ensemble average energy is calculated as shown in Figs. 11, 12 and 13. A good agreement between responses by MLSR-SEA and validation models can be observed. In fact, the inclination angle lead to the in-plane motion which corresponds to much larger stiffness. Hence, the z-direction motion of beam is restricted, and this cause the decrease in energy level of the beam in comparison with the case 1. However, the energy response of plate 2 seems to be less different for case 1 (Fig. 3) and case 2 (Fig. 12). This is because the torsion of the beam contributes much to the energy transfer between plate 1 and plate 2, and the torsional motion is barely influenced by inclination angle. It can be concluded from this case study that both the bending and torsional stiffness of the beam affect the energy transfer effectiveness. The realization of vibratory energy barrier requires the consideration of both bending and torsional motions in high frequency range. The final case study (Case 3) looks into the structure containing the plate as the joint transmitting the energy from one subsystem to the other. The schematic geometry of this case is shown in Fig. 14. This case represents the power transmission model for coupled plates. The parameters of the plates 1a, 1b and 1c are the  Table 1, and those of plates 2a and 2b are the same as that of plate 2 in Table 1. The plate 1a is seen as the main source plate on which the excitation P is introduced, and all other plates are the receiving plates. When the energy flows from plate 1a to plate 1b and than from plate 1b to plate 1c through the joints (plate 2a and plate 2b), out-of-plane motion mainly excites in-plane (transverse) motion of the vertical plates (2a and 2b). The in-plane motions in plate 2a and 2b corresponding to large transverse stiffness can hinder the dominant energy transmission between plate 1a, 1b and 1c. The ensemble average energy for plates 1a, 1b and 1c is shown in Figs. 15, 16 and 17. It can be seen that all these three methods provide very close average energy response. Especially for the results by MLSR-SEA and FE-SEA, the responses are almost coincident, for the reason that within the calculation frequency range, the deterministic modelings by MLSR and FE methods for in-plane motion are nearly the same. Even though the in-plane motion of the vertical plates is considered, in fact, the energy transfer is mostly realized by the out-of-plane motion of the vertical plates. It has been calculated that the first resonance frequency of the in-plane motion of vertical plates is very large and beyond the frequency range of analysis. Hence, in this case study, the motion along the vertical direction at the coupling area can be assumed to be zero.

Conclusion
The present investigation derives a new hybrid meshless-SEA formulation as an extension to the FE-SEA framework. The rational behind the new formulation is associated with an increase in the energy response accuracy through a meshless technique which could substitute the FE method in the modeling of the deterministic components within a complex dynamical systems. In this paper, a meshless method MLSR is used to construct the model for deterministic components. Because the MLS method does not enforce boundary conditions, penalty method and Lagrange multiplier method are employed for the enforcement. The MLSR-SEA formulation regarding the penalty method shares a very similar form with FE-SEA since the penalty terms do not change the dimensions of the dynamic stiffness matrix, but requires a specified penalty factor. However, the formulation by Lagrange multiplier method requires expansion of the dynamic stiffness matrix. This expansion makes the formulation slightly more complex. For validation purpose, three cases considering both beams (one-dimensional structure) and plates (twodimensional structure) as the deterministic component were investigated to take a step to real-life structures. A benchmark model based on the Ritz method and MCS was employed, and FE-SEA and FE-MCS were also applied to enhance the validation. In those three case studies, the hybrid MLSR-SEA model yields highly consistent results to the other analysis results. The case studies could provide a solid validation for the proposed MLSR-SEA formulations. This work lays the foundation for the application of the meshless-SEA formulation to more complex structures. As described in the Introduction, MMs has advantages over FE in certain engineering fields. Hence, applying hybrid MLSR-SEA model to these fields, e.g., crack propagation within deterministic components, can be a further investigation.

Appendix A: MLS-Ritz method
This section of appendix introduces the meshless MLSR method. In the first subsection, MLS is employed to construct the shape functions to approximate the displacement. One feature of the construction is that the nodes, along with their corresponding compact support domain, need to be specified in the process of building the approximations. The displacement at the nodes provides basis for the approximation and the compact support domain defines the localization of the shape functions. In the present section, the shape functions are generated by interpolating the displacement of nodes through polynomials in the least-squares sense. The coefficients of the polynomials are the functions of the spatial coordinates at the nodes, i.e., the coefficients are 'moving' along with the coordinates. The second subsection introduces the Ritz method for the computational procedure. Since the shape functions by MLS lack information of essential boundary conditions, implementing the constraints by boundary condition becomes necessary in the procedure. Both penalty method and Lagrange multipliers method allow to enforce the boundary conditions, but lead to different forms of the governing equations. Since they may yield different mathematical formulations when coupled with SEA subsystem, both derivations are presented below.

A.1 Moving least square (MLS)
Within the framework of MMs, the moving least square (MLS) technique employed hereinafter is based on the interpolation on the moving nodes in the least square sense. The approximated displacements at the spatial coordinate x can be written as a form of finite function sequence [32], where p(x) is the polynomial form vector of the N m basis functions; p i (x) represents ith basis function; a i (x) denotes the coefficient of ith basis function. The weighted square residual between exact displacement and the approximated one can be constructed as follows: where N n is the number of nodes within the compact support domain; W (x − x i ) denoted in the rest of the paper as W i (x) is the weight function varying with the distance from x to the node x i . It can be denoted that closer distance to the node tends to encompass more weight. The minimization of the weighted square residual J in Eq. (A2) with respect to a(x) is achieved by solving the following equation: Then a(x) can be derived as where The node parameter vector q in Eq. (A4) can be seen as the generalized displacement; A(x) is referred as the weighted moment matrix. Substituting Eq. (A4) into Eq. (A1) gives the estimation of displacement at x as and where (x) is the vector of shape functions. It is noted that MLS can realize high-order differentiable shape functions by introducing only a few basis functions, which contributes to a high accuracy of the Ritz method.

A.2 MLS-Ritz formulation by penalty method
The following considers the Ritz method which is based on the principle of minimum potential energy and applied to structural vibration. Two steps are included in this procedure: One is to derive the total potential energy by considering the elastic potential and external input potential, and the other is to minimize the energy functional by variational principle. Notably, because shape functions by MLS lack of information regarding boundary condition, the energy functional should include an extra term in least square sense to penalize the local deviation from boundary condition according to penalty method. Consider a dynamic structural model defined in domain with the following stress-strain and strain-displacement relationships: where σ denotes the stress tensor; ε is the strain tensor; D is matrix of material constants related to the constitutive law; L is a differential operator regarding the relationship between strain and displacement. The essential boundary conditions in boundary b are assumed as where the vectorū contains the prescribed boundary constraints. The energy functional is expressed as the sum of the elastic potential energy E e , potential energy from the external input E ext and the penalized residual in terms of the boundary condition, and is written as where α is the penalty factor and the last term denotes the penalty to the energy functional. The elastic potential energy can be expressed as Consider the inertial force as a type of body force. Then, the energy from external loading can be written as where ρ is the mass density; f b and f t are the external body force and the external surface force in the domain b and t , respectively. Substituting Eqs. (A5), (A7), (A8), (A11), (A12) into Eq. (A10) and minimizing the functional according to the Ritz procedure as give the governing equation as follows: where M is referred as the global mass matrix; K is the global stiffness matrix; F ext is the external force vector; K α is the global penalty stiffness matrix; F α is the penalty force vector.

MLSR formulation by Lagrange multiplier
When applying the Lagrange multipliers, the term regarding the constraints is different from that by penalty method in Eq. (A10). But the elastic potential energy and the energy from external force have the same expressions as in Eqs. (A11) and (A12). The functional of total potential energy with respect to Lagrange multipliers can be written as where λ is the Lagrange multiplier vector, the size of which is determined by the number of constraints. Similarly, substituting Eqs. (A5), (A7), (A8), (A11) and (A12) into Eq. (A20) and then minimizing the energy functional as δ δq = 0 and δ δλ = 0, yield the governing equation in the form of Lagrange multipliers which can be given as where M, K and F ext have the same form as Eq. (A15)-(A17).