Numerical Analysis of Propagation of Nonlinear Waves in Prestressed Solids

The details of numerical algorithms implemented in CAE FIDESYS for the analysis of the propagation of nonlinear waves in elastic and viscoelastic bodies are discussed. It’s taken into account that waves propagation lead to new strains which superimpose on existing stresses (induced anisotropy) in the media. For the formulation of problem we used the theory of repeated superposition of large strains. The details of numerical algorithms for the analysis of the propagation of nonlinear waves in elastic and viscoelastic bodies are discussed. The implementation of the spectral element method for the nonlinear dynamic problems of elasticity under finite strains is considered. Some details of parallel computing on multicore and multiprocessor systems for the problems of nonlinear dynamic elasticity are presented. The results of numerical experiments obtained in CAE Fidesys are shown, in particular: the analysis of propagation of nonlinear shock wave; the analysis of propagation of surface waves; the dynamic processes related with the origination of a hole in a weakly compressible nonlinear-viscoelastic material.


Introduction
The dynamic effects are important in nonlinear solid mechanics (Blend, 1969).The subject of special interest is the propagation of waves in prestressed solids.The preliminary deformation leads to some important effects such as the dependence of wave velocity on the initial strains, the stress-induced anisotropy, and others.These effects are analysed in the present paper.For the formulation of problems we used the theory of repeated superposition of large strains (Levin, 1998).It's taken into account that waves propagation leads to new strains which superimpose on existing strains (induced anisotropy) in the media.The details of numerical algorithms implemented in CAE FIDESYS for the analysis of the propagation of nonlinear waves in elastic and viscoelastic bodies are discussed (Levin & Vershinin, 2008;Levin, Zingerman, Vershinin, Freiman, & Yangirova, 2013).The implementation of the spectral element method for the nonlinear dynamic problems of elasticity under finite strains is considered.Some details of parallel computing on multicore and multiprocessor systems for the problems of nonlinear dynamic elasticity are presented.The results of numerical experiments obtained in CAE Fidesys (Levin, Zingerman, Vershinin, & Nikiforov, 2012; the official website of Fidesys Limited, 2015) are shown, in particular: the analysis of propagation of nonlinear shock wave; the analysis of propagation of surface waves; the dynamic processes related with the origination of a hole in a weakly compressible nonlinear-viscoelastic material.

The Mathematical Statement of Problem
As elastic waves propagation occurs in the body, already having finite strains, it leads to the additional strains.The theory of repeated superposition of large strains is used for the problem solution (Levin, 1998;Levin & Vershinin, 2008;Levin et al., 2013).Let us consider the basic equations, initial conditions, and boundary conditions of this theory.Multiple states of a body are distinguished: the initial (undeformed) state; n-1 intermediate states, for which body goes step by step by successively applied external effects or, in the case of the viscoelastic material,processes taking place in it; final or current (n+1)-th state, to which body goes after the application of all external loads to it in the predetermined order.
The problem is formulated for the (n+1)-th state in coordinate basis of the n-th state.
The equation of motion Boundary conditions Constitutive equations for Murnaghan's material is the deformation gradient in transition from the initial state to the current (n+1)-th state; is the displacement vector defining the transition from the n-th state to the (n+1)-th state; n r is the radius-vector of a particle in the n-th state; is the generalized stress tensor in the base of the k-th state under transition from the initial state to the (n+1)-th state; at k=0, this tensor is the second Piola-Kirchhoff stress tensor; is the generalized strain tensor in the base of the initial state that describes strains accumulated in transition from the initial state to the (n+1)-th state; n ν is the normal vector to the boundary n Γ of a body in the n-th state.The superscript T denotes the transposition, and the superscript -T denotes the superposition of transposition and inversion of the second-rank tensors.

Numerical Solution
We will seek a solution of a presented system of partial differential equations using spectral element method (SEM) (Komatitsch & Violette, 1998;Komatitsch & Tromp, 2002;Komatitsch, Goddeke, Erlebacher & Michea, 2010), and Galerkin's weak formulation (Zienkiewicz & Taylor, 1971).The distinctive features of SEM, as applied to the numerical solution of discussed problems, are: 1) A fast computational process due to a diagonal mass matrix (Komatitsch & Violette, 1998;Komatitsch & Tromp, 2002), thanks to a special kind of cubature formula for the domain integration.In addition, an explicit time integration scheme is used.
2) A high accuracy of the solution's approximation for a small number of mesh elements.A numerical approximation error is estimated as [ ] , where for SEM, h is a typical size of the mesh elements, m is the order of approximation inside a spectral element, and h u is a numerical solution.
3) A convenient parallelization of computations for multithreading and multiprocessing systems, for example, using OpenMP and MPI technologies.High order of spectral element and data locality inside an element make this method well-parallelizable on hybrid massively parallel architectures with accelerators (Komatitsch et al., 2010;Charara, Vershinin, Sabitov, & Pekar, 2011).
Let us rewrite an original differential formulation in a weak form (Zienkiewicz & Taylor, 1971).Equation ( 1) is multiplied on a shape function i N and integrated over the region Ω occupied by a body.
Taking into account the boundary conditions (3) and using Green's formula we get: by the shape functions is used: where m K is the number of nodes in a spectral element mesh, and ) ( are unknown nodal displacements. Substituting equation ( 12) in ( 11) we get the following Galerkin system of ordinary differential equations: here U is a vector of unknown nodal displacements of Let us remark that nonlinearity in ( 13) is presented only for displacements and for their time derivatives.This observation allows us to apply well-known time integration schemes that were developed for solving dynamic problems of linear elasticity.We will use a time discretization scheme SS22 (Zienkiewicz & Taylor, 1971) for the solving the obtained system (13).This method is a second-order-accuracy method in time.Another similar method GN22 is considered also by Zienkiewicz & Taylor, 1971.
Let us suppose that we know a solution at the j-th time stepj U .Then at the (j+1)-th step we will seek the solution in the following way: where , and α is an unknown correction vector.
Substituting the expression ( 14) in equation ( 13) for some time moment between n t and 1 + n t , we get: Having solved the given system of nonlinear algebraic equations (15) using the Newton iterations (Zienkiewicz & Taylor, 1971) one can find α and obtain + by means of ( 14).The absolute stability of that method is guaranteed while (Zienkiewicz & Taylor, 1971).It is recommended to take as a time step a minimum characteristic size of the spectral element divided by the maximum velocity of elastic waves (Zienkiewicz & Taylor, 1971).At 0 2 = θ the time integration scheme is fully explicit.This algorithm is implemented in CAE Fidesys [5] based on the spectral element method.

Results
Let us present results of solving model problems.Let us consider square elastic plate with height-to-width aspect ratio 1:20, fixed by joints, as is shown in Figure 1.Plate consists of two parts.Left part is not deformed; while the right one is loaded by tensile load along the horizontal axis (as a result there are induced strains in that part).The load is removed instantaneously at the initial moment.The system comes out from the equilibrium state and two waves start to propagate in the plate from its center: to the left -a loading wave, to the right -an unloading wave.The dependency of loading and unloading waves' velocities on the value of the initial loading is shown in Figure 2. One can see that the loading and unloading waves' velocities differ essentially -and that difference increases with the initial loading.This result is due to the preliminary loading of the body's part, where the loading wave propagates.

Figure 2. The dependency of waves' velocities on the initial stress level
We present further a series of tests on numerical convergence in order to verify the results of numerical solution.Five calculations were performed with decreasing time-step and element size.The parameters of these calculations are shown in Table 1.The presented graphs show that black and red lines, which correspond to the maximum mesh refinement, almost coincide, which proves the mesh convergence in space.Let us consider further a problem of modelling of cylindrical surface waves propagation in the elastic body that is pre-loaded by the gravity force.For the convenience, all calculations will be performed in relative (dimensionless) quantities.A cylindrical sample (Figure 4) with the height 10, external radius 1 and internal radius 0.1 made of the Murnaghan material (Lurie, 1990) with , and the initial density 2200 0 = ρ , is loaded by the gravity force directed down along the cylindrical axis of symmetry.The force produces initial strains inside the solid.

A point force acting along the direction
. The force's amplitude varies in time in the following way: is the Cauchy-Green deformation tensor, ) ( The origination of the hole is modelled by the instantaneous removal of loadings from its boundary at time moment t=1 sec.We note that posing and solving problems of holes' origination under finite strains are only possible using the theory of repeated superposition of finite (large) strains (Levin, 1998;Levin et al., 2013).

Conclusion
An approach to the numerical solution of the dynamic problems of solid mechanics for the pre-loaded nonlinear elastic and viscoelastic bodies under the finite strains is considered in the article.The problem statement is done using the theory of repeated superposition of large strains.The numerical algorithms for problems solving are based on the spectral element method.This method is applied for solving the considered problems for the first time.Developed algorithms have a possibility for an efficient parallelization of the computational process.The presented solution of a 3D dynamic problem for the hollow cylinder which requires significant computing resources, demonstrates an efficiency of developed algorithms from the computational point of view.
An approach, described in the article, allows numerically solving the problems of propagation of 3D acoustic waves in a nonlinear elastic or viscoelastic loaded body with initial, in general terms, finite strains.This approach develops the mechanical models proposed by Levin (1998).
The results of numerical solutions of the problems presented in the article demonstrate the significant impact of pre-loading on the stress-strain states of solids.For example, in case of the problem of a wave propagation in the partially loaded plate a velocity of an unloading wave was 40% more for the initial loading p=2.5G than for the initial loading p=1G.For the problem of wave propagation caused by the action of a point force in the pre-loaded hollow cylinder, a wave front has the shape which differs significantly from the spherical one.This fact is explained by an induced anisotropy caused by an impact of pre-loading along the axis of the cylinder.
In the authors' opinion while using dynamic strength criterions under finite strains, the results of solving such problems and problems of a non-stationary origination of inclusions in loaded solids could be used for the dynamic problems of nonlinear strength theory.

Figure 1 ..
Figure 1.The geometry of the problem and boundary conditions Table 1.Mesh parameters for the test on numerical convergence The step of mesh refinement The number of nodes along the longest side of the plate The time step x1 the true stress tensor component 11 σ along the middle line of the plate at different time moments is shown in Figure 3 for different mesh sizes (colored).

Figure 3 .
Figure 3.The distribution of 11 σ along the plate at different time moments

Figure 4 .
Figure 4.The geometry of the problem of cylindrical waves' propagation

Figure 5 .
Figure 5.The displacement field at different time moments.The elliptic wave fronts clearly prove the presented pre-stressed solid state with an induced anisotropy

Figure 6
Figure 6 shows the distribution of true stress tensor component 22 σ in the time moments t =1.03 sec (left) and t