Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and Education" (GRID'2021), Dubna, Russia, July 5-9, 2021 THE GRID-CHARACTERISTIC METHOD FOR APPLIED DYNAMIC PROBLEMS OF FRACTURED AND ANISOTROPIC MEDIA Vasily Golubev Laboratory of Applied Numerical Geophysics, Moscow Institute of Physics and Technology, 9 Institutskiy per., Dolgoprudny, Moscow Region, 141701, Russian Federation E-mail: w.golubev@mail.ru The wave processes are occurred in different technological applications: the seismic survey process, the non-destructive material quality control, the ultrasound medical technique. In this research the dynamic loading problem of complex media is investigated. The grid-characteristic approach is extended for this case. The general approach on the curvilinear structured mesh is considered. The standard splitting technique is used for reducing the initial multidimensional system to a set of one- dimensional transport equations. This technique also allows to fulfill automatically physically correct linear contact conditions. The thin plate loading process is numerically simulated in the full wave three-dimensional case. The whole spectrum of elastic waves initiated on the fracture is observed. The two-dimensional mode conversion experiment is simulated for the anisotropic inclusion. Keywords: anisotropic medium; grid-characteristic method; mathematical modeling; waves Vasily Golubev Copyright Β© 2021 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). 32 Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and Education" (GRID'2021), Dubna, Russia, July 5-9, 2021 1. Introduction The wave processes are occurred in different technological applications: the seismic survey process, the non-destructive material quality control, the ultrasound medical technique. This phenomenon is described by the hyperbolic equation system. It should be noted that the analytical solution can’t be found for general heterogeneous materials or complex domain geometries. Recently, numerous numerical methods were developed for this mathematical problem: the discontinuous Galerkin method [1], the staggered-grid method [2], the boundary-element method [3] and the grid- characteristic method [4]. The latest one is being actively developed last years. In the paper [5] the numerical algorithm for the simulation of periodic non-linear deformable media was proposed. In [6] its application to the medicine ultrasound problem was demonstrated. In [7] the domain covered by materials with different rheological equations was considered. The proposed method was thoroughly tested in [8, 9]. Recently, an effective approach for the simulation of arbitrary oriented fractured inclusions was published [10]. This work is directed to the extension of the grid-characteristic approach to the general anisotropic media. To illustrate the general concept the problem of the isotropic fractured thin plate is considered initially. The final algorithm was successfully applied to the wave mode conversion problem. 2. Mathematical Model and Numerical Method The dynamic behavior of elastic media is described by the linear elastic system πœŒπ‘£βƒ—Μ‡ = βˆ‡ βˆ™ 𝜎, (1) where 𝑣⃗ is the velocity vector, 𝜌 is the medium density, 𝜎 is the stress tensor. For an isotropic case, the relationship between stresses and deformations are (πœ† and πœ‡ are Lame parameters) πœŽΜ‡ = πœ†(βˆ‡ βˆ™ 𝑣⃗)𝐼 + πœ‡(βˆ‡ βŠ— 𝑣⃗ + (βˆ‡ βŠ— 𝑣⃗)𝑇 ). (2) In the two-dimensional case, for anisotropic media there is more general connection between stress and strain tensors 𝜎π‘₯π‘₯ 𝐢11 𝐢12 𝐢16 πœ€π‘₯π‘₯ 𝜎 ( 𝑦𝑦 ) = (𝐢21 𝐢22 𝐢26 ) (πœ€π‘¦π‘¦ ), 𝜎π‘₯𝑦 (3) 𝐢61 𝐢62 𝐢66 πœ€π‘₯𝑦 where πœ€ is the strain tensor. 𝑇 Let’s take into consideration the vector of unknown functions 𝑒 βƒ—βƒ— = (𝑣π‘₯ , 𝑣𝑦 , 𝜎π‘₯π‘₯ , πœŽπ‘¦π‘¦ , 𝜎π‘₯𝑦 ) . In both cases, the govern system of equations can be represented in the canonical form βˆ‚π‘’ βƒ—βƒ— βˆ‚π‘’ βƒ—βƒ— βˆ‚π‘’ βƒ—βƒ— (4) = 𝐴π‘₯ + 𝐴𝑦 . βˆ‚π‘‘ βˆ‚x βˆ‚y For the anisotropic case, the matrixes are 33 Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and Education" (GRID'2021), Dubna, Russia, July 5-9, 2021 1 1 0 0 0 0 0 0 0 0 𝜌 𝜌 1 1 0 0 0 0 0 0 0 0 𝜌 𝜌 𝐴π‘₯ = 𝐢16 𝐢16 (5) 𝐢11 0 0 0 , 𝐴𝑦 = 𝐢12 0 0 0 . 2 2 𝐢26 𝐢26 𝐢12 0 0 0 𝐢22 0 0 0 2 2 𝐢66 𝐢66 (𝐢16 2 0 0 0) ( 2 𝐢26 0 0 0) The grid-characteristic method uses characteristic properties of the hyperbolic equation system. Here, its application is shown for the anisotropic medium model. To solve numerically the hyperbolic system of equations, we need to split it by two directions. For rectangular grids, it would be axes 𝑂π‘₯, 𝑂𝑦; for curvilinear grids, it must be an arbitrary direction π‘‚πœ‰. After rotating counterclockwise, the given coordinate system 𝑂π‘₯𝑦 by an angle 𝛽 for a point (π‘₯, 𝑦) and stretching new coordinates (πœ‰, πœ‚) by coefficients (not equal to 1, in the general case) π‘™πœ‰ , π‘™πœ‚ , respectively, the equations will change to βˆ‚π‘’ βƒ—βƒ— βˆ‚π‘’ βƒ—βƒ— βˆ‚π‘’ βƒ—βƒ— (6) = π΄πœ‰ + π΄πœ‚ , βˆ‚π‘‘ βˆ‚ΞΎ βˆ‚Ξ· π΄πœ‰ = π‘™πœ‰ (𝐴π‘₯ cos 𝛽 + 𝐴𝑦 sin 𝛽), (7) π΄πœ‚ = π‘™πœ‚ (βˆ’π΄π‘₯ sin 𝛽 + 𝐴𝑦 cos 𝛽). (8) Let us consider the homogenous equation along the direction π‘‚πœ‰: βˆ‚π‘’ βƒ—βƒ— βˆ‚π‘’ βƒ—βƒ— (9) βˆ’ π΄πœ‰ βƒ—βƒ—. =0 βˆ‚π‘‘ βˆ‚ΞΎ βƒ—βƒ—βƒ—βƒ—βƒ—0 = (cos 𝛽 , sin 𝛽)𝑇 . Orthogonal vector βƒ—βƒ—βƒ—βƒ—βƒ— π‘‚πœ‰ is given by a unit direction vector 𝑛 𝑛1 is given by 𝑇 βƒ—βƒ—βƒ—βƒ—βƒ—0 = (βˆ’sin 𝛽 , cos 𝛽) . The hyperbolicity of the equations means that matrix π΄πœ‰ has a full set of 𝑛 eigenvectors. This allows us to represent the matrix using its spectral decomposition: π΄πœ‰ = Ξ©βˆ’1 ΛΩ. (10) Here, the rows of Ξ© are left eigenvectors of π΄πœ‰ , the columns of Ξ©βˆ’1 are (right) eigenvectors of π΄πœ‰ ; Ξ› is the diagonal matrix consisting of eigenvalues, absolute values of which have the physical meaning of wave velocities. Substituting it into the initial system and multiplying both sides by the matrix Ξ© on the left, we can introduce a substitution πœ” βƒ—βƒ— = Ω𝑒 βƒ—βƒ—, (11) where πœ” βƒ—βƒ— consists of the Riemann invariants, and we can finally obtain βˆ‚πœ” βƒ—βƒ— βˆ‚πœ” βƒ—βƒ— (12) βˆ’Ξ› βƒ—βƒ—, =0 βˆ‚π‘‘ βˆ‚ΞΎ which is a system of independent (as Ξ› is a diagonal matrix) transport equations. For each equation, the value on the next time step is defined by the following expression: πœ”π‘– (𝑑 + Δ𝑑, πœ‰, πœ‚) = πœ”π‘– (𝑑, πœ‰ + πœ†π‘– Δ𝑑, πœ‚), (13) and is calculated using the approximation of the necessary order. On each time step, the increment to the original variables is computed and added. Since some Riemann invariant corresponding to zero eigenvalues results in zero additive part, it does not need to be calculated. Finally, the inverse transformation from the additive parts of the Riemann invariants to the given unknowns is done using the appropriate formulae. 34 Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and Education" (GRID'2021), Dubna, Russia, July 5-9, 2021 3. Simulation Results In this work, the problem of the elastic wave propagation in a thin plate 1 mm thick was considered. A through crack with a length of 10 cm was set on the symmetry axis. A computational grid with a step of 0.1 mm was constructed, covering a thin plate in all three dimensions. The horizontal dimensions of the computational area were equal to 30 Γ— 20 cm. A non-reflecting boundary condition was used on the lateral faces. The elastic medium properties were described through elastic wave velocities: 6153 m/s (longitudinal wave), 3099 m/s (shear wave). The material density was 2700 kg/m3. To specify the source, on the upper and lower plate surfaces, nodes were selected that lie on circles with a radius of 2 mm, shifted 10 cm to the right along the OX axis, and 5 cm up along the OY axis. A force directed along the radius to the center of the disturbance sources was set as the signal source. The time dependence was chosen to be limited in the time with the periodic amplitude modulation of the form 3-cycle Hann-window-modulated sinusoidal tone burst with a frequency of 600 kHz. To solve one-dimensional transport equations, a third-order accurate scheme was used on an extended template; the time step was chosen based on the Courant condition. Figure 1a shows the distribution of the horizontal 𝑣𝑦 component of the velocity arising on the upper surface of the plate at time 3 * 10βˆ’5 s. The figure clearly shows the formation of diffracted SH0 and S0 waves, reflected S0 and SH0 waves and Rayleigh waves propagating along the crack surface. The possibility of using the grid-characteristic method for calculating dynamic processes in thin plates with defects was confirmed. Π°) b) Figure 1. (a) The spatial distribution of the 𝑣𝑦 on the top surface at 𝑑 = 3 βˆ— 10βˆ’5 s. The fracture is represented by the white color. (b) The computational domain of the anisotropic problem. In the paper [11] the effect of the energy transfer from P-waves to S-waves was experimentally obtained. Our extension of the grid-characteristic approach to anisotropic media allowed us to carry out the same experiment numerically. The computational domain is represented at Figure 1b. The total size was 1.5 x 1 m. It consisted of four independent rectangular meshes with the spatial step equals to 0.5 mm. In the middle, the medium was described by the anisotropic model with following parameters: 𝜌 = 1920 kg/m3, 𝐢11 = 12.98 βˆ— 109 Pa, 𝐢12 = 2.77 βˆ— 109 Pa, 𝐢16 = 5 βˆ— 109 Pa, 𝐢22 = 79.2 βˆ— 109 Pa, 𝐢26 = 4.88 βˆ— 109 Pa, 𝐢66 = 13.05 βˆ— 109 Pa. The other parts were described by the isotropic model with parameters: 𝑉𝑃 = 6242 m/s, 𝑉𝑆 = 3144 m/s, 𝜌 = 2700 kg/m3. The time step 1.6 βˆ— 10βˆ’8 s was chosen based on the Courant condition. The explicit solution of the contact problem between isotopic and anisotropic media was used. It was constructed based on the glue condition involving four independent equalities on the contact boundary. The grid-characteristic method allowed to fulfill them automatically. In the numerical experiment the incidence P-wave with the main frequency of 90 kHz was used. The propagated and reflected waves are depicted at the Figure 2. The initiation of the shear wave right after the anisotropic inclusion is clearly seen. 35 Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and Education" (GRID'2021), Dubna, Russia, July 5-9, 2021 Figure 2. The spatial distributions of the 𝑣π‘₯ (left) and 𝑣𝑦 (right) are represented by colors 4. Conclusion In this research the dynamic loading problem of fractured and anisotropic media was investigated. The grid-characteristic approach was extended for this case. The general approach on the curvilinear structured mesh was considered. The standard splitting technique was used for reducing the initial multidimensional system to a set of one-dimensional transport equations. This technique also allowed us to fulfill automatically any physically correct linear contact conditions. The thin plate loading process was numerically simulated in the full wave three-dimensional statement. The whole spectrum of elastic waves initiated on the fracture was observed. The two- dimensional mode conversion experiment was successfully simulated based on described theoretical investigations for anisotropic media. The obtained results can be used for the computer simulation of complex applied dynamical problems. 5. Acknowledgement We thank our colleague Nikolay Khokhlov for fruitful discussions. This work has been carried out using computing resources of the federal collective usage center Complex for Simulation and Data Processing for Mega-science Facilities at NRC β€œKurchatov Institute”, http://ckp.nrcki.ru/. References [1] Zhang Z., Li D. M., Cheng Y. M. and Liew K. M. The improved element-free Galerkin method for three-dimensional wave equation // Acta Mechanica Sinica. 2012. V. 28 (3). P. 808-818 [2] Lisitsa V. V. and Vishnevsky D. M. On peculiarities of the Lebedev scheme for simulation of elastic wave propagation in anisotropic media // Sib. Zh. Vychisl. Mat. 2011. V. 14 (2). P. 155-167 [3] Zhaoand X. G., Rose J. L. Boundary element modeling for defect characterization potential in a wave guide // International Journal of Solids and Structures. 2003. V. 40. P. 2645-2658 [4] Kholodov A. S. The construction of difference schemes of increased order of accuracy for equations of hyperbolic type // USSR Computational Mathematics and Mathematical Physics. 1980. V. 20 (6). P. 234-253 [5] Nikitin I. S., Burago N. G., Golubev V. I., Nikitin A. D. Methods for Calculating the Dynamics of Layered and Block Media with Nonlinear Contact Conditions // Smart Innovation, Systems and Technologies. 2020. V. 173. P. 171-183 [6] Favorskaya A. V., Golubev V. I. Elastic and acoustic approximations for solving direct problems of human head ultrasonic study // Procedia Computer Science. 2020. V. 176. P. 2566-2575 36 Proceedings of the 9th International Conference "Distributed Computing and Grid Technologies in Science and Education" (GRID'2021), Dubna, Russia, July 5-9, 2021 [7] Golubev V., Shevchenko A., Petrov I. Simulation of Seismic Wave Propagation in a Multicomponent Oil Deposit Model // International Journal of Applied Mechanics. 2020. V. 12 (8). 2050084 [8] Favorskaya A., Golubev V. Study of Anisotropy of Seismic Response from Fractured Media // Smart Innovation, Systems and Technologies. 2021. V. 238. P. 231–240 [9] Beklemysheva K., Golubev V., Petrov I., Vasyukov A. Determining effects of impact loading on residual strength of fiber-metal laminates with grid-characteristic numerical method // Chinese Journal of Aeronautics. 2021. V. 34 (7). P. 1-12 [10] Khokhlov N., Favorskaya A., Stetsyuk V., Mitskovets I. Grid-characteristic method using Chimera meshes for simulation of elastic waves scattering on geological fractured zones // Journal of Computational Physics. 2021. V. 446. 110637 [11] Yang X., Kweun J. M. and Kim Y. Y. Theory for Perfect Transmodal Fabry-Perot Interferometer // Sci. Rep. 2018. V. 8, P. 69 37