Harmonic-Modal Hybrid Reduced Order Model for the Efﬁcient Integration of Non-Linear Soil Dynamics

: Nonlinear behavior of soils during a seismic event has a predominant role in current site response analysis. Soil response analysis, and more concretely laboratory data, indicate that the stress-strain relationship of soils is nonlinear and exhibits hysteresis. An equivalent linearization method, in which non-linear characteristics of shear modulus and damping factor of soils are modeled as equivalent linear relations of the shear strain is usually applied, but this assumption, however, may lead to a conservative approach of the seismic design. In this paper, we propose an alternative analysis formulation, able to address forced response simulation of soils exhibiting their characteristic nonlinear behavior. The proposed approach combines ingredients of modal and harmonic analyses enabling efﬁcient time-integration of nonlinear soil behaviors based on the ofﬂine construction of a dynamic response parametric solution by using Proper Generalized Decomposition (PGD)-based model order reduction technique.


Introduction
Structural solid dynamics is usually formulated either in the time or in the frequency domains [1]. When considering a general linear viscoelastic behavior, momentum balance leads to the semi-discretized form where M, C, and K are, respectively, the mass, damping, and stiffness matrices; U the vector that contains the nodal displacements; and F the nodal excitations (forces). From now on, we assume the mechanical system described by N nodal degrees of freedom, which gives the size of the different matrices and vectors involved in Equation (1). The direct integration of the previous discrete form, which consists of N second order coupled ordinary differential equations, can be performed using either, explicit or implicit time integration schemes.
Addressing fast transient dynamics can be usually accomplished using explicit integrations that require satisfying stability conditions affecting the largest integration time-step ∆t max , closely related to the size of the elements involved in the mesh M covering the domain Ω. Many robust integration schemas exist and are widely employed, from the very popular Newmark [2] to more advanced schemes able to preserve inherent mechanical properties [3].
In many cases, explicit integrations are combined with mass-lumping to recover a diagonal mass matrix, in which inverse becomes trivial, fact that improves significantly the integration efficiency. In the same spirit, modal analysis offers an appealing alternative route. The last is based on the fact that the properties (positivity and symmetry) of mass and stiffness matrices in Equation (1), imply the existence of a basis {φ 1 , φ 2 , . . . , φ N } (in which construction will be detailed later), such that the two following relations apply: with δ the Kroenecker delta, and By defining matrix P, in which columns correspond to vectors φ i , the two previous relations can be rewritten in their more compact matrix counterparts and P T KP = K, (5) I is the identity matrix, and K is the diagonal matrix with entries K ii = κ i . Thus, if C = 0 (undamped dynamics) or C = αM + βK (proportional damping), using the coordinates ϕ, U = Pϕ, the problem (1) can be rewritten as with C diagonal, i.e., C ij = c i δ ij (C = 0 in the undamped case). Equation (6) constitutes a system of N uncoupled second order ordinary differential equations that can be integrated very efficiently for providing transient responses.
Another alternative route consists in considering the frequency instead of the time as coordinate, which is applying the Fourier transform to problem (1), which constitutes a choice particularly appealing in the case of forced regimes and linear models. In the case of nonlinear models, such a route needs appropriate procedures detailed later.
The frequency-based counterpart of Equation (1), constituting the so-called harmonic formulation, reads −ω 2 M + iωC where U and F refer to the Fourier transform of the nodal displacement and forces, respectively, U and F. As indicated, the main limitations of such a formulation are: (i) first, the treatment of nonlinear models that requires specific procedures; and second (ii) the necessity of solving the problem for each angular frequency ω involved by the loading or present in the induced displacement (when nonlinearities are addressed). However, the use of harmonic analysis avoids the necessity of constructing the basis P that needs the solution of the large-size eigenproblem related to −ω 2 M + K φ = 0, (8) in which eigenvectors φ i diagonalize the mass and stiffness matrix, as well as the one related to the damping as soon as a proportional damping is assumed. The main limitation of modal analysis is not only the necessity of solving such an eigenproblem for extracting the normal modes but the lack of validity of such basis in the case of parametric models in which the different matrices change when modifying the parameters grouped in vector µ, i.e., M(µ) and K(µ). The solution of parametric eigenproblems remains a recurrent and still open problem.
Another advantage of the harmonic formulation is the possibility of considering more general damping. In mechanics of materials, sometimes damping is assumed scaling with the inverse of frequency [4]. The interested reader can refer to Reference [5] that analyzed the theoretical consequences of assuming a frequency dependent dashpot parameter. In Reference [5], it was proved that, even if such a choice succeeded to fit experimental data, when coming back to the time-domain, causality is lost, and then the resulting expressions in the time-domain were called non-equations. Anyway, it is important to note that, even when considering complex nonlinear frequency dependent damping C(ω) in Equation (7), the problem in the frequency domain remains linear because, here, the frequency is a model parameter (or a model extra-coordinate within the Proper Generalized Decomposition (PGD) framework addressed later).
In the context of model order reduction (MOR), in Reference [6], authors proposed a Proper Orthogonal Decomposition (POD)-based reduced order modeling operating in the time domain. Ladeveze and coworkers proposed an extension of their radial approximation [7] for addressing mid-frequency dynamics, the so called variational theory of complex rays [8].
In our former works, we considered a PGD (Proper Generalized Decomposition) formulation for constructing a parametric transfer function [9] that allowed efficient solutions of transient dynamics operating in the time-domain. On the other hand, the separation of variables, at the heart of PGD [10,11], was extensively employed in the harmonic domain for solving multi-parametric dynamics, and was successfully extended to the non-linear case, and then combined with modal analysis [12].
In soil dynamics, the strain-stress relationship indicates that soil behavior during a seismic event is highly nonlinear and exhibits hysteresis [13,14]. Nonlinear effects, are defined in terms of the shear modulus reduction and an increase of damping ratio as a function of the strain level. Various approaches are available to model this dependence, such as equivalent-linear (EL) or nonlinear (NL) analysis. Equivalent linear analysis is by far the most commonly utilized procedure in practice, however, may lead to a conservative approach [15]. Therefore, a better characterization of soil behavior requires investigating the influences of non-linearity in the response analysis of the site using an efficient numerical model to address more precise behavior models (nonlinear) without compromising the computational solution efficiency.
Even if, as indicated above, nonlinear visco-elasto-dynamics can be addressed under the stringent real-time constraint by combining modal and harmonic analysis [16], which procedure fails when addressing a stratified media with different material properties to be considered parametrically. On the other hand, a fully parametric PGD-based solution could be envisaged [17], but its complexity scales with the size of the problem, making difficult the construction of that parametric solution. This paper aimed at proposing, testing, and validating a new procedure embracing the two just referred, where the harmonic formulation is combined with a truncated (reduced) modal basis, for calculating in real-time the response of the structural system to any loading in any soil.
For that purpose, the present work considers the hybrid framework, consisting of an harmonic formulation operating in the modal representation. To circumvent the issue related to the parametric dependence of the modal basis, as well as the issue related to the nonlinear behaviors, we propose a linearization and a way of addressing the parametric behavior, which guarantee the invariability of the modal basis. Moreover, in order to alleviate the issue of addressing all the frequencies that the solutions could involve, the angular frequency ω is introduced as an extra-coordinate in the solution representation.
This new numerical technique allows obtaining responses of nonlinear models defined in stratified media (as it is the case in soils) in real-time. As the response to any loading in any parametrized medium can be explored very fast, in almost real-time, as just indicated, its application in optimization and reliability analyses, as well as in structural monitoring and control, could be envisaged.
In the present paper, after describing the main concepts of the proposed methodology in Section 2, Section 3 introduces the nonlinear soil-dynamics problem that is solved in Section 4 under the stringent real-time constraints and discussed accordingly.

Extended Harmonic-Modal Hybrid (XHMH) Approach
This section revisits the hybrid formulation [12,16] that will be then applied for solving efficiently nonlinear dynamics of soils. For that purpose, we consider the harmonic formulation with a proportional damping that using the associated modal basis, with the normal modes grouped in matrix P, using variables U = Pξ and taking into account relations (4) and (5) reads that results in a system of N decoupled algebraic equations with i = 1, 2, . . . , N, in which the solution reads From ξ(ω) calculated with (12), the nodal amplitudes results U (ω) = Pξ(ω), which allows computing the nodal displacements U(t) by using the inverse Fourier transform.

Nonlinear Dynamics
Many engineering applications cannot be modeled as a linear system. In soil dynamics, for example, the relation between stress and strain is highly nonlinear and hysteretic; thus, an accurate nonlinear modeling, addressing the soil stress-strain behavior is required.
In the nonlinear case, the general semi-discretized momentum balance can be written as where R(U) is a nonlinear contribution that depends on U(t).
The simplest linearization compatible with the just described hybrid framework consists in calculating the nonlinear term R(U) from the solution at the previous iteration of the nonlinear solution procedure and considering it as an extra-loading, i.e., where U − (t) correspond to the time-dependent solution at the previous iteration. Now, all the rationale previously introduced applies, with F being now the Fourier transform of R(U − (t)) + F(t).
To improve the numerical efficiency, we consider the expression of the loading in the frequency domain, expressed in the modal basis, f(ω) = P T F (ω) and consider an approximation of it in the frequency domain with N l (ω) the approximation functions in the frequency domain. Thus, the problem solution reads The expression above allows the offline calculation of the inverse transform of the approximation functions N l (ω), making possible an online reconstruction in almost real-time [16].

Parametric Damping in Stratified Media
As just discussed, by assuming a proportional damping, the dynamical problem can be fully diagonalized, and, in presence of nonlinearities, the last are linearized and transformed into an effective extra-loading.
However, when addressing the dynamics of stratified soils, the damping of each layer is different and can depend on the mechanical solution itself in the nonlinear case.
For facilitating the computation, one is tempted to consider a parametric damping by decomposing its associated global matrix C as where L is the number of layers composing the stratified media, and ζ 1 , . . . , ζ L the associated damping parameters. The main consequence of such a decomposition is that a fully diagonalization is not possible anymore. Thus, the most immediate route consists in considering the whole damping as an effective loading by evaluating it at the previous iteration solution of the nonlinear solution procedure as discussed in the previous section. However, such a route compromised the convergence in most of the numerically analyzed cases.
For that reason, we consider in what follows the usual harmonic formulation: that as indicated cannot be fully diagonalized. However, in order to compute efficiently the parametric solution U (ω, ζ 1 , . . . , ζ L ) by invoking the PGD [10,11], an appealing alternative consists in using a reduced modal basis, which are the modes φ i related to the highest eigenvectors associated with the eigenproblem (8), grouped in the truncated matrixP, which allows defining the amplitude approximation U ≈Pξ, which introduced into Equation (18) and premultiplying byP T reads

Nonlinear Soil Behavior
It is well established in geotechnical engineering that soil response is nonlinear beyond a certain level of deformation. Stress-strain relationships for the levels of shear deformation produced by large earthquakes are nonlinear and hysteretic, as has been confirmed by numerous results of vibratory and cyclic loading tests on soil samples.
The performance of cyclic nonlinear models can be illustrated by a very simple example in which the shape of the backbone curve is described by τ = f (γ) and depends on the parameters describing the nonlinear behavior of soils [18,19].
Among the hyperbolic type models, the most famous and commonly used model is the one initially proposed by Kondner and Zelasko (KZ) [20], lately redefined by Hardin and Drnevich [21], where τ is the shear stress, γ is the shear strain, G is the undisturbed shear modulus (taken at the origin), and τ max is the shear strength (the maximum stress that material can support in the initial state, τ max = Gγ r ). Then, (21) results, In the numerical tests here addressed, the above equation is rewritten as where γ r is the reference deformation, and c can be considered as a curve fitting parameter (let us note that, when c = 0, the nonlinear case reduces to the linear one). The other parameter related with this nonlinearity is the damping characteristics of the soil (represented by the damping ratio ζ), which is calculated by the ratio of the energy dissipated W D and the one associated to the elastic deformation W E , from the hysteretic loop, i.e., Ishihara proposed, in Reference [22], the following expression of the damping ratio: that using (22) leads to with γ a = max(γ). Now, considering the linearization introduced in the previous section, it results where R(t) account the nonlinear contribution of the constitutive Equation (23), i.e., with τ and γ evaluated at the previous iteration of the nonlinear solver. Equation (27) can be expressed in the form of Equation (20) to compute efficiently the parametric solution. When considering a proportional damping, the damping ratio (26) will naturally appear explicitly in the left-hand side of the equation, which will represent a sort of parametric transfer function. In the proposed approach, the iterative procedure is initialized from the small-strain shear modulus and damping. Then, damping is updated at each iteration step, in the left-hand side, from Equation (26). This method fits perfectly with any damping model proposed in the literature. As soon as such parametric solutions are available, the non-linear problem can be solved in real time because no new calculation is needed; the non-linear solver only needs to particularize online the parametric solution calculated offline.

Numerical Results
In this section, a numerical test for illustrating the potentialities of the technique just proposed is addressed. Here, we consider a soil deposit consisting of 5 strata. The PGD method was used to calculate the parametric solution for the displacement fieldξ(ω, ζ 1 , ζ 2 , ζ 3 , ζ 4 , ζ 5 ), as detailed in Appendix A.
The different domains were discretized by considering 100 nodes for the spatial coordinate, 1023 nodes required for Fast Fourier Transform (FFT), and 100 nodes for the damping coordinates. Even if 100 nodes for discretizing the parametric domain seem too much, as the calculation of functions depending on the parameters does not imply the solution of any linear system, it is preferable to consider a rich-enough discretization to be sure of representing accurately the parametric solution.
Different methodologies were taken into account in the numerical simulations: (i) the Equivalent Linear Response Analysis (EERA); (ii) the Extended Harmonic Modal Hybrid Approach (XHMH) proposed before; and (iii) the nonlinear procedure used in DEEPSOIL software. In the numerical simulations, it was considered 5 layers and a half-space, with the properties listed in Table 1. The Equivalent Linear analysis was carried out using the appropriate material curves given by Reference [23], and, for DEEPSOIL analysis, it was considered the modified hyperbolic model MKZ proposed by Reference [20], as with β and s the parameters that adjust the shape of the backbone curve. The different approaches considered an initial estimation of damping ratio ζ 0 (assuming a linear regime, which are small strains); it was considered for XHMH analysis c = 1, and ,for DEEPSOIL analysis, β = 0.9, s = 0.85, and σ re f = 0.18 for each stratum, respectively. Table 1. Soil properties, with Gmax the initial shear modulus, ζ 0 (%) the damping ratio at small strains, γ r (%) the reference strain, and σ v the effective vertical stress.

Stratum
Material Thickness (m) Density (kg/m 3 ) Gmax (MPa) ζ 0 (%) γ r (%) σ v (kPa) In soil analysis, input motion is defined from the response spectrum, or its corresponding time history at bedrock surface (outcrop), known as rock outcropping motion. From the outcropping motion, the objective is to predict the bedrock motion covered by the soil deposit. Thus, the bedrock half-space can be substituted with boundary condition whereU s the rock outcropping velocity (assumed measurable), andU * o the velocity at the base of the soil column [17,24], which coincides with the soil-bedrock half-space interface and τ * o the shear stress at that position. The solution is obtained assuming the acceleration time history in the outcropping is known; in our simulation, this is the one shown in Figure 1 (recorded at rock site).    Figures 3 and 4 illustrate the time evolution of the acceleration and shear strain, respectively, at 3 different layers (1, 3, and 5), taking into account different methodologies. In Figure 3, the acceleration predictions from non-linear methods show similar behaviors, but the maximum acceleration predictions, especially at the surface, are significantly different. In this case, the proposed approach is closer to the EERA solution. Figure 4 the shear strain time history from nonlinear methods are reasonably close for most cases compared to EERA analysis, but the peak strain predictions are significantly different. As previously discussed the XMHM addresses the real nonlinear problem, while EERA performs on an equivalent linear behavior. The differences between both nonlinear solutions are due to the different soil properties considered and also to the different mesh resolution (being finer the one associated with the XMHM). The XMHM solution was validated by performing a standard time-integration of the considered nonlinear model. Figure 5 shows the modulus reduction, damping, and shear stress curves for the three different methodologies. It can be noticed that the results for the modulus reduction curve are similar. The damping and shear stress curves exhibit differences, particularly at high strains, depending on the considered model, which involve slightly different damping modeling. Figure 6 shows PGA, the maximum displacements at the top of each layer and maximum strain at the mid-depth of each layer for the three considered solution techniques. Certain differences can be observed for the PGA values between the DEEPSOIL nonlinear model and the XHMH model, however, the maximum displacements are nearly exact matches. For the maximum values of shear strain, the three approaches show variations. These variations may be due to the different damping models for each method or the selection of curve fit values.  Finally, it can be observed in Figures 3-6, which despite the differences that the three methods exhibit, the global behavior remains similar. The computational methodology here proposed, which is the truncated XHMH technique, proceeds extremely fast in the nonlinear case, even when addressing a stratified soil with different material properties, enabling the response evaluation instantaneously (real-time). The difference between the methods is the way to address the nonlinear behaviors. The proposed approach does not pretend to discuss on the different models' accuracy. It simply aims at proving that general nonlinear soil models can be solved with the same computational performances than linear models.

Conclusions
This paper proposed a new dynamic calculation of nonlinear soil behavior based on a hybrid technique able to compute very fast solutions. The approach combines a modal analysis, a harmonic space-frequency description of the dynamic problem, the introduction of material parameters as model extra-coordinates, and an online integration that proceeds by particularizing the parametric solution for the damping parameters, and then updating the level of deformation from the just calculated solution.
The performances of this new approach enable real-time computations because of the fact that it only needs particularizing an already computed parametric solution instead of performing a new resolution as usual solution procedures.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The Proper Generalized Decomposition Applied to the Calculation of the Parametric Solution
In this section, we describe the separated representation construction of, for the sake of simplicity, a single-stratum parametric displacement. For that purpose, we consider the generic formulation −ω 2 M + iωζC + K ξ =f, and look for a separated representatioñ Now, considering the weighted residual form where Ω = Ω ω × Ω ζ , the m-rank PGD approximation reads with the test functionξ * expressed from where X k is the vector of nodal displacement, and W k (ω) and M k (ζ) are the functions related with the frequency (ω) and the damping (ζ). The introduction of Equations (A4) and (A5) into (A3) results in a non-linear problem. We proceed by considering the simplest linearization strategy, an alternated directions fixed point algorithm. The three problems to be solved for calculating the functions X m , W m (ω), and M m (ζ) are described below: Now, integrals in Ω ω × Ω ζ can be calculated: Substituting Equation (A8) into Equation (A7), the last reduces to