The requirement for hydrostatic initialisation in LS-DYNA / USA finite element models

The LS-DYNA/USA (Underwater Shock Analysis) coupled finite element codes are being investigated as a tool for predicting the local response of compliant plate structures subjected to far-field underwater explosion. It had previously been observed in LS-DYNA/USA models that extraneous pressure build-ups emanating from the DAA (doubly asymptotic approximation) boundaries may occur in the surrounding fluid region of the model, which inevitably lead to erroneous modelling of fluid-structure interaction and inaccurate structural responses. These instabilities typically result in divergence of the solution and eventually premature termination of the simulation. After a comprehensive investigation, it was found that the instabilities did not arise if the finite element model was hydrostatically initialised before conducting the LS-DYNA/USA simulation. The purpose of this study is to investigate the need for achieving hydrostatic equilibrium prior to the modelling of the shock wave propagation through the fluid-structure media. The method for achieving static equilibrium with the current version of the LS-DYNA/USA software is presented. The example simulations presented show that the hydrostatic initialisation procedure is effective in removing instabilities occurring at the DAA-fluid boundary, associated with the USA ambient hydrostatic pressure condition.


Introduction
The LS-DYNA explicit, non-linear, finite element analysis code [3], coupled with the Underwater Shock Analysis (USA) boundary element code [1,2,6] is being used by DSTO to predict local and global structural responses for floating and submerged naval vessels to underwater explosions.This capability is being developed with support and advice from Livermore Software Technology Corporation (LSTC) and Engineering Technology Corporation1 (ETC).This capability is highly desirable and of value to naval architects and to designers of onboard equipment and for the development of shock requirements for future RAN vessels.
The USA code uses the doubly asymptotic approximation (DAA) method, which is an approximation to the exact boundary-element representation of the surrounding fluid, which can be either an infinite or semiinfinite acoustic medium.It is called doubly asymptotic because it approaches the exact equations, and hence the exact answers, governing the fluid in steady state problems at the asymptotic limits of zero and infinite frequencies.Similarly, for transient response problems the DAA method approaches the exact equations at high frequency (early time, shock phase) and at low frequencies (late time, bubble phase).A smooth transition is assumed to take place between the asymptotic limits.
The requirement for hydrostatic initialisation arose following the appearance of extraneous pressures at DAA boundaries in the results from LS-DYNA/USA runs.In a series of simulations modelling the shock response of submerged air-backed plates to far-field underwater explosions [4], it was noticed that pressure instabilities occurred at the DAA boundaries surrounding the fluid region.These instabilities increased in magnitude with time, significantly affecting the pressure field in the fluid and consequently affecting the structural response of the target.After an extensive investigation of the LS-DYNA/USA modelling procedure, involving parametric studies using a range of input variables, it was discovered that the instabilities only occurred when the hydrostatic pressure parameter, HYDPRE, in the TIMINT input deck was set at a non-zero value.HYDPRE is a static pressure, which is added to the dynamic load pressure, which is subsequently passed back to LSDYNA from the USA code.
It was hypothesised that the instabilities at the DAA boundary (involving the HYDPRE parameter) could be avoided if the finite element model was in static equilibrium prior to performing the LS-DYNA/USA simulation.When the hydrostatic pressure is applied by the USA code, it is assumed that the system is already in the preloaded state when it is passed to the USA code.Consequently, pressure discontinuities at the DAA boundaries are avoided and the likelihood of instabilities in the DAA staggered solution procedure for solving the structural and fluid equations [5] are substantially reduced.
The development of the pressure field in the fluid adjacent to the responding target is contingent on the initial condition of the finite element model.For the fluid-structure interaction to be simulated accurately, the initial boundary conditions of the submerged structure and surrounding fluid, prior to the interaction of the shock wave, must be incorporated into the finite element model.This necessitates that the model be pre-loaded to achieve a state of static equilibrium, prior to the excitation by the shock wave.The procedure adopted by DSTO utilises LS-DYNA to initialise the structure and fluid to a static equilibrium state as a result of the ambient hydrostatic pressure, although alternative methods would probably be more CPU efficient (eg.by using an implicit finite element analysis code).
The hydrostatic initialisation procedure involves subjecting the structural and surrounding fluid meshes to ambient hydrostatic pressure conditions and allowing sufficient run time for the pressure field in the fluid to converge to a constant, equilibrium value across the fluid mesh.Similarly, the structural parameters, such as displacements and strains, must also converge to a static value.The hydrostatic initialisation procedure is achieved by using the dynamic relaxation feature of LS-DYNA.This feature has only recently become available for problems solved using LS-DYNA/USA.The nodal coordinates, shell element stresses and fluid pressures generated during the dynamic relaxation phase are then used by LS-DYNA as initial conditions for the finite element model, prior to the shock simulation.
The procedure for applying hydrostatic initialisation is presented in relation to an example problem involv-ing a fully submerged, air-backed square plate fully clamped to a rigid box-like structure.A series of finite element modelling results, highlighting the effect of hydrostatic initialisation, are presented below.The effect of hydrostatic initialisation on the structural response is investigated via a series of finite element simulations and provide an insight into it's effect on the overall model response.
Two LS-DYNA fluid materials are used in simulations in this investigation.These are: -Material-1 (MAT ELASTIC) -which is used for modelling inviscid, irrotational flow.It is a fluid option of the elastic material #1 and is the most suitable material for modelling fluid-structure interactions emanating from far-field underwater explosions.Material-1 has a numerical, rather than a physically meaningful, viscosity coefficient and consequently is only useful where viscosity is not of interest.As Material-1 does not use an equation-of-state, the acoustic speed is constant and therefore it is not suitable for modelling the development of a shock wave (ie.close to an explosion source).Material-1 has the flexibility of being able to be modelled with either Euler or Lagrangian meshes.Furthermore, the ALE capability can also be applied to this material.These capabilities would permit the modelling of relatively large structural motions (ie. at the fluid-structure interface) as compared to Material-90, which is limited to small structural displacements.-Material-90 (MAT ACOUSTIC) -is a linear, Eulerian, acoustic element that was formulated specifically for modelling low-pressure, acoustic shock wave propagation.Only pressure information is available to the user with respect to Material-90 meshes.However, Material-90 is a very cheap element in terms of CPU time and as such is useful for modelling large fluid volumes.Material-90 is limited to small structural displacements, particularly since the formulation is limited to a linear, Eulerian mesh.This is better understood when considering that the Eulerian fluid nodes at the fluid-structure boundary are coincident with Lagrangian structural nodes.Hence for any significant structural motions, the fluid elements immediately adjacent to the structural mesh will become stretched or distorted.Consequently, Material-90 is only suitable for modelling small amplitude, linear structural behaviour.

Finite element model
A quarter-symmetry finite element model of the shock rig with 99,825 degrees of freedom (dof) was produced using the MSC/Patran (v7.5) pre-processing software (refer Fig. 1).The finite element model consisted of 2500 Hughes-Lui formulation shell elements.Specific boundary conditions were required for nodes coincident with the xz-and yz-symmetry planes shown in Fig. 2. No particular boundary conditions were applied to the plate edges, since the plate and supporting shock rig was being modelled as a free body in a fluid medium.The 3 mm steel target plate and the mild steel supporting box structure were modelled using the LS-DYNA Material 3, which incorporates a bi-linear stress-strain model for modelling the elastic/plastic nature of a material.

Hydrostatic initialisation procedure
Hydrostatic initialisation was achieved for finite element models using Material-1 using the dynamic relaxation feature of LS-DYNA.The dynamic relaxation feature is not available for models using Material-90.The dynamic relaxation capability of LS-DYNA is a means of approximating a solution to a quasi-static process that would otherwise have been performed by an implicit code.Dynamic relaxation can be used to obtain the initial stress and nodal displacement field prior to the commencement of the transient analysis.This is in contrast to system damping, which is applied during the analysis phase to all or part of a finite element model.Dynamic relaxation is essentially a critically, or near-critically, damped dynamic system.Damping is applied in the form: where: µ is the input damping factor (0 < µ < 1; µ = 1 represents no damping); v is the nodal velocity at time 'n'; ∆t is the time step increment; a n is the nodal acceleration at time 'n'.
The dynamic relaxation phase continues until a convergence criterion is met.The convergence criterion is satisfied when the kinetic energy, E k , is reduced to some fraction of the maximum kinetic energy, E max k , referred to as the convergence tolerance, drtol, where: The dynamic relaxation parameters used were: Tolerance (drtol): 10 −6 relaxation factor (µ): 0.9999 max.termination time: 1000 ms actual termination time: 520 ms Although this procedure is relatively expensive with respect to CPU time, the dynamic relaxation procedure need only be run once for a particular fluid-structure model.A range of explosion scenarios can subsequently be modelled by running only the final time integration phase2 of the LS-DYNA/USA software.Alternatively, an implicit structural code could be used to generate the static equilibrium condition.For this investigation it was decided that it would be advantageous in terms of software compatibility to use the same software package for the hydrostatic initialisation as well as the structural deformation analyses.Hydrostatic initialisation was achieved by applying a static load of 0.022 MPa, equivalent to a hydrostatic pressure experienced at approximately 2.3 m depth, 3to the external free surfaces of the fluid material for the duration of the dynamic relaxation phase.During this period, the fluid pressures eventually converged to the hydrostatic pressure value 4 and the structure settled under the hydrostatic pressure load.
For this experiment, it was found that the dynamic relaxation phase of the LS-DYNA run terminated after 520 ms of real time.At this time, the required tolerance factor (ie. 10 −6 ) had been achieved.Figures 3, 4 and 5 illustrate the convergence of the displacement-time history of the plate centre, the pressure-time history in the fluid adjacent to the plate centre and the kinetic energytime history of the finite element model, respectively.Figure 3 indicates that the 3 mm 'test plate' centre had converged to a deflection of 4.2 mm from its unloaded position.Also, Fig. 4 shows that the fluid pressures had converged to the hydrostatic pressure value over the entire fluid region.
To determine the effect of hydrostatic initialisation on the finite element model, finite element simulations were produced with and without dynamic relaxation for comparison purposes.For each finite element run, the load described by equation ( 1) was applied.
where the peak pressure P 0 = 5.68 MPa and the exponential decay constant θ = 0.08 ms.

LS-DYNA/USA simulations
Four LS-DYNA/USA simulations were performed as listed in Table 1.The quarter symmetry finite element model shown in Fig. 1 has 4480 DAA 1 elements coincident with all exterior fluid segments (excluding symmetry boundaries).The acoustic pressure load, described by equation ( 1), was applied 25 cm along the z-axis from the centre of the plate at the DAA boundary for all Material-1 simulations.For the two Material-90 simulations, the load was initiated one fluid element back from the structure within the fluid.The pressure and decay constants are adjusted accordingly.
The spherical shock wave emanates from the nominal charge source, which is 3.8 m along the z-axis, perpendicular to the target plate centre5 as shown in Fig. 6.50 ms of structural response were generated for each of the five simulations, as this was sufficient time to capture the first peak positive and negative structural motions of the 3 mm plate.

LS-DYNA/USA Model #1
The first simulation listed in Table 1 was performed as a baseline measurement for comparison with subsequent runs.No hydrostatic pressure was applied at the DAA boundaries (i.e. the hydrostatic pressure parameter, HYDPRE, was set to zero in the USA TIMINT input deck) nor was the dynamic relaxation feature of LS-DYNA utilised in this run.The exponentially decaying pressure load described by equation 1 produced an inward (i.e. in the direction of propagation of the applied load) plate deflection to a displacement of 7.4 mm.The plate cycled back through its rest position and continued outwards into the fluid to a maximum displacement of 1.2 mm.

LS-DYNA/USA Model #2
The hydrostatic pressure parameter, HYDPRE, was set at 0.022 MPa in the USA TIMINT input deck (corresponding to a depth of 2.3 m).As in Model #1, dynamic relaxation was not utilised.It was observed from the results that extraneous pressures arising at the DAA boundaries gradually increase in magnitude with time, and even during the first several milliseconds of response, lead to the formation of high-pressure regions, as shown in Fig. 7.These spurious pressures increase in magnitude and volume with time and ultimately significantly affect the model's response, hence the gross deflection of the 3 mm plate.
It was hypothesised that the pressure instabilities arose because the finite element model had not been hydrostatically initialised prior to the application of the dynamic pressure load; see case #3 below.

LS-DYNA/USA Model #3
The hydrostatic pressure parameter, HYDPRE, was set at 0.022 MPa in the USA TIMINT input deck and the dynamic relaxation capability in LS-DYNA was utilised using the parameters listed in Section 3.This model was successfully terminated after 50 ms.The hydrostatic initialisation procedure was considered successful since no DAA instabilities had been observed in the fluid region.

LS-DYNA/USA Model #4
As a useful comparison of fluid materials, Material-1 was replaced by the acoustic element formulation Material-90.As mentioned in the introduction, Material-90 is inappropriate for use in models where large structural deformations are likely to occur.The hydrostatic pressure parameter, HYDPRE, was set to zero in the USA TIMINT input deck.The dynamic relaxation feature in LS-DYNA is not available with version 940.2 of LS-DYNA for Material-90.The pressure field in Material-90 was initialised according to gravitational and depth parameters provided in the LS-DYNA deck.However, although the fluid pressures were initialised, the model was not hydrostatically initialised, since the structure requires a finite period of time to achieve static equilibrium under ambient pressure conditions.6   * * Hydrostatic initialisation feature is not available for Material-90.However, the fluid pressure is initialised before the load application.
The load was applied one element layer back from the structure's wet surface and produced an inward plate deflection to a displacement of 7.4 mm.The plate cycled back through its rest position and continued outwards into the fluid to a maximum displacement of 1.2 mm.No pressure instabilities were observed in the fluid region.

LS-DYNA/USA Model #5
This run is the same as Model #4, except that the hydrostatic pressure parameter, HYDPRE, was set at 0.022 MPa.The load was applied in the same manner as Model #4 and produced a similar result.The plate deflected inwards to a displacement of 7.7 mm, then cycled back through its rest position and continued outwards into the fluid to a maximum displacement of 1.1 mm.No pressure instabilities were observed in the fluid region.

Discussion
The application of hydrostatic pressure (via the HY-DPRE parameter in the USA TIMINT deck) causes spurious pressures to arise at the DAA boundaries when using Material-1 as a fluid.With HYDPRE 0 applied at the DAA boundary elements and fluid pressures within these boundaries initially set to zero, it is reasonable to assume that this condition will initially lead to a pressure discontinuity at the DAA boundary.
The appearance of a local cavitation zone adjacent to the plate, caused by the motion of the responding target, effectively reduces the overall load applied to the plate.The existence and the extent of the local cavitation zone are dependent upon the plate's velocity, with respect to the adjacent water velocity, shortly after the interaction with the incident pressure wave.As a consequence, the accuracy of the finite element model is dependent upon the initial pressure condition of the fluid in the model, prior to application of the pressure load at the DAA boundary.Hence the authors propose that the accuracy of the finite element modelling results would be improved if the model were hydrostatically initialised prior to the application of the shock load.Also, it has been shown that build up of spurious pressures in the Material-1 fluid will not occur once the model has been hydrostatically initialised.

DAA boundaries gross plate deformation
The effect of hydrostatic initialisation on the fluid pressure field is illustrated in pressure-time profiles plotted in Fig. 8 for a fluid element halfway between the centre of the test plate and the outer surface of the fluid region.The two graphs plotted in Figure 8 represent models with and without hydrostatic initialisation, respectively.For the latter model, spurious pressures are initiated at the DAA boundaries and eventually lead to an erroneously high magnitude pressure field in the vicinity of the test plate as shown by the pressure fringe plot in Fig. 7.
The runs reported use the LS-DYNA elastic-fluid Material-1 and the acoustic fluid Material-90 in models #1 to #3 and models #4 and #5, respectively.The instability problems were not observed when Material-90 was used as the fluid material although, as explained in the introduction, material-90 is only suitable for modelling small amplitude, linear structural behaviour.The Lagrangian structural displacement amplitudes toward the fluid are also limited by the thickness of the Eulerian fluid elements adjacent to the structure.

Conclusions
-Spurious pressures may arise in fluid regions in the vicinity of DAA boundaries for LS-DYNA/USA models, particularly those using Material-1, which have not been hydrostatically initialised.The spurious pressure regions may inevitably lead to less accurate structural responses being predicted and furthermore, may cause premature termination of the LS-DYNA/USA simulation.-DAA boundary instabilities were observed for the LS-DYNA fluid Material-1, but not for Material-90 (although the latter is only suitable for modelling small amplitude, linear structural displacements).The pressures in the Material-90 elements are initialised according to physical parameters supplied in the LS-DYNA input deck.-The ability of LS-DYNA/USA finite element models to predict structural responses involving fluidstructure interactions may provide a more accurate representation if they are hydrostatically initialised prior to load application since this would better reflect the condition at the time of dynamic loading.

Fig. 1 .Fig. 2 .
Fig. 1.The quarter symmetry models of the (a) fluid volume surrounding the box structure and (b) the box structure and target plate.

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3.The displacement-time history for the centre of the 3 mm plate indicating displacement converges after 520 ms of real time.

Fig. 6 .
Fig. 6.The charge-structure displacement geometry used in the finite element model for this study.

Fig. 7 .
Fig. 7. Finite element pressure-fringe plot for Model #2 (not hydrostatically initialised, HYDPRE = 0) at 5.39 ms showing the erroneous pressure build-up in the fluid region (structural mesh removed for illustrative purposes) and the subsequent gross distortion of the fluid mesh at the fluid-structure interface.Light and dark regions indicate low and high pressures, respectively over an approximate range of 0 to 10 MPa.

Fig. 8 .
Fig. 8. Finite element pressure-time profiles for the fluid element adjacent to the centre of the 3 mm test plate for (a) Model #2 (i.e.not hydrostatically initialised) and (b) Model #3 (i.e.hydrostatically initialised) indicated by the dashed and solid lines, respectively.
The Material-3 parameters used in the finite element model were:

Table 1 LS
-DYNA/USA model parameters for dynamic relaxation and choice of fluid material.The resulting peak deflections of the 3 mm-plate centre are also listed * Hydrostatic deflections following dynamic relaxation process.