Numerical solution to the problem of a stress wave reflecting in a spherical thick-walled reservoir loaded with a time depending internal pressure

In the paper the problem of a stress wave reflecting in the wall of a spherical thick-walled reservoir loaded with a time depending internal pressure was investigated. The material of reservoir wall was assumed to be linearly elastic and isotropic. The stress wave generated on the internal surface of reservoir propagates towards outside and reflects on the external surface, generating the reflected component of wave which superimposes on the incident one due to the linearity of the problem. Consecutive reflections of the wave front generate the successive components of wave. The Laplace transforms determining solution independently for every reflected component of the stress wave were presented. The case of the step pressure, which describes the simplest model of immediate detonation, was considered. Due to the complex form of functions determining in the domain of time the analytical solutions for further reflections, the problem was integrated numerically. It was shown that for the selected finite number of reflections the describing equations can be solved together as one ordinary differential equations system. From the introductory analysis of results it follows that the summed up displacements of particular spherical layers of reservoir wall oscillate around their static values, despite the fact that the amplitude of particular components of wave approaches infinity.


Introduction
In the military technology there are present numerous constructions containing components in the shape of a spherical or a cylindrical thick-walled reservoir. They are often, from the nature of their operation, loaded internally by surge-pressure. This group contains not only the parts of armament but also the laboratory devices like manometric bombs and casings used in the ring method of dynamical testing of materials or explosion containment vessels.
Cylindrical and spherical symmetry of above mentioned objects facilitates the analysis of their dynamics, making the problem spatially one-dimensional [1]. Nowadays the basic tool used for modeling and analysis of the problems of dynamically loaded constructions are FEM/CFD based software packages. Nevertheless simple one-dimensional engineering theories are still not out-ofdate. They give the straight insight into interactions between the physical parameters of the problem, in contrast with the results of computer simulations, where some relations can be hard to catch.
In the case of cylindrical symmetry emerges the problem of necessary assumption concerning parameters behavior along axial direction. Moreover for a cylindrical stress wave in general there is no direct analytical solution. It is determined by Bessel functions usually represented with series [1]. Spherical symmetry lacks this problems and therefore analysis of this case was presented.

Formulation of the problem
Let us consider a thick-walled spherical reservoir of internal radius a and external b made of linearly-elastic a e-mail: ZielenkiewiczM@witu.mil.pl isotropic material. Its internal surface is loaded with an arbitrary pressure p(t). The problem is solved within the scope of linear elasticity theory. Due to its spherical symmetry it can be assumed as the spatially one dimensional boundary value problem depending in the spherical coordinate system only on the Lagrangian coordinate r and time t.
Therefore the states of stress and strain in the reservoir material can be represented by the following components: σ r -radial stress, σ ϕ = σ θ -circumferential stresses, ε r -radial strain, ε ϕ = ε θ -circumferential strains. The rest of the components of the stress and strain tensors equal to zero in the considered coordinate system.
According to the linear elasticity theory and generalized Hooke's law [2] we have where u is the radial displacement, E and ν denote Young's modulus and Poisson's ratio, respectively. For real materials the range of elasticity is limited. In ductile metals reaching by Huber-Mises-Hencky (HMH) reduced stress the yield point is assumed to be the limit. Therefore additionally in order to enable the assessment of limitation of practical application of solution HMH reduced stress was included in analysis. In case of a spherical symmetry it overlaps with Tresca hypothesis of maximum shear stress [3]  For an infinitesimal element of the linearly-elastic medium, the equation of motion can be written in the form [4] ∂σ r ∂r where ρ 0 is the initial density of medium. Eliminating the stresses σ r and σ ϕ from equation (3) by means of expressions (1) we obtain ∂ ∂r where , The quantity c e denotes the velocity of spherical stress wave propagation in linearly elastic medium and c 0 is the thin-rod wave velocity. The internal pressure p(t) acting on the internal surface generates a spherical stress wave propagating towards outside. As its front reaches the external surface it reflects from it, generating the reflected component of wave which superimposes on the incident one. When the reflected component reaches the internal surface the similar situation takes place, it reflects and generates the third wave component superimposing on the others and so on. The trajectory of the wave front together with the set of initialboundary conditions was presented in figure 1. As can be seen the wave trajectory divides the area between straight lines r = a and r = b into the series of triangles marked with Arabic numerals. Summed up parameters of solution in k. area were distinguished with subscript in parentheses (k) and the consecutive components of fields of parameters with subscript k. For example, for the radial stress we can write The general boundary conditions on the boundary surfaces for the equation (4) that have to be fulfilled for all areas are However none of the incident wave components generating in successive reservoir spherical sections non zero radial stress can alter the boundary conditions (7) and (8). Therefore the courses of radial stresses generated with the reflected and the incident wave have to sum up to zero on the corresponding boundary surfaces. It would be fulfilled only if the source of reflected wave were the radial stress of corresponding incident wave taken with minus sign. This reasoning leads to the following series of additional boundary conditions for the consecutive reflected wave components where i = 1, 2, 3, . . . Furthermore, it should be on the corresponding segment of wave trajectory.

General analytical solution
Because due to the spherical symmetry the medium movement is irrotational, the problem can be solved with the use of the method of scalar potential [4]. The field of displacement can be determined by the scalar potential ϕ, in case of the spherical symmetry being Substituting it into equation (4) gives as a result the equation of plain wave determining the quantity rϕ with general solution where twice differentiable function F represents the outgoing part of solution and G the incoming part, respectively. However, the only real source of wave lies on the internal surface, so in order to keep the uniform notation the wave functions determining the incoming wave components were also denoted with F. This way from equation (13) we have where τ k is the argument of wave function equal zero on the corresponding segment of wave trajectory, namely for the outgoing wave components k = 2i − 1 and for the incoming wave components k = 2i, respectively.

DYMAT 2012
Function F, together with its first and second derivatives with respect to its argument, through equations (1), (12) and (14) determines all mechanical parameters of reservoir Function F itself can be determined with the use of boundary conditions. For the first wave component the solution overlaps with the solution of the problem of a spherical cavity in an isotropic elastic medium [1]. Boundary condition (8) together with relation (21) yields For subsequent components taking into account the following facts and using boundary conditions (9) and (10) again together with equation (21) allow to write on the external surface and on the internal reservoir surface, respectively.
From the form of equations (24), (28) and (29) it results that they can be solved sequentially, obtaining the solution of every equation on the basis of the solution of preceding one. But as the domain of all arguments τ k is the same -0, + ∞) -they in fact constitute the ordinary differential equations (ODEs) system of infinite number of equations with unknown functions F k (τ) and F k (τ) with one independent variable τ. From initial-boundary condition (11) and relation (18) follow the homogeneous initial conditions for the system The shape of differential equations (24) and (29) determining the outgoing wave components indicates that they describe the harmonic oscillator of underdamped angular eigenfrequency [5] with excitation. On the other hand it is worth to notice that for the incoming components the equation (28) has similar form, but with the negative damping constant. The equations can be solved recursively using Laplace transform together with its properties [6]. After introducing the following quantities the transforms of the wave functions can be written in the following forms, respectively, for the outgoing wave component (including the first) and for the incoming component It will suffice to find the inverse transforms to determine the solution. The sought derivatives can be found either by common differentiation in the domain of time or through their transforms obtained by multiplying expressions (33) and (34) by s and s 2 . The most of the analytical functions suitable for determining pressure course p(t) have Laplace transform in the form of rational function. Therefore to find the inverse transforms of sought functions the method of partial fraction expansion can be applied.

Numerical solution to the problem
As can be easily shown, even for the simplest case of loading like the step pressure the analytical solution of the 04042-p.3

EPJ Web of Conferences
problem is rather useless for analysis due to its growing complexity with the consecutive wave reflections. On the other hand the considered equations are easy to solve numerically together, as one ODEs system, in the domain of time.
The differential equations are the second order. However, to integrate ODEs numerically they have to be converted to the system of ODEs of first order, which can be written in the compact vector form [7,8], After assigning (36) the system of equations (24), (28) and (29) will get the following form with homogeneous initial condition Again the structure of the system allows to use recursive formulae, this time for obtaining consecutive components of f vector. Obviously, the finite number of equations has to be chosen for computations, depending on the number of areas of interest. The problem formulated above was coded with the use of MathWorks MATLAB 6.1. The system was integrated by means of implemented solvers and then the sought parameters were calculated algebraically using equations (18-22) and through yet corresponding wave arguments (15) and (16).

Preliminary analysis of numerical solution for the step pressure
The step pressure p 0 was considered as the input function. In order to facilitate the analysis of calculations results the following dimensionless quantities were introduced  It was possible, because the problem is linear. However, the calculations itself were conducted on dimensional parameters for steel reservoir: E = 2.1 · 10 5 MPa, ν = 0.3, ρ 0 = 7800 kg/m 3 , a = 1 m loaded with pressure p 0 = 400 MPa. After making use of relations (39) the main influence on the shape of investigated parameters has Poisson's ratio ν [9,10]. But, as was shown in these papers, in wide range of ν parameter values the changes of solution are only quantitative and relatively small. Therefore 0.3 can be treated as some mean value. The spatial graph of U/P quantity for ν = 0.3, β = 10 and the step pressure was presented in figure 2. The considerable thickness of reservoir wall was selected intentionally in order to highlight some characteristic features of solution. The outline of the wave front reflecting from the boundary surfaces is well visible. Every pass of wave front through particular spherical section of reservoir wall generates damped oscillations, approximately around the static value. It can also be seen that the consecutive disturbances generated by passing wave front last longer with every pass.
It is worth to notice that although the incoming wave components are determined by the differential equations of oscillator with negative damping, the amplitudes of summed up courses of displacement in particular spherical sections are of the same order and do not increase with time, as should be expected.
The analogical planar plots for β = 1.5 were presented in figure 3. The courses were plotted for internal (ξ = 1), external (ξ = β) and mean (ξ = (β + 1)/2) surface. As can be seen, due to the relatively small wall thickness the interval between wave reflections, marked with visible 04042-p.4  discontinuities of first derivatives, is significantly shorter than the natural period of surface damped oscillations. The summed up courses of displacement in the studied spherical sections are close to harmonic oscillations around the static values. The oscillations have common phase and comparable amplitudes close to the static values. As should be expected, this way the system is closing to the case of single degree of freedom thin-walled shell.

DYMAT 2012
In figure 4 the analogical graphs for β = 2 were shown. In this case the oscillations are more irregular. Moreover, the amplitude on the internal surface is significantly higher than on the two others, which in turn are more close. On the other hand the internal surface amplitude became significantly lower than static value, opposite to the amplitude on the external and on the mean surface, which became higher in this comparison. The courses for β = 2.5 presented in figure 5 confirm all tendencies described above. Additionally it can be observed that the displacements on the external and on the mean surface attain considerable negative values. Nevertheless, the phase of main frequency oscillations in studied sections seem to remain common.
In an analogous way the plots of reduced stress S z were shown in figure 6. The main visible difference is the presence of their significant discontinuities on the mean surface accompanying every pass of the wave front. On the internal surface the only discontinuity occurs at the initial instant, while the external surface lacks any. Moreover the decrease of internal surface amplitude is visibly larger in comparison with the static value than for corresponding course of displacement.

Conclusions
From the analysis presented above the following conclusions can be drawn: 1. Although the nature of phenomena suggests the subsequent solving, the equations determining the problem of spherical stress wave reflections in the thick-walled reservoir can be integrated both sequentially or together as ODEs with one independent argument. The latter is the most convenient for numerical integration. However, in this case the former additionally allows the recursive computation of right-sides of the system for the given number of reflections of interest. 2. With the decrease of reservoir wall thickness the summed up courses of displacement are closing to regular harmonic oscillations characteristic for the case of single degree of freedom thin-walled shell. 3. With the increase of reservoir wall thickness the interval between wave reflections becomes larger than the natural period of internal surface damped oscillations. Therefore, the shape of disturbances generated by passing wave front in particular spherical sections of the reservoir wall can be easily distinguished. 4. Despite the fact that the incoming components of wave are determined with the differential equation of oscillator with negative damping and therefore their amplitudes approach infinity, the summed up variations of parameters in spherical sections remain of the same order. 5. The discontinuities in the course of reduced stress on the mean surface result from its dependence on the second derivatives of wave functions, which have non-zero values at the initial instants, i. e. on the corresponding segments of wave trajectory. On the boundary surfaces the discontinuities compensate due to the boundary conditions.