- Email: [email protected]
An adaptive grid procedure for high intensity pulses propagating in air
Applied Acoustics, Vol. 48, No. 3, pp. 205-235, 1996 Copyright a 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0003-682X/96/%1 5.00 + 0.00 0003-682X(95)00060-7
ELSEVIER
An Adaptive Grid Procedure for High Intensity Pulses Propagating in Air P. Caine & M. West* The Telford Institute of Acoustics, University of Salford, Salford, UK, M5 4WT (Received 4 April 1995; revised version received 21 September 1995; accepted 29 September 1995)
ABSTRACT This paper describes a procedure for accurately modelling the evolution of a short duration high intensity acoustic impulse propagating over a large distance. An algorithm is described which achieves this objective by using a number of moving and non-overlapping grids. Copyright 0 1996 Elsevier Science Ltd Keywords: High intensity
acoustic
pulse, adaptive
grid algorithm.
1 INTRODUCTION Many existing numerical procedures for solving the problem of a large amplitude pulse propagating in one, two or three dimensions normally use a uniformly discretized spatial domain. ’ This demands substantial computing resources especially for two and three dimensions if the spatial domain is large compared with the smallest wavelength of interest. The size of the domain measured in the uniformly discretized mesh interval units depends on the starting scale of the phenomenon generating the pulse. For the practical application of interest considered here the phenomena are usually small explosive charges (for example, 0.5-8 kg of PE4 explosive) where the starting scale in the region where the pulse emerges into the air is typically centimetres or even millimetres to adequately resolve the waveform. The domain of interest extends typically to 1 km so that the number of points in *To whom correspondence
should be addressed 205
206
P. Cuinc, M. Wc’xt
the domain would be of the order 105- lo6 per dimension! A similar number of timesteps would be needed for the pulse to reach the end of the spatial range. The complete calculation therefore requires an order of 10’0-10’2 units of computation per dimension, each local unit of computation involves a solution of the full non-linear fluid dynamic equations. It is clear that this computation would be impractical on modest desktop computer hardware. We are mainly interested in problems which have one predominant direction. In this paper we investigate only the one-dimensional case. Recently adaptive grid techniques have been developed which can provide solutions in realistic times using modest computer hardware for problems in one and two dimensions.’ The organizational aspects of these techniques are complex and difficult to understand and implement because of their generality. In this paper by focusing on a very specific problem, namely a propagating single pulse. we have reduced the complexity suficiently to permit the development of an adaptive grid algorithm (AGA) which can be described more easily than the general algorithms. The adaptive grid procedure described here reduces the number of points to be computed from order ns2 to order n, log n,, where n,= rangej(finest mesh size). In this paper. although we concentrate on a homogeneous one-dimensional system of equations for the algorithm, we can extend this to the more practically useful case of one modelling the non-homogeneous system of equations governing a problem with spherical symmetry without serious difficulty. The model described here includes a novel method of maintaining conservation at the block (level) interfaces. This feature avoids interpolation in space and time at these level interfaces. Such interpolations are commonly found in other adaptive schemes. The method described here uses loss of one mesh point on the left-hand side of a block (level) and a corresponding gain of one point on the right-hand side of a mesh as the basis of the refinement procedure. This procedure for the treatment of fine-coarse interfaces, which requires shedding and gaining mesh points. moves the level along without the need for mesh refinement. The frame following capability of this system is obvious when low amplitude constant sound speed linear pulses are considered. The level configuration is designed so that the finest level mesh always adequately covers the parts of the pulse containing the most rapid variation. The other coarser mesh levels deal with less ‘severe’ gradients in the waveform and give us the capability of handling large spatial regions at economic computational speeds. Triggering of refinement is not necessary in the model in the linear case where sound speed does not change. However in the non-linear case we aim to configure the refinement so that:
High intensity pulses propagating
in air
201
(a) when a level is full the next coarser level is activated; (b) the pulse front is always close to the right-hand side of the finest mesh. If the pulse starts to drift too far to the left within the finest level then we interpolate from the level below to gather new mesh points. This interpolation is the ‘weakest’ point of the calculation and is most liable to numerical error but is carried out at the farthest possible point from the pulse front. As part of the interpolation procedure we used a first-order interpolation of the fluxes which maintains conservation. Our primary interest is in one-dimensional Cartesian cylindrically and spherically symmetric propagation problems. Generalization of the method to two space dimensions may be possible but would require a much more complex adaptive algorithm which would have to include a mesh control system triggered by wavefront slope. This will be addressed in a future publication. For many problems of interest the propagation in one direction is dominant and an approximate treatment of the components in the other directions is sufficient. In this case a one-dimensional treatment with the same methodology as that given in this paper will provide adequate threedimensional solutions.
2 OVERVIEW We wish to construct an adaptive multigrid scheme which will track an evolving pulse. In order to achieve this objective we employ a system of nonoverlapping grids or levels. Each level L has its own, unique, uniform mesh spacing. The levels are labelled 1, 2, .. .. LMAX. The convention adopted here is the higher the level number, the finer is the level’s mesh spacing. We impose the constraint for any given level L that its neighbouring (adjacent) levels can only be L + 1 and L-l (except for the first and the last levels of course!). There is a fixed mesh size ratio between levels and the total number of available mesh points in any level is a pre-assigned constant. When the algorithm is initiated only the finest level is switched on. This growing finest level mesh will cover the whole birth process of the emerging pulse. Subsequently, once all the points in any given level have been used, including the finest of course, then the coarser, adjacent, level is activated. When we advance At, at level L we would require m (an example with m = 3 is shown in Fig. la) time steps at level L + 1 to reach the same time. (Note that these individual time steps at level L + 1 are not generally equal except in the linear case.) The algorithm is arranged so that the finest mesh is used first and is the dominant level for the whole calculation.
P. Caine, M. West
208
(b) t
‘-1,
/ hAx(LMAX)
&MAX
x
kMAx(LMAX- 1) ~X,LMAX~I,
Fig. 1. (a) Schematic diagram of the time/space computational interface between two neighbouring levels: L and L + I, (b) Schematic diagram of the evolution of the leading edges of various levels for the linear problem. (c) Schematic diagram of the evolution of the leading edges for various levels for a non-linear problem.
The above requirement for maintaining time synchronicity is easily understood in the linear case where the sound speed, co, is fixed. On the x--t diagram the leading edge of the ‘frame’ for any level L is a characteristic line with slope l/co (Fig. 1b). The ‘frame’ is considered here to be all points in level L. When all the points in level L MAP have been used up we start a new characteristic line for level L MAX-l and so on. Because the behaviour is linear, from the above figure, each time we advance one spatial step Ax(L) we must advance one corresponding time step Al(L)= Ax(L)/co. For two levels, say L and L- 1, the coarser level, L- 1, has spatial steps which are by definition m times the size of those of the finer level, L. Similarly the
High intensity pulses propagating in air
(cl
LMAX-I
Fig. 1. -
209
LMAX
contd.
coarser level must have time steps which are m times those of the finer level. Therefore to maintain synchronicity, m steps at level L will be equivalent to one at level L-l. In the non-linear case we have curved characteristics and the time steps are no longer constant but increase with time until the linear conditions apply (Fig. lc). However, the number m of time steps at level L to be equivalent to one time step at level L- 1, is still equal to the ratio Ax(L-1)/0x(L). In this non-linear case, the frame speed is no longer co and varies as the calculation proceeds. In order to ensure that the frame contains the pulse we use a purpose built code to adjust the allocation of points. We want to take, as near as possible (allowing for stability constraintssee section 8.4) the maximum time step appropriate for each level. 2.1 Conservation law form We obtain the generic form of the solution procedure in Section 3 containing flux terms which are separately determined using the Roe algorithm.3 2.2 Dynamic pointer database
We consider our one-dimensional region to be made up from a set of grids (or levels) each of which contains a number of equidistant mesh points. Each
level, L, has a mesh spacing &C(L) which is an integer multiple of the finest mesh spacing, LLx(L MAX). Dyadic meshes are used typically. In Section 4 we develop a dynamic pointer database which contains a full specification of each level and its relation to the other levels. Included in the database is a structured set of real elements which monitor the absolute time and the time steps of each level. 2.3 Generation of level pattern sequences We require a processing order for each level in our system. This stems from the different time steps used at each level and the requirement for time synchronicity. If. for example, we have level LM*X initially and then level L MAX-l starts up then the number of time steps we can advance at level L MAX-l is governed by the progress at level LMAX. We cannot switch on another level unless (a) the finer level is full and (b) the time steps are synchronous This sets constraints on the level processing order. A special algorithm has been developed for generating a sequence of level indices in the order to be processed starting at the finest level, LM*x. This is described in section 5. 2.4 Treatment of level interfaces Each level must have an interface with the adjacent finer and/or coarser level. The accessibility of data across these interfaces is essential for solutions to be generated at the next time step. We describe a procedure for handling the solutions at these interfaces in section 6. 2.5 Adaptive grid algorithm The complete algorithm bringing together the imposed structures is presented in flow chart terms in section 7 with the master routine for advancing one time step at any level L treated initially as a ‘black box’. This master routine is then described in detail in relation to its components in section 8. 2.6 Prototype code-a Preliminary
first appraisal
tests with a multilevel
code are described
in section 9.
High
intensity
pulses propagating
in air
211
3 CONSERVATION LAW FORM OF THE FLOW EQUATIONS USING INTEGRAL FORMULATION The one-dimensional
form of the equations governing the flow are: (3.1)
where the flow vector u is given by
0 P pu
Q=
(3.2)
e
p being the density, u the particle velocity and e is given by e
=-----+ipu* P
(3.3)
Y-l
The pressure is p. The flux function is
( 1 PU
_F=
p+pu*
(3.4)
(e + P)U
All of the variables p, u, p are ‘total’* quantities. We integrate eqn (3-l) over a rectangular domain on a space (x)-time (t) diagram (see Fig. 2 below) .\-, XI 11 11 lJ(rl, x)dx I(to, x)dx + _F, xl)dt - _F(t,xo)dt = 9 (3.5) .I s s s ‘0
50
fo
Iu
We take small x and t intervals, Ax and At, respectively, as shown in Fig. 3. Applying eqn (3.5) to the shaded area in Fig. 3 with xo=xk--AX/~, x1 = xk + Ax/2 and to = t(“), t 1= tfn+ ‘) then
where we define \r+A.\/*
U” = -
-!T
1
A.Y
u( 6”’ , x)dx s Yl-A.Y/Z
*Total means background
(atmosphere)
value plus acoustic
or over-value.
(3.7)
P. Caine, M. West
212
Fig. 2. Schematic
n
diagram
k-
Fig. 3. Region of integration
I
x
XI
YI
of a grid cell in the x, r computational
k
I
k+
I
used in the local construction
domain
I
of the finite difference
scheme.
and (3.8)
We are using data at a time t w to compute a new solution with eqn (3.6) at a time 6”“). The first-order calculation assumes that Q is constant for x between xk-Ax/2 and xk + AX/~. In our implementation the flux functions _Fin eqn (3.6) are obtained using Roe’s method3 where
4 DYNAMIC
POINTER
DATABASE
4.1 Definition of one mesh at level L The level L is divided into Knox intervals and has Knox+ shown in Fig. 4. The coarsest mesh is at level 1 and the finest at level LM*~.
1 points as
213
High intensity pulses propagating in air
4.2 Integer control array
An integer array IC has been set up with a batch or grouped structure. Each batch of integers has a separate function. The array index is
kL=O
I
I
I
I kL =
&i’
kMA;lL’
L
Fig. 4. Mesh point labelling system for a level L. TABLE 1
Meaning of the batches Batch number (nb)
Meaning
_ 0 1 2 3 4 5 6 7 8
9 10 11 12 13 14 15 16 17 18
Starting k value in units of that level (usually zero for all levels) No. of spatial steps at level L equivalent to one step at level L-l
lZ
Total number of mesh points at level Final k value at each level. (Note this can be a small number initially to cut down computation for an evolutionary flow) Boundary type at locations kL=O at left-hand side. We have two boundary types: 1. source; 2. fine-coarse interface (coarse on left, fine on right). Boundary type at locations kl(max) at right-hand side. We have two boundary types: 3. outgoing wave (background values); 4. coarse-fine interface (coarse on left, fine on right). Starting address (i.e. index) for level L data in master array XA Switch to indicate processing status at level L: 1 for ON, -1 for OFF, 0 for WAITING STATE lll -1, shift backward one unit; 0, no shift; 1, shift forward one unit; m, shift by m units (positive or negative) l12 Local counter for level L It counts number of applications of the integrator state at L before level L- 1 can be updated ‘13 1, level L full; 0, level L not full. k index, at point, at level L where the maximum derivative (steepest waveform slope) is located Gives the amount of ‘free space’ in level L, i.e. the number of points which can be downgraded such that they can be allocated to level L- 1 l16 Record of starting values in batch 4 Colour index for waveform plot at each level L
214
P. Cuine, M. West
where nb is the batch number and L the level. We can summarize the ‘batches’ to be used by our code in Table 1. We note that IC values apply mainly to spatial features of each level L. From Table 1 note the following: *2: for level 1 the number of steps per lower level is set to 1. “11: level LMAX is always switched ON. When level L is full (all available points allocated to that level have been used up) we set the switch for level L- 1 to 0 which indicates that level L-l can now be CREATED from level L data. Note, we must be careful to delay the actual creation of data in level L-l until sufficient integrations have been done at level L to complete the first level L-l time step. “12: this value controls the frame position such that effectively it follows or tracks the propagating pulse in space. This is achieved by moving (typically) ALL points at level L one spatial step left reassigning the spatial index values accordingly (see Fig. 5). This means one point at the right of the frame is lost and a new point at the left is gained. This procedure is applicable strictly only when we can use a constant frame speed cO. In the non-linear case we require an adjustment to the number of points transferred to guarantee that the LMAx frame covers the pulse. “13: each application of the integrator will move us forward one time step at level L. We have to move a pre-set number of steps forward at level L (given in batch 2) before we can move forward one time step at level L- 1. This relates to time synchronicity to be discussed later. *16: in note 12 above the description refers to a single step shift. This current batch refers to the number of steps by which the frame can be shifted. 4.3 Temporal and spatial dynamic real data array This array, RDN, holds real data values for step sizes in space and time, and positions in space and time. The space and time steps are related through the CFL stability condition applicable to the integrator. As with the IC array, the RDN array is broken down into batches with batch index KQ,so the array has the form
Fig. 5. Resulting
frame position
after shift of data points from right to left.
215
High intensity pulses propagating in air
RDN(nb* LMAX + L) where J&AX > L 2 1. The array description is given in Table 2.
5 GENERATION
OF LEVEL PATTERN
SEQUENCE
Our objective here is to configure our adaptive grid algorithm using a pattern of level indices to specify their processing order. The following examples illustrate how the patterns are set up. 5.1 Example 1 Suppose we have only two levels (indices 1 and 2) and the finest level (index 2) has to be moved forward in time by five steps to say correspond to a single time step at level 1. The pattern or sequence of level indices to be processed after level 2 is full would be: 2, 2, etc.
2, 2,
2, 2,
2, 2,
2, 2,
1 1
TABLE 2
RDN array Batch number (nh)
Description
0
Current time at level L t= y
At,
where nL is number of time steps so far in level L. Note At,. varies (Fig. 5) with time Time increment for level L, At, (Fig. 5) Left-hand side value of x at level L, XL(O),in absolute units (m) (Fig. 4) Right-hand side value of x at level L, xl(max) (in m) Spatial mesh size for level L, AxL (Fig. 4). Note: (c + u),,~ AtL/AxL < 1 which is the CFL condition, where c is the local sound speed and u the local particle velocity. AxL is fixed, AtL is updated and computed Safety factor for level L replacing 1 in the CFL condition by a factor < 1 X= At/Ax values for each level Maximum time step for each level
P. Caine. M. West
216
The first line is repeated during the calculation process and is defined here as the ‘fundamental’ pattern which we try to find each time. 5.2 Example 2 This example considers three levels present with the following number of time steps in a given level being equivalent to one time step in the level below. Level 1
Level 2
Level 3
1
3
2
For example, we need two time steps at level 3 to be equivalent to one time step in level 2. We process in level 3 first until it is full. The fundamental pattern then is: 3, 3, 3.
3, 3, 3,
2 2 2,
1
We can generalize now for any number of levels with an arbitrary number of steps relating two adjacent levels. The procedure (implemented in our code) is as follows: (i)
Set PI(e) =
(LMAX)
(ii) Set fPI
= (PZ(e),
N(1) PI(e), .. ... PI(e))
Here N(1) is the number of time steps at level LMA~ to reach one time step at level LM~x-- I. (iii) Set P2(e) = (Pl,
LMAX -
1)
We regard this as an elemental pattern (e), just as Pi(e) elemental. P2(e) is repeated in pattern P2 below. (iv) Set a pattern P2 P2 = (P2(e), P2(e), ...... P2(e)) cN(2) -
was
Here N(2) is the number of time steps at level &*x-l to reach the same time step at the level below, that is level L~,4x--2, and so on.
High intensity pulses propagating in air
217
Let us examine example 2 using the above representation Pi(e) Pl
= =
P2(e)
=
P2
=
P3(e)
=
(3) (3, (3, (3, (3,
3) 3; 2) 3, 2; 3, 3, 2; 3, 3, 2) 3, 2, 3, 3, 2, 3, 3, 2; 1)
We can represent the actual processing order for each space-time cell in the table below.
T t
37
-_)
level 1
level 2
level 3 ---------)
x
In this example we have four spatial cells available at levels 2 and 3.
6 TREATMENT
OF LEVEL INTERFACES
6.1 Overview of requirements
Consider an interface, P, between two levels L and L-l (see Fig. 6). We take P as a point belonging to both levels. In order to advance our solution at any point by one time step we will need data at a point on either side of it (see eqn (6.1)). This is straightforward within a level but at the interface, say P, we do not have a level L point on its left-hand side. If P is considered to be in level L- 1, however, we do have a suitable point QL_, in level L, which provides the required data. In order to find QL we need to interpolate. One solution to the interfacing problem then is to obtain data at QL and QL-i (see Fig. 6). In the above procedure we have had to use an interpolation, which we would like to avoid if possible, and also we have not considered flux
P. Caine, M. West
218
L-I
Fig. 6. Schematic
diagram
Interface
L
of the spatial distribution of mesh between levels L- 1 and L.
points
near
the interface
conservation. In order to address these two issues we have devised an alternative procedure for interfacing, described below. 6.2 Non-interpolation
and conservation maintaining interface treatment
We can only advance a point at any level (say L) by one time step if we have data at points on the left and right of the point in question. It is clear, therefore, that we can only obtain solutions at points labelled as solid circles in Fig. 7. This means that at both levels L and L-l we can obtain all the required new solutions at all points except the point marked 0 (Q). We can obtain this unknown solution using the procedure described below. The solutions within any level L are obtained using a discrete flow integration (see eqn (3.6)) (6.1) where kL is a mesh point at level L; Ax= is mesh spacing (constant); AtL is a time step, XL= AtL/AxL; _Fis a flux vector (taken as an average value and assigned to a point in the centre of the t-x cell), Where _Uand _Fare given, respectively, in eqns (3-2) and (3-4). The three points used to obtain the solution at the next time step are indexed kL-& kL, k,+i. This flow integration can thus operate within the level L or level L-l. However, it cannot obtain the solution at point Q. 6.3 Treatment of an interface which guarantees conservation is maintained 6.3. I Solution procedure at the interface Consider the points at level L to the left of the interface at x0. All such points
can be advanced one time step using eqn (6-l). The point x0 could also be advanced provided we can extract the point x0+ AxL from level L+ 1. This point is identical to the point x0 + m AxL + 1 (m = 4 in Fig. 8).
High intensity pulses propagating
219
in air
.-.-. .-.-.-. F.-.-.-.-.
LQ
-.
.-.
r
I
I
I
I
I
I
I
I
I
P
Fig. 7. Scematic
diagram of spatial/temporal distribution of mesh points hood of the interface between two adjacent levels.
Fig. 8 Sequence of left-hand
in the neighbour-
mesh points at level L + 1 needed before level L can be updated.
Consider all points in level L+ I to the right of x0, and including x0. These points can all be advanced forward one time step At=+ I with the exception of the leftmost point at x o. Successive time solutions leads to the triangular ‘hole’ in the solutions possible at level L+ 1. Provided we can accept that the coarser level L can be extended AxL to the right we can manage without needing to fill this triangular hole. The above procedure for extracting the solution at x~“+’ from (x0-Ax=) and (xo+m Ax L+ i)n does not guarantee conservation nor would an interpolation done at t @J+‘1 between points at (x0-AxL)n+’ and (xo+m x,5+ i)n+‘. The procedure below does maintain conservation and first-order accuracy whilst avoiding interpolation. 6.3.2 A conservation maintaining procedure Let us examine the region around point x0 in Fig. 9. In what .follows, the integrations are of the form given in eqns (3.7) and (3.8), respectively. In order to maintain conservation we have: I!H~AW
- &,4B
+(&y
-&~f)=!l
(6.2)
P. Cuine. M. West
220
YII
h
h-l
Fig. 9. Schematic
diagram
used
ht
to describe the conservation.
treatment
I
of the
interface
employing
The capital letters which are subscripts of _Udelineate the integration path in x. Similarly the capital letters which are subscripts of _F delineate the integration path in t. We must split up the path HAB into HA at level L and AB at level L+ 1 (we regard all of H’A’B’ as lying in level L). Hence eqn (6.2) becomes l,&‘& - (r/H* + c;,,)
The flux along BB’ can be broken
+ (EBR’ - &d
= 0
(6.3)
down into two parts using eqn (3.8):
--FBB’
=
&j[, + &
(6.4)
The flux _FDrr is not known. _FD, can be found provided level L + 1 has been processed. Let us examine the two triangles ABD and DB’C’ in Fig. 9. Applying the conservation law to triangle ABD we have llAI> - YAB + ~BD - --FAD=
O
(6.5)
We note that path AD represents a time and space variation. This demands a change in our definition of &, and I&, making it more general but encompassing the definitions in eqns (3.7) and (3.8).
(6.6)
F
-AD
=
.I’
_F(t. s(t))dt
flZi
(6.7)
High intensity pulses propagating in air
221
Similarly for triangle DB’C’
u _B'C'
u + @DC’- EDB’) =0 - -DC
(6.8)
Substituting for Ena, from eqn (6.4) in eqn (6.3) and using eqn (6.5) for Z&a and eqn (6.8) for &,a,, then eqn (6.3) becomes U”‘tj’= &A + &c’ - _Urjc - @AC’- &iH’)
(6.9)
It can be shown using a Taylor series expansion that _ILAC,= $1;
+ O(Ax2, At’, 2AxAt)
(6.10)
I& = El;; f 0(Ax2, At’, 2AxAt)
(6.11)
and 2
We can write out each of the terms in eqn (6.9) using the n, k notation as follows
Here _UE’,lis the value of _Utaken to be that given by the normal solution procedure at level L + 1. The remaining terms, gAcj, and &, are given in eqns (6.10) and (6.11). We can then write eqn (6.9) as follows
(6.13)
Comparing eqn (6.13) with eqn (3.6) we can see that we have acquired a correction term on the last line of eqn (6.13). Using a Taylor series expansion of this term we can show that it is of second order. Thus, using eqn (6.13) would give first-order accuracy. The solution in eqn (6.13) can be shown to maintain conservation by checking that the difference between the integrals of _U with respect to x, taken at time indices n and n + 1, is equal to the flux difference taken between the left and right edges of the rectangular computational region in I-X space. The integrals are approximated by summing the _Uvalues at each point and multiplying by the appropriate Ax. We can incorporate the ‘awkward’ region around the interface in this integration process.
P. Caine, M. Wesl
222
7 DESCRIPTION
OF THE ALGORITHM
In what follows we use two large arrays XA and XB which both contain flow data for all existing levels L, where L belongs to the set { 1, 2,..., LMAX}. The data for the set of mesh points belonging to a level L is contained in the array XA and when a solution at a subsequent time is computed it is stored in a corresponding location in the array XB. After this procedure is completed, the appropriate part of XB is written into the XA array ready to repeat the process. 7.1 Flow chart of the scheme
The chart (see the Appendix) describes the section which we have called the ‘MASTER TIME LOOP’ in our code. We start by planting our initial data in the XA array for all levels. These initial data can either be the initial startup data or the last run’s data dumped to disk. 7.2 Lock into ‘pattern’ structure The pattern structure described in Section 5 is constructed here from an array NCD which holds the level processing pattern. The array NCD has two indices LC and J, where LC is the row index and J is the column index. The pattern sequence stored in NCD has a level index in each element. These level indices are obtained with the same algorithm as that described in Section 5. In the first instance we simply pick up the level index L and process at this level. However, we must also carry information about other levels which will be required to: (a) maintain time synchronicity between levels; (b) flag when we can switch processing to a lower* level; (c) ensure spatial matching of different level outputs. Each of the level indices stored in array NCD corresponds to a set of points with a mesh interval Ax= (fixed) and a time interval At, (variable). The number of points, in general, differs between the levels. 7.3 Calculation for a single time step forwards, AtL (routine MOVITIN
in
the code)
This routine contains the main fluid dynamical calculation for all points (at level L) which are advanced in the calculation by one time step from t to ‘Lower refers to a coarser
mesh
High intensity pulses propagating in air
223
t+ AtL. This procedure is based on the Roe algorithm which has been described in Ref. 3. The output solutions are held in array XB which has the same format as XA. Once all the solutions held in XB are known we write XB back into XA ready for the next time step calculation. The routine returns a time increment for level L MAX which maintains numerical stability and is stored in variable DTIME. 7.4 Updating the time step and the current time 7.4.1 Time step We pick up the time step AtLMAX(DTIME from MOVlTIN) which is the time step in thefinest mesh. ALL other mesh levels have their time steps AtL based on this value. We assume here that the finest level contains the most rapid variation of the solution. If L= LMAX the time step in LM*x is At LMAX = DTIME *
safety factor
(7.1)
In the code, this is written RDN(2 * LMAX) = DTIME * RDN(5 * LM,A,X + L)
(7.2)
(The safety factor is close to unity and is present to avoid code failure due to numerical instabilities.) The time step at the other levels is AtL = 2 AtLMAX(m) m=l
(7.3)
where m is the time step index from 1 to n, where n is given by n=
(7.4)
AxLIAxL~~~
Here n=n(L, L MAX) represents the number of intervals at the finest level (LM,& required to span or cover one level L interval. Note that the AxLs are constant throughout the calculation. We must use a spatial step ratio here rather than a time step ratio because the spatial steps are invariant. The above procedure ensures all levels are stable. Equation (7.3) is designed to maintain time synchronicity of the various levels. We let the finest level LM~X dictate the choice of temporal mesh spacing because that level is configured to contain the most rapidly varying part of the solution. Provided the finest level is stable we can show that all other levels must also be stable. The CFL criterion for stability of the level LMAX is ~(LMAx,
4
%.4AXW< ] Ax
hAX
1
(7.5)
224
P. Caine, M. West
for all m. Here X(L fix, m) represents from eqns (7.3) and (7.4)
the largest eigenvalue.
For level L
(7.6) Each term on the right-hand side is < 1 so that the stability criterion is met at level L. Note that this is applicable only if the pulse’s rapidly varying parts (or the parts with the largest X) are contained in level LMAX otherwise > 1 and instability ensues. X(L)IXLMAX) 7.4.2 Current time We need a master clock simply to inform us of the current time for labelling output. The current time is updated just after routine MOVlTIN by adding the current time step at level L to the last time value prior to the fluid flow calculation. If this level L is ON RDN(L)
We also update
= RDN(L) + RDN(LMAx + L)
the corresponding
time index
NT(L) = NT(L) + 1
again if the level is ON. In the case where L#L MAX we bypass eqn (7.1) onwards At, to zero
and set time step
Atl, = 0
In the code, this is written RDN(LMAx +
L) = 0
This statement is made just after we have performed a level L calculation. Since we use the level LMAX as the dominant level, each time we compute at level LMAX, we update the value of At, using eqn (7-2) in preparation for the next call to level L. 7.5 Restarting the calculation Having got to this point we have completed processing at one level L which has advanced one time step At,. Note, we may have processed many time steps at other levels greater than L. We hold the value of the last processed level L in the integer variable LASTL and then increment our level-pattern column counter, J, ready to go back to label 2001 (see flow chart) to process the next level obtained from the pattern, i.e. to advance all points in that level by one time increment. Depending on the pattern configuration this could be the same level we have just processed; then we would be advancing that level by an additional time step.
High intensity pulses propagating in air
225
This process continues until a negative NCD element is encountered which starts the procedure described below. 7.6 Restart pattern pointer to a column 1 If a negative L value is picked up from NCD we move our pointer back to column J = 1 on the same row we have already been working on. We then check to see if the level (LASTL- 1) is ON using our IC value IC( 11 * LM+,x+ (LASTL - 1)) (-1 OFF, 0 WAITING, 1 ON). If a value of - 1 is returned then although we must have processed level LASTL (to have reached this point) it cannot have been filled to capacity otherwise the above IC value would have returned a 0 or 1. The IC value is set in the routine MOVlTIN according to whether the level’s storage is full or not. In the case where the above IC element contains a 0 or 1 we know level LASTL is full and we must now move our pattern pointer to the next row of array NCD. 7.7 Restart pattern pointer at the next row We perform the trivial check to ensure that we have not passed the last row when we would jump back into the last row. This last row contains the full pattern sequence so that every level is then being processed. Provided our maximum allowed time (number of time steps) has not been exceeded we return to the beginning of the last row, column 1 of the NCD array at label 2003 (see flow chart).
8 MASTER
ROUTINE FOR ADVANCING ONE TIME STEP AT LEVEL L (ROUTINE MOVlTIN)
The flow chart for this routine is appended (see Figs A.1 and A.2 in the Appendix). The routine organizes the collection of data. It then computes solutions within one level which advance one time step at that level. Referring to the flow chart the routine first checks whether the level is to be created. In this case we create the level L by simply collecting data from level L+ 1 which has already been processed and plant it as a solution.
226
P. Caine, M. West
If the level L exists and is switched on then data are picked up at the current time according to whether it is: (a) a left-hand boundary point or (b) an interior point or (c) a right-hand boundary point Part (a) is shown on the flow chart between labels 9992 and 9993. The interior points are treated by the code from label 9993 to the end
Levels I IO
IO09070 flO-
40 -
20 -IO __--
0 I0 ii)
1 35
Distance.
__-.
_*--
..:
.,:
,..’
JO
m (x IO”)
x0 -L
70 -
5
ho-
”
51)~~
$ 40 2 2 30a zIO0-IO-
_--1 5.5
__--
___..‘.
..... .
I 60 Distance. m ix IO”)
Fig. 10. A series of ‘snap shots’ showing
I 65
the pulse position
in space.
221
High intensity pulses propagating in air
of the 3001 labelled DO loop. The last part of the code treats the case (c). Note that in section (a) we SAVE fluxes ready for use at level L-l and in section (c) we COLLECT and use pre-stored fluxes from level L+ 1. In advancing one time step we incorporate the moving frame concept by adding one new spatial point at the right of the spatial dataset and losing one point on its left at the new time in the case of a ‘full’ level. If we are working at level LM~X the ‘grabbed’ point is a background value, otherwise the value is taken from the next higher level, L + 1.
Levels
I IO 10090 x0 70 60 50 403020 IO ___----
___--’
. . . . ..“L
o-10
I 70
I 75
Distance,
I 80
tn (x 10”)
Levels
IlO-
IOO9080 -2
70-
x
ho-
z*
50-
2 z M 2 L
403020IO-
___---
___---.
,......
.;
o70
7s
Distance, Fig. 10. -
m (x 10”)
cod.
80
P. Caine. M. West
228
IO ~
__---
______.
,’
,:
0.r
Fig. IO.
9 A TEST
CASE
USING
cwntd.
THE
PROTOTYPE
CODE
For a first test of our prototype adaptive grid algorithm we have used a starting square pulse with a duration r of I ms and an amplitude of 1 atmosphere over-pressure. The Rankine-Hugoniot (RH) shock relations were used to compute the corresponding starting density and particle velocity values. At the end of the pulse duration, after T ms, we assume that all quantities drop back to background atmospheric values. (The RH relations are strictly only exact at the initial discontinuity.) Over the 1 ms time the pressure amplitude emerging from the source placed at the origin is 1 atmosphere and drops back down when the front edge of the pulse has travelled forward a distance cr, where c is averaged over the time r and is much greater than co, the low amplitude sound speed.
TABLE3 Level index
As(L)(m)
5
0.01
4 3 2 I
0.02 0.04 0.08 0.16
Number of’points
200 200 200 400 1600
Length
2 4 8 32 254
High intensity pulses propagating in air
229
We want to predict the field here for a maximum range of 300 m using say five mesh levels. We take a dyadic mesh ratio between adjacent meshes. For our finest mesh LMAXwe can set the mesh spacing such that over the distance CTwe have at least 50 points which will adequately resolve the pulse. Taking an average c of 500 m/s (say) would require a starting mesh spacing, Ax(LM,&, of 1 cm. At the outset we must prescribe the size of each of our levels in terms of number of points. In doing this we make allowance for expected spreading of the pulse. In Table 3 the values used in the test are listed.
10 CONCLUDING
REMARKS
A comparison of the test case predictions with those obtained using only the finest grid show that the use of the AGA with a five-level diadic grid system incurs negligible errors. This excellent result has been achieved because the sharpest part of the pulse is always well inside the finest mesh level LMAX, and because of the dyadic ratio. Use of larger mesh ratios which would give much increased speed causes deterioration in accuracy around the mesh interfaces. The number of mesh levels possible depends on available computer memory. We have managed to use a lo-level diadic mesh to provide a solution on a PC with 16 Mbytes of RAM. The model has been tailored for outgoing pulse propagation, and the test case presented conforms to this description. However the model has also been tested in cases where outgoing waves were present.
ACKNOWLEDGEMENTS The authors wish to thank the Ministry of Defence, Directorate of Health and Safety for their financial support for this work and Professor T. Toro of the Manchester Metropolitan University for his helpful suggestions and encouragement.
REFERENCES 1. Caine, P., West, M. & Walkden, F., A new non-linear solution procedure for one, two and three dimensional sound propagation. Applied Acoustics, 35 (1992) 2. 3.
297-3 10. Quirk, J. J., PhD thesis, Cranfield Institute of Technology, UK, 1990. Roe, P. L., Characteristic-based schemes for the Euler equations. Reviews in Fluid Mechanics,
18 (1986)
337-365.
Annual
P. Caine. M. Wesl
230
APPENDIX: DEMONSTRATION PERFORMANCE OF THE
OF THE IMPROVEMENT ADAPTIVE ALGORITHM
IN
If we make the simplifying assumptions below we can obtain an expression for the total number of computed solutions required to reach a given remote location distant X=mAx from the source at the origin, where Ax is the finest mesh interval, The assumptions are:
(1) the mesh system (2) each mesh level
is dyadic; has n mesh points; (3) the pulse is adequately covered by the n mesh steps at the finest level; and (4) processing in any given level ceases when all mesh points at the coarsest level have been used. of m spatial steps such
If we have r levels then we require a total number that 2’-‘) = m
n( 1 + 2 +
(Al)
and r=
ln(m/n + 1) In 2
642)
If we were working with the finest level only then the total number of mesh points processed (product of time steps and spatial steps) is given by P-t =im2
643)
When all levels are in use for each level LMA~-S have to process is 1
of points
(m - n[l + 2 + . ..2”])
~~~~~~~~ __s = - n2 +
2
Summing
the number
n
(A4)
2”
all the (PT)~~~~_’ for s from 0 to r-l,
we
r= LMAX
we
obtain r-2
P TOTAL
=
(PT)[_MAX_S
c
=
nm(r
-
1) + g(4
-
+)
+
i
n2
(A5)
0
s
From eqn (A2) for m>>n, r z ln(m/n)/ln
PTOTAL
=
nm
2 so that eqn (A5) becomes
ln(mln)
’
In 2
- 1 + 0(n2) I
(W
231
High intensity pulses propagating in air
Plant initial data in
LC is the row counter for the level control pattern matrix NCD. J is the column counter for matrix NCD.
Start at 1st column of NCD Pick up element (LC,J)
of NCD. This represents the level L being processed.
lYes MOVlTlN Calculate the flow soln for level L for one time step Pick up data from array XA and plant soln in array XB.
t If level = LMAX find DTIME t Update
approp
IC and RDN
DTIME is the step allowable which ensures (Bigger values cause failure.)
t COPY XA to XB 1 RDN(L) = RDN(L) + RDN(LMAX + L) Updates the time RDN(L) for level L to current value if it is switched on.
max time (computed stability. would
If level L is switched on IC( 11*LMAX+L).NE.I RDN(LMAX+L) contains the time step for level L
Update the number of time stepsbyl forlevelL
m
RDN(2*LMAX) = DTlME * RDN(S*LMAX+L) Update the time step using DTIME from
Fig. Al.
Flow chart
for adaptive
grid
program.
232
P. Caine, hf. West
RDN(LMAX
Set the time step for the level which has been processed to zero
+ Li = 0.0
t Plot results if destred
1 Store last level processed Increase column number by I, so can move horizontally Set column counter to I again
Is the level below LASTL switched on for processing? (Level
c
LASTL-
I IS“on”)
GRAPH plot if required
LC = LC + I
NB LC is the row counter Look at the next row
(While LC is less than or equal to the max level indicator)
233
High intensity pulses propagating in air
NO
Create level L
EXIT YYYI
Level L exists I X(L) = At(L) / Ax(L)
From array XA pick up source data for t = I,,(L). Plant in array UP
k,,(L)
1
I
Find source data at t = t,,(L) + At(L). Plant in array UPN.
Get the data for the 1st point and plant in UP
t Send UPN into XB array. Use KSTB to shunt data if necessary. (It Get the data from XA for point k,,(L) + I Plant in UC
t KK = K,(L) + 2
I KK=KK+
I
DO 3001
k Return from 300 t
I
Get the data from XA point KK. Plant in UN.
I Fig. A2. Flow chart for procedure MOVITIN.
P. Caine, M. West
234
1 Calculate flux FKMH from iJP and UC 1 If the level 13 “full” bake the flux FKMH along the diagonal of left hand \Ide of level FKMH -> FLUXT 4 Calculate the flux FKPH from UC and UN i FDIFF = FKPH
FKMH
# UCN = IJC (FKPH - FKMH) Send UCN to the XB array. U\e KSTB Set up DTIME.
This i\ the maxunum
copy UC -> UP t DO 300
1
I
Copy UN -> UC
(3001) NO We now treat the la\t two pomtr
We know that background value\ cx,\t ahead of the last point considered. Calculate new solution. Plant in UCN EXIT of the level
Fig. A2. -
if nece\wry
c~mtd
w
High intensity pulses propagating in air
235
The first term is dominant here. As an example take the following typical values: n = 500 m = 106. From eqn (A3) the non-adaptive calculation would require 5x 10” calculation steps. (From eqn (A2), r = 11.) From eqn (A6) the adaptive calculation would require approximately 5 x lo9 calculation steps-less than l/lOOth of the non-adaptive calculation.
ELSEVIER
An Adaptive Grid Procedure for High Intensity Pulses Propagating in Air P. Caine & M. West* The Telford Institute of Acoustics, University of Salford, Salford, UK, M5 4WT (Received 4 April 1995; revised version received 21 September 1995; accepted 29 September 1995)
ABSTRACT This paper describes a procedure for accurately modelling the evolution of a short duration high intensity acoustic impulse propagating over a large distance. An algorithm is described which achieves this objective by using a number of moving and non-overlapping grids. Copyright 0 1996 Elsevier Science Ltd Keywords: High intensity
acoustic
pulse, adaptive
grid algorithm.
1 INTRODUCTION Many existing numerical procedures for solving the problem of a large amplitude pulse propagating in one, two or three dimensions normally use a uniformly discretized spatial domain. ’ This demands substantial computing resources especially for two and three dimensions if the spatial domain is large compared with the smallest wavelength of interest. The size of the domain measured in the uniformly discretized mesh interval units depends on the starting scale of the phenomenon generating the pulse. For the practical application of interest considered here the phenomena are usually small explosive charges (for example, 0.5-8 kg of PE4 explosive) where the starting scale in the region where the pulse emerges into the air is typically centimetres or even millimetres to adequately resolve the waveform. The domain of interest extends typically to 1 km so that the number of points in *To whom correspondence
should be addressed 205
206
P. Cuinc, M. Wc’xt
the domain would be of the order 105- lo6 per dimension! A similar number of timesteps would be needed for the pulse to reach the end of the spatial range. The complete calculation therefore requires an order of 10’0-10’2 units of computation per dimension, each local unit of computation involves a solution of the full non-linear fluid dynamic equations. It is clear that this computation would be impractical on modest desktop computer hardware. We are mainly interested in problems which have one predominant direction. In this paper we investigate only the one-dimensional case. Recently adaptive grid techniques have been developed which can provide solutions in realistic times using modest computer hardware for problems in one and two dimensions.’ The organizational aspects of these techniques are complex and difficult to understand and implement because of their generality. In this paper by focusing on a very specific problem, namely a propagating single pulse. we have reduced the complexity suficiently to permit the development of an adaptive grid algorithm (AGA) which can be described more easily than the general algorithms. The adaptive grid procedure described here reduces the number of points to be computed from order ns2 to order n, log n,, where n,= rangej(finest mesh size). In this paper. although we concentrate on a homogeneous one-dimensional system of equations for the algorithm, we can extend this to the more practically useful case of one modelling the non-homogeneous system of equations governing a problem with spherical symmetry without serious difficulty. The model described here includes a novel method of maintaining conservation at the block (level) interfaces. This feature avoids interpolation in space and time at these level interfaces. Such interpolations are commonly found in other adaptive schemes. The method described here uses loss of one mesh point on the left-hand side of a block (level) and a corresponding gain of one point on the right-hand side of a mesh as the basis of the refinement procedure. This procedure for the treatment of fine-coarse interfaces, which requires shedding and gaining mesh points. moves the level along without the need for mesh refinement. The frame following capability of this system is obvious when low amplitude constant sound speed linear pulses are considered. The level configuration is designed so that the finest level mesh always adequately covers the parts of the pulse containing the most rapid variation. The other coarser mesh levels deal with less ‘severe’ gradients in the waveform and give us the capability of handling large spatial regions at economic computational speeds. Triggering of refinement is not necessary in the model in the linear case where sound speed does not change. However in the non-linear case we aim to configure the refinement so that:
High intensity pulses propagating
in air
201
(a) when a level is full the next coarser level is activated; (b) the pulse front is always close to the right-hand side of the finest mesh. If the pulse starts to drift too far to the left within the finest level then we interpolate from the level below to gather new mesh points. This interpolation is the ‘weakest’ point of the calculation and is most liable to numerical error but is carried out at the farthest possible point from the pulse front. As part of the interpolation procedure we used a first-order interpolation of the fluxes which maintains conservation. Our primary interest is in one-dimensional Cartesian cylindrically and spherically symmetric propagation problems. Generalization of the method to two space dimensions may be possible but would require a much more complex adaptive algorithm which would have to include a mesh control system triggered by wavefront slope. This will be addressed in a future publication. For many problems of interest the propagation in one direction is dominant and an approximate treatment of the components in the other directions is sufficient. In this case a one-dimensional treatment with the same methodology as that given in this paper will provide adequate threedimensional solutions.
2 OVERVIEW We wish to construct an adaptive multigrid scheme which will track an evolving pulse. In order to achieve this objective we employ a system of nonoverlapping grids or levels. Each level L has its own, unique, uniform mesh spacing. The levels are labelled 1, 2, .. .. LMAX. The convention adopted here is the higher the level number, the finer is the level’s mesh spacing. We impose the constraint for any given level L that its neighbouring (adjacent) levels can only be L + 1 and L-l (except for the first and the last levels of course!). There is a fixed mesh size ratio between levels and the total number of available mesh points in any level is a pre-assigned constant. When the algorithm is initiated only the finest level is switched on. This growing finest level mesh will cover the whole birth process of the emerging pulse. Subsequently, once all the points in any given level have been used, including the finest of course, then the coarser, adjacent, level is activated. When we advance At, at level L we would require m (an example with m = 3 is shown in Fig. la) time steps at level L + 1 to reach the same time. (Note that these individual time steps at level L + 1 are not generally equal except in the linear case.) The algorithm is arranged so that the finest mesh is used first and is the dominant level for the whole calculation.
P. Caine, M. West
208
(b) t
‘-1,
/ hAx(LMAX)
&MAX
x
kMAx(LMAX- 1) ~X,LMAX~I,
Fig. 1. (a) Schematic diagram of the time/space computational interface between two neighbouring levels: L and L + I, (b) Schematic diagram of the evolution of the leading edges of various levels for the linear problem. (c) Schematic diagram of the evolution of the leading edges for various levels for a non-linear problem.
The above requirement for maintaining time synchronicity is easily understood in the linear case where the sound speed, co, is fixed. On the x--t diagram the leading edge of the ‘frame’ for any level L is a characteristic line with slope l/co (Fig. 1b). The ‘frame’ is considered here to be all points in level L. When all the points in level L MAP have been used up we start a new characteristic line for level L MAX-l and so on. Because the behaviour is linear, from the above figure, each time we advance one spatial step Ax(L) we must advance one corresponding time step Al(L)= Ax(L)/co. For two levels, say L and L- 1, the coarser level, L- 1, has spatial steps which are by definition m times the size of those of the finer level, L. Similarly the
High intensity pulses propagating in air
(cl
LMAX-I
Fig. 1. -
209
LMAX
contd.
coarser level must have time steps which are m times those of the finer level. Therefore to maintain synchronicity, m steps at level L will be equivalent to one at level L-l. In the non-linear case we have curved characteristics and the time steps are no longer constant but increase with time until the linear conditions apply (Fig. lc). However, the number m of time steps at level L to be equivalent to one time step at level L- 1, is still equal to the ratio Ax(L-1)/0x(L). In this non-linear case, the frame speed is no longer co and varies as the calculation proceeds. In order to ensure that the frame contains the pulse we use a purpose built code to adjust the allocation of points. We want to take, as near as possible (allowing for stability constraintssee section 8.4) the maximum time step appropriate for each level. 2.1 Conservation law form We obtain the generic form of the solution procedure in Section 3 containing flux terms which are separately determined using the Roe algorithm.3 2.2 Dynamic pointer database
We consider our one-dimensional region to be made up from a set of grids (or levels) each of which contains a number of equidistant mesh points. Each
level, L, has a mesh spacing &C(L) which is an integer multiple of the finest mesh spacing, LLx(L MAX). Dyadic meshes are used typically. In Section 4 we develop a dynamic pointer database which contains a full specification of each level and its relation to the other levels. Included in the database is a structured set of real elements which monitor the absolute time and the time steps of each level. 2.3 Generation of level pattern sequences We require a processing order for each level in our system. This stems from the different time steps used at each level and the requirement for time synchronicity. If. for example, we have level LM*X initially and then level L MAX-l starts up then the number of time steps we can advance at level L MAX-l is governed by the progress at level LMAX. We cannot switch on another level unless (a) the finer level is full and (b) the time steps are synchronous This sets constraints on the level processing order. A special algorithm has been developed for generating a sequence of level indices in the order to be processed starting at the finest level, LM*x. This is described in section 5. 2.4 Treatment of level interfaces Each level must have an interface with the adjacent finer and/or coarser level. The accessibility of data across these interfaces is essential for solutions to be generated at the next time step. We describe a procedure for handling the solutions at these interfaces in section 6. 2.5 Adaptive grid algorithm The complete algorithm bringing together the imposed structures is presented in flow chart terms in section 7 with the master routine for advancing one time step at any level L treated initially as a ‘black box’. This master routine is then described in detail in relation to its components in section 8. 2.6 Prototype code-a Preliminary
first appraisal
tests with a multilevel
code are described
in section 9.
High
intensity
pulses propagating
in air
211
3 CONSERVATION LAW FORM OF THE FLOW EQUATIONS USING INTEGRAL FORMULATION The one-dimensional
form of the equations governing the flow are: (3.1)
where the flow vector u is given by
0 P pu
Q=
(3.2)
e
p being the density, u the particle velocity and e is given by e
=-----+ipu* P
(3.3)
Y-l
The pressure is p. The flux function is
( 1 PU
_F=
p+pu*
(3.4)
(e + P)U
All of the variables p, u, p are ‘total’* quantities. We integrate eqn (3-l) over a rectangular domain on a space (x)-time (t) diagram (see Fig. 2 below) .\-, XI 11 11 lJ(rl, x)dx I(to, x)dx + _F, xl)dt - _F(t,xo)dt = 9 (3.5) .I s s s ‘0
50
fo
Iu
We take small x and t intervals, Ax and At, respectively, as shown in Fig. 3. Applying eqn (3.5) to the shaded area in Fig. 3 with xo=xk--AX/~, x1 = xk + Ax/2 and to = t(“), t 1= tfn+ ‘) then
where we define \r+A.\/*
U” = -
-!T
1
A.Y
u( 6”’ , x)dx s Yl-A.Y/Z
*Total means background
(atmosphere)
value plus acoustic
or over-value.
(3.7)
P. Caine, M. West
212
Fig. 2. Schematic
n
diagram
k-
Fig. 3. Region of integration
I
x
XI
YI
of a grid cell in the x, r computational
k
I
k+
I
used in the local construction
domain
I
of the finite difference
scheme.
and (3.8)
We are using data at a time t w to compute a new solution with eqn (3.6) at a time 6”“). The first-order calculation assumes that Q is constant for x between xk-Ax/2 and xk + AX/~. In our implementation the flux functions _Fin eqn (3.6) are obtained using Roe’s method3 where
4 DYNAMIC
POINTER
DATABASE
4.1 Definition of one mesh at level L The level L is divided into Knox intervals and has Knox+ shown in Fig. 4. The coarsest mesh is at level 1 and the finest at level LM*~.
1 points as
213
High intensity pulses propagating in air
4.2 Integer control array
An integer array IC has been set up with a batch or grouped structure. Each batch of integers has a separate function. The array index is
kL=O
I
I
I
I kL =
&i’
kMA;lL’
L
Fig. 4. Mesh point labelling system for a level L. TABLE 1
Meaning of the batches Batch number (nb)
Meaning
_ 0 1 2 3 4 5 6 7 8
9 10 11 12 13 14 15 16 17 18
Starting k value in units of that level (usually zero for all levels) No. of spatial steps at level L equivalent to one step at level L-l
lZ
Total number of mesh points at level Final k value at each level. (Note this can be a small number initially to cut down computation for an evolutionary flow) Boundary type at locations kL=O at left-hand side. We have two boundary types: 1. source; 2. fine-coarse interface (coarse on left, fine on right). Boundary type at locations kl(max) at right-hand side. We have two boundary types: 3. outgoing wave (background values); 4. coarse-fine interface (coarse on left, fine on right). Starting address (i.e. index) for level L data in master array XA Switch to indicate processing status at level L: 1 for ON, -1 for OFF, 0 for WAITING STATE lll -1, shift backward one unit; 0, no shift; 1, shift forward one unit; m, shift by m units (positive or negative) l12 Local counter for level L It counts number of applications of the integrator state at L before level L- 1 can be updated ‘13 1, level L full; 0, level L not full. k index, at point, at level L where the maximum derivative (steepest waveform slope) is located Gives the amount of ‘free space’ in level L, i.e. the number of points which can be downgraded such that they can be allocated to level L- 1 l16 Record of starting values in batch 4 Colour index for waveform plot at each level L
214
P. Cuine, M. West
where nb is the batch number and L the level. We can summarize the ‘batches’ to be used by our code in Table 1. We note that IC values apply mainly to spatial features of each level L. From Table 1 note the following: *2: for level 1 the number of steps per lower level is set to 1. “11: level LMAX is always switched ON. When level L is full (all available points allocated to that level have been used up) we set the switch for level L- 1 to 0 which indicates that level L-l can now be CREATED from level L data. Note, we must be careful to delay the actual creation of data in level L-l until sufficient integrations have been done at level L to complete the first level L-l time step. “12: this value controls the frame position such that effectively it follows or tracks the propagating pulse in space. This is achieved by moving (typically) ALL points at level L one spatial step left reassigning the spatial index values accordingly (see Fig. 5). This means one point at the right of the frame is lost and a new point at the left is gained. This procedure is applicable strictly only when we can use a constant frame speed cO. In the non-linear case we require an adjustment to the number of points transferred to guarantee that the LMAx frame covers the pulse. “13: each application of the integrator will move us forward one time step at level L. We have to move a pre-set number of steps forward at level L (given in batch 2) before we can move forward one time step at level L- 1. This relates to time synchronicity to be discussed later. *16: in note 12 above the description refers to a single step shift. This current batch refers to the number of steps by which the frame can be shifted. 4.3 Temporal and spatial dynamic real data array This array, RDN, holds real data values for step sizes in space and time, and positions in space and time. The space and time steps are related through the CFL stability condition applicable to the integrator. As with the IC array, the RDN array is broken down into batches with batch index KQ,so the array has the form
Fig. 5. Resulting
frame position
after shift of data points from right to left.
215
High intensity pulses propagating in air
RDN(nb* LMAX + L) where J&AX > L 2 1. The array description is given in Table 2.
5 GENERATION
OF LEVEL PATTERN
SEQUENCE
Our objective here is to configure our adaptive grid algorithm using a pattern of level indices to specify their processing order. The following examples illustrate how the patterns are set up. 5.1 Example 1 Suppose we have only two levels (indices 1 and 2) and the finest level (index 2) has to be moved forward in time by five steps to say correspond to a single time step at level 1. The pattern or sequence of level indices to be processed after level 2 is full would be: 2, 2, etc.
2, 2,
2, 2,
2, 2,
2, 2,
1 1
TABLE 2
RDN array Batch number (nh)
Description
0
Current time at level L t= y
At,
where nL is number of time steps so far in level L. Note At,. varies (Fig. 5) with time Time increment for level L, At, (Fig. 5) Left-hand side value of x at level L, XL(O),in absolute units (m) (Fig. 4) Right-hand side value of x at level L, xl(max) (in m) Spatial mesh size for level L, AxL (Fig. 4). Note: (c + u),,~ AtL/AxL < 1 which is the CFL condition, where c is the local sound speed and u the local particle velocity. AxL is fixed, AtL is updated and computed Safety factor for level L replacing 1 in the CFL condition by a factor < 1 X= At/Ax values for each level Maximum time step for each level
P. Caine. M. West
216
The first line is repeated during the calculation process and is defined here as the ‘fundamental’ pattern which we try to find each time. 5.2 Example 2 This example considers three levels present with the following number of time steps in a given level being equivalent to one time step in the level below. Level 1
Level 2
Level 3
1
3
2
For example, we need two time steps at level 3 to be equivalent to one time step in level 2. We process in level 3 first until it is full. The fundamental pattern then is: 3, 3, 3.
3, 3, 3,
2 2 2,
1
We can generalize now for any number of levels with an arbitrary number of steps relating two adjacent levels. The procedure (implemented in our code) is as follows: (i)
Set PI(e) =
(LMAX)
(ii) Set fPI
= (PZ(e),
N(1) PI(e), .. ... PI(e))
Here N(1) is the number of time steps at level LMA~ to reach one time step at level LM~x-- I. (iii) Set P2(e) = (Pl,
LMAX -
1)
We regard this as an elemental pattern (e), just as Pi(e) elemental. P2(e) is repeated in pattern P2 below. (iv) Set a pattern P2 P2 = (P2(e), P2(e), ...... P2(e)) cN(2) -
was
Here N(2) is the number of time steps at level &*x-l to reach the same time step at the level below, that is level L~,4x--2, and so on.
High intensity pulses propagating in air
217
Let us examine example 2 using the above representation Pi(e) Pl
= =
P2(e)
=
P2
=
P3(e)
=
(3) (3, (3, (3, (3,
3) 3; 2) 3, 2; 3, 3, 2; 3, 3, 2) 3, 2, 3, 3, 2, 3, 3, 2; 1)
We can represent the actual processing order for each space-time cell in the table below.
T t
37
-_)
level 1
level 2
level 3 ---------)
x
In this example we have four spatial cells available at levels 2 and 3.
6 TREATMENT
OF LEVEL INTERFACES
6.1 Overview of requirements
Consider an interface, P, between two levels L and L-l (see Fig. 6). We take P as a point belonging to both levels. In order to advance our solution at any point by one time step we will need data at a point on either side of it (see eqn (6.1)). This is straightforward within a level but at the interface, say P, we do not have a level L point on its left-hand side. If P is considered to be in level L- 1, however, we do have a suitable point QL_, in level L, which provides the required data. In order to find QL we need to interpolate. One solution to the interfacing problem then is to obtain data at QL and QL-i (see Fig. 6). In the above procedure we have had to use an interpolation, which we would like to avoid if possible, and also we have not considered flux
P. Caine, M. West
218
L-I
Fig. 6. Schematic
diagram
Interface
L
of the spatial distribution of mesh between levels L- 1 and L.
points
near
the interface
conservation. In order to address these two issues we have devised an alternative procedure for interfacing, described below. 6.2 Non-interpolation
and conservation maintaining interface treatment
We can only advance a point at any level (say L) by one time step if we have data at points on the left and right of the point in question. It is clear, therefore, that we can only obtain solutions at points labelled as solid circles in Fig. 7. This means that at both levels L and L-l we can obtain all the required new solutions at all points except the point marked 0 (Q). We can obtain this unknown solution using the procedure described below. The solutions within any level L are obtained using a discrete flow integration (see eqn (3.6)) (6.1) where kL is a mesh point at level L; Ax= is mesh spacing (constant); AtL is a time step, XL= AtL/AxL; _Fis a flux vector (taken as an average value and assigned to a point in the centre of the t-x cell), Where _Uand _Fare given, respectively, in eqns (3-2) and (3-4). The three points used to obtain the solution at the next time step are indexed kL-& kL, k,+i. This flow integration can thus operate within the level L or level L-l. However, it cannot obtain the solution at point Q. 6.3 Treatment of an interface which guarantees conservation is maintained 6.3. I Solution procedure at the interface Consider the points at level L to the left of the interface at x0. All such points
can be advanced one time step using eqn (6-l). The point x0 could also be advanced provided we can extract the point x0+ AxL from level L+ 1. This point is identical to the point x0 + m AxL + 1 (m = 4 in Fig. 8).
High intensity pulses propagating
219
in air
.-.-. .-.-.-. F.-.-.-.-.
LQ
-.
.-.
r
I
I
I
I
I
I
I
I
I
P
Fig. 7. Scematic
diagram of spatial/temporal distribution of mesh points hood of the interface between two adjacent levels.
Fig. 8 Sequence of left-hand
in the neighbour-
mesh points at level L + 1 needed before level L can be updated.
Consider all points in level L+ I to the right of x0, and including x0. These points can all be advanced forward one time step At=+ I with the exception of the leftmost point at x o. Successive time solutions leads to the triangular ‘hole’ in the solutions possible at level L+ 1. Provided we can accept that the coarser level L can be extended AxL to the right we can manage without needing to fill this triangular hole. The above procedure for extracting the solution at x~“+’ from (x0-Ax=) and (xo+m Ax L+ i)n does not guarantee conservation nor would an interpolation done at t @J+‘1 between points at (x0-AxL)n+’ and (xo+m x,5+ i)n+‘. The procedure below does maintain conservation and first-order accuracy whilst avoiding interpolation. 6.3.2 A conservation maintaining procedure Let us examine the region around point x0 in Fig. 9. In what .follows, the integrations are of the form given in eqns (3.7) and (3.8), respectively. In order to maintain conservation we have: I!H~AW
- &,4B
+(&y
-&~f)=!l
(6.2)
P. Cuine. M. West
220
YII
h
h-l
Fig. 9. Schematic
diagram
used
ht
to describe the conservation.
treatment
I
of the
interface
employing
The capital letters which are subscripts of _Udelineate the integration path in x. Similarly the capital letters which are subscripts of _F delineate the integration path in t. We must split up the path HAB into HA at level L and AB at level L+ 1 (we regard all of H’A’B’ as lying in level L). Hence eqn (6.2) becomes l,&‘& - (r/H* + c;,,)
The flux along BB’ can be broken
+ (EBR’ - &d
= 0
(6.3)
down into two parts using eqn (3.8):
--FBB’
=
&j[, + &
(6.4)
The flux _FDrr is not known. _FD, can be found provided level L + 1 has been processed. Let us examine the two triangles ABD and DB’C’ in Fig. 9. Applying the conservation law to triangle ABD we have llAI> - YAB + ~BD - --FAD=
O
(6.5)
We note that path AD represents a time and space variation. This demands a change in our definition of &, and I&, making it more general but encompassing the definitions in eqns (3.7) and (3.8).
(6.6)
F
-AD
=
.I’
_F(t. s(t))dt
flZi
(6.7)
High intensity pulses propagating in air
221
Similarly for triangle DB’C’
u _B'C'
u + @DC’- EDB’) =0 - -DC
(6.8)
Substituting for Ena, from eqn (6.4) in eqn (6.3) and using eqn (6.5) for Z&a and eqn (6.8) for &,a,, then eqn (6.3) becomes U”‘tj’= &A + &c’ - _Urjc - @AC’- &iH’)
(6.9)
It can be shown using a Taylor series expansion that _ILAC,= $1;
+ O(Ax2, At’, 2AxAt)
(6.10)
I& = El;; f 0(Ax2, At’, 2AxAt)
(6.11)
and 2
We can write out each of the terms in eqn (6.9) using the n, k notation as follows
Here _UE’,lis the value of _Utaken to be that given by the normal solution procedure at level L + 1. The remaining terms, gAcj, and &, are given in eqns (6.10) and (6.11). We can then write eqn (6.9) as follows
(6.13)
Comparing eqn (6.13) with eqn (3.6) we can see that we have acquired a correction term on the last line of eqn (6.13). Using a Taylor series expansion of this term we can show that it is of second order. Thus, using eqn (6.13) would give first-order accuracy. The solution in eqn (6.13) can be shown to maintain conservation by checking that the difference between the integrals of _U with respect to x, taken at time indices n and n + 1, is equal to the flux difference taken between the left and right edges of the rectangular computational region in I-X space. The integrals are approximated by summing the _Uvalues at each point and multiplying by the appropriate Ax. We can incorporate the ‘awkward’ region around the interface in this integration process.
P. Caine, M. Wesl
222
7 DESCRIPTION
OF THE ALGORITHM
In what follows we use two large arrays XA and XB which both contain flow data for all existing levels L, where L belongs to the set { 1, 2,..., LMAX}. The data for the set of mesh points belonging to a level L is contained in the array XA and when a solution at a subsequent time is computed it is stored in a corresponding location in the array XB. After this procedure is completed, the appropriate part of XB is written into the XA array ready to repeat the process. 7.1 Flow chart of the scheme
The chart (see the Appendix) describes the section which we have called the ‘MASTER TIME LOOP’ in our code. We start by planting our initial data in the XA array for all levels. These initial data can either be the initial startup data or the last run’s data dumped to disk. 7.2 Lock into ‘pattern’ structure The pattern structure described in Section 5 is constructed here from an array NCD which holds the level processing pattern. The array NCD has two indices LC and J, where LC is the row index and J is the column index. The pattern sequence stored in NCD has a level index in each element. These level indices are obtained with the same algorithm as that described in Section 5. In the first instance we simply pick up the level index L and process at this level. However, we must also carry information about other levels which will be required to: (a) maintain time synchronicity between levels; (b) flag when we can switch processing to a lower* level; (c) ensure spatial matching of different level outputs. Each of the level indices stored in array NCD corresponds to a set of points with a mesh interval Ax= (fixed) and a time interval At, (variable). The number of points, in general, differs between the levels. 7.3 Calculation for a single time step forwards, AtL (routine MOVITIN
in
the code)
This routine contains the main fluid dynamical calculation for all points (at level L) which are advanced in the calculation by one time step from t to ‘Lower refers to a coarser
mesh
High intensity pulses propagating in air
223
t+ AtL. This procedure is based on the Roe algorithm which has been described in Ref. 3. The output solutions are held in array XB which has the same format as XA. Once all the solutions held in XB are known we write XB back into XA ready for the next time step calculation. The routine returns a time increment for level L MAX which maintains numerical stability and is stored in variable DTIME. 7.4 Updating the time step and the current time 7.4.1 Time step We pick up the time step AtLMAX(DTIME from MOVlTIN) which is the time step in thefinest mesh. ALL other mesh levels have their time steps AtL based on this value. We assume here that the finest level contains the most rapid variation of the solution. If L= LMAX the time step in LM*x is At LMAX = DTIME *
safety factor
(7.1)
In the code, this is written RDN(2 * LMAX) = DTIME * RDN(5 * LM,A,X + L)
(7.2)
(The safety factor is close to unity and is present to avoid code failure due to numerical instabilities.) The time step at the other levels is AtL = 2 AtLMAX(m) m=l
(7.3)
where m is the time step index from 1 to n, where n is given by n=
(7.4)
AxLIAxL~~~
Here n=n(L, L MAX) represents the number of intervals at the finest level (LM,& required to span or cover one level L interval. Note that the AxLs are constant throughout the calculation. We must use a spatial step ratio here rather than a time step ratio because the spatial steps are invariant. The above procedure ensures all levels are stable. Equation (7.3) is designed to maintain time synchronicity of the various levels. We let the finest level LM~X dictate the choice of temporal mesh spacing because that level is configured to contain the most rapidly varying part of the solution. Provided the finest level is stable we can show that all other levels must also be stable. The CFL criterion for stability of the level LMAX is ~(LMAx,
4
%.4AXW< ] Ax
hAX
1
(7.5)
224
P. Caine, M. West
for all m. Here X(L fix, m) represents from eqns (7.3) and (7.4)
the largest eigenvalue.
For level L
(7.6) Each term on the right-hand side is < 1 so that the stability criterion is met at level L. Note that this is applicable only if the pulse’s rapidly varying parts (or the parts with the largest X) are contained in level LMAX otherwise > 1 and instability ensues. X(L)IXLMAX) 7.4.2 Current time We need a master clock simply to inform us of the current time for labelling output. The current time is updated just after routine MOVlTIN by adding the current time step at level L to the last time value prior to the fluid flow calculation. If this level L is ON RDN(L)
We also update
= RDN(L) + RDN(LMAx + L)
the corresponding
time index
NT(L) = NT(L) + 1
again if the level is ON. In the case where L#L MAX we bypass eqn (7.1) onwards At, to zero
and set time step
Atl, = 0
In the code, this is written RDN(LMAx +
L) = 0
This statement is made just after we have performed a level L calculation. Since we use the level LMAX as the dominant level, each time we compute at level LMAX, we update the value of At, using eqn (7-2) in preparation for the next call to level L. 7.5 Restarting the calculation Having got to this point we have completed processing at one level L which has advanced one time step At,. Note, we may have processed many time steps at other levels greater than L. We hold the value of the last processed level L in the integer variable LASTL and then increment our level-pattern column counter, J, ready to go back to label 2001 (see flow chart) to process the next level obtained from the pattern, i.e. to advance all points in that level by one time increment. Depending on the pattern configuration this could be the same level we have just processed; then we would be advancing that level by an additional time step.
High intensity pulses propagating in air
225
This process continues until a negative NCD element is encountered which starts the procedure described below. 7.6 Restart pattern pointer to a column 1 If a negative L value is picked up from NCD we move our pointer back to column J = 1 on the same row we have already been working on. We then check to see if the level (LASTL- 1) is ON using our IC value IC( 11 * LM+,x+ (LASTL - 1)) (-1 OFF, 0 WAITING, 1 ON). If a value of - 1 is returned then although we must have processed level LASTL (to have reached this point) it cannot have been filled to capacity otherwise the above IC value would have returned a 0 or 1. The IC value is set in the routine MOVlTIN according to whether the level’s storage is full or not. In the case where the above IC element contains a 0 or 1 we know level LASTL is full and we must now move our pattern pointer to the next row of array NCD. 7.7 Restart pattern pointer at the next row We perform the trivial check to ensure that we have not passed the last row when we would jump back into the last row. This last row contains the full pattern sequence so that every level is then being processed. Provided our maximum allowed time (number of time steps) has not been exceeded we return to the beginning of the last row, column 1 of the NCD array at label 2003 (see flow chart).
8 MASTER
ROUTINE FOR ADVANCING ONE TIME STEP AT LEVEL L (ROUTINE MOVlTIN)
The flow chart for this routine is appended (see Figs A.1 and A.2 in the Appendix). The routine organizes the collection of data. It then computes solutions within one level which advance one time step at that level. Referring to the flow chart the routine first checks whether the level is to be created. In this case we create the level L by simply collecting data from level L+ 1 which has already been processed and plant it as a solution.
226
P. Caine, M. West
If the level L exists and is switched on then data are picked up at the current time according to whether it is: (a) a left-hand boundary point or (b) an interior point or (c) a right-hand boundary point Part (a) is shown on the flow chart between labels 9992 and 9993. The interior points are treated by the code from label 9993 to the end
Levels I IO
IO09070 flO-
40 -
20 -IO __--
0 I0 ii)
1 35
Distance.
__-.
_*--
..:
.,:
,..’
JO
m (x IO”)
x0 -L
70 -
5
ho-
”
51)~~
$ 40 2 2 30a zIO0-IO-
_--1 5.5
__--
___..‘.
..... .
I 60 Distance. m ix IO”)
Fig. 10. A series of ‘snap shots’ showing
I 65
the pulse position
in space.
221
High intensity pulses propagating in air
of the 3001 labelled DO loop. The last part of the code treats the case (c). Note that in section (a) we SAVE fluxes ready for use at level L-l and in section (c) we COLLECT and use pre-stored fluxes from level L+ 1. In advancing one time step we incorporate the moving frame concept by adding one new spatial point at the right of the spatial dataset and losing one point on its left at the new time in the case of a ‘full’ level. If we are working at level LM~X the ‘grabbed’ point is a background value, otherwise the value is taken from the next higher level, L + 1.
Levels
I IO 10090 x0 70 60 50 403020 IO ___----
___--’
. . . . ..“L
o-10
I 70
I 75
Distance,
I 80
tn (x 10”)
Levels
IlO-
IOO9080 -2
70-
x
ho-
z*
50-
2 z M 2 L
403020IO-
___---
___---.
,......
.;
o70
7s
Distance, Fig. 10. -
m (x 10”)
cod.
80
P. Caine. M. West
228
IO ~
__---
______.
,’
,:
0.r
Fig. IO.
9 A TEST
CASE
USING
cwntd.
THE
PROTOTYPE
CODE
For a first test of our prototype adaptive grid algorithm we have used a starting square pulse with a duration r of I ms and an amplitude of 1 atmosphere over-pressure. The Rankine-Hugoniot (RH) shock relations were used to compute the corresponding starting density and particle velocity values. At the end of the pulse duration, after T ms, we assume that all quantities drop back to background atmospheric values. (The RH relations are strictly only exact at the initial discontinuity.) Over the 1 ms time the pressure amplitude emerging from the source placed at the origin is 1 atmosphere and drops back down when the front edge of the pulse has travelled forward a distance cr, where c is averaged over the time r and is much greater than co, the low amplitude sound speed.
TABLE3 Level index
As(L)(m)
5
0.01
4 3 2 I
0.02 0.04 0.08 0.16
Number of’points
200 200 200 400 1600
Length
2 4 8 32 254
High intensity pulses propagating in air
229
We want to predict the field here for a maximum range of 300 m using say five mesh levels. We take a dyadic mesh ratio between adjacent meshes. For our finest mesh LMAXwe can set the mesh spacing such that over the distance CTwe have at least 50 points which will adequately resolve the pulse. Taking an average c of 500 m/s (say) would require a starting mesh spacing, Ax(LM,&, of 1 cm. At the outset we must prescribe the size of each of our levels in terms of number of points. In doing this we make allowance for expected spreading of the pulse. In Table 3 the values used in the test are listed.
10 CONCLUDING
REMARKS
A comparison of the test case predictions with those obtained using only the finest grid show that the use of the AGA with a five-level diadic grid system incurs negligible errors. This excellent result has been achieved because the sharpest part of the pulse is always well inside the finest mesh level LMAX, and because of the dyadic ratio. Use of larger mesh ratios which would give much increased speed causes deterioration in accuracy around the mesh interfaces. The number of mesh levels possible depends on available computer memory. We have managed to use a lo-level diadic mesh to provide a solution on a PC with 16 Mbytes of RAM. The model has been tailored for outgoing pulse propagation, and the test case presented conforms to this description. However the model has also been tested in cases where outgoing waves were present.
ACKNOWLEDGEMENTS The authors wish to thank the Ministry of Defence, Directorate of Health and Safety for their financial support for this work and Professor T. Toro of the Manchester Metropolitan University for his helpful suggestions and encouragement.
REFERENCES 1. Caine, P., West, M. & Walkden, F., A new non-linear solution procedure for one, two and three dimensional sound propagation. Applied Acoustics, 35 (1992) 2. 3.
297-3 10. Quirk, J. J., PhD thesis, Cranfield Institute of Technology, UK, 1990. Roe, P. L., Characteristic-based schemes for the Euler equations. Reviews in Fluid Mechanics,
18 (1986)
337-365.
Annual
P. Caine. M. Wesl
230
APPENDIX: DEMONSTRATION PERFORMANCE OF THE
OF THE IMPROVEMENT ADAPTIVE ALGORITHM
IN
If we make the simplifying assumptions below we can obtain an expression for the total number of computed solutions required to reach a given remote location distant X=mAx from the source at the origin, where Ax is the finest mesh interval, The assumptions are:
(1) the mesh system (2) each mesh level
is dyadic; has n mesh points; (3) the pulse is adequately covered by the n mesh steps at the finest level; and (4) processing in any given level ceases when all mesh points at the coarsest level have been used. of m spatial steps such
If we have r levels then we require a total number that 2’-‘) = m
n( 1 + 2 +
(Al)
and r=
ln(m/n + 1) In 2
642)
If we were working with the finest level only then the total number of mesh points processed (product of time steps and spatial steps) is given by P-t =im2
643)
When all levels are in use for each level LMA~-S have to process is 1
of points
(m - n[l + 2 + . ..2”])
~~~~~~~~ __s = - n2 +
2
Summing
the number
n
(A4)
2”
all the (PT)~~~~_’ for s from 0 to r-l,
we
r= LMAX
we
obtain r-2
P TOTAL
=
(PT)[_MAX_S
c
=
nm(r
-
1) + g(4
-
+)
+
i
n2
(A5)
0
s
From eqn (A2) for m>>n, r z ln(m/n)/ln
PTOTAL
=
nm
2 so that eqn (A5) becomes
ln(mln)
’
In 2
- 1 + 0(n2) I
(W
231
High intensity pulses propagating in air
Plant initial data in
LC is the row counter for the level control pattern matrix NCD. J is the column counter for matrix NCD.
Start at 1st column of NCD Pick up element (LC,J)
of NCD. This represents the level L being processed.
lYes MOVlTlN Calculate the flow soln for level L for one time step Pick up data from array XA and plant soln in array XB.
t If level = LMAX find DTIME t Update
approp
IC and RDN
DTIME is the step allowable which ensures (Bigger values cause failure.)
t COPY XA to XB 1 RDN(L) = RDN(L) + RDN(LMAX + L) Updates the time RDN(L) for level L to current value if it is switched on.
max time (computed stability. would
If level L is switched on IC( 11*LMAX+L).NE.I RDN(LMAX+L) contains the time step for level L
Update the number of time stepsbyl forlevelL
m
RDN(2*LMAX) = DTlME * RDN(S*LMAX+L) Update the time step using DTIME from
Fig. Al.
Flow chart
for adaptive
grid
program.
232
P. Caine, hf. West
RDN(LMAX
Set the time step for the level which has been processed to zero
+ Li = 0.0
t Plot results if destred
1 Store last level processed Increase column number by I, so can move horizontally Set column counter to I again
Is the level below LASTL switched on for processing? (Level
c
LASTL-
I IS“on”)
GRAPH plot if required
LC = LC + I
NB LC is the row counter Look at the next row
(While LC is less than or equal to the max level indicator)
233
High intensity pulses propagating in air
NO
Create level L
EXIT YYYI
Level L exists I X(L) = At(L) / Ax(L)
From array XA pick up source data for t = I,,(L). Plant in array UP
k,,(L)
1
I
Find source data at t = t,,(L) + At(L). Plant in array UPN.
Get the data for the 1st point and plant in UP
t Send UPN into XB array. Use KSTB to shunt data if necessary. (It Get the data from XA for point k,,(L) + I Plant in UC
t KK = K,(L) + 2
I KK=KK+
I
DO 3001
k Return from 300 t
I
Get the data from XA point KK. Plant in UN.
I Fig. A2. Flow chart for procedure MOVITIN.
P. Caine, M. West
234
1 Calculate flux FKMH from iJP and UC 1 If the level 13 “full” bake the flux FKMH along the diagonal of left hand \Ide of level FKMH -> FLUXT 4 Calculate the flux FKPH from UC and UN i FDIFF = FKPH
FKMH
# UCN = IJC (FKPH - FKMH) Send UCN to the XB array. U\e KSTB Set up DTIME.
This i\ the maxunum
copy UC -> UP t DO 300
1
I
Copy UN -> UC
(3001) NO We now treat the la\t two pomtr
We know that background value\ cx,\t ahead of the last point considered. Calculate new solution. Plant in UCN EXIT of the level
Fig. A2. -
if nece\wry
c~mtd
w
High intensity pulses propagating in air
235
The first term is dominant here. As an example take the following typical values: n = 500 m = 106. From eqn (A3) the non-adaptive calculation would require 5x 10” calculation steps. (From eqn (A2), r = 11.) From eqn (A6) the adaptive calculation would require approximately 5 x lo9 calculation steps-less than l/lOOth of the non-adaptive calculation.