Mathematical modeling of the dynamics of layered and block media with nonlinear contact conditions on supercomputers

Continual models of solid media with a discrete set of slip planes (layered, block media) and with nonlinear slip conditions at the contact boundaries of structural elements are constructed. The constitutive equations of the resulting system of equations contain a small viscosity parameter in the denominator of nonlinear free terms. For a stable numerical solution of a system of differential equations, an explicit-implicit method is proposed with an explicit approximation of the equations of motion and an implicit approximation of the constitutive relations containing a small parameter. From implicit nonlinear difference approximations analytically, using the perturbation method, various effective formulas for the correction of stress components after an elastic time step are obtained. To calculate the elastic step, we used a grid-characteristic method on hexahedral grids, which allowed us to significantly increase the speed of calculations and simulate a non-stationary three-dimensional problem of generating a response from an oriented layered or block cracked cluster located in a homogeneous medium.


Introduction
Continual models of deformable solid media with a discrete set of slip planes (layered, block media) and with nonlinear slip conditions at contact boundaries can be obtained with the method of the asymptotic homogenization [1] or with the discrete variant of the slippage theory [2]. In all of these cases, in the constitutive system of equations with the non-linear free term and the small stress relaxation time are included. For the stable numerical solution of differential equations, the explicit-implicit method with the explicit approximation of motion equations and the implicit approximation of constitutive equations containing small parameter in the denominator of the free term was proposed. A set of effective formulas for the stress tensor correcting at the elastic step was obtained. Numerical simulations of the dynamic scattering process for subsurface layered and block objects in elastic 2D and 3D media were carried out using modern high-performance computing systems.
Here k γ k ω -coefficients of weak elastic tangential and normal bond of layers, η the viscosity coefficient, q -the Coulomb friction coefficient, ⟨F (y)⟩ = F (y)H(y) , H(y) the Heaviside function, H(y) = 0 if y < 0 , H(y) = 1 if y ≥ 0 . The contact plane with the interaction condition defined is called the slip-delamination plane.
It should be noticed, that previously in paper [2] for the delamination regime the full delamination condition was used :Ω ≥ 0 : τ = σ n = 0, that is the asymptotic case of the weak delamination when k ω → 0, k γ → 0 ; Ω = [u n ]/ε the normalized discontinuity of normal displacements at the contact boundary, defining by the equationΩ = ω . Also, for the case of the weak delamination ω =σ n /k ω inequalities σ n ≥ 0 and Ω ≥ 0 are equivalent.
From the numerical point of view, small parameters k γ k ω are regularizations that allow us to prevent the oscillations occurrence when sharp changing the compress boundary to the full delamination condition.
To construct the continuum model with a set of these slip-delamination planes we are going to deal with γ and ω as discontinuous functions of space and time. Also we will use main relationships from the slippage theory as many other authors. It allows us to take into account contributions from γ and ω into non-elastic strain tensors e γ and e ω respectively: e γ = (n ⊗ γ + γ ⊗ n)/2, γ · n = 0, e ω = (n ⊗ ω + ω ⊗ n)/2 = ωn ⊗ n, ω = ωn The full strain tensor equals to the sum of elastic and non-elastic parts: e = e e + e γ + e ω , e = (∇v + ∇v T )/2 Here v is macroscopic velocity of medium particles, e e is the elastic strain tensor accordingly to the Hookes law:σ = λ(e e : I)I + 2µe e The final equation is the motion equation:

System of equations for a layered medium
In layered medium, containing a set of elastic layers, the only system of slip-delamination planes are possible with the normal n. If the normal to contact boundaries n is oriented along the coordinate x 2 then its components are n j = δ 2 j . Through conditions for γ and ω, corresponding to the local contact conditions 1), are: Through conditions for γ and ω, corresponding to the local contact conditions 2), are: Based on the chosen normal direction the final system of equations for this model is:

System of equations for a block medium
The block medium consists of parallelepiped elastic elements with three possible slipdelamination planes. These planes are defined with normals n (s) , s = 1, 2, 3. In this case the non-elastic strain tensors are: If three normals to slip-delamination planes are oriented along coordinate axis of the Cartesian system then n (s) j = δ s j , where δ s j is Kronecker's symbol. Through conditions for γ (i) and ω (i) , corresponding to the local contact conditions 1), are: Through conditions for γ (i) and ω (i) , corresponding to the local contact conditions 2), are: As for the layered medium, we can rewrite the main system in the suitable form:

Numerical method
Both formulated systems are semi-linear hyperbolic systems and the numerical solution can be obtained with different explicit schemes. However, the slippage process switches on the nonlinear free term with small viscosity parameter in the denominator. The system transforms into the form with small parameter and ordinary explicit schemes will not be stable. To overcome this problem, the usage of explicit-implicit method is proposed. The implicit approximation is used only for equations that contain small term in the denominator. All other equations are approximated with the explicit scheme.
Lets describe this approach forσ 2j for compressed contact case σ 22 < 0 for the layered medium:σ Implicit approximation with first order time approximation is: Here indices n + 1 and n mean current and next time layers, ∆t is the time step. We assume, that values v n+1 i σ n+1 jj were calculated with the explicit schemes for elastic equations. Solving this equation and also the same equation for normal tensions we can obtain correcting formulas. When σ n+1 When σ n+1 22e ≥ 0 (delamination process) Coefficients for weak shear α = k γ /µ and stretching β = k ω /(λ + 2µ), α, β << 1, ∆t is the stress value after the elastic step. These formulas, in fact, is the adjustment of elastic stresses for the friction cone with viscous corrections. We used here the simplest formula for σ 2j n+1 e only to illustrate the derivation process.
In this research we used the grid-characteristic method on hexahedral meshes with the third order of the spatial approximation [3,4,5]. The monotonization was done based on the gridcharacteristic monotonicity criterion, based on the internal structure of the linear transport equation. It decreases the approximation of the method to the second order on discontinuities. It allows us to increase significantly the precision and the robustness of the numerical simulation. The seismic response from the cluster of cracks and blocks in the homogenous medium was calculated in 2D and 3D cases.

Simulation results
On Fig. 1 two 2D wave fields of scattered seismic waves on the layered structure with different orientations for contact condition 2) are shown. The influence of the layer orientation on the scattered wave field is clearly seen. Seismic response is non-symmetric for non-vertical layers orientation. Based on this fact it is possible to make some assumptions about the main direction in the subsurface fractured medium. On Fig. 2-3 3D wave fields of scattered seismic waves on layered and block structures for contact condition 1) are shown. The whole medium was a parallelepiped with sizes 10 x 10 x 3 km. P-wave velocity was 4500 m/s, S-wave velocity 2250 m/s, density 2500 kg/m3. The day surface condition was used on the upper side. In the center of the medium along OX and OY with the depth 50 m the point source with the 30 Hz Ricker time dependence function was applied. At the 2 km depth the object was set with sizes 3 x 3 x 0.2 km. Two different models were compared: block and layered with the normal vector along the axis OX. The parallelepiped mesh with 5 m cells, containing 2.4 billion of nodes was constructed. The time step was 1 ms, and totally 2000 steps were done. To visualize results we save each 40 time steps: 3 slices: 2 main vertical slices and one horizontal equals to the day surface. The total computation time was 5.5 hours on 100 computation cores for both mathematical models. As expected, the amplitude of the seismic response for the block medium is significantly higher due to the existence of horizontal contact planes, that reflect the incident wave (see left Fig. 2, 3). Also, the 3D asymmetry of the response from the layered medium is clearly seen (see right Fig. 2, 3). It is explained by the presence of the preferred direction along the normal vector.  Figure 3. 3D -scattered seismic waves, t= 0.8 sec: in block (left), in layered cluster (right) .

Conclusion
Continual models of solid media with a discrete set of slip planes (layered, block media) and with nonlinear slip conditions at the contact boundaries of structural elements were constructed. For a stable numerical solution of a system of differential equations, an explicitimplicit method is proposed with an explicit approximation of the equations of motion and an implicit approximation of the constitutive relations containing a small parameter in the denominator. A set of effective formulas for the stress tensor correcting at the elastic step was obtained. Numerical simulations of the dynamic scattering process for subsurface layered and block objects in elastic 2D and 3D media were carried out using modern high-performance computing systems.