Formulation for wave propagation in dissipative media and its application to absorbing layers in elastoplastic analysis using mathematical programming

In this article, we propose a new solution scheme for modeling elastoplastic problems with stress wave propagation in dissipative media. The scheme is founded on a generalized Hellinger–Reissner (HR) variational principle. The principle renders the discretized boundary‐value problem into an equivalent second‐order cone programming (SOCP) problem that can be resolved in mathematical programming using the advanced optimization algorithm—the interior point method. In such a way, the developed method not only inherits admirable features of the SOCP‐based finite element method in solving elastoplastic problems but also enables the enforcement of absorbing layers (i.e., Caughey absorbing layer), which is essential in modeling stress wave propagation problems, to absorb wave energy. The proposed scheme is validated via the comparison between analytical and numerical results for seismic wave propagation in dissipative media. Its application to elastoplastic dynamic problems with stress wave propagation is also illustrated to demonstrate its efficiency.


INTRODUCTION
Stress wave propagation in the ground is related to many geotechnical and geological problems. Under dynamic loadings, geomaterials and/or structures undergo different degrees of yielding, leading to very complicated deformation mechanisms. Typical examples are the earthquake-induced hazards such as landslides 1 and damage of infrastructure. 2 Additionally, dynamic analysis of geotechnical problems and design in engineering often require the accurate description of the stress wave propagation. The penetration of gravity anchors during installation is an example that the stress wave generated by the impact disturbs the surrounding soils whose mechanical properties change sequentially. 3 The analytical solutions of these problems are usually not available due to the high nonlinearity of the associated governing equations and therefore they are commonly resolved via numerical methods such as the finite element method (FEM). [4][5][6][7] When developing nonlinear finite element programs, the nested scheme is usually adopted that the solution is sought out via iterations between global structural levels, where out-of-balance forces are minimized, and local material points, where stresses are integrated implicitly or explicitly based on constitutive laws. The iteration between the two levels is conducted in the manner of elastic prediction and plastic correction. 8 While the nested scheme is popular, the FEM in mathematical programming (MP) has long been recognized as a possibility. 9 Notably, the FEM in second-order cone programming (SOCP) has been shown as an efficient and robust strategy. In SOCP-based FEM, the governing equations are cast as a standard SOCP problem which is then resolved using the interior point method (IPM). This solution strategy possesses many admirable features, such as the forthright treatment of singularity in some yield criteria (i.e., Mohr-Coulomb and Drucker-Prager criteria), the natural extension from single-surface to multi-surface plasticity, among others. 10 Moreover, the convergence of modern IPM is super linear and the number of iterations for desirable convergence is almost unaffected by the number of unknowns. 11,12 In light of its advantages, this solution strategy has been extended for solving many challenging problems in geotechnical and geological engineering such as limit analysis, 13,14 elastoplastic analysis, 10,[15][16][17] viscoplastic analysis, 18,19 contact analysis 20,21 , analyses of fracture propagation 22,23 , etc. Many recent publications indicate that this MP strategy is very efficient when analyzing engineering problems with plasticity. 13,18,[24][25][26] Despite its development, the possibility of the SOCP-based FEM for modeling problems with stress wave propagation, to the authors' best knowledge, has not been explored. A key challenge associated is the enforcement of artificial boundary conditions in the SOCP-based FEM framework. By absorbing wave energy, the artificial boundary condition ensures that the finite domain used in numerical modeling can mimic an unbounded domain.
The most powerful artificial boundary condition is probably the perfectly matched layer (PML). The PML introduces a rigorously designed absorbing boundary layer to attenuate the incident wave. [27][28][29] The implementation of PML in finite element modeling, however, is nontrivial [28][29][30][31][32][33][34] and, so far, it has been mostly used for elastic or viscoelastic problems. The absorbing boundary layer is an alternative in which waves are absorbed by gradually increasing the damping coefficients. Pioneered by the works in the 1970s and 1980s, [35][36][37] several versions of the absorbing layer method have been developed. [38][39][40][41] Among them, the Caughey absorbing layer method (CALM), based on the classical concept of the Caughey damping 42 and with the rheological interpretation and the optimization of Rayleigh coefficients, 40 is an efficient multi-direction method. More recently, CALM was developed by a rigorous multi-layer strategy 43 in which the strong form of the elastic wave propagation in a damping Rayleigh medium was implemented. The basic idea of CALM is to attenuate the waves in an absorbing layer of finite thickness using the Caughey damping. 40 Compared to the PML, the implementation of the CALM is simpler while acceptable absorbing efficiency is guaranteed. 40,43 In these simulations for wave propagation in elastic and viscoelastic media, a hybrid time-integration scheme (explicit time integration for the domain of interest and implicit time integration for the absorbing subdomains) is developed to avoid the decrease of critical time step due to damping effects. 43,44 To further extend the capability of the SOCP-based approaches for analyzing engineering problems, in this study we explore the possibility of the SOCP-based FEM for elastoplastic analysis involving stress wave propagation. Specifically, a new generalized Hellinger-Reissner (HR) variational principle is proposed for the incremental analysis of elastoplastic dynamic problems with absorbing layers. The proposed generalized HR variational principle leads to a fully implicit finite element formulation in the form of standard SOCP that can be resolved using the advanced IPM. To validate the developed formulation and solution scheme, a series of seismic wave propagation problems are simulated, in which the performance of the CALM is also investigated. Finally, the absorbing boundary layer is used for an elastoplastic problem to demonstrate its efficiency in energy absorption.

GOVERNING EQUATIONS
In this section, we introduce the governing equations of elastic wave propagation in a dissipative medium characterized by Rayleigh damping and the constitutive relationship for elastoplastic materials.

Medium without damping
Considering a 2D domain V delimited by a boundary S, the set of equations for dynamic analyses are as below.
a. Equilibrium equations where , b, and a denote the stress, body force, density, and acceleration, respectively. The acceleration a is the time derivative of velocity v.

c. Boundary conditions
where N is the unit normal vector, t is the tractions and u p is the prescribed displacements. Equations (3) and (4) are the traction boundary condition and displacement boundary condition, respectively.
where S is the compliance matrix, is the plastic multiplier and G is the plastic potential. The incremental form of Equation (5) is given as follows: For elastoplastic materials, the yield function F( ) is used to determine the elastic or plastic state and thus leads to the compact form for elastoplastic media

Dissipative medium with the Rayleigh damping
Rayleigh damping is a second-order Caughey damping formulation, 42 in which M and K are the Rayleigh coefficients of the mass and stiffness matrices, respectively. Therefore, the governing equations of seismic wave propagation in a dissipative elastic medium characterized by Rayleigh damping can be expressed as follows 43 : in which, K is a stress tensor related to the stiffness component of the damping force and M is a force vector related to the mass component of the damping force. In the above expression, the Rayleigh damping is equivalent to the combination of a viscous Kelvin-Voigt model and a dashpot 43 as illustrated in Figure 1. As seen, to construct the Kelvin-Voigt model, a dashpot is added to the right-hand side of Equation (9) with another one included parallelly on the left-hand side. The Rayleigh damping coefficients M and K are relevant to the predominant frequency 0 , and are assumed to satisfy the relation 40,43 : where is the damping ratio.

Time discretization
By means of the -method, 46 the stresses and velocities become leading to the time discretized equation of Equation (9): Hereafter, the subscripts n and n+1 denote the known and unknown states and Δt is the time step. By introducing the new variable , the equations are rearranged as below in which F I G U R E 1 Rheological interpretation of the adopted Rayleigh damping model. 43 E is the elastic modulus.
The velocity v n+1 is updated using the incremental displacement Δu at the new state according to the time discretization of velocities shown in Equation (14): The constants 1 and 2 are both set as 0.5 in the simulations to ensure an implicit scheme of a second-order accuracy in time.

GENERALIZED HR VARIATIONAL PRINCIPLE
While numerical schemes for finite element analysis are usually derived from the principle of virtual work where displacements are the sole master field, the multi-field variational principles can be adopted for mixed finite element analysis. 47 Based on the Hellinger-Reissner (HR) variational principle, 48 the MP formulations have been developed for elastoplastic, elasto-viscoplastic, and poro-elastoplastic problems. 10,16,19,49,50 It is shown that the dynamic responses of various media can be simulated efficiently and precisely based on the HR variational principle using the adopted implicit time integration scheme. 16,19,50 Moreover, the previously established HR variational principles for various problems can be integrated together to account for problems involving more features, for example, elastoplasticity with contact. 16 To link this work with previous studies, in this section we first introduce the verified HR functionals for static elasticity 8,48 and dynamic elastoplasticity. 10, 16 We then propose a modified HR variational principle for dissipative medium of the Rayleigh damping in dynamics. The proposed HR variational principle is further verified against the governing equations in Section 2 by checking functional derivatives.

Static elasticity
For an elastic solid subject to quasi-static loads, the HR functional is where both the displacement u and the stress are master fields. The variation of the functional, L( , u) = 0, results in a saddle point that renders neither the maximum nor the minimum of the functional. The problem of static elasticity then can be cast as an equivalent min-max optimization problem 8 : whose incremental form is

Dynamic elastoplasticity
It was demonstrated in References 10,16 that the generalized HR functional for dynamic elastoplastic analysis is where is an arbitrarily small constant and is a slack variable. The logarithmic barrier function includes the condition that > 0. This functional can be verified against the governing equations for elastoplastic material illustrated in Section 2.1 by checking functional derivatives. 10,16 The associated min-max optimization problem is

Dynamic elastoplasticity with the Rayleigh damping
To include the damping terms in the governing equations, for example, Equation (9), we construct an HR type functional considering dynamic elastoplasticity with the Rayleigh damping based on (24).
In the above functional, K n+1 and M n+1 are two additional master fields. The equivalence between this functional and the governing equations shown in Section 2 can be proven by taking functional derivatives Equation (27) corresponds to Equations (16) and (17). Equations (28)(29)(30) are related to Equation (7). Note that in Equation (28), the associated-flow rule (F = G) is implicitly included. The details of proof can be found in Reference 10. Equation (31) is the expression of n+1 as shown in Equation (18). Equation (32) . With the expressions v = Δu/Δt and K n+1 = 1/ 1 K from Equations (14) and (19), we obtain SK = K ∇(v) which is Equation (10). In a similar manner, Equation (33) proved to be Equation (11). If the material is elastic (Δλ = 0), Equations (27)(28)(29)(30)(31)(32)(33) correspond to the governing equations for dynamic elastic analysis of dissipative medium characterized by Rayleigh damping (see Section 2.2). If the damping force is negligible ( M and K are very small), then Equations (27)-(33) represent the dynamic elastoplastic analysis (see Section 2.1). 10 Hence, the time-discretized governing equations with the Rayleigh damping shown in Section 2 are recovered.

Mixed triangular element
The finite element approximations are carried out using shape functions for master fields. Following the standard finite element notation, the state variables within each element can be related to those on element nodes via: wherê,û,̂,K, andM are stresses, displacements, dynamic forces, stiffness-damping forces and mass-damping forces at nodes, and N σ , N u , N γ , N K , and N M are the matrix forms of the shape functions. The relationship between strains and nodal displacements is where B u are the partial derivatives of the displacement shape functions with respect to the global coordinates x and y. A mixed isoparametric triangular element (see Figure 2) similar to that proposed in Reference 10 is adopted in this study. Even though the use of quadratic interpolation nodes introduces additional degrees of freedom (compared to linear interpolation nodes), the adopted mixed triangular element has proven to be robust in overcoming the locking issue in elastoplastic computations of challenging problems. 14,19,50 In the dynamic elastoplastic analyses treated in the previous works 14,19,50 , linear interpolation nodes were used for stress , while quadratic interpolation nodes were used for displacement u and inertial force . Given that the two damping terms K and M are related to stiffness and mass components, respectively, in this study six nodes are used for N u , N γ and N M and three integration nodes are used for N σ and N K .

Discretized variational formulation
Using Equations (34) and (35), the functional (26) is discretized as follows: where F I G U R E 2 Quadratic/linear isotropic triangular element utilized in the simulation.
Following Zhang et al., 49 the displacement boundary condition can be imposed by introducing a new variable, known as the nodal reaction forcer n+1 . Thereby, the expression ofL in Equation (36) is modified as follows: where two additional terms (i.e., the underlined) are included for the enforcement of prescribed displacement. The selection matrix, E u , consisting of 0 and 1 is used to indicate the nodes for which the displacementû p is prescribed.
Furthermore, the functional derivatives of Δû andr n+1 are correspondingly obtained as follows: One sees that Equation (44) corresponds to the displacement boundary condition shown in Equation (4).

Reformulated mathematical program
The proposed variational formulation (42) can be remolded as a standard SOCP problem in the form: where x = (x 1 , x 2 , … , x n ) T consists of the field variables, T and y are the matrix and vector related to the linear transformation, p is the vector of coefficients for the objective function, ℵ is the tensor product of the second-order cones such that ℵ = ℵ 1 × · · · × ℵ l . The second-order cones can be in the type of • Rotated quadratic cone The standard SOCP (45) is a minimization of a linear objective function subject to linear constraints and the second-order cones. The functional (42) can be first reformulated as a maximization problem (by solving the minimization part analytically) following the trace of the previous works, 10,16,51 that is: where the equilibrium condition is now a linear equality constraint. The yield function F (̂n +1 ) ≤ 0 can be formulated as a conic constraint 13,18,19,52 according to the expression of different types of yield criteria. Specifically, here we present the implementation of the plane strain Mohr-Coulomb criterion that will be used in the following numerical example. Similar steps can be applied to other types of yield functions (e.g., Drucker-Prager, Von-Mises, and multi-surface criteria) as shown in References 13,18,19,52 and references therein.
The plane strain Mohr-Coulomb yield criterion at each stress interpolation point is where is the friction angle and c is the cohesion. Equation (49) is an inequality related to the stress components and material strength. By introducing a new variable and linear equality between and , the yield function can be cast in terms of a quadratic cone 10 : One can prove that Equation (50) where N G is the number of Gauss points. By doing so, the quadratic objective function in (48) is converted into a linear function as shown in (51), with the help of four auxiliary scalars (s 1 , s 2 , s 3 , and s 4 ). The minimization problem (51) is now exactly the standard form shown in (45) and can be solved using the interior-point algorithms available in optimization engines. Algorithm 1 illustrates the solution steps of the proposed numerical scheme. A reference of the implementation including the construction of the optimization problem is introduced in Reference 53 in detail.

Algorithm 1. Solution scheme
For the ith incremental analysis: 1. Update history variables such as the velocity, the acceleration and the stress; 2. Form global matrices A, B, C, D, D 2 , H, and d according to Equations (37)- (41) and Equation (50); 3. Define the index matrix E u and prescribed displacementû p ; 4. Prescribe the traction force t and calculate f via Equation (41); 5. Construct the corresponding T and y and p shown in (45) according to (51); 6. Define the conic constraints in (51) using the standard cones (46) and (47); 7. Submit the SOCP program to the optimization solver; 8. Extract/update the variables using the solution returned from the solver, for instance: (a) find the stresŝn +1 in the solution of x; (b) calculate the velocityv n+1 using Equation (20) with Δû that is contained in the dual variables of x; (c) update the strain using Equations (2) and (6). Return to step 1 with i = i + 1.

NUMERICAL EXAMPLES
In this section, we apply the proposed approach to dynamic elastic analysis with damping (Sections 5.1-5.3), stress wave propagation in a cylinder (Section 5.4) and dynamic elastoplastic analysis with a damping layer (Section 5.5). All simulations are performed on a ThinkPad laptop with 2.90 GHz CPU (i7-7820 HQ) and 16.0 GB memory in Windows 10 system. We construct the optimization problem (51) using MATLAB and the problem is solved by MOSEK Optimization Toolbox for MATLAB.

1D P-wave propagation in Rayleigh damping material
The 1D P-wave propagation in a homogeneous Rayleigh damping material is investigated here to verify the proposed variational formulation and its implementation. The benchmark was studied numerically and analytically in Zafati et al., 43 where the Ricker source wavelet is prescribed at the left end of a bar shown in Figure 3. In Equation (52), A p is the amplitude, t is time, f p is the predominant frequency and t s is the time delay. The adopted Ric parameters are illustrated in Figure 3. The damping coefficient M and K are computed according to Equation (12) and the values of are: 0 = 0, 1 = 0.5, and 2 = 1. For the case of 0 = 0 (undamped condition), M and K are set as sufficiently small constants to mimic an elastic medium. The comparison between numerical results (element edge length h = 1 m) and the analytical solutions are presented in terms of the normalized displacement u in Figure 4, where great agreements are obtained.
To investigate the convergence of the method with respect to mesh size, the mixed triangular elements with edge length h being 1, 10, and 16 m are used, leading to grids counting as many as of 1861, 256, and 118, respectively. As shown in Figure 5, the simulation results agree with each other very well. Note that the element size is 1/240, 1/24, and 1/15 of the wavelength computed from the predominant frequency f p and the velocity of the P-waves. For each simulation, a total of 1500 incremental steps are used to simulate the wave propagation in 15 s. The averaged computational costs of simulations with different damping ratios for one time step analysis are 3.98 s, 0.49 s, and 0.2 s for element size of 1, 10 and 16 m, respectively.

1D P-wave propagation in a continuously absorbing medium
A purely transparent open boundary, where no reflection can take place and the wave is totally transmitted, is hard to implement in finite domains since residual reflection usually remains. A strategy to reduce the undesired reflected wave is  Figure 3 using 0 = 0, 1 = 0.5, and 2 = 1, compared to analytical solutions. 43 to extend the computational domain beyond the boundary where one places a medium of continuously changing damping coefficients able to absorb the incoming wave, as proposed in References 40,41 To illustrate this strategy, we take an elastic bar with the same geometry and subjected to the same propagating wave as in the previous example. However, the numerical simulations differ since we add an absorbing layer to the right end of the elastic medium. The original boundary becomes therefore an inner interface separating media with elastic and inelastic properties. The additional absorbing medium is described through a continuously damping ratio, i.e. (x). Making reference to Figure 6, it is assumed that the function (x) increases linearly from the interface, where it is zero, to the right end, that is the new right end boundary, where displacement is imposed to be zero, as proposed in Rajagopal et al. 41 The expression for the damping ratio can be given as follows:

F I G U R E 4 Time history of displacement at point C of
The parameter max is computed according to Equation (12). In this example three different values of max are adopted: max = 0.5, max = 1 and max = 2. The reflection wave is quantified by ref , which is the percentage between the amplitudes of the reflected and the incident waves. Numerical results obtained by the devised approach are displayed in Figure 7, in which the time is normalized by the period of Ric and the displacement is normalized by the maximum displacement of Ric. The edge length of the element is set as 10 m, and a time step of Δt = 0.01 s is adopted in the simulations of 40 s wave propagation.
In general, a one-wavelength long inelastic layer is needed to keep ref close to or under 5%. Comparing Figure 7B, C, one sees that the use of a strong damping ratio increases the reflection wave due to the increased change of damping property. The choice of the thickness layer, the damping coefficient and k(x) has also been investigated in F I G U R E 5 Time history of displacement at point C of Figure 3 with three element sizes.

F I G U R E 6 1D P-wave propagation in an elastic bar, where the rhs boundary is transformed into an inner interface. Beyond it, a
continuously absorbing layer is set: k(x) is the linear function computed from the interface (k = 0) to the new right end (k = 1) and A 0 is the coefficient for the length of the absorbing layer with A 0 being 0.5, 1 and 1.5. Point B is the monitoring point in the domain.
References 39-41, where some optimized manners are suggested. For this example, the best results are the ones displayed in Figure 7B. This example serves as a further validation of the proposed mathematical formulation given the demonstrated capability of handling a 1D continuously absorbing medium and mimicking a purely transmitting boundary. The averaged computational costs of different damping ratios for one time step analysis are 0.27 s, 0.37 s, and 0.49 s for A 0 = 0.5, 1, and 1.5, respectively.
A more sophisticated way of absorbing a wave across an interface is the multi-layer strategy, 43 where the CALM can be further modified by using prescribed layer-dependent expressions for the elastic modulus E. The absorbing layer consists of a series of N e absorbing sublayers (i) that are characterized by the damping ratio i . The elastic modulus E i of each layer i is designed by the following expressions: in which E 0 is the elastic modulus of the elastic medium closest to the interface of . In the simulation, each element inside the damping region is chosen as a sublayer with the designed elastic modulus and damping ratio. The length of the absorbing layer is chosen as one wavelength A 0 = 1. We plot the time history of displacement in the left panel of Figure 8, where the results from optimized and unoptimized simulations are nearly the same. We then display the reflection in the right panel of Figure 8, which shows the improvements with the use of such a multi-layer strategy with specified elastic moduli under different values of max . The averaged computational cost of different damping ratios for one time step analysis is 0.38 s.

2D Lamb's problem
The Lamb's problem is a typical analytical benchmark in seismic problems. It has been treated as a classical problem in the validation of numerical methods of wave propagation. Lamb tests study the response of a given point of a homogeneous elastic half-space that is excited by a prescribed surface vertical and/or horizontal load. The proposed numerical approach with the CALM strategy applied in 1D cases is used here in 2D to investigate its performance in a more general case, where P-wave, S-wave, and surface waves are involved. A vertical load of the same Ricker wavelet is applied on the surface with a monitoring point C located on its right at a 0.5 distance. As displayed in Figure 9, the concerned area is a square with an edge length of 2 and is embedded in an absorbing layer ( max = 1) with thickness set to λ. The parameters are: V P = 80 m/s, V S = 56.6 m/s and = 1700 kg/m 3 . We consider two mesh configurations: (i) a dense mesh (MESH1) consists of 67,050 triangular elements with an edge length is defined for all points in the damping region, with x i being the smallest Euclidean distance of the given point to the interface. We assume that the maximum value of the adopted damping ratio is equal to max , so we have the expression for k(x i ): with the distributions of E i and i obtained from Equation (54). The value of k(x i ) is illustrated in Figure 9A with colors of the damping layer. The vertical Ric displacement is applied at the source, and the time period of the simulation is 25 s with Δt = 0.1 s. As shown in Figure 10, the reflection wave is under 2%, and the accuracy of the presented method with embedding damping region can be regarded as satisfactory for seismic wave modeling in a bounded domain. The obtained numerical results (red and green dashed lines in Figure 10) agree well with each other and also with the one (black dashed line) from the finite difference method. 54 The scattered fields of the normalized displacement at four instants (MESH1) are displayed in Figure 11 which confirms the good performance of the method. The average computational costs for one incremental step analysis of MESH1 and MESH2 are 3.3 min and 21 s, respectively, and no obvious differences are observed using different time steps.

A thick-walled cylinder subjected to internal pressure
To further demonstrate the capability of the proposed approach for the dynamic elastoplastic analysis of engineering problems, we consider a thick-walled cylinder subjected to suddenly applied internal pressure. 55,56 The material is assumed to be elastoplastic with von Mises yield criterion under plane-strain conditions. The material properties are as follows: the elastic modulus E = 21 × 10 4 MPa, Poisson's ratio = 0.3, density = 7850 kg/m 3 , internal pressure p = 185 N/mm 2 and yield stress Y = 355 MPa. Due to the symmetry, only a quarter of the cylinder is simulated. Two mesh configurations MESH1 (186 elements) and MESH2 (911 elements) are considered, and the setup of MESH1 is shown in Figure 12A. The time step is Δt = 1 × 10 −5 s, the same as the one adopted in Reference 56. In Yang at al., 56 a scaled boundary finite element method (SBFEM) is proposed to solve this problem using 158 polygons and the results were compared with the ABAQUS simulation result using a very fine mesh (i.e., 24,313 linear elements). As shown in Figure 12B, the numerical results obtained from MESH1 and MESH2 in this study agree slightly better with the ABAQUS result using a very fine mesh, compared to the SBFEM results. The increase of the number of elements, F I G U R E 10 Vertical and horizontal displacements at Point C of Figure 9. The normalized displacements u x and u y are calculated by displacements normalized to the corresponding maximum displacements at C, respectively. ref x and ref y denote the reflection percentage of horizontal and vertical displacements, respectively. The resulting reflection wave is confined within the 2% shadow band. For the FDM with PML, the reflection wave is at a lower magnitude and it is confined within 1%.

F I G U R E 11
Normalized displacement u field at four instants (t = 5, 10, 15 and 25 s). Normalization is computed with respect to the amplitude of the incident vertical Ric wavelet. The area of study is delimited by the black dashed line. that is, from MESH1 to MESH2, also slightly improves the results. Therefore, we show that the adopted method is an efficient approach for elastoplastic analysis. The average computational costs for one time step analysis of MESH1 and MESH2 are 0.74 and 4.43 s, respectively.

Strip footing
In this part, we apply the proposed numerical scheme to the dynamic analysis of elastoplastic soils. Its performance is investigated through a classical soil mechanics problem, i.e. strip footing. Similar to the Lamb's problem in Section 5.3, the problem domain is subjected to a vertical load. We choose the same damping layer thickness with optimized elastic modulus (see Equations (54) and (55)), elastic parameters, and Ricker wavelet as in Section 5.3, while the computational domain consists of two parts, that is, elastoplastic material and damping layer. The plasticity of the elastoplastic material is described by the Tresca criterion ( = 0 • in Equation (49)) with cohesion being 20 kPa. The damping layer is modeled as an elastic material characterized by Rayleigh damping. As shown in Figure 13A, the computational domain is discretized using 2297 triangular elements and the elastoplastic region consists of 1710 elements. The width of the vertical loading is 10 m, and we refine a square region with side length being 10 m. The amplitude of vertical loading is 0.05 m, and we run 30 s simulation using a fixed time step Δt = 0.1 s. The x-direction displacement along the left boundary is constrained, while the right and the bottom boundaries are fixed. We consider two cases (i.e., with damping and without damping) to present the performance of damping layer (coarsely discretized using 587 elements) in the study. The numerical solution obtained from a large physical domain without damping (4600 m × 3800 m) is used here as a reference solution shown in Figure 13B. We consider five values for the maximum damping ratio and max = 0 represents the case without damping. We plot the normalized reaction force (sum of the vertical nodal forces at the loading nodes, normalized with respect to the maximum absolute value of the reference solution) in Figure 13B. As seen, the reflection wave induces large fluctuations in the time history of the force if damping is not applied. With a large damping ratio ( max = 4 and max = 8), the final reaction force is nearly the same as the reference solution. Note that these two values also increase the wave reflection during a time period from t = 4 s to t = 5.5 s, as can be seen in the zoomed area. Even though max = 1 also provides good performance in absorbing waves, the optimal value here is max = 4. This has also been investigated in Reference 39 that the generic rules can be applied to design good absorbing layers, while the optimization of parameters of absorbing layers depends on the application and the property of signals. In the most recent studies, [39][40][41]43,44 only the elastic wave propagation is considered. Our simulation results show that the value of damping ratio may need to be further optimized for some plasticity problems. The guideline of choosing the maximum damping ratio max is beyond the scope of this study, while we will try to investigate it in future works when we apply this model to different applications. To further present the performance of the absorbing layer, we also plot the incremental displacement field in Figure 13C, comparing the cases with and without damping. The distribution of the equivalent plastic strain (PEEQ) is shown in Figure 13D, which is in line with the one from the typical elastoplastic analysis of strip footing. The average computational cost for one time step analysis is 12.1 s.

CONCLUSIONS
An MP formulation for dynamic elastoplastic analysis problems with absorbing layers is proposed to mimic an unbounded domain. The stress wave propagation in the absorbing layers (i.e., dissipative media) is characterized by Rayleigh damping. Specifically, based on the strong form of the governing equations, the mass and stiffness components of the damping terms are incorporated into a generalized Hellinger-Reissner (HR) functional by introducing two additional master variables. The proposed functional falls in a min-max optimization problem that can be reformulated as a standard SOCP problem and solved using the interior-point method available in optimization solvers, which is of high efficiency and accuracy. The validation of the proposed mathematical formulation is carried out by comparing our numerical results to the analytical solution for 1D P-wave propagation in a dissipative medium. Its efficiency for modeling seismic wave propagation with absorbing boundary conditions by using increasing-damping layers to mimic an unbounded medium in both 1D and 2D cases is also verified. Further, the stress wave propagation in elastoplastic medium is performed to show the efficiency and accuracy of this method by comparing the results to other published numerical data. Finally, the dynamic elastoplastic analysis of a strip footing problem with the damping layer in absorbing the reflection waves is simulated to show the robustness of the proposed method and also to investigate the performance of the absorbing layer. We show that the developed approach can be used not only for seismic wave propagation in geological engineering but also for elastoplastic analysis of geotechnical problems with stress wave propagation and unbounded domain. Moreover, the developed mathematical formulation can be further integrated into such MP-based approaches, facilitating the use of interior-point algorithms in analyzing different problems with Rayleigh damping.

ACKNOWLEDGMENTS
The author Liang Wang greatly appreciates the financial support from the cooperation agreement between the University of Bologna and the China Scholarship Council (No. 201606060161). The author Liang Wang thanks the members of the computational seismology group led by Prof. Peter Moczo for the discussions and suggestions during his visiting period. The authors thank Dr. David Gregor for providing FDM results shown in Figure 10.