Convergence study of 1 – D CFD 1 cell size for shockwave parameters using Autodyn hydrocode

The use of CFD codes has become widespread to solve the shockwave propagation problem, as they are able to successfully model explosions. Using the structured grid spherically and cylindrical charges cannot be modelled in 3D, thus remapping techniques have been developed to map spherical shockwaves from 1D to 2D and then to 3D. This method highlights the importance of the accuracy of the 1D model, which is mainly dependent on the mesh size. This study aims at finding the required division ratio for the calculation domain.


Introduction
Apart from the empirical formulas for shockwave parameters mentioned in TM5-1300 (De-partment of the Army, Navy and the Air Force, 1990), there is a growing use of more ad-vanced techniques in blast effects engineering. Since CFD codes provide the possibility of calculation in 3D only on structured grids, one has to implement mapping techniques to get fine results. Remapping is a procedure where we calculate the explosion of a spherical charge only in 1D, and when it hits the ground, or any obstacle -and the problem turns into a 2D or 3D problemwe switch to the spatial solver, and copy the available 1D domain into the 3D model with all specific transport parameters. This way the spherical or cylindrical propa-gation of a shockwave can easily be handled in the 3D structured grid.
This procedure though requires caution, because care has to be taken about the accuracy of the 1D calculation. The accuracy of the 1D model will determine the accuracy of the whole time consumingcalculation. Accuracy mainly depends on the mesh size, so the aim of this study is to carry out a convergence study for this problem, in order to determine the necessary mesh size for accurate results. At first, the model setup is discussed, then the results are presented, and then, the convergence study is carried out based on the analysis results, and finally a conclusion is made and the required cell size is described.

Model setup and analysis settings
The most widespread, up-to-date method to investigate effects of shockwaves is using the Jones-Wilkins-Lee equation of state of explosives linked with a numerical Euler flow solv-er code. (Nagy, 2012) The equation of the Jones-Wilkins-Lee equation of state for TNT is based on Lee, Finger, Collins (1973). 3 Air is modelled as ideal gas. Initial condition of air is 101.33 kPa ambient pressure and 15 °C temperature. An Euler-Multimaterial solver is used to run the calculation. (Century Dynamics, 2006) The convergence study has been carried out on a domain, which is 6 m long with a spher-ical TNT charge, having a radius of 60 mm which corresponds to a charge weight of 1.47 kg. The one dimensional domain is visualized in Autodyn by a wedge, the angle of which is determined by the software. The spherical charge is located in the center with a given radius, the rest of the domain is air. Gauge points are set up every 100 mm measured from the center. Figure 1. shows the layout of the model, while Figure 2. shows a shaded pressure-distance diagram and a pressure-time function for one of the gauge points.  Covered Z values range from 0.439 m/kg 1/3 to 3.953 m/kg 1/3 , and it has to be pointed out that these values describe blast scenarios in a general way, so the convergence results to be determined and deductions to be made are independent of the actual charge size.   Figure  3. It can be seen that impulse values show a little scattering about the converged value, there is no significant difference between results obtained by different mesh sizes. These findings are also valid for the other gauge points. So we can state without further investigation, that impulse values can be trusted no matter which domain division was used to calculate it, while pressure values require further investigation.

Convergence study
The idea of a convergence study is to declare a specific value for the examined variable, which is considered accurate. This specific value should be an analytical solution, but in this problem, no analytical valueeither for overpressure or for impulseis possible to determine. Some authors carried out a convergence study for CFD mesh cell size using Kingery-Bulmash blast parameters as a comparison. (Yanchao et al., 2008) The validity of this method is questionable, since it is known, that Kingery-Bulmash curves contain in-accuracies as well, and field-test data is dependent on more variables showing scattering when measured. Ambient pressure, temperature, wind, point of detonation, shape of charge, homogeneity of charge are all variables with scattering, while all these are handled as fix parameters in CFD-calculation, and empirical constants only play role in the EOS of explosive material. Because of this difference in the number of variables which they are depending on, comparing results of the Kingery-Bulmash curves with CFD-calculation results is not a reliable method to measure the accuracy of mesh size. The value which we consider accurate is determined based on a relative error limit. We use the domain division ratio number "N" to describe the number of cells. As "N" increases, the number of cells increases, and results are getting more accurate, they converge to a specific yet unknownvalue. Introducing the relative error ε for the n-th division step, P ε n= P n n-1 we can declare a P pressure value to be accurate if its relative error compared to the pre-vious N division step is below a convergence criteria, which we take here as 2%. This theory requires • equal increments in N, • which is chosen in a way so that calculation time remains within reasonable limits.
As discussed previously, the first step was to determine the relative error of overpressure values for each n>1 division increment step. Calculated ε values can be found in Table 3. It can be stated that the N=5000 division meets the convergence criteria at all distances. Based on this statement, results obtained with N=5000 are declared accurate. Absolute error of a certain step is calculated using the formula: P ε n.abs = N=5000 n The absolute error, relative error and pressure values for gauge at 2 m can be seen on Figure 4. It can be clearly seen how the absolute error decreases as N increases, and how pressure values converge. However, relative error is nonmonotonic (shows more increasing and decreasing regions), which is because of the unequal increments in N, and it can be seen here, in addition to why the relative error with unequal N increments can not be used for declaring the accurate value.      Table 4. show that the farther the gauge point is from the center, the smaller the N number is required for an absolute error smaller than 5%. However if we only calculate the required N for the given R value, where the gauge is located, we get a non-monotonic function, which is plotted on Figure 5.   Figure 5. Required N function corresponding to 5% pressure error Using Figure 5. one can determine the required division of domain based on Z and R input variables. It is worth pointing out that function shown in Figure 5. has a peak at about Z=1, and it is a non-monotonic function. As a conservative approach it can be stated that dividing the scaled distance into 400 cells is enough to get accurate results for all scaled distances.

Conclusion
A convergence study has been carried out for blast parameters using 1-D Autodyn hydrocode. The investigation pointed out that impulse shows no dependency on the mesh size, while pressure values change considerably with mesh size. With the help of a parametric approach it has been shown that it's possible to create a "required N" function for domain division corresponding to a desired error limit. In the study, the 5% absolute error criteria was chosen and the required N function was presented in terms of scaled distance, so this function can be used for other blast scenarios to estimate the required cell size for domain division. Further investigation should include the effect and influence of charge division, and possibilities of the use of cell size zoning.