C W Hirt's current address is Science Applications Inc, La Jolta, California
A numerical fluid-dynamics computing technique is presented that combines the Implicit Continuous-fluid Eulerian (ICE) and the Arbitrary Lagrangian-Eulerian. (ALE) methods. An implicit treatment of the pressure equation similar to that in ICE enables the calculation of flows at all speeds from supersonic to far subsonic. In addition, the vertices of the computing grid may be moved with the fluid in normal Lagrangian fashion or be held fixed in a Eulerian manner, or be moved in some arbitrary way to give a continuous rezoning capability, as in the ALE method. Greater distortions in the fluid motion can be handled than would be allowed by a purely Lagrangian method, and with more resolution than is afforded by a purely Eulerian method. The report describes the combined (ICED-ALE) technique in the framework of a computer program called YAQUI, for which the complete flow diagram and FORTRAN index listing are provided. Representative calculations illustrate some of the features of YAQUI, and include both computer-generated plots and numerical listings.
Over the past decade, there has been considerable progress in development of computer techniques for solution of multidimensional problems in fluid dynamics. A number of basic techniques have become well established, and useful and practical applications are being made to an ever-increasing range of problems in many fields. Because of computer storage and time limitations, numerical methods obviously cannot afford the luxury of following the dynamics of each and every molecule of the fluid at hand, but must, instead, depend upon following the dynamics of a finite, discrete set of fluid elements. Therefore, the region of interest is usually subdivided into a finite grid or mesh of computing zones, associating with each zones or vertex the local values of the quantities of interest, such as mass, energy, and velocity. The governing differential equations are approximated by finite-difference forms in relation to the grid, and this set of equations is then solved repeatedly over the domain to advance the solution through finite intervals of time, analogous to the frames of a motion picture.
Given this basic description, however, there are two fundamentally important considerations beyond which the various techniques differ. The first of these considerations is the flow-speed regime of interest, and the second is the interrelationship of the grid and the fluid, These two points will be discussed separately and then brought together.
The type of fluid flows that have been most amenable to calculation are generally those that can be characterized as either compressible or incompresible. Compressible, or high-speed, flows are those in which the fluid speed is comparable to or faster than the local material sound speed, and they are therefore governed only by local influences. In the incompressible or low-speed regime, however, fluid speeds are much less than material sound speeds, and disturbances at any point must, for all practical purposes, be felt instantaneously through-out the entire domain. As a result, the numerical stability restrictions for high-speed flows produce intolerably small time steps at low flow speeds, On the other hand, low-speed methods cannot sense compressibility effects produced by increased flow speeds, as no equation of state is used. Unfortuately, many fluid-dynamics problems of interest do not fall at either of these two extremes, and they are therefore not accurately calculated by either high- or low-speed methods. Examples are flows that are initially supersonic but rapidly become subsonic, or flows that are supersonic in one region or direction and far subsonic in another. Consequently, much effort is being placed on developing techniques to calculate in this intermediate regime.
The second point concerns the relationship between the fluid and the coordinate grid. Traditionally, there have been two basic viewpoints for both high- and low-speed flows. The first is Lagrangian, in which the mesh of grid points is embedded in the fluid and moves with it. Clear delineation of fluid interfaces and well-resolved details of the flow are afforded, but the approach is limited by its inability to cope easily with strong distortions, which so often characterize flows of interest. The second basic viewpoint, know as Eulerian, treats the mesh as a fixed reference frame through which the fluid moves. Strong distortions can be handled with relative ease, but generally at the expense of precise interface definition and resolution of detail.
Because of the obvious shortcomings of purely high-speed and purely low-speed methods, coupled with the shortcomings of purely Lagrangian vs purely Eulerian approaches, increasing enphasis is being placed on development of ever more sophisticated hybrid techniques.
Presently, the most successful method for calculating flows at all speeds is the Implicit Continuous-fluid Eulerian (ICE) technique [1], in which the flow may vary from supersonic to far subsonic. This is enabled by an implicit treatment of the pressure calculation. The method is extremely versatile and can be used for calculations in one, two, or three space dimensions, allowing for arbitrary equation of state.
Simultaneously, techniques have been developed that succeed to a great extent in combining the best features of both the Lagrangian and Eulerian approaches, In some methods. Lagranian particles are used to define fluid interfaces or free surfaces, or to define the fluid itself within a Eulerian mesh. There are other approaches, however, that have no basic dependence on particles. One such is the Arbitrary Lagrangian-Eulerian (ALE) [2] method, a low-speed technique that allows the vertices to move with the fluid in normal Lagrangian fashion or be held fixed in a Eulerian manner, or to move on some arbitrarily specified way to give a continuous rezoning capability. Greater distortions in the fluid motion can be handled than would be allowed by a purely Lagrangian method, with more resolution than is afforded by a purely Eulerian method.
This report describes a combination of the ICE and ALE schemes (ICED-ALE) in a computer program called YAQUI. It is based on the most recent improvements available in both the ICE and ALE methods, together with other improvements made possible by the marriage of the two schemes.
Although much more work remains to be done in the development of such hybrid techniques, YAQUI has established itself as a versatile tool for studying flows at all speeds and it has the capability of continuous rezoning.
The ICE technique was originally developed by F H Harlow, and the ALE technique by C W Hirt, who originated the ICED-ALE combination as it is represented in YAQUI, The code as it stands, however, represents the efforts of a number of people who have experimented with many alternatives along the way, and who have provided valuable contributions to its development. We are grateful for the help of our colleagues F H Harlow, T D Butler, K M Ruppel, J U Brackbill, R A Gentry. and W E Pracht of Group T-3, and for that of E M Jones and R C Anderson of Group J-10.
Inasmuch as the underlying technique has been discussed in detail elsewhere [2], this report will start from that point. First, the finite-difference equations will be presented as they appear in YAQUI, then the code itself will be discussed in detail.
The basic hydrodynamic part of each cycle of the ICED-ALE method is divided into three distinct subsections or phases. The first phase is a typical, explicit Lagrangian calculation. The second is an iteraticn that provides advanced pressures for the momentum equations and advanced compression for the mass equation. These ensure the stability of the method with respect to sound-signal propagation. Finally, the third phase, called the rezone section, performs all the convective-flux calculations, which must be included if the mesh is not purely Lagrangian.
The computing mesh consists of a two-dimenensional network of quadrilateral cells, and it will handle calculations in either cylindrical or plane (Cartesian) coordinates. Calculations in cylindrical coordinates are scaled to unit azimuthal angle, thus allowing the equations to be written without any π factors The radial coordinate is denoted by r or x, and the axial coordinate by z or y, with the origin located at the lower left corner of the mesh. The coordinate names in the equations in this report are x and y. The coordinate named r is used to detemine the geometry: r is always set equal to x for cylindrical coordinates, but the expressions automatically reduce to Cartesian expressions if all r's are set to unity. The vertices of the cells are labeled with the indices i and j, which increase in the radial and axial directions, respectively, Cell centers are denoted by half-integer indices i + l/2 and j + 1/2. The mesh of cells is I (IBAR) cells wide by J (JBAR) cells high.
I and J have a bar above in the original.
The mesh illustrated in Fig. 1 is in cylindrical coordinates, where the cells are sections of toroids of revolution about the cylinder.
The variables in an ICED-ALE grid are of two types: those defined at vertices, and those defined at cell centers. The principal variables are shown in Fig. 2, where coordinates (x and y), velocities (u and v), and masses (M) are defined at vertices, and the densities (ρ), pressures (p), volumes (V), and energies are defined at the cell centers. E is the specific total energy, and I is the specific internal energy.
In the equations that follow, the superscript n denotes the beginning-of-cycle values, The advancement of the solution through a time step, of duration δt, provides values at the beginning of the next (n+1) cycle. Intermediate values are typically labeled with a tilde for the results of Phase 1, or with a subscript L for the results of Phase 2.
Input Quantities: The input data supply the initial values of x, y, u, and v at the vertices and ρ and I for cells.
Preliminary Calculations (each cycle):
(1) The radius of each vertex is calculated as
r = x in cylindrical coordinates, or r = 1 in plane coordinates,
(2) Cell pressures are calculated at the beginning of each cycle using an equation of state
p = p (ρ,I)
although the equation of state may be bypassed for small Mach numbers. This is discussed further in Sec. III F, Incompressible Flow Calculations.
(3) Cell volumes are given by
The subscript notation for vertex quantities has been simplified to that shown in Fig. 2. It is used throughout this report and in the YAQUI code.
(4) With the cell volumes defined, the masses at cell centers can be computed from
but because most references are to the vertex masses, it is convenient to replace the cell masses immediately by vertex masses:
To maintain energy conservation throughout the entire calculational cycle, it is necessary to calculate and store E, the total specific energy per cell. However, the pressure iteration in Phase 2 requires a set of internal energies, I. One could get by with only a computer storage matrix of internal energies, by updating I during the iteration so that the total energy was conserved. The extra calculation required to do this, however, especially within an iteration, makes it seem reasonable to keep E and I separately. Therefore, we maintain a field of E's throughout the cycle, where initially
to which we will later add the various work and dissipation terms.
In this section we carry out a typical fully explicit Lagrangian calculation. with no grid motion, to obtain vertex values of the tilde velocities, ũ and ṽ and the change in total energy per unit mass, Q.
These three quantities are calculated in several steps. The following formulas show how these values accumulate from the contributions of each step. The appropriate initial values are, for each vertex,
Where Ax and Ay are body accelerations or accelerations from other forces applied at the vertices, and the superscript n denotes the beginning-of-cycle values. In most cases of interest, Ax and Ay are set equal to the gravity components
It is also helpful to insert a small, controlled, artificial diffusive acceleration into Ax and Ay at this point, To see the reason for this, consider the integration area that will be used for updating the velocity components at vertex (ji) in Fig. 3. The region is surrounded by dashed lines connecting the vertices (ji+1), (j+1i), (ji-1), and (j-1i) that will influence the accelerations at vertex (ji), but in the equations the acceleration computed from the surface stresses is independent of the vertex's location within the integration area, Although proper rezoning will tend to keep the vertex near the center of the region, and aid in obtaining the most accurate results, consider the integration area for the next cell (ji+1), indicated by the dotted lines in Fig. 3.
This indicates that values at four different vertices, (ji+2), (j+1i+1), (ji), and (j-1i+1), will enter. Although this definition of an integration area provides flexibility, there is a definite lack of communication between neighboring vertices, which can allow slight relative oscillations to arise in the velocity field. Introduction of a small restoring acceleration at each vertex, based upon the local velocity field, can prevent any vertex from deviating too strongly from its neighbors and couple the alternate nodes more strongly, This is done in YAQUI by introducing a weighted average of the neighboring vertex velocities, We can write
and anc is a coefficient that governs the amount of coupling, and upon which there is a stringent stability requirement, discussed in Sec. II F. It is possible to interpret anc as the number of cycles required for the vertex velocity to nearly equal the average of the neighboring velocities. The effect of this formulation becomes more apparent if one considers the initial tilde velocities that result from it:
which show the effective interpolation among neighbors that is added in. Note that for anc = 1.0, the technique becomes identical to a procedure that Lax introduced many years ago. To avoid the difficulty of that procedure as δt →0, it would be appropriate to take anc = a'nc/δt in which case a'nc is the actual relaxation time, rather than the number of cycles for relaxation.
Next, the appropriate initial vertex energy change is calculated as
but this is inappropriate because of the way we calculate Phase 1. The initial tilde velocities we have established contain only the body accelerations. and inserting this part into the initial Q's, before the pressure forces have been calculated, can cause the Q's, and hence the E's, to depart steadily from the correct value.
This effect would be manifested, for example, in a simple hydrostatic equilibrium, in which the velocities at time n are zero, To maintain equilibrium, we know that the final tilde velocities must also equal zero, but the gravitational accelerations in the initial tilde velocities would be repeatedly added into the Q's every cycle. This condition would not arise if we were to hold off the Q calculation until the final, complete tilde velocities were available, We choose, however, to form the Q's simultaneously in the same next step that will adjust for the various forces applied through pressure and viscous stresses over the control volumes surrounding each vertex. The changes in the initial ũ, ṽ, and Q values are computed by sweeping through the cells and suitably adjusting the four vertices of each cell, Thue the net result at each vertex is the cumulative contribution from each of the four surrounding cells. This technique of initializing vertices and then accumulating contributions from the cells is preferable to sweeping the vertices themselves as it is less dependent on boundary conditions and requires calculation of auxiliary cell quantities only once per cycle.
Thus, the set of energy changes that we obtain at corner vertices is subsequently assigned to the adjacent cells in a manner that preserves total energy while updating the vertex velocities.
This second step for each vertex proceeds as follows. First, for each cell, we calculate the divergence,
and the components of the viscous stress tensor: [3]
Here μ is the shear stress, and λ=ξ=⅔u, where ξ is the coefficient of dilatational viscosity. The corresponding finite-difference equations for these quantities are:
In the above equations, V is the cell volume and all the velocities on the right are the beginning-of-cycle values at time n, not tilde velocities, The coefficient CYL appearing in Πxx, Πxy and Πθ equals 1.0 when used in cylindrical coordinates, or 0.0 when used in plane coordinates. Note also the cyclic increase in index values in each term.
Next, with the D and Π terms calculated for a cell, the resulting changes in the ũ and ṽ velocities can be calculated at the four cell vertices as follows. Start by defining
where p is the cell pressure previoualy calculated from the equation of state. The energy changes are similarly calculated for each of the four vertices, but the p's are handled in a special mass-weighted fashion to improve accuracy when contact surfaces are present.. First, calculate the Q contributions without the work terms:
When all cells have been so treated, we can distribute the vertex energy changes, Q, into the stored cell-center energies E, to form an <E~>, which is denoted by brackets < > to identify that the pressures were omitted from the Q terms:
Next, convert <E~> values throughout the mesh to E's by sweeping all cells, Define the following mass-weighted ratios for each cell:
Although cell-center masses are no longer available, they can be approximated easily here by averaging four vertex masses. Finally, we can write
We have observed tbat this technique of calculating E~ in two steps is useful in enhancing the sharpness of shock fronts as well as contact surfaces.
The E formulation does, indeed, conserve total energy, and this can be shown as follows. If we sum over all cells
where the last sum has been changed from cells to vertex l, and the coefficient of Ql is precisely the mass of vertex l. Energy conservation is ensured, because the MlQl cancel in pairs when summed.
This completes the calculations associated with Phase 1 of the ICED-ALE cycle, If, at this point, one were to move coordinates with the ũ and ṽ velocities and calculate new densities, the result would be a typical, explicit. Lagrangian calculation.
We now need an implicit treatment to eliminate the Courant-like restriction on high sound speed that usually is required to ensure computational stability. This is accomplished by iterating the tilde- quantities from Phase 1 so to provide an advanced-time set of pressures for use in the momentum equations. These pressures, in turn, must reflect the new densities that will be calculated with the new velocities. In other words, the new densities are computed from coordinates obtained using accelerations that are functions of the new densities. Such an implicit treatment can, indeed, prevent instabilities at high sound speeds. For a completely incompressible flow, for example, the iteration tends to keep the ρ of each cell constant as the sound speed approaches infinity. The implicit coupling of p and ρ forces the cell to return to its initial ρ value, as ρ changes force corresponding pressure changes.
The implicit Phase-2 calculation proceeds as follows. First, we initialize velocities, densities, and pressures, where
uL = ũ vL = ṽ ρL = nρ pL = np
The subscript L identifies those quantities to be updated during the iteration. (In Phase 3, uL, vL and ρL will be further changed to their final values, n+1u, n+1v, and n+1ρ). As the tilde quantities ũ, ṽ, and np need not be saved for any other purpose, one can simply rename the tilde velocity arrays and the pressure array without any actual storage transfers. The quantity np, however, appears again in the Phase-3 convective flux equation, thus requiring separate storage for the ρL values.
In addition to the starting values of uL, vL, ρL, and pL, one must keep the nI values available for each cell in order to compute cell pressures. The Q values are no longer needed after Phase 1.
Second, we sweep the mesh systematically in i and j and make the following calculations for each cell.
(Note that this is the same divergence equation that appeared in Phase 1, except that the velocities are at step L instead of time n.) From the mass equation, we define
This prescription for A is exact only when the cells are rectangular, but it is much simpler and quicker to calculate than the full general value, which may be preferable when the zoning deviates strongly from the rectangular, as errors in A alter the rate of convergence.
In the above, C2 = (∂p/∂ρ)s is the square of the adiabatic sound speed, and Δr and Δz represent the average δr and δz of the cell:
Δr = ½ (x2 - x4 + x1 - x3) Δz = ½ (y2 - y4 + y3 - y1)
If the mesh is srongly rotated or distorted, more sophisticated Δr and Δz expressions may be required.
Note that the adiabatic sound speed should be used here, because the Lagrangian representation in this phase is accomplishing the changes in pressure through simultaneous changes in ρ and I (even though the latter is not being calculated at this point). This is in contrast to the purely Eulerian calculation described in previous papers on ICE methodology [1], in which the change of pressure through the interaction phase results from changes in density only, the full change in internal energy being calculated separately and incorporated into the pressure-change calculation as a separate step. In this purely Eulerian technique it is in the isothermal sound speed that is accordingly required in the implicit pressure-calculation phase.
(4) With D, S, and A defined, we can calculate the necessary pressure change for the cell.
δp = ωAS
where ω is a relaxation factor. Straight relaxation is given by ω=1. An optimum overrelaxation in many cases is ω=1.5 to 1.7, whereas ω>2 will lead to an unstable iteration.
(5) The convergence test is
|δp/pmax| < ε
where pmax is a current maximum pressure in the system. If the test fails for any cell, a flag is set to indicate that another iteration path through the mesh will be necessary.
(6) With δp calculated, we can update the ρL and pL values for the cell.
ρL = ρ(pL+nI) or ρL = ρL + (δp / c2) pL = pL + δp
(7) Now we can adjust uL and vL values for the four vertices of the cell:
uL1 = uL1 + δt/2M1 r1 (y2 − y4) δp uL2 = uL2 + δt/2M2 r2 (y3 − y1) δp uL3 = uL3 + δt/2M3 r3 (y4 − y2) δp uL4 = uL4 + δt/2M4 r4 (y1 − y3) δp vL1 = vL1 + δt/2M1 #189;(r2 + r4) (x4 − x2) δp vL2 = vL2 + δt/2M2 #189;(r3 + r1) (x1 − x3) δp vL3 = vL3 + δt/2M3 #189;(r4 + r2) (x2 − x4) δp vL4 = vL4 + δt/2M4 #189;(r1 + r3) (x3 − x1) δp
When all cells satisfy the convergence test, the iteration, using the above seven steps, is terminated. At this point, the quantities uL, vL, ρL, and pL describe the results of an implicit Lagrangian calculation that is not subject to the Courant condition.One could now move coordinates to complete the calculation if no rezoning were necessary. Note that ρL was calculated in terms of nx, not n+1x. The neglect of high order terms causes n+1ρ to differ slightly from ρL, but the approximation has not caused any difficulty. The ρL is used only in the pressure iteration, whereas in Phase 3, n+1ρ will be calculated from nρ by means of conservative fluxing in the mass equation.
In summary, at the end of Phase 2, we have in storage the n-time values of x, y, r, ρ, V, M, and I, as well as E ~ and the iterated values of uL, vL, pL and ρL.
The final phase of the ICED-ALE cycle computes the necessary rezoning changes i.e. convective and diffusive fluxes.
Assume at this point that a field of grid vertex velocities, uG and vG, have been assigned in some appropriate fashion with respect to a fixed, Eulerian reference frame. Thus for a purely Eulerian calculation,
uG = vG = 0
At the other extreme, a purely Lagrangian calculation would use:
uG = uL vG = vL
In general, the grid velocities may be any designated functions, and as such they are neither purely Eulerian nor purely Lagrangian.
There are two types of quantities to be updated in the rezone: cell quantities M (or ρ) and E (or I), and vertex quantities u and v, The procedure is to compute the cell quantities first, then change n+1Mcell to n+1Mvertex and compute the vertex quantities. Finally, n+1I can be calculated.
The rezoning can be accomplished using either the old nx, ny coordinates or the new n+1x, n+1y coordinates. The differences in rezoned quantities that results from these difference coordinate coordinates are of the order δt2, and they can be neglected for most purposes. Our procedure is to move the coordinates before the rezone calculations, as numerical methods are usually slightly more stable when time-advanced quantities are used. The new coordinates for all vertices are given by
n+1x = nx + uG δt n+1y = ny + vG δt and n+1r = n+1x for cylindrical coordinates or n+1r = 1.0 for plane coordinates
The mass and energy are rezoned on a cell-by-cell basis. For every cell, we must first calculate flux coefficients for each of the four faces, using the new coordinates:
Note that the flux conditions are zeo in Lagrangian calculations, as uG=uL and vG=vL, and that FR for cell (i+½,j+½) is equal to −FL for cell (i+3/2,j+½), and FT for cell (i+½,j+½) is equal to −FE for cell (i+½,j+3/2).
Recall that the momentum equations in the tilde calculations already contains diffuse terms, through a general stress-tensor deviator, which are used to represent true viscosity or to ensure computational stability. A slight instability results, however, if old-time values are used in the convective flux terms in the mass and energy rezone, although the stability of the mass equation is enhanced by the use of the partially advanced density, ρL. In general, it seems preferable to prevent the instability at its source, rather than to add a separate diffusion process. (The truncation errors responsible for instability are not really in a full diffusion form when more than one dimension is considered.) Therefore we will use flux expressions that can be adjusted towards a partial donor-cell treatment. It is convenient to embody the flux coefficiants, FR, FT, FL, and FB, within expressions that allow various differencing forms determined from input constants α0 and β0:
where sign FR, for example, is +1 if FR ≥ 0 and -1 if FR < 0, and the input constants allow these combinations:
α0 = 0 and β0 = 0 = centered α0 = 1 and β0 = 0 = donor cell α0 = 0 and β0 = 2 = interpolated donor cell
Note, however, that α0 must be sufficiently positive [2] fpr the mass equation o be stable. As full donor-cell differencing is too diffusive for most circumstances, generally 0 < α0 < 1.
The new mass and energy for a cell (i+½,j+½) are given by
Before updating the vertex quantities, we next calculate n+1V (as per the equation in Section I C, item (3), which in turn allows us to calculate
The new volume and density replace nV and nρ. The new vertex masses are then calculated using the Mij equation given in Sec. I C, item (4).
To adjust the vertex values of uL, and vL, for rezoning, we set initial values at all vertices,where
The second sweep adjusts the four corner vertex values of each cell in a manner analogous to that in the first two phases. For each cell in turn, we define several quantities:
Here, the sign has the same meaning as it did in the mass- and energy-flux expressions given above, and the α0 and β0 are the same quantities as in those expressions or may be chosen independently. Because a greater proportion of donor-cell differencing is required to stabilize the mass equation than the momentum equations, it is well to use a different (smaller) α0 in the momentum equations, Given the values of F13, F24, α13, and α24 for a given cell, the vertex contributions are given by:
Finally, the nev velocity field allows us to calculate the new specific internal energies for all cells:
where the u and v values are the n+1 values just calculated.
Various boundary treatments can be used in an ICED-ALE program [2], but we discuss here only the simple case of straight, rectangular reflective boundaries on all four sides of the mesh. (For ease of understanding, the version of the code presented in the following sections is limited to this one case.) The reflective boundaries considered are free-slip walls, the left boundary becoming the axis of symmetry for calculations in cylindrical coordinates. The criterion for any boundary condition is that velocities on boundary vertices be set in a suitable fashion. For free-slip walls, this means that normal wall velocities must be kept zero throughout the calculation. In the three phases of the calculational cycle, particular attention must therefore be given to the following:
(1) After the Phase-1 tilde velocity calculations, normal wall velocities must be reset to zero, i.e., ũ=0 on the left and right boundaries, and ṽ=0 on the top and bottom boundaries.
(2) During the pressure iteration in Phase 2, the normal wall velocities must be kept zero. Therefore, the appropriate uL and vL component(s) must be set to zero in boundary cells before proceeding to the next cell in the iteration.
(3) During the rezoning of cell quantities, cells adjacent to boundaries do refer to ρ and E values outside the walls, but these terms have zero coefficients, as they may be left unspecified.
(4) After the n+1u and n+1v calculations, the normal wall velocities must be zeroed again. in a manner analogous to that used for the Phase-1 tilde velocities. As described here, the normal velocities are set assuming that the boundary ia truly horizontal or vertical. Generally. however, any boundary (except the axis) may be curvilinear; then the normal velocity becomes a function of both the u and v components, requiring more careful treatment.
It is important to note that no pressure boundary conditions are required in YAQUI, This is a direct benefit of the Phase-2 iteration procedure.
It is useful to surround the mesh with a band of fictitious cells (described in Sec, II B) to aid in the treatment of the boundary conditions. Generally, ρ and E should simply be set forever to zero in the fictitious cells. This allows calculation of appropriate zero fluxes at the boundaries in Phase 3. In many applications, however, it is useful to allow the rigid walls of the mesh to expand in the rezone. Then fluid is swept up, and appropriate ambient values of ρ and E must be maintained in the fictitious cells. An example of this is shown in the sample code version included in this report, Here a uniform exterior E is generated in the setup and is allowed to remain constant for all time. The rezone calculates appropriate exterior ρL values to maintain atmospheric equilibrium. These new exterior ρL values subsequently become the final exterior n+1ρ values for the cycle. Rezoning is discussed further in Sec. III B.
Here we describe the principal structural details of the LASL ICED-ALE computing program, called YAQUI, whose flow diagram and listing appear in Appendix A and Appendix B, respectively. YAQUI was written aa a CDC-7600 production code for specific contractual purposes. As such, it embodies a number of features to make efficient use of computer storage and time. As was anticipated. however, the same basic code has been developed in several directions by a number of investigators, so it was purposely constructed in a modular form. The physical arrangement of these modules corresponds to their logical sequence in the computing cycle to the greatest degree practicable. The loss of efficiency in certain regions that results from having the entire code in FORTRAN IV, rather than machine language, is hopefully counterbalanced by increased readability for most users and the simplification of adapting it for use on computers other than the CDC-6600/7600 series.
As depicted in Fig. 4, YAQUI is built in an overlay fashion to minimize the use of small core memory (SCM), which is the fast memory on the CDC-7600, The main overlay, (0,0), which always resides in SCM, contains the main controlling program, YAQUI. Subservient to it are the programs in the two primary overlays, (1,0) and (2,0), which reside on disk storage. YASET is the setup program, and YAQUI1 performs the three-phase ICED-ALE caculations.
The structure within each of these three programs is further detailed in Fig. 5, which shows the UPDATE notation used in the actual code.
In addition to the main progrm, YAQUI, the (0,0) overlay contains the subroutine LOOP and FILMCO. LOOP handles the three-row buffering scheme that shuttles cell data between large core memory (LCM) and the SCM common YSC1. The details of cell-data storage snd the buffering scheme are given in Sec. II C and Sec. II-D. FILMCO (for film coordinates) computes the scaling for the microfilm plots. Because these two subroutines are required by both of the primary overlays, it is expedient to place them in the main overlay. Because they are thus always resident, the primary overlays can access them at will. Also in the main overlay is the common YSC2, which contains all the SCM data that must be maintained from cycle to cycle, and which is the SCM portion of the information written on tape for restarting purposes.
Two LCM blocks are initially defined in the main program: YLC1 is the storage block for cell data, and YLC2 is the storage block for the optional particles, described in Sec. II E.
To set up a calculation from initial input data. the main program calls YASET, the (1,0) overlay program, from the disk, and surrenders control to it. This overlay is placed in SCM immediately following the (0,0) overlay. YASET itself is only a two-instruction program: it prints YASET CALLED, so that the user can monitor his path through the overlays, and then immediately calls subroutine YASET1. YASET1 performs the actual setup, and in turn calls upon PARTGEN, to generate particles if specified, and MESHMKR, which creates the computing mesh with its initial cell and vertex quantities. When the problem setup is complete, YASET1 returns control to the (0,0) main overlay program.
To calculate after setting up, the main program calls YAQUI1, the (2,0) primary overlay, from the disk, and surrenders control to it. Because this is an overlay of tha same level as (1,0), it covers the image of (1,0) in SCM, being read in also to the locations immediately following the (0,0) overlay and thus making efficient use of SCM space. Like YASET, YAQUI1 is a two-instruction program: it prints YAQUI2 CALLED, and immediately calls subroutine YAQUI2. Should the prgram abort because of an unexpected error, the user can quickly ascertain which program he is in, inasmuch as the range of instruction addresses for both the (1,0) and (2,0) overlays is the same.
YAQUI2 is the largest section of code in the entire computer program. It contains the three-phase ICED-ALE, whose calculational cycles are repeated continuously under the direction of the Control Region. This region is strategically placed immediately after the np calculation, at which point in the cycle all the quantities that represent the complete solution at a given instant in problem time are available. The control region provides all microfilm plots and cell data prints. Also, it updates the problem time, t, by the current δt, performs tape dumps and tape restarts, and senses problem completion or an impending operating-system time limit. In the latter two events, it returns control to the main program, which, in turn, searches the input queue for further tasks, If there are none, the job is ended.
YAQUI2 makes use of two subroutines: PARTMOV moves and plots particles, and REZONE calculates new grid vertex velocities, uG, and vG, and new vertex coordinates, x, y, and r, for Phase 3 if the flow is neither pure Lagrangian nor pure Eulerian. If, however, the flow is pure Lagrangian or pure Eulerian, these velocities and the new coordinates are directly calculated in YAQUI2 in a simple, straightforward fashion. The REZONE package is really a roll-your-own section in which the user creates rezoning logic appropriate to his particular needs. (In the version of YAQUI presented here, the REZONE subroutine is an example of a possible way to follow the rise of debris from an atmospheric burst.)
To restart a calculation from a tape dump, the main program bypasses the (1,0) overlay and calls (2,0) directly, The restart condition is sensed immediately by YAQUI2, and the control region reads in the tape dump, placing the data in SCM and LCM as required, and turns control over to the point in the calculational cycle that will continue the problem from where it left off when the tape dump was made.
An examination of Fig. 2 shows some variables centered at vertices and some at cell centers, a common occurrence in Lagrangian computing methods, In FORTRAN, one can reference xji, simply as X(I,J), but pj+½i+½ cannot be referenced by a half-integer index, so the convention has evolved that P(I,J) refers to this pressure. Thus the indices I and J refer to a quantity lying at the lower left vertex of a cell, or at the cell center, depending upon where the quantity is defined. In YAQUI, (I,J) is replaced simply by (IJ), as only single subscripts are used for computer efficiency. In the YAQUI subscript notation, the letter P stands for +, and M stands for -. Thus, we write
IJ = (i,j) IPJ = (i+1,j) IJM = (i,j-1) IPJP = (i+1,j+1) etc
Such a notation permits easy readability of programmed difference equations in the code, Figure 6 shows the single subscripts typically seen in reference to vertex quantities, and Fig. 7 shows subscripts referring to cell quantities.
As the number of vertices in either direction is one greater than the number of cells, it is apparent that the grid in computer storage must be at least (I+1) by (J+1) in size. Because our indexing refers to cell centers and lower left vertices, we must allow one extra column of storage on the right and one extra row along the top. YAQUI includes one extra row along the bottom in addition, giving a mesh that is (I+1) by (J+2) in extent. These exterior zones are known as fictitious cells, and having them on three sides helps in the treatment of expanding meshes and certain boundary conditions, Note that fictitious cells are not used on the left, however. The code was basically intended for calculations in cylindrical coordinates, in which the left boundary is an axis of symmetry. In plane coordinates, it becomes a rigid free-slip wall, or plane of symmetry. The omission of fictitious cells on the left implies that no fluxing will ever be desired on that side, and the code would have to be modified to allow such a feature. The actual YAQUI mesh for the conceptual mesh of Fig. 1 is shown in Fig. 8. Coordinates are not calculated for fictitious cell vertices.
Obviously, double DO 1oops in FORTRAN to cover all vertices would have the limits J= 2 to JP2 and I= 1 to IP1. Similarly, loops to cover all cell centers would have the limits J= 2 to JP1 and I= 1 to IBAR.
The YAQUI code waa designed for running finely resolved calculations, implying several thousand computing cells. In addition to the basic fluid dynamics, space has been left in SCM for the later inclusion of code to deal with other physical phenomena, such as magnetohydrodynamics, turbulence, and chemistry effects, Therefore, all cell data are maintained on LCM, and the code deals with only three rows of cells at a time in SCM processing. Clearly, the optimum procedure is that which requires the minimum number of read/write references to LCM. Accordingly, the cell variables are stored in an interleaved fashion, in which all the variables for a given cell are stored contiguously, followed by all the variables for the next cell, etc. (Contrast this with the traditional method of storing cell variables in individual I by J blocks for each variable, This scheme is appropriate when the computing code is designed for smaller meshes that will always fit in SCM.) In the version of YAQUI presented here, the full calculational cycle, including the optional particles, requires 35 different cell variables, but we are able to get by with using only 14 storage words per cell. This is made possible by retaining quantities during a cycle only as long as they are needed, and then using their storage words for other quantities. Figure 9 shows the allocation of the 14 storage words for a YAQUI cell in the (1,0) and (2,0) overlays. The ordering from left to right corresponds to the actual order in which quantities are calculated in the code. A black dot implies that the quantity currently in the given storage word is referenced to calculate the quantity specified at the top of the column. The open dots (not clear which dots in the original diagram are open) in the rezone imply that x, r, y, and V may be referenced, depending upon the particular rezone. Note that the vertex masses, Mv and the cell volumes, V, are stored as reciprocals for increased computer efficiency. Because most references to M, and V are in denominators of equations, the time-consuming divide operation is thus avoided much of the time.
The quantities 2δt (1/Δr2 + 1/Δz2), known in the code as DELSM, and 1/c2, known as RCSQ, are invariant through the Phase-2 iteration. It is, therefore, expedient to compute their values throughout the mesh beforehand to avoid needless and repetitive calculation within the iteration itself.
In the convective-flux part of Phase 3, the n+1 values of Mc (the cell mass), 1/Mv (the vertex mass), and u and v are initially stored in vacant slots. Their n values are still required through the calculation of the momentum equations, after which the new masses and velocities can be transferred to their ordinary storage words. This places them in their proper locations as the n values going into Phase 1 of the next cycle.
The contour quantity (CQ) in the control region denotes the field of some chosen cell variable for which a contour plot is drawn on microfilm. The quantities referred to in the PARTMOV subroutine are described in Sec. II E.
Charts such as Fig. 9 have proven extremely useful in initially planning the storage before a code is written. but they are equally useful thereafter as an aid in visualizing what quantities are available at a given point in the calculational cycle, and where storage vacancies exist.
YLC1, the storage block for cell data on LCM, contains a single array, AA1, dimensioned at 13100010 words in the version of the code presented here. Because 14 words per cell are used in this version, a maximum of 9357 cells (the product of I+1 and J+2) are available.
Subroutine LOOP, in the (0,0) main overlay, shuttles the cell data between the large LCM array and a small buffer in SCM where it is operated on. Generally, LOOP maintains three rows of the grid in SCM at a time: the row being processed and the rows above and below. All calculations affecting cell data are actually performed directly on the current contents of the buffer. Because the cell data are interleaved in storage, all quantities pertaining to the three rows of cells are instantly available. The schematic flow diagram and sample FORTRAN double DO-loop in Fig. 10 show how the buffering takes place.
(1) Before the DO statements are entered, a CALL is made to the START entry of LOOP, START reads in the entire contents of the bottommost three rows of the grid from LCM to the SCM buffer, placing row j=1 in the buffer section designated row 1/3; likewise row j=2 is read into row 2/3, and row j=3 is read into row 3/3. Rows 1/3, 2/3, and 3/3 are contiguous in SCM, and like their counterparts in LCM each contains NQI=NQ * IP1 words, where NQ is the number of quantities, or storage words, per cell. With the three rows read in, the calling program needs to know how to access data in the buffer. Thia information is provided by the setting of the indices IJM, IJ, and IJP to point to the first words for the i=1 column or cells in each row. Thus, IJM is set to the first-word address (f.w.a.) of SCM row 1/3; similarly, IJ points to the f.w.a. of 2/3, and IJP points to the f.w.a. of 3/3. Note the indicator IBUF which is set to 1; it will control the subsequent reading and writing of individual rows and the resetting of the three indices. With the first three rows of cells read in and the basic indices set, control is returned to the calling program.
(2) The double DO loops are initiated. Secondary indices are needed for cells not lying immediately above or below cell IJ, so IPJ (= 1+1,J) and IPJP (= i+1,j+1), which initially refer to the i=2 column of cells, are easily obtained by applying increments of the number-of-storage-words-per-cell, NQ, to the primary indices IJ and IJP. In the example shown in Fig. 10, we are able to calculate the radius of a cell as the simple average of the radii of its four vertices. The terminal statement of the inner DO loop, which counts columns within each row, is statement No. 89. Note how the primary indices, IJ and IJP, are first advanced to the next column in the row. The inner loop is repeated until the row is completed, at which time control passes to the CALL LOOP statement.
(3) The LOOP entry immediately writes row IJM back onto LCM, and depending on the setting of IBUF, goes to statement No. 10, 20, or 30. Because IBUF was initially set to 1, control passes to statement No. 10 in our example. Now the indices IJP, IJ, and IJM are reset to point to different SCM rows - IJP to the vacated row 1/3, IJ to 3/3, and IJM to 2/3. IBUF is reset to 2 to control the next entry to LOOP and control passes to statement No. 40 which will read row j=4, the new IJP row, into SCM row 1/3. Note that no unnecessary shuffling of data in SCM has taken place; row j-1 was read out and replaced by row j+1, and the three indices were reset to point to where the rows j+1, j, and j-1 are located. As depicted at the bottom of Fig. 10, the grid rows in SCM are in their actual logical order only every third row.
(4) LP returns to statement No, 99 In the sample calling program, and rows are processed similarly until all the rows specified by the J DO-loop have been processed, at which time control passes to the CALL D0NE statement.
(5) The DONE entry is really only a cleaning-up operation: because further reading is unnecessary, it merely writes the final two rows, j and j+1 (JP1 and JP2, respectively) back out onto LCM.
Not indicated in the flow of Fig. 10 is the incrementing of the relative address indices for reading and writing LCM. These are initially set to 0, and incremented by NQI as processing progresses up the mesh.
Given this three-row buffering subroutine, the user needs only to include the CALLs to START, LOOP, and DONE at the appropriate points in the DO-loops and to increment by NQ words for each cell within a row. Other than that, the logic he must know is no more complex than if the data were entirely in SCM.
The number of storage words per cell in the YAQUI version presented here is seen to be NQ=14, as per Fig. 9. This may be increased very simply by adding the new variables to the EQUIVALENCE and DIMENSION statements in the Comdecks EQVREAL and DIMEN, respectively, and redefining NQ in the (0,0) main program at one place only.
The SCM buffer is in common YSC1 which contains a single array, AASC, dimensioned at 422410 words in this YAQUI version. Because AASC must be able to hold three complete rows at once, NQ=14 means I<100.
Other entry points in the LOOP subroutine, not shown in Fig. 10, allow the user to access any cell at random easily, and to perform DO loops with the J index reversed so as to sweep from top to bottom. Examples of this flexible routine are seen in numerous places throughout YAQUI.
The basic ICED-ALE scheme has no dependence upon particles, but deals entirely with field variables related to the computing grid. It is useful, however, to include particles in many problems of interest. These may serve either as true markers, which are carried along by the flow but have no influence upon it, or they may act upon the surrounding fluid. In the latter case, we must store velocity components, a mass, and a drag coefficient for each particle, in addition to the usual coordinates. This permits calculation of momentum changes experienced by the particles owing to fluid forces, these changes in turn being subtracted from the fluid momentum on a cell-by-cell basis.
For simplicity, we shall first discuss the basic particle-moving scheme used in YAQUI, and then describe the inclusion of the momentum-exchange feature.
Our technique for moving particles in a general, quadrilateral grid is based on use of a temporary, uniform rectangular cell grid superimposed on the YAQUI grid, Given a velocity field related to a unifom grid, particles are easily moved by an ordinary interpolated area-weighting scheme. The problem, then, is to define a velocity field on the superimposed grid (hereafter called the particle grid) so that it reasonably approximates the velocity field of the current YAQUI grid. This is done as follows.
(a) First, we define the particle grid by specifying its cell-edge lengths, Δx and Δy, and its overall dimensions PXR and PYT. Because we use a vacant part of the YAQUI cell storage to store particle grid quantities, we do not allow the number of zones in either the r or z direction to exceed that of the YAQUI grid. We generally, however, run at this maximum for the best resolution. In most calculations, the dimensions Δx and Δy are chosen so that the particle grid just encompasses the region covered by the regular rectilinear or curvilinear grid, although in some cases when particles are used solely in some specific region, the particle grid may be placed only over the region of interest.
(b) The second step, after definition and location of the particle grid, is to sweep the YAQUI vertices systematically and do the following for each.
If the vertex lies outside the particle grid, skip to the next vertex. To be included, the vertex must have x<PXR and PYB<y<PYT, where PYB is the y coordinate of the bottom of the particle grid.
Determine (i,j) of the particle-grid cell that contains the vertex.
Assign to each of the four corners of the particle grid cell (i,j) x and y moments (Mx and My) and mass (M0 according to
where w and h are defined as shown in Fig, 11. Use similar expressions with the same weighting factors for My and M0.
(c) Finally, after moments and mass have been assigned from all YAQUI vertices onto the appropriate particle-grid vertices, we calculate the u and v velocities of the particle-grid vertices as PU=Mx/M0 and PV=My/M0 which are the velocities to be used in moving the particles.
(d) To move a particle, we first determine in which cell (i,j) of the particle grid it is located. It is then moved according to an interpolated velocity based on the corner velocities of the particle-grid cell. If the particle lies at position (w,h) as shown in Fig. 12, then
and similarly for the vp velocity and the yp coordinates.
In particle-fluid momentum exchange the particles do not necessarily move with the fluid velocity. The basic particle mover is extended to include reaction with the surrounding fluid as follows.
(a) In addition to storing PU and PV for all particle-grid vertices, we store the quantities ΔMx and ΔMy, which are the x and y momentum changes made by fluid forces suffered by the particles near each vertex. These quantities must therefore be subtracted from the fluid momentum on a cell-by-cell basis. In YAQUI, we felt that to within the accuracy of the particle mover itself, we could split the word storage used for particle-grid velocities, combine FU and ΔMx at one-half word each, and similarly combine PV and ΔMy
(b) The calculation of ΔMx and ΔMy proceeds as follows. The velocity of each particle is governed by the equation of motion
in which uf1 is the up of the previous equation - the area-weighted value of the particle-grid velocities at the particle location, urand is the velocity contribution from turbilaent fluctuations [4], and np is a drag coefficient. The x-momentum change of the particle is
where mp is the particle mass. This momentum is distributed to the particle-grid vertices in much the same manner that uf1 was calculated. Thus, if the particle is in cell (i,j), the corresponding changes at the vertices are given by
A similar distribution is performed for ΔMy, where
Because ΔMx and ΔMy are calculated through a summation, their values must be initialized at zero each cycle.
Note that whereas the basic particle mover required that the particle coordinates (xp and yp) be stored from cycle to cycle, the momentum exchange requires in addition up, vp, mp, and np. In YAQUI, particle-storage words are split like particle-grid storage. Thus, the six quantities per particle are kept in three words, where xp is combined with up, yp is combined with vp, and np is combined with mp.
(c) After all particles have been moved and their momentum changes recorded at the particle-grid vertices, these momentum changes must be inserted into the fluid momentum field. This is done by sweeping the YAQUI vertices in the same manner as that used to set up the particle-grid velocities: Determine in which particle-grid cell (i,j) each Lagrangian vertex, is located. Because the mass, M0, associated with each particle-grid vertex is still in storage, the change in velocity components of the Lagrangian vertices can be calculated easily, The adjusted velocity component, u, of Fig. 11 is given by
and v is given similarly, with ΔMx replaced by ΔMy. These expressions conserve momentum.
.The YAQUI particle mover has been written with the momentum-exchange feature built in. To calculate with true marker particles only, however, we merely set all mp=0, np=0 and urand=vrand=0 and bypass all ΔMx and all ΔMy calculations.
Two-fluid dynamics can be performed without using particles in a purely Lagranian manner when the fluid distortion are not severe [2], whereas for incompressible flows involving large distortions, two-fluid dynamics can be calculated by tying the particle motions strongly to the fluid in which the particles are embedded.
The particle passes are chosen so as to supplement the density already contributed by the background fluid, the sum of that density and the particle density being the total density of the second fluid. (More generally, the presence of a spatially varying density in the second fluid can likewise be represented by appropriate choice of particle masses.) The masses can be negative or positive.
In the absence of a free surface, the effects of gravity are most efficiently represented by separating the pressure of the background fluid into two parts, the uniform gradient in equilibrium with gravity and the departure from this. As a result, only the departure pressure is obtained by iteration, its boundary condition being zero gradient on the top and bottom walls. The gravitational acceleration on the particles (i.e,, on the difference between the densities of the fluids) then remains as the only exterior force field. To allow for this, one must accordingly supply a separate specification for the gravitational acceleration on the particles, designated by gzp.
Particle storage in YAQUI is maintained in the LCM block named YLC2, which contains a single array, AA2, dimensioned at 13100010 words. Because the particle data are stored using three words per particle, a maximum of 43,666 particles may be used in the version of the code presented here.
The automatic calculation of the time step, δt, is included as an option in YAQUI, primarily on the basis of two stability conditions, one of which is imposed by the viscous stresses, with coefficients λ and μ, and the other of which is associated with the convective fluxes for Eulerian calculations, or with the prevention of negative volumes for Lagrangian calculations,
The viscous-stress stability conditions are tested in the Phase-1 calculations in conjunction with the calculation of the viscosity coefficients for the stress-tensor terms. As described below, the code creates effective values of λ and μ on a cell-by-cell basis, as determined by a combination of input quantities and local flow conditions. From these considerations, it then calculates a tentative δt, labeled δtv for use in the next cycle of calculation.
The convective-flux limitation can be imposed during the rezoning of M and E in Phase 3, From an examination of the flux coefficients FR, FT, FL, and FB, calculated for each cell as defied in Sec. I F, along with the divergence, D, and α0, the code obtains another competing tentative δt, labeled δtc.
The δt actually chosen for the next calculational cycle, then, is the smaller of δtv and δtc, δt=min(δtv, δtc). Subsequently, the initial values of δtv, and δtc, that go into the next cycle's tests differ by some factor, δtfac, which is usually slightly larger than unity, times the new δt just chosen. This permits the t to increase when conditions become more stable. Because the δt is always chosen for the next cycle of calculation, it can be argued that it is always a cycle behind. Ideally, the δt chosen should be for use in the present cycle, as it is based upon present conditions, but this would be more difficult to accomplish. The one-cycle lag, however, presents no problems, as the δt is always small enough that significant changes in the flow field occur only over a number of cycles of calculation. Accuracy considerations alone demand this, in addition to the requirements of numerical stability.
There is a great deal of latitude in how the viscosity coefficients may be determined for Phase 1. Governing the use of the input values of λ and μ is the input quantity ξ, an integer exponent used in conjunction with ρij. Three possible forms of viscosity are allowed, depending upon the definition of ξ:
(1) ξ=1 will allow a read-in value of artificial kinematic viscosity. The input values of λ and μ must be chosen with regard to the numerical stability requirements for expected flow conditions [2].
(2) ξ=0 is used when the input values of λ and μ represent the real, physical coefficients of viscosity.
(3) ξ = -1 forces the code to seek its ovn viscosity on the basis of local numerical-stability conditions in the flow. Note that the actual numerical values of the input λ and μ are immaterial when ξ = -1. Only the ratio λ/μ will be considered for dividing the total viscosity between λ and μ.
The effective λ and μ used in the viscous terms of the equations are, in all three cases, given by
when ξ=1 or 0. When ξ=−1, k is determined directly from the numerical-stability requirement:
(λ + 2 μ) / ρ < > ½u2 δt + ½u' δx2
We define
where ε is a coefficient, u12 is the square of a representative velocity at vertex (ij) of the cell
u12 = (u2 + v2 )ij
and u'δx2 is approximated by the maximum uδx of the cell times a factor (1/n),
in which #916;r and Δz have the usual definitions of average δr and δz:
Δr = ½(x2 − x4 + x1 − x3) Δz = ½(y2 − y4 + y3 − y1)
Use of ξ has removed the restriction for infinitesimal δt's that would otherwise be required in very low-density regions, as when ξ=1 or −1, λ=ρ times some quantity, and μ=ρ times some quantity, so the condition we must satisfy for stability,
and the dependence upon ρmin has been entirely removed. Moreover, for ξ=1, λ has been converted to a kinematic form more convenient for artificial viscosity in problems involving large density variations.
The λD term appearing in the Πxx, Πyy, and Πθ equations is calculated as
as D is applied only in compressive regions, that is, when it is negative.
With the viscous effects included through the stress tensors, the crucial question for determining δtv is the stability condition
which roughly states that momentum must diffuse less than one cell width per time step. Because (λ + 2μ)=4/3 μ + ξ, the right side of the above expression is always positive. Further, the alternate node coupler in Phase 1 introduces another stability condition, which can be shown to always be ANC(≡ 1/anc) < 1.
Combining these two conditions, we obtain
from which δtv may be reset as
δtv = min (δtv,Quantityij) ,
thus allowing every cell a chance to participate in the selection.
As mentioned earlier, several criteria in Phase 3 influence the choice of δtc. One requirement is that material cannot be fluxed more than one cell width per time step, as the flux approximations are based on the implicit assumption that exchanges occur only between adjacent cells. Therefore, δtc must be based in part upon the quantity
|max(|FR|,|FT|,|FL|,|FB|)/Volume|ij
if the flow has any Eulerian features. In Lagrangian cases, Dδt can provide the same measure that flux/volume does for Eulerian cases, as both expressions have the appropriate δu/δx δt dimensions,
Besides monitoring these two quantities, δtc must take into account the differencing scheme itself. It can be shown that for stability we must require
Uδt/δx < (2α0) / (1+α02)
in which U=ufluid−ugrid, and α0 is a measure of the donor-cell proportion in the mass equation, as described in Sec. 1F. (The right side of the above condition has its maximum when α0= 1.) For accuracy, we restrict the limit to only a quarter of this amount:
α0/2(1+α02)
Combining these three conditions, then, yields the crucial quantity for determining δtc:
according to which we reset δtc, if necessary, by
δtc = min (δtc, Quantityij )>
on a cell-to-cell basis, as we did to calculate δt (For computational purposes, the denominators of both the above and the previous Quantityij expressions should contain an added constant (on the order of 10-10) to ensure that they do not vanish.)
The generation of the initial mesh and fluid configuration in YAQUI is the responsibility of the (1,0) subroutine, MESHMKR. This subroutine must provide the starting x, y, r, u, and v values for all vertices, and ρ and I values for all zones. Data punched on input cards, described fully in Sec. C below, provide MESHMKR the information necessary to perform this task for a variety of circumstances.
The first consideration is to define the coordinates for all vertices. The input quantities I, J, δr and δz (= IBAR, JBAR, DR,and DZ in the program) are the four fundamental quantities in grid generation. They permit creation of a grid of uniform δr by δz zones whose origin, vertex (1,2), lies at coordinates (0.0, 0.0), The addition of a fifth quantity yB (= YB), the y coordinate of vertex (1,2), allows the entire mesh to be displaced upward. This is useful for calculations involving expanding meshes. The initial, basic part of MESHMKR generates exactly this uniform grid.
The version of MESHMKR presented here further allows the option of nonunifom zoning. As depicted in Fig. 13, the previously generated grid lines cay be shifted vertically and horizontally, with zone size increasing continuously outside of some remaining inner area of uniform zones, The region of uniform zones occupies Iuniform (=IUNF) by Juniform (=JUNF) zones, centered at Jcenter (JCEN) zones up fro the j=2 bottom boundary line. IUNF and JUNF may range from values of 1 and 2, respectively, implying variable zones throughout, up to values of IBAR and JBAR, implying uniform zones throughout. The input coefficient FREZ provides the expansion ratio for the zones lying outside the IUNF by JUNF region. A relationship of the form xi=xi-1 + FREZ (xi-1-xi-2) is used to locate grid lines lying to the right of IUNF, above JCEN+JUNF/2, and below JCEN-JUNF/2. For accuracy, FREZ generally should not exceed about 1.1. The above expression will retain uniform zoning throughout if FREZ=1.0. A simple program modification would allow for different expansion rates in the two directions.
If variable zoning is employed (FREZ > 1.0), user calculation of YB would be inconvenient, so, instead, the input quantity REZYO, which is the y coordinate of the center of expansion, y0, determines the vertical placement of the mesh. REZYO refers to the YAQUI vertex (1, JCEN + 2), and allows YAQUI to calculate the actual value of YB.
Although the variable zoning shown in Fig. 13 for this version of MESHMKR is of a simple rectilinear form, we emphasize that neither the technique nor the code is by any means limited to this. MESHMKR may be modified easily to create any curvilinear grid, and, indeed, simple iterative techniques [5] have been used in MESHHKR to define a variety of more complex grid shapes for special applications.
With the basic grid (x, y, and r values) defined, the second consideration is the initial u, v, ρ, and I values to define the fluid. Input data cards are read which define regions containing integral numbers of zones, and specify the initial values of the four variables assigned to all zones lying within each region. Use of these data cards is fully described in Sec. C below.
This version of MESHMKR includes the special capability of setting up atmospheric-explosion calculations. In this case, it creates an entire background of ambient atmosphere through use of one of the fluid-region data cards, In addition, ρ and E values are initialized in the external or fictitious zones, as the grid will later be expanded in the REZONE. MESHMKR then adjusts the uniform ρ field of this ambient atmosphere to a state of gravitational equilibrium by use of an algorithm that accounts for nonuniform zoning.
When combining variable zoning with an equilibrium atnosphere, the user must consider zone height in relation to scale height. It is calculationally inaccurate to allow a zone height to exceed one scale height, and for zones larger than this, a change in sign can be introduced into the ρ field, Therefore. it is wise to ensure that the condition
δz < (γ-1)Iamb / |gz|
is satisfied throughout the mesh. In the above expression, Iamb is the ambient specific internal energy. As coded in this report, Iamb (given by REZSIE) is a constant, but when atmospheric conditions allow it to increase with increasing altitude, larger zones may be employed in the region of increased Iamb. However, the above condition must still be satisfied.
Upon the ambient background, the spherical burst may be defined in any manner that the user wishes. We usually employ a special set of data cards to define the upper right quadrant of the burst. These cards are provided by a one-dimensional spherical code whose purpose is to calculate the early-time dynamics. The cards are arranged in relation to a set of uniform Eulerian zones, each data card in the set specifying a pair of relative i and j cell indices and the associated ρ, I, v, and u values. With a j index specified in YAQUI to correspond to the center of the burst, MESHMKR creates only the upper right and the mirror-image lower right quadrants, taking advantage of the cylindrical symmetry of the burst. The data are superimposed over the previously defined ambient atmosphere, overwriting a part of it that is restricted to lie within the IUNF by JUNF area, whose uniform zones are identical in size to those of the one-dimensional code. (This restriction would be unnecessary if one were to interpolate the input data separately.) The velocity components specified on the data cards are located at cell edges, creating a minor complication, as YAQUI velocities must be centered at the vertices. As a result, MESHMKR. must store these velocity values in temporary locstions as they are read in. After the entire set of data cards has been processed, MESHMKR can transform the field through appropriate averaging to form a vertex-centered velocity field.
Whatever logic the user chooses to employ in grid generation, MESHMKR's work is finished when all initial x, y, t, u, v, ρ, and I values have been defined and appropriately stored, This information enables YASET1 to calculate the initial values of the remaining basic cell and vertex quantities (Mc, V, E, and Mv) in a straightforward manner.
Rezoning, which is grid motion relative to fluid motion, occurs in any flow that is not purely Lagrangian. Indeed, purely Eulerian flow is a rezone flow, and is unique only in that the grid motion is such as to maintain the grid in a fixed location, When rezoning occurs, there is a convective flux of mass, momentum, and energy from one zone to its neighbors, which must be properly accounted for. For fluxing accuracy, the grid velocities or the time step must be restricted so that fluid is fluxed less than one cell per cycle, as it is assumed in the equations that exchanges take place only between neighboring zones.
These considerations are dealt with in Phase 3, in which grid velocities, uG and vG, and from them the resulting new x, y, and r coordinates, are determined. For the extremes of purely Lagrangian and purely Eulerian flows, these grid velocities may be specified quite simply. In the Lagrangian limit, the grid velocities are identical to the Lagrangian velocities resulting from Phase 2, uG=uL and vG-vL. In the Eulerian limit, the grid velocities are identically zero. In YAQUI, these two cases are treated in Phase 3 in YAQUI2 itself. The (2,0) subroutine REZONE is called to define the grid velocities and the resulting coordinates for any flow that is neither of these two extremes. This roll your own subroutine allows the user to force the grid to follow one or more features of the flow continuously.
The sample REZONE included here shows just one of a variety of possible schemes for following the dynamics of an atmospheric-explosion calculation. It is a good example because it shows the three basic objectives that must be met by a rezone for such a calculation.
(1) The mesh must be expanded so as to maintain the boundaries ahead of the strong radially expanding shock.
(2) The mesh must rise in the atmosphere at the rate of fireball rise, and
(3) The zoning must resolve the details of the torus in the central region finely, but may be much coarser in the outlying regions, for computer efficiency.
(The variable zoning discussed in Sec. III A and shown in Fig. 13 can provide a good beginning for thie last aspect.) The following briefly explains how REZONE meets these three objectives in this version.
(1) The mesh expansion is controlled by monitoring the largest absolute uL, or vL fluid velocity (umax) along a column or (vmax) along a row, several cells in from each rigid boundary, thereby allowing signals to be sensed before they can reach the boundaries. The normal grid velocity assigned to the boundary vertices is then calculated to be the square root of the product of this maximum velocity times the largest absolute uL or vL velocity (Vmax) in the entire grid:
(uG)i=1 = 0 for all j (uG)i=IP1 = [Vmax (umax)i=I-6]½ for all j (vG)j=2 = [Vmax (vmax)j=9]½ for all i (vG)j=JP2 = [Vmax (vmax)j=J-14]½ for all i
(2) The overall upward rise or translational velocity (vT ) of the mesh can be determined by tracking the rising maximum point of some representative feature of the flow. We have found that the vorticity will serve this purpose if care is taken. Although the rising fireball torus will soon develop a strong vorticity field, the vorticity profile flattens with time, developing a vertically elongated plateau of the larger values, Upon this plateau, the maximumn point itself may move around significantly from cycle to cycle. If the grid translation is tied to such a shifting point, the result is discontinuous up and down translation, perhaps moving the grid several zone heights all at once. A smoother and more reliable quantity to follow than the pure vorticity would be some weighted average vorticity. One possible form that we have used successfully is based upon the quantity
which is summed over all cells except those near the rigid boundaries. Then vT is calculated as
where ycen is the y coordinate about which the fireball should be kept centered, and ΩG is an under-relaxation factor used to ensure smooth rezoning.
(3) The technique for rezoning the interior grid lines also makes use of the center of maximum vorticity, requiring the radial center,
as well as the vertical yc, to move. The interior vertices are then made to satisfy the relations
for all j, and
where the coefficient BG determines how tightly the vertices are drawn in towards the center of vorticity (xc,yc), and therefore governs the level of resolution in the fireball torus.
In terms of grid velocities, these results are obtained by setting
The minimum zone size is related not only to the value of βG, but also to the value of I. For a given I, it is helpful to be able to estimate a priori an appropriate value for BG to achieve some desired level of resolution. The relationship can be shown to be
where δrmin is the minimum δr after relaxation of the grid, xR is the x coordinate of the right boundary after all grid expansion has taken place, and βF=(2βG+βG2) ½. The procedure is to obtain solutions for various values of βG, but the final βG chosen as optimum for a given problem will probably differ slightly from the value suggested by the relationship, and it is generally obtained empirically.
The rest of the REZONE subroutine, from statement No. 1200 to the end, is somewhat more general than the preceding part. New values of x, y, and r are calculated for all vertices, using the values of uG and vG in expressions identical to those in Phase 3 of YAQUI2, This is done in REZONE, however, to enable the following adjustments to be made before RETURNing.
(1) The particle grid parameters are adjusted to fit the new grid.
(2) Subroutine FILMCO is called to adjust all the film-plot scaling parameters, and finally,
(3) The ρL values in the exterior zones are recalculated using the new coordinates.
Formatted input data cards provide the information necessary to specify a problem setup. The number of cards required varies according to the problem. However, the following cards must always appear.
Card No. 1: IBAR, JBAR, IUNF, JUNF, JCEN, DR, DZ, CYL, GRDVEL, A0, A0M, B0, KXI (Format 5I4, 7F8.3, I4), where:
IBAR= I, the number of real zones in the r or x direction.
JBAR= J, the number of real zones in the z or y direction.
IUNF, JUNF, JCEN, and FREZ (see Card No. 4 below) allow one form of variable orthogonal zoning in the initial grid generation. Refer to Sec. III A and to Fig. 13.
DR= δr, the cell size in the r or x direction in the uniform region.
DZ= δz, the cell size in the z or y direction in the uniform region.
(Note: The user may wish to completely override the specifications of IUNF, JUNF, JCEN, DR, and DZ in MESHMKR.)
CYL=1.0 for cylindrical geometry, or =0.3 for plane geometry.
GRDVEL= grid velocity, 0.0 = pure Eulerian, 1.0 = pure Lagrangian, 2.0 = REZONE.
A0=α0 coefficient in the Phase-3 momentum equations.
A0M= α0 coefficient In the Phase-3 mass and energy equations.
B0= β0 coefficient in the Phase-3 mass, energy, and momentum equations.
FXI=ξ, the exponent of ρ that determines the form of viscosity in the problem. Refer to Sec. II F.
Card No. 2: NAME (Format 10A8), where columns 2-80 of this card are used for problem identification on prints and plots. Column 1 should not be used because it is treated as a carriage control. If desired, the card may be entirely blank, but it must always be included.
Card No. 3: MU, LAM, OM, EPS, GR, GZ, ASQ, RON, GM1 (Format 9F8.3), where:
MU= μinput, LAM= λinput Viscosity coefficients. Refer to Secs. 1 D and II F.
OM=ω, the Phase-2 iteration relaxation parameter. The value ω>1 provides overrelaxation, whereas ω>2 is unstable. Refer to Sec. I E.
EPS= ε, the Phase-2 iteration convergence criterion, typically on the order of 10-5 (specifying convergence to within 10-5of the maximum pressure in the system at a given instant), but ε may be greater or smaller depending upon the problem. Refer to Sec. IE.
GR=gr, gravity felt by the fluid in the r or x direction, which may be + or − to pull rightward or leftward, respectively.
GZ=gz, gravity felt by the fluid in the z or y direction, which may be + or − to pull upward or downward, respectively.
The quantities ASQ, RON, and GM1 are applicable to the stiffened gas equation of state, which appears in the code version of this report.
ASQ= a2 the zero-temperature sound speed.
RON= ρo, normal material density.
GM1= (γ-1), where γ is a constant characteristic of the gas, becoming the ratio of specific heat at constant pressure to the specific heat at constant volume, γ= cp/cv, if the gas is truly polytropic.
Card No. 4: PREZ, YB, REZYO, REZUE, REZVE, REZVT, REZRN, REZSIE (Format 8F8.3), where: FREZ, YB, and REZYO are parameters relating to the zoning and grid location in the initial grid generation, Refer to Sec. III A.
REZUE,REZVE, REZVT Available for use in REZONE to specify grid expansion (ue,ve)) and translation (vT velocities, if these velocities are constant.REZRON= the ρ0 of the ambient atmosphere at altitude REZY0 at t=t0.
REZSIE= the specific internal energy of the ambient atmosphere. In the code listing in this report, REZSIE is a value that remains constant in space and time.
Card No. 5: IBP, JPX, PDR, PDZ, PYB, GZP, IMOMX (Format 214, 4F8.3, I4). This card supplies the parameters for the optional particles deocribed in Sec. II E.
IBP= Iparticles, JBP=Jparticles If no particles are to be used, set IBP=0. Then the rest of Card No. 5 is unused, so proceed to Card No. 6. For particle usage, IBP and JHP are the I and J of the particle overlay grid. IBP≤IBAR, and JBP≤JBAR.
PDR=Δx,PDZ=Δy, the (uniform Δx and Δy of the particle-grid cells. See Fig. 11. In variable-zoned meshes, these values are calculated automatically in the setup. Similarly, PDR and PDZ are recalculated by REZONE. In both places, the code version presented here stretches the particle grid to just cover the farthest points on the bottom, top, and right edges of the YAQUI grid.
PYB= YBparticles, the displacement of the particle-grid lower edge measured relative to YB. To superimpose exactly, set PYB=0.0 and allow the code to adjust the particle grid automatically, as described above.
GZP=gz felt by the particles, which may or may not be equal to GZ (see Card No. 3).
IMOMX-1.0 for the momentum-exchange option, =0.0 otherwise.
Card No. 6: T, DT, T20MD, TLIMD, TFIN, LPR, ICOLOR (Format 5F8.3, 214), where:
T= t0, the problem starting time, usually zero.
DT=δt0, the initial δt. The first cycle is automatically run with δt=δt0/10, then the second and third cycles are run with δt=δt0. From cycle 4 on, the δt is chosen automatically as described in Sec. II F.
T20MD=1.0 to force tape dumps every 20 min of central processor (CP) time for restarting purposes, or = 0.0 to bypass this option.
TLIMD=1.0 to force a tape dump and RETURN to the (0,0) overlay just before reaching the CP time limit specified on the JOB card; > 1.0 to force tape dump and RETURN immediately after cycle 0 output; =0.0 to run out to a full time limit with no tape dump.
TWFIN= problem finish time. When this is reached (t > TWFIN), control will RETURN to the (0,0) overlay.
(Upon RETURN to (0,0) for either the TLIMD or TWFIN condition, the (0,0) main program YAQUI searches the input queue for further tasks.) 1= zone prints on microfilm only, 2= zone print on both film and printer, 3 = zone prints on printer only. These are described more fully in Sec. III D.
IC0LOR > 0 plots particles in red, and anything else on film in white, obviously effective only with color processing. ICOLOR < 0 implies normal black-and-white processing,
Card No. 7: (DTO(N), N=1,10) is used in conjunction with
Card No. 8: (DTOC(N), N=1,10) (both are Fornat 10F8.3), where DTOn specifies the problem-time output interval for both plots and prints. DTOCn, specifies the time at which to change to DTOn+1. As an example, assume that t is in seconds, and that output is wanted every 1/4 sec for the first second, then every l/2 sec up to 4 sec of problem time, then every 1/2 sec until t = 10, then every 2 sec until t = 50, and every 10 sec until t = 200. One would use
DTO(1-10) = 0.25, 0.5, 1.0, 2.0, 10.0, DTOC(1-10) = 1.0, 4.0, 10.0, 50.0, 200.0
To keep the output time interval fixed throughout run, specify DTO(1)=(interval) and DTOC (1) > TWFIN.
(Note: When an output time is being approached, the automatic δt routine will choose a special δt for one cycle so that the output occurs at the precise time desired).
The above eight cards pertain to all YAQUI setups. They have defined a basic grid and provided the parameters for its use. What remains to be defined is the contents of this grid - particle regions and fluid regions. Because tbese regions vary with the problem geometry, the number of cards in the rest of the input deck varies widely. The procedure beyond Card No. 8 is to define the particle regions first, if any exist, then finally to define fluid regions.
Card No. 9: DRPAR, DZPAR, XC, YC, XD, YD, UPAR, VPAR, MTE, DRAG (Format 10F8.3). This is a particle-region card, to be expected only if IBP> u on Card No. 5. One card of the above format must be provided for each discrete particle region in the mesh. In the present version of PARTGEN, a particle region may be one of two shapes - cylindrical or spherical (rectangular or circular in plane geometry). These two general shapes are shown in Fig. 14 with the named dimensions that specify them, The four dimensions (XC, YC, XD, and YD) are input in true distance units because the particle regions are not constrained to follow some edges. For a cylinder, XC and YC specify the coordinates of the lower left corner, and XD and YD specify those of the upper right corner. For a sphere, YD must always be identically zero to enable PARTGEN to distinguish it fram a cylinder, YC specifies the position of the center, measured up the axis, and XD specifies the sphere radius. Note that the y dimensions are defined relative to y=0, not relative to the bottom of the mesh, (This was allowed so that particles might origiaally be placed outside an expanding mesh, but the user should not try to move any particle while it is still outside the mesh, as the present logic in PARTMOV assumes that all particles lie within the mesh.)
DRPAR, DZPAR: Particle spacing in the r or x and z or y directions, in problem units.
XC, YC, XD, YD: Particle-region dimensions in problem units, relative to x = 0 and y = 0. See Fig. 14.
UPAR,VPAR: Initial u and v velocity components for the particles in this region; will be 0.0 for true marker particles.
MTS: Mass per particle =MTE*rparticle. Use MTE=0.0 for markers.
DRAG: Drag coefficient η for these particles. Use DRAG=1010 for markers.
Particle-region cards are processed individually, and the number of particle regions is unlimited. If particles are used at all, the set of particle-region cards terminates with the final card having DRPAR<0.0 and the rest of the card is unused. Therefore, the number of particle-region cards in a YAQUI deck is either zero, or two or more.
Card No. 9: if no particles are used. If particles are used, however, then this card follows the DRPAR=0.0 card: NB, MR, NT, NL, UI, VI, RO1, SIEI (Format 4I4, 4F8.3). This is a fluid-region card, one card of this format being provided for each discrete fluid region in the mesh. The allowed fluid region covers some specified number of zones, as shown in Fig. 15 with the named dimensions that define it. The four dimensions (NB, NR, NT, NL) are given in integer numbers of cells to emphasize that the four corners of the region must coincide with cell vertices.
Thus, NL and NB specify how many cells in from the left and up to the vertex where the lower left corner of the region is located, and NR and NT similarly locate the vertex of the upper right corner. Even if the grid is not originally orthogonal, specifying two diagonal corners uniquely specifies the zones that will be included in the region. To use a single fluid region as an entire ambient background, set NL=NB=0, NR=I, NT=J.
NB, NR,NT,NL: number of zones (integers only). Refer to Fig.15.
UI=uI, VI=vI: the initial velocity components to be assigned to all vertices in the fluid region.
ROI = ρI, SIEI=SIEI: the initial density and specific internal energy to be assigned to all zones in the fluid region.
Fluid-region cards, like particle-region cards, are processed individually, ard the number of fluid regions is also unlimited. The version of MESHMKR presented in the code listing in this report expects that at least one fluid region will be defined by this type of card. The set of fluid-region cards terminates with the final card having NR = 0 and the rest of the card unused. Therefore, a minimum of two fluid-region cards must be present in a YAQUI deck,
To set up an atmospheric-explosion calculation, as described in Sec. III A, the final NR card (following the ambient fluid NR card) has NR=1000, instead of NR =0, with NT = the j index of the burst point in the full YAQUI grid. Generally, NT=JCEN+2, where JCEN was defined on Card No. 1. The set of special data cards follows the NR= 1000 card, and contains I, JJ, ROI, SIEI, VI, and UI (format 215, 4(4X,E11.5)), one card per Eulerian zone. The j index is called JJ here to emphasize that it is relative to the definition j=1 at the burst point on the data cards. VI is the v velocity centered on the top edge of the Eulerian zone, and UI is the u velocity centered along the right edge. These cards are read and processed individually, the set terminating with a card having I=0.
This completes the discussion of the input data cards. The final card normally placed at the end of the input deck is in reality the first card for the next problem. The first quantity on Card No. 1 is IBAR, and its value determines the action to be taken by YAQUI. If IBAR > 0 it is valid for use as I, and YASET is called. The value IBAR = 0 indicates a tape restart, and IBAR < 0 indicates that the end of data has been reached. Thus a negative IBAR card is the appropriate way to terminate a deck, and hence, the run.
The particle plots are made by plotting the xp and yp coordinates of all particles, and are provided automatically when particles are used.
The zone plot is included for all Lagrangian or REZONE runs (GRDVEL > 1.). For purely Eulerian calculations (GRDVEL = 0.), the zone plot is provided only at cycle 0. The labels of minimum and maximum δr and δz on the zone plot are really unambiguous only for orthogonal grids. The general form used in their calculation was intended to make the labels meaningful for slightly distorted grids.
The velocity-vector plot shows the direction of fluid flow and the relative magnitude of the velocities. Vectors are plotted originating at each vertex, denoted by a +, and have a length and direction proportional to the vertex velocity components. If (x1,y1) are the coordinates of vertex (i,j), the coordinates of the vector end point (x2,y2) are given by
x2 = x1 + (uij) (DROU) and y2 = y1 + (vij) (DROU) where DROU is a scaling coefficient defined as DROU = (0.9) (VELmax) (xi-IP1/I)
This coefficient is recalculated whenever a velocity-vector plot is drawn, and it scales the length of a vector drawn for the largest u or v velocity in the system at that instant (VELmax) to be 9/10th the length of the average zone (xi-IP1/I). This method ensures that the vectors are always of reasonable length, regardless of velocity magnitude, The plot is deleted if there are no significant velocities in the entire systeml.
The contour plots are drawn by a routine that creates plots for any cell-centered quantity stored in CQ, and they are composed of connected vector segments joining points of equal quantity, just as the lines on a contour map join points of equal altitude. The plots may be either linear or logarithmic in contour increment. Logarithmic plots are are more useful for atmospheric studies and are provided for the isopycnics and isotherms, whereas the vorticity plot is linear.
The printed information consists primarily of a listing of the principal field variables over the entire grid, One line is printed for each zone, giving the 12 quantities i, j, xij, yij, uij, vij, Iij, ρij, Vij, Dij, Mij, and pij for the zone of its lower left vertex.
A two-line short print is provided every cycle, and it contains the following 13 quantities:
For either δtv or δtc, if the printout indicates that the limiting zone is zone (1,2), the tentative next time step, δtv=δtc= (δt) (δtac), is small enough to satisfy the stability and accuracy requirements at every point in the mesh.
The short print is provided on fanfold paper regardless of the LPR setting, and on microfilm if LPR=1 or LPR=2. LPR primarily controls the destination of the full zone prints, where:
LPR = 1 gives zone prints on microfilm only, LPR = 2 gives zone prints on film and paper, LPR = 3 gives zone prints on paper only.
If IPR=0, no information is printed on microfilm. This case is intended for motion picture use. and the only microfilm output is a particle plot, For movies, the user should hold the δt constant, and set DTO = δt or 2δt. The code is easily altered to provide some plot other than a particle plot for the movie if desired, or to have a frame shared by several different types of plots.
Tape dumps are staged out as T20MD and/or TLIMD, as described in Sec. III C. The quantities dumped are the contents of the SCM common YSC2, the LCM Block YLC1, and, if particles are used, the LCM block YLC2.
A tape restart is performed by staging in the dump tape as Fileset 7. The input deck consists of an IBAR= 0 data card. where JBAR = the dump number on the tape and is used as a check.
Conceptually, the YAQUI code in this report should be able to calculate a truly incompressible flow, defined as a flow in which the sound speed is vastly greater than the fluid speed. Practically, however, the code should be slightly modified to render it suitable for handling such calculations. The equations we use are intended for flows containing compressibility effects, and they, indeed, differ from those we would choose for a fully incopressible flow technique.
In incompressible flow, variations in I can be neglected unless buoyancy effects are important, and as ρ is essentially a constant for each fluid element, the mass equation reduces to the requirement for vanishing velocity divergence, Using an equation of state is therefore unnecessary, because the changes in pressure arise as a direct consequence of the dynamics. In YAQUI, however, the equation of state is inherent: Phase 2 assumes it througn the appearance or c2, and the equation of state is used directly to update pL into the newn+1p, to account for ρ changes that occurred in the Phase-3 convection. Nevertheless, the implicit treatment [1] should still enable YAQUI to handle incompressible flows. In practice, we see this to be true for Mach numbers down to about 0.01. For lower Mach numbers, however, there are three features in YAQUI that introduce difficulties. The first difficulty occurs in the (Lagrangian) iterative phase, where we compute ρL from an equation that leaves δt2 errors. The second arises in the convective flux calculation, which is treated explicitly in Phase 3. This introduces nonzero values into the velocity divergence and, consequently, allow the densities to change. The third arises in the internal energy calculation, in which nonvanishing values introduce fluctuations into the internal energy field. When the overall level of internal energy is high, these fluctuations are reflected in large variations of n+1p, which cannot be efficiently corrected in the subsequent pressure iteration. As a solution to the problem, we have bypassed the equation-of-state calculation after cycle 0, and instead use a "n+1p=pL, in incompressible flows.
Yet another choice would be to iterate Phase 2 to much greater accuracy, which would not be very economical, especially in view of the vastly increased computer time requirements. We cannot run with the limit of ε = 0 in Phase 2, but rather use a value more on the order of, say, 10-5 which leaves relative errors of that order in pL.
These considerations can be illustrated with the stiffened gas equation of state, appearing in the code version of this report. For this, n is given by
np = a2 (nρ − ρ0) + (γ − 1)n ρnI
The incompressible limit can be described by a2→ ∞ forcing the (γ − 1)ρI part of the equation to be negligible, or by I→ ∞. Because true ∞ cannot be used on the computer, we might choose, say, a2 = 1010. Even this less-than-infinity a2 is large enough to magnify any slight ρ erors into appreciable variations in the p field.
To implement the n+1p = PL logic in YAQUI, the storage requirements must be considered. Examination of Fig. 9 reveals that PL, in storage word 11, is not saved in Phase 3. The simplest way to preserve PL, throughout the cycle is to create a 15th word of storage and store PL in it after Phase 2. Then, at the beginning of the next cycle, n+1p can be set from it quite easily, Note that then one must set NQ = 15.
The standard Phase-2 treatment is to bypass the updating of the vertices of any cell whose δp satisfies the convergence test, the argument being that the slightly improved accuracy is generally not worth the extra computer time required to obtain it. When using n+1p = pL, however, it becomes more appropriate to update the vertices of all cells, whether or not the convergence test was satisfied.
The following list provides the names, descriptions, and sources of all quantities in the SCM COMMON YSC2 in the (0,0) overlay. This COMMON is of fundamental importance in communication between the various overlays and their subroutines, It contains all the SCM-based information that must be maintained from cycle to cycle, and it is the SCM portion of the tape-dump data.
The sources in the list are keyed to the following symbols:
I = Supplied as part of the standard input data. The parenthetical symbol that follows I specifies where this quantity is read.
O = (0,0) Main Overlay
L = (0,0) Subroutine LOOP
F = (0,0) Subroutine FILMCO
S = (1,0) Subroutine YASET1
P = (1,0) Subroutine PARTGEN
2 = (2,0) Subroutine YAQUI2
R = (2,0) Subroutine REZONE
Multiple sources indicate that the quantity is re-calculated.
| NAME | DESCRIPTION | SOURCE |
|---|---|---|
| AA | Dummy word, always the first word in COMMON. | -- |
| ANC | aNC, alternate node-coupler coefficient. | S |
| ASQ | a2, the zero-temperature sound speed for the stiffened gas equation of state. | I(S) |
| A0 | α0, determines Phase-3 momentum differencing form. | I(0) |
| A0FAC | α0M/[2(1+α0M2)],used in calculating δtc | S |
| A0M | α0M, the α0 used in Phase-3 M and E calculations | I(0) |
| B0 | >β0, determines Phase-3 differencing form, used with α0 and α0M. | I(0) |
| COLAMU | (1+ε)(λ+2μ), used in Phase-1 viscocity-coefficient calculation. | S |
| CYL | =1, if cylindrical coordinates, =0, if plane coordinates. | I(0) |
| DR | δr, the cell size in the radial direction if uniformly zoned. | I(0) |
| DT | δt, the time set, subject | I(S),2 |
| DTC | δtc, competing δt based on Phase-3 convective flux and divergence considerations. | 2 |
| DTFAC | Initial δtv and δtc each cycle are given by δtv=δtc=(δt)(δtfac). | 2 |
| DTGR | δt*gr | 2 |
| DTGZ | δt*gz | 2 |
| DTGZP | δt*gzp | 2 |
| DTO | Problem time interval between outputs (plots and prints). | I(S) |
| DTOC | Problem time at which to change to next DTO in the set. | I(S) |
| DTO16 | δt/16 | 2 |
| DTO2 | δt/2 | 2 |
| DTO4 | δt/4 | 2 |
| DTO8 | δt/8 | 2 |
| DTPOS | δt possible for the cycle, but actual δt may be reduced to adjust to output time. | 2 |
| DTV | δtv, competing δt based on Phase-1 viscous stress considerations. | 2 |
| DT8 | δt*8. | 2 |
| DZ | δz, the cell size in the axial direction if uniforml zoned. | I(0) |
| EM10 | 10-10, epsilon added to terms to ensure that they do not vanish. | S |
| EPS | ε, convergence criterion for the Phase-2 iteration. | I(S) |
| FIBP | Floating-point equivalent of Ip | P |
| FIPXL | Floating-point frame coordinate for left edge of particle plot. | F |
| FIPXR | Floating-point frame coordinate for right edge of particle plot. | F |
| FIPYB | Floating-point frame coordinate for bottom edge of particle plot. | F |
| FIPYT | Floating-point frame coordinate for top edge of particle plot. | F |
| FIXL | Floating-point frame coordinate for left edge of regular plots. | F |
| FIXR | Floating-point frame coordinate for right edge of regular plots. | F |
| FIYB | Floating-point frame coordinate for bottom edge of regular plots. | F |
| FIYT | Floating-point frame coordinate for top edge of regular plots. | F |
| FJBP | Floating-point equivalent of Jp. | P |
| FREZ | Expansion coefficient for zoning; = 1.0 if uniform throughout. | I(S) |
| GGM1 | γ(γ-1), in which γ is the equation-of-state specific heat ratio if the gas is truly polytropic. | S |
| GM1 | (γ-1). | I(S) |
| GR | gr, gravity component in the r direction, ±. | I(S) |
| GRDVEL | =0. if pure Eulerian, =1. if Lagrangian, =2. if REZONE. | I(0) |
| GZ | gz, gravity component in the z direction, ±. | I(S) |
| GZP | gzp, gz felt by the particles. May be equal to GZ. | I(S) |
| I | Index i. In COMMON because of ENTRY SETIJ in LOOP. | -- |
| IBAR | I, the number of interior fluid zones in the r direction. | I(0) |
| IBP | Ip, the number of particle-grid zones in the r direction. | I(S) |
| IBP1 | Ip+1, index of right-most column of particle-grid vertices. | P |
| ICOLOR | =1 for color movie, =0 for black and white processing. | I(S) |
| IDTO | Index for DTO and DTOC tables. | S,2 |
| IJ | Index for cell (i,j), initialized by LOOP. | L |
| IJM | Index for cell (i,j-1), initialized by LOOP. | L |
| IJP | Index for cell (i,j+1), initialized by LOOP. | L |
| IJPS | Index for cell (i,j+1), saved for later reference to cell (i,j+1). | L |
| IMOME3 | IMOMX+1000. forces resetting of J in Statement No. 2020 in PARTMOV if IMOMX=1. | P |
| IMOMX | =1 if particle-fluid momentum exchange, = 0 otherwise. | I(S) |
| IM1 | I=1,index of next-to-last zone or vertex in column. | S |
| IM6 | I-6, in usual large grids, this column is in somewhat from the right. | S |
| IPAR | LOCF(AA2), the address of LCM block AA2, for tape dump. | P |
| IPXL | Integer frame coordinate for left edge of particle plot. | F |
| IPXR | Integer frame coordinate for right edge of particle plot. | F |
| IPYB | Integer frame coordinate for bottom edge of particle plot. | F |
| IPYT | Integer frame coordinate for top edge of particle plot. | F |
| IP1 | I+1, index of rightmost column of grid vertices. | S |
| IP2 | I+2, index used in reversed DO loops. | S |
| ISCF1 | ISC2-NQ, the relative first word address (f.w.a) of i=I+1 zone in SCM buffer row 1/3. | S |
| ISCF2 | ISCF1+NQI, the relative f.w.a of i=I+1 zone in SCM buffer row 2/3. | S |
| ISC2 | NQI+1, the relative f.w.a of i=1 zone in SCM buffer row 2/3. | S |
| ISC3 | ISC2+NQI, the relative f.w.a of i=1 zone in SCM buffer row 3/3. | S |
| ITV | JP1*NQI, the relative f.w.a of the J+2 row in LCM storage. | S |
| IUNF | IUNF, the number of zones with uniform initial δr (DR). | I(0) |
| IXL | Integer frame coordinate for left edge of regular plots. | F |
| IXR | Integer frame coordinate for right edge of regular plots. | F |
| IYB | Integer frame coordinate for bottom edge of regular plots. | F |
| IYT | Integer frame coordinate for top edge of regular plots. | F |
| J | Index j. In COMMON because of ENTRY R1ROW and W1ROW in LOOP. | -- |
| JBAR | J, the number of interior fluid zones in the z-direction. | I(0) |
| JBP | Jp, the number of particle-grid zones in the z direction. | I(S) |
| JBP2 | Jp+2, index of the topmost row of particle-grid vertices. | P |
| JCEN | Number of zones up to center of uniform-grid region. | I(0) |
| JM10 | J-10. In usual large grids, this row is down from the top. | S |
| JM14 | J-14. In usual large grids, this row is down from the top. | S |
| JP1 | J+1, index of topmost row of inerior zones. | S |
| JP2 | J+2, index of topmost row of grid vertices. | S |
| JP4 | J+4, index used in reversed DO loops and in LCM clearing. | S |
| JP4O2 | (J+4)/2, j index at midpoint of full YAQUI grid. | S |
| JUNF | JUNF, the number of zones with uniform initial δz (DZ). | I(0) |
| JUNFO2 | JUNF/2 uniform zones lie above JCEN, and JUNF/2 lie below. | S |
| KXI | ε, the ρ exponent that determines the viscosity form. | I(0) |
| LAM | λinput viscosity coefficient. A real number. | I(S) |
| LJP2 | First word address of last zone in row JP2 when in usual SCM buffer row. | S |
| LPB | Number of words, truncated to a multiple of 3, that will fit in NQI-wd. | P |
| LPR | Determines output options of film and printer. | I(S) |
| MU | μinput viscosity coefficient. A real number. | I(S) |
| NAME | The problem identification from input card No. 2, up to 79 characters. | I(S) |
| NCYC | Number of calculational cycles completed. | S,2 |
| NLC | Number of words of LCM block AA1 actually in use, for tape dump. | S |
| NPS | Number of words of LCM particle block AA2 actually in use, for tape dump. | P |
| NPT | Total number of particles generated | P |
| NQ | Number of quantities, or storage words, per cell. | 0 |
| NQI | NQ*IP1, the number of words for a full row of zones. | S |
| NQIB | NQ*IBAR, the number of words back to zone i=1 when at i=I+1 in SCM. | S |
| NQI2 | NQI*2, the number of words in two full rows of zones, for PARTMOV. | P |
| NSC | Number of words in this SCM COMMON, for tape dump. | S,2 |
| NUMIT | Number of iterations required for PHASE-2 convergence. | 2 |
| NUMTD | Number of the next tape dump. | S,2 |
| OM | ω, the Phase-2 iteration relaxation factor. | I(S) |
| OMANC | (1-aNC, used in δtv calculation. | S |
| OMCYL | (1-CYL), used in calculating r from x. | S |
| OMEM10 | (1-10-10). | S |
| OPEM10 | (1+10-10). | S |
| PDR | The uniform Δx of the particle grid. | S,P,R |
| PDZ | The uniform Δy of the particle grid. | S,P,R |
| PXCONV | Frame-conversion coefficient for particle-plot x direction | F |
| PXL | X coordinate of left edge of particle grid, in problem units. | F |
| PXR | X coordinate of left edge of particle grid, in problem units. | F |
| PXRP | PXR*OPEM10, test comparand in particle-grid stepping. | F |
| PYB | Y coordinate of bottom edge of particle grid in problem units. | S,F,R |
| PYBM | PYB*OPEM10, test comparand in particle-grid mapping. | F |
| PYCONV | Frame-conversion coefficient for particle-plot y direction. | F |
| PYT | Y coordinate of top edge of particle grid in problem units. | F |
| PYTP | PYT*OMEM10, test comparand in particle grid mapping. | F |
| RDT | 1/δ. | 2 |
| REZRON | ρ0 of the ambient atmosphere at altitude REZYO at t=t0. | I(S) |
| REZSIE | The (constant) specific internal energy of the ambient atmosphere. | I(S) |
| REZUE | Grid-expansion u velocity, available for REZONE use. | I(S),R |
| REZVE | Grid-expansionv velocity, available for REZONE use. | I(S),R |
| REZVT | Grid-translation velocity, available for REZONE use. | I(S),R |
| REZY0 | Y coordinate of center of expansion, refers to YAQUI vertex (1,JCEN+2). | I(S) |
| RIBAR | Reciprocal of I. | S |
| RIBJB | Reciprocal of I*J used in control region grind calculation. | S |
| RIBP | Reciprocal of Ip. | P |
| RJBP | Reciprocal of Jp. | P |
| ROMFR | Reciprocal of (1.-FREZ). | S |
| RON | ρ0, normal density for equation-of-state use. | I(S) |
| RPDR | Reciprocal of Δx. | F |
| RPDRDZ | Reciprocal of (Δx*Δy). | F |
| RPDZ | Reciprocal of Δy. | F |
| T | t, the problem time. | I(S),2 |
| THIRD | 1./3. | S |
| TLIMD | = 1.0 to force a tape dump and RETURN before time limit. | I(S) |
| TOUT | The next output time for plots/prints. | S,2 |
| TWFIN | Time-When-To-Finish: calculation completed when t≥TWFIN. | I(S) |
| T20MD | =1.0 to force tape dumps every 20' CP time. | I(S) |
| VV | Velocity-vector plot-scaling coefficient, = 9/10 average δr. | F |
| XCONV | Frame-conversion coefficient for regular plots, x direction. | F |
| XL | X coordinate of leftmost vertex of the grid, in problem units. | F |
| XR | X coordinate of rightmost vertex of the grid, in problem units. | F |
| YB | Y coordinate of bottommost vertex of the grid, in problem units. | F |
| YCONV | Frame-conversion coefficient for regular plots, y direction. | F |
| YT | Y coordinate of topmost vertex of the grid, in problem units. | F |
| ZZ | Dummy word, always the final word in the COMMON | -- |
Here we present results from several YAQUI calculations. Emphasis is on the method's versatility in handling a given problem, rather than on presenting a wide variety of different examples [2].
The flexibility of the Arbitrary Lagrangian-Eulerian approach is illustrated in the calculation of a one-dimensional shock tube, performed first in a Lagrangian fashion, and then with a full Eulerian rezone. This example is followed by sequences at very early times from three calculations of a low-altitude explosion, first Lagrangian, next Eulerian, then with the REZONE subroutine as presented in this report.
The versatility of the YAQUI particle technique is illustrated at one extreme by the marker particles carried along with the fluid in the low-altitude explosion calculations, where the particles have no influence on the flow, and at the other extreme by calculations in which the particles govern the fluid dynamics through the momentum-exchange feature.
Finally, we present listed results from a particle-fluid momentum-exchange calculation, for those readers who may find a benchmark calculation useful.
Detailed discussions of various YAQUI calculations will be presented elsewhere, and no attempt is made here to describe a variety of late-time results .
The two examples in Fig. l6 and Fig. 17 are selected from a series of one-dimensional shock-tube test cases, although they do not necessarily represent the best that YAQUI can do for this problem, they clearly demonstrate that satisfactory results can be obtained in both the Lagrangian and Eulerian limits. The figures show the profiles (heavy lines) of velocity, pressure, specific internal energy, and density from a pure Lagrangian (GRDVEL = 1.0) calculation and then a pure Eulerian (GRDVBL = 0.0) calculation of a 2:1 density-ratio shock tube, along with the theoretical solution (lighter lines) to the problem [3].
The calculations were performed in a plane mesh 60 cells long by 1 cell high, allowing 30 cells for each fluid region. The initial ρ was 0.2 on the left and 0.1 on the right, and the initial specific internal energy was 0.18. The Initial cell size was δr = δz = 1/3, the viscosity coefficients were λ=0.002 and μ=0.0, and, in addition, the gas was polytropic with γ=5/3. The Eulerian shock tube was run with full donor-cell differencing (α0=α0M=1.0, α0=0.0). At t=0, the diaphram separating the two fluid regions was instantaneously removed, causing a shock to advance into the lower density region, and a rarefaction to propagate back from the contact surface into the higher density region. In both calculations, δt was held constant at 0.1, and the profiles shown in Fig. l6 and Fig. 17 are at t=10.0, Such calculations typically require 20 to 30 sec of CDC-7600 time to run to t=15.0, producing plots and prints every unit of time.
These examples demonstrate three distinct approaches to the treatment of grid motion in a typical low-altitude explosion calculation. The sets of six plots in each of Fig. 18, Fig. l9, Fig. 21, Fig. 23, and Fig. 24 represent the marker particles, computing mesh, and velocity vectors (top) and isopycnic, isotherm, and vorticity contour plots (bottom).
Figure 18 shows the various plots at time t=0, immediately after superposing the explosion density, energy, and velocity data, which were provided by a one-dimensional spherical code, onto a unifom 26 by 52 cell YAQUI computing grid that already contained an appropriate ambient background. This procedure was described in Sec. III A. In the particle plot, the explosion debris is represented by a hemisphere of particles, surrounded by more widely spaced particles in an adjacent region of the ambient atmosphere. These marker particles do not enter directly into the calculation, but are used solely as an aid to flow visualization. Note that the velocity, density, and energy fields are well developed, but that the vorticity field is not, and indeed, will not be well established for about the first two seconds of problem time.
Figure l9 shows the same six plots at t=1 sec from a calculation of this problem in which the interior vertices were allowed to move in a purely Lagrangian fashion (GRDVEL=1.0). The rigid walls of the computing mesh are held fixed, causing the strong radially expanding shock to reflect back into the central region shortly after t=2 sec. This effect is visible in Fig. 20, which shows the appearance of tbe velocity vectors at t=2 and t=3 sec.
Figure 21, shows the six basic plots from a pure Eulerian calculation (GRDVEL=0.0) of the same problem at t=1 sec. More resolution is available in the central region, and less resolution is given to the shock front, than in the Lagrangian calculation, As in the Lagrangian calculation, the walls are rigid, and the t=2 and t=3 sec velocity-vector plots of Fig. 22 show the same strong wall reflection as did Fig. 20.
In reality, the edges of an atmospheric region are not rigid walls, so to calculate such an atmospheric explosion beyond the first two or so seconds, with this degree of resolution, would require one of several possible alternatives:
(1) A vastly larger computing mesh could be used, but this obviously is not economical in terms of computer storage aad time requirements.
(2) Continuative outflow boundaries would allow the strong radially expanding shock to leave the system with a minimum of upstream disturbance, but the subsequent rise of the explosion debris sucks material up behind it in the central column, causing the bottom and right walls of the mesh to become inflow boundaries. Appropriate inflow conditions are difficult to define, suggesting again a larger computing mesh to avoid this difficulty.
(3) A third choice, which we exploit in YAQUI, is to allow the entire mesh to expand at a rate that will keep the reflective boundaries out ahead of the radial shock while it has significant strength. At the same time, we vary the sizes of the interior zones to provide high resolution in the central region and much coarser resolution in the outlying regions, which still allows the use of the same number of cells.
Figure 23 and Figure 24 show such a calculation, using the REZONE subroutine exactly as provided in tba code version of this report. The problem input is identical to the preceding cases except that GRDVEL=2.0. As the problem proceeds, the mesh is continuously enlarged at a rate that depends upon the magnitude of the velocities approaching the boundaries. This expansion leaves a region without particles around the outer regions of the mesh, which is already evident by t=1 sec (Fig. 23). By t=5 sec (Fig. 24), the initial mesh radius has already increased by 50%, allowing the calculation to run to much later times without boundary interference than do either the Lagrangian or pure Eulerian approaches. Note in Fig. 24 that the velocities near the rigid walls are negligible, and that the vorticity field has become well established,
Because the computer is programmed to draw pictures of a fixed size, the frame scales of Fig. 23 and Fig. 24 differ and are further quite different from the scale in the preceding figures. (Information printed on the film below each plot provides the necessary specifications to properly interpret the plot.)
Figure 23 and Figure 24 represent only the early stages of a calculation that has been made feasible through the use of continuous rezoning and mesh expansion. These techniques, combined with an appropriate mesh translation that follow the debris rise, allow the dynamics to be followed for several hundred seconds of problem time. A wide variety of REZONE subroutines have been used with success, each tailored to provide optimum reaulta for a particular problem.
We generally enhance this approach by combining it with an initial grid containing variable cell sizes, aa described in Sec. III A. Figure 25 shows a setup configuration for the same problem in which the cells are expanded (FREZ=1.1) beyond uniformly zoned 16 by 32 cell central region. This affords high resolution where required, at the same time allowing the continuous rezoning and expansion to take place much more gradually, as in this particular case the initial mesh encompasses a much larger volume.
The CDC-7600 CP time per cell per cycle (grinds) averages approximately 0.50 msec at two iterations per cycle, increasing by about 0.03 to 0.04 msec for each additional iteration required for convergence in Phase 2. Calculations such as this, with 1500- to 1600-cell meshes can be followed to over 200 sec of problem time in well under 1hr of CDC-7600 time, with generous amounts of output along the way.
An example of the particle-fluid momentum-exchange feature described in Sec. II E 2 is illustrated in the particle-drag problem of Fig. 25. The first set of three frames show the initial particles, velocity vectors, and the (Eulerian) computing grid with cylindrical symmetry and rigid free-slip boundaries. In this calculation, a sphere of particles, each of which has a finite mass and drag coefficient, is immersed in a fluid of uniform density and energy, representing a two-fluid configuration in which the density of the heavier spherical part is given by the sum of the background fluid and particle densities. Initially, there are no velocities in the system; the entire dynamics of the calculation result from a gravitational force upon the particles but not upon the fluid. This causes the sphere of particles to fall and deform, producing a pronounced circulation pattern within the fluid.
The evolution of this process is shown in the remaining seven sets of plots in Fig. 26, at times of 9, 12, 15, 18, 21, 24, and 27. Each set of three frames consists of a particle plot and the velocity vectors and vorticity contours for the fluid. Note that the effects of drag soon retard the leading edge of the sphere relative to the shielded trailing edge. The sphere is deformed into a cup, with a vortex ring around the rim. At time t=21, the cup collides with the bottom wall of the mesh and is seen gradually settling into place thereafter. By time t=27, only the rolled rim retains any definition, but it, too, will soon collapse into the rest of the particles. The circulation pattern will persist for some time, until viscous effects gradually damp it out.
The following pages are abstracted from the microfilm output of a particle-fluid momentum-exchange test calculation, They are included as an aid to the reader who uses YAQUI. allowing him to set up the same problem and compare results.
The input data are listed in their entirety, and include all information necessary to specify the problem. Subsequent pages show the initial particle and zone plot configurations at time 0, along with a sample frame abstracted from the cell print. This includes a row across the mesh halfway up, cutting through the initial position of the sphere of particles. This same frame of print output is included for cycle 1, to show the initial changes in the fluid variables.
For t=1.0 (cycle 7). we present six frames. These include plots of the particles, and for the fluid, the velocity vectors and contours of density, specific internal energy, and vorticity, followed by the sample listing. The normal velocities on the symmetry axis are, of course, nonphysical. They result from the momentum carried by the particles and distributed to the cell vertices. After each cycle. these velocities are reset to zero, so no buildup can occur. This will, however, act as a sink for momentum, which would be easy to correct if it became a problem.
Finally. we present the same six frames at t=9.0 (cycle 232). Note in the listing that the circulation pattern is quite evident in the wake of the particles. The u velocities at this height on the axis are now zero, as the particles are no longer present here to contribute momentum changes.
The CDC-7600 CP time for this calculation was 305 sec for 265 calculational cycles (to time t=9.54181). After the first 100 cycles, the number of iterations required for convergence in each cycle stabilized at 4, for which the grinds (CP time per cell per cycle) averaged about 0.637 msec. Comparison with the grinds for the low-altitude explosion calculation (0.56 msec) indicates that slightly more time is required for the momentum-exchange option.
The following listings have most of the values as Fortran E format with a fraction followed by the power of 10. Where appropriate these have been shortened to the standard F format with an appropriate number of decimal places.
PARTICLE-FLUID MOMENTUM EXCHANGE TEST
IBAR = 26
JBAR = 52
IUNF = 26
JUNF = 52
JCEN = 26
DR = 1.0
DZ = 1.0
CYL = 1.0
GRDVEL = 0.0
A0 = 0.75
A0M = 0.75
B0 = 0.0
KXI = -1
MU = 1.0
LAM = 10.0
OM = 1.0
EPS = 0.0001
GR = 0.0
GZ = 0.0
ASQ = 100.0
RON = 1.0
GM1 = 0.0
FREZ = 1.0
YB = 0.0
REZY0 = 0.0
REZUE = 0.0
REZVE = 0.0
REZVT = 0.0
REZRON = 0.0
REZSIE = 0.0
IBP = 26
JBP = 52
PDR = 1.0
PDZ = 1.0
PYB = 0.0
GZP = 1.0
IMOMX = 1
T = 0.0
DT = 0.1
T20MD = 0.0
TLIMD = 0.0
TWFIN = 10.0
LPR = 1
ICOLOR = 0
DTO(1-10) = 1.0 -0.0 -0.0 -0.0 -0.0
-0.0 -0.0 -0.0 -0.0 -0.0
DTOC(1-10) = 100.0 -0.0 -0.0 -0.0 -0.0
-0.0 -0.0 -0.0 -0.0 -0.0
DRPAR = 0.5
DZPAR = 0.5
XC = 0.0
YC = 26.0
XD = 0.0
YD = 0.0
UPAR = 0.0
VPAR = 0.0
MTE = 0.25
DRAG = 1.0
406 PARTICLES GENERATED, WITH TOTAL MASS = 346.375
NB = 0
NR = 26
NT = 52
NL = 0
UI = 0.0
VI = 0.0
RO1 = 1.0
SIEI = 1.0
PDR=1.0 PDZ=1.0 PXR=26.0 PYB=0.0 PYT=52.0 YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=0. CYCLE = 0
ZONES DRMIN=1.0 DRMAX=1.0 DZMIN=1.0 DZMAX=1.0 XR=26.0 YB=0. T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=0.0 CYCLE = 0
IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=0. CYCLE = 0 I J X Y U V SIE RHO VOL D M P 8 26 7.0 24.0 0.0 0.0 1.0 1.0 7.5 0.0 7.0 0.0 9 26 8.0 24.0 0.0 0.0 1.0 1.0 8.5 0.0 8.0 0.0 10 26 9.0 24.0 0.0 0.0 1.0 1.0 9.5 0.0 9.0 0.0 11 26 10.0 24.0 0.0 0.0 1.0 1.0 10.5 0.0 10.0 0.0 12 26 11.0 24.0 0.0 0.0 1.0 1.0 11.5 0.0 11.0 0.0 13 26 12.0 24.0 0.0 0.0 1.0 1.0 12.5 0.0 12.0 0.0 14 26 13.0 24.0 0.0 0.0 1.0 1.0 13.5 0.0 13.0 0.0 15 26 14.0 24.0 0.0 0.0 1.0 1.0 14.5 0.0 14.0 0.0 16 26 15.0 24.0 0.0 0.0 1.0 1.0 16.5 0.0 15.0 0.0 17 26 16.0 24.0 0.0 0.0 1.0 1.0 16.5 0.0 16.0 0.0 18 26 17.0 24.0 0.0 0.0 1.0 1.0 17.5 0.0 17.0 0.0 19 26 18.0 24.0 0.0 0.0 1.0 1.0 18.5 0.0 18.0 0.0 20 26 19.0 24.0 0.0 0.0 1.0 1.0 19.5 0.0 19.0 0.0 21 26 20.0 24.0 0.0 0.0 1.0 1.0 20.5 0.0 20.0 0.0 22 26 21.0 24.0 0.0 0.0 1.0 1.0 21.5 0.0 21.0 0.0 23 26 22.0 24.0 0.0 0.0 1.0 1.0 22.5 0.0 22.0 0.0 24 26 23.0 24.0 0.0 0.0 1.0 1.0 23.5 0.0 23.0 0.0 25 26 24.0 24.0 0.0 0.0 1.0 1.0 24.5 0.0 24.0 0.0 26 26 25.0 24.0 0.0 0.0 1.0 1.0 25.5 0.0 25.0 0.0 27 26 26.0 24.0 0.0 0.0 0.0 0.0 0.0 0.0 12.75 0.0 1 27 0.0 25.0 0.0 0.0 1.0 1.0 0.5 0.0 0.25 0.0 2 27 1.0 25.0 0.0 0.0 1.0 1.0 1.5 0.0 1.0 0.0 3 27 2.0 25.0 0.0 0.0 1.0 1.0 2.5 0.0 2.0 0.0 4 27 3.0 25.0 0.0 0.0 1.0 1.0 3.5 0.0 3.0 0.0 5 27 4.0 25.0 0.0 0.0 1.0 1.0 4.5 0.0 4.0 0.0 6 27 5.0 25.0 0.0 0.0 1.0 1.0 5.5 0.0 5.0 0.0 7 27 6.0 25.0 0.0 0.0 1.0 1.0 6.5 0.0 6.0 0.0 8 27 7.0 25.0 0.0 0.0 1.0 1.0 7.5 0.0 7.0 0.0 9 27 8.0 25.0 0.0 0.0 1.0 1.0 8.5 0.0 8.0 0.0 10 27 9.0 25.0 0.0 0.0 1.0 1.0 9.5 0.0 9.0 0.0 11 27 10.0 25.0 0.0 0.0 1.0 1.0 10.5 0.0 10.0 0.0 12 27 11.0 25.0 0.0 0.0 1.0 1.0 11.5 0.0 11.0 0.0 13 27 12.0 25.0 0.0 0.0 1.0 1.0 12.5 0.0 12.0 0.0 14 27 13.0 25.0 0.0 0.0 1.0 1.0 13.5 0.0 13.0 0.0 15 27 14.0 25.0 0.0 0.0 1.0 1.0 14.5 0.0 14.0 0.0 16 27 15.0 25.0 0.0 0.0 1.0 1.0 16.5 0.0 15.0 0.0 17 27 16.0 25.0 0.0 0.0 1.0 1.0 16.5 0.0 16.0 0.0 18 27 17.0 25.0 0.0 0.0 1.0 1.0 17.5 0.0 17.0 0.0 19 27 18.0 25.0 0.0 0.0 1.0 1.0 18.5 0.0 18.0 0.0 20 27 19.0 25.0 0.0 0.0 1.0 1.0 19.5 0.0 19.0 0.0 21 27 20.0 25.0 0.0 0.0 1.0 1.0 20.5 0.0 20.0 0.0 22 27 21.0 25.0 0.0 0.0 1.0 1.0 21.5 0.0 21.0 0.0 23 27 22.0 25.0 0.0 0.0 1.0 1.0 22.5 0.0 22.0 0.0 24 27 23.0 25.0 0.0 0.0 1.0 1.0 23.5 0.0 23.0 0.0 25 27 24.0 25.0 0.0 0.0 1.0 1.0 24.5 0.0 24.0 0.0 26 27 25.0 25.0 0.0 0.0 1.0 1.0 25.5 0.0 25.0 0.0 27 27 26.0 25.0 0.0 0.0 0.0 0.0 0.0 0.0 12.75 0.0 1 28 0.0 26.0 0.0 -0.000074 1.0 1.0 0.5 0.0 0.25 0.0 2 28 1.0 26.0 0.0 -0.000099 1.0 1.0 1.5 0.0 1.0 0.0 3 28 2.0 26.0 0.0 -0.000099 1.0 1.0 2.5 0.0 2.0 0.0 4 28 3.0 26.0 0.0 -0.000099 1.0 1.0 3.5 0.0 3.0 0.0 5 28 4.0 26.0 0.0 -0.000099 1.0 1.0 4.5 0.0 4.0 0.0 6 28 5.0 26.0 0.0 -0.000099 1.0 1.0 5.5 0.0 5.0 0.0 7 28 6.0 26.0 0.0 -0.000099 1.0 1.0 6.5 0.0 6.0 0.0 8 28 7.0 26.0 0.0 -0.000099 1.0 1.0 7.5 0.0 7.0 0.0 9 28 8.0 26.0 0.0 -0.000047 1.0 1.0 8.5 0.0 8.0 0.0 10 28 9.0 26.0 0.0 0.0 1.0 1.0 9.5 0.0 9.0 0.0 11 28 10.0 26.0 0.0 0.0 1.0 1.0 10.5 0.0 10.0 0.0 12 28 11.0 26.0 0.0 0.0 1.0 1.0 11.5 0.0 11.0 0.0 13 28 12.0 26.0 0.0 0.0 1.0 1.0 12.5 0.0 12.0 0.0 14 28 13.0 26.0 0.0 0.0 1.0 1.0 13.5 0.0 13.0 0.0 15 28 14.0 26.0 0.0 0.0 1.0 1.0 14.5 0.0 14.0 0.0
IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=1. CYCLE = 1 I J X Y U V SIE RHO VOL D M P 8 26 7.0 24.0 0.0 -0.000092 1.0 1.0 7.5 -0.000012 7.0 0.0 9 26 8.0 24.0 0.0 -0.000029 1.0 1.0 8.5 -8.992 8.0 0.0 10 26 9.0 24.0 0.0 0.0 1.0 1.0 9.5 0.0 9.0 0.0 11 26 10.0 24.0 0.0 0.0 1.0 1.0 10.5 0.0 10.0 0.0 12 26 11.0 24.0 0.0 0.0 1.0 1.0 11.5 0.0 11.0 0.0 13 26 12.0 24.0 0.0 0.0 1.0 1.0 12.5 0.0 12.0 0.0 14 26 13.0 24.0 0.0 0.0 1.0 1.0 13.5 0.0 13.0 0.0 15 26 14.0 24.0 0.0 0.0 1.0 1.0 14.5 0.0 14.0 0.0 16 26 15.0 24.0 0.0 0.0 1.0 1.0 16.5 0.0 15.0 0.0 17 26 16.0 24.0 0.0 0.0 1.0 1.0 16.5 0.0 16.0 0.0 18 26 17.0 24.0 0.0 0.0 1.0 1.0 17.5 0.0 17.0 0.0 19 26 18.0 24.0 0.0 0.0 1.0 1.0 18.5 0.0 18.0 0.0 20 26 19.0 24.0 0.0 0.0 1.0 1.0 19.5 0.0 19.0 0.0 21 26 20.0 24.0 0.0 0.0 1.0 1.0 20.5 0.0 20.0 0.00000001 22 26 21.0 24.0 0.0 0.0 1.0 1.0 21.5 0.0 21.0 0.0 23 26 22.0 24.0 0.0 0.0 1.0 1.0 22.5 0.0 22.0 0.0 24 26 23.0 24.0 0.0 0.0 1.0 1.0 23.5 0.0 23.0 0.00000001 25 26 24.0 24.0 0.0 0.0 1.0 1.0 24.5 0.0 24.0 0.00000001 26 26 25.0 24.0 0.0 0.0 1.0 1.0 25.5 0.0 25.0 0.0 27 26 26.0 24.0 0.0 0.0 0.0 0.0 0.0 0.0 12.75 0.078 1 27 0.0 25.0 0.0 -0.000074 1.0 1.0 0.5 0.0 0.25 0.0 2 27 1.0 25.0 0.0 -0.000099 1.0 1.0 1.5 0.0 1.0 0.0 3 27 2.0 25.0 0.0 -0.000099 1.0 1.0 2.5 0.0 2.0 0.0 4 27 3.0 25.0 0.0 -0.000099 1.0 1.0 3.5 0.0 3.0 0.0 5 27 4.0 25.0 0.0 -0.000099 1.0 1.0 4.5 0.0 4.0 0.0 6 27 5.0 25.0 0.0 -0.000099 1.0 1.0 5.5 0.0 5.0 0.0 7 27 6.0 25.0 0.0 -0.000099 1.0 1.0 6.5 0.0 6.0 0.0 8 27 7.0 25.0 0.0 -0.000099 1.0 1.0 7.5 0.0 7.0 0.0 9 27 8.0 25.0 0.0 -0.000047 1.0 1.0 8.5 0.0 8.0 0.0 10 27 9.0 25.0 0.0 0.0 1.0 1.0 9.5 0.0 9.0 0.0 11 27 10.0 25.0 0.0 0.0 1.0 1.0 10.5 0.0 10.0 0.0 12 27 11.0 25.0 0.0 0.0 1.0 1.0 11.5 0.0 11.0 0.0 13 27 12.0 25.0 0.0 0.0 1.0 1.0 12.5 0.0 12.0 0.0 14 27 13.0 25.0 0.0 0.0 1.0 1.0 13.5 0.0 13.0 0.0 15 27 14.0 25.0 0.0 0.0 1.0 1.0 14.5 0.0 14.0 0.0 16 27 15.0 25.0 0.0 0.0 1.0 1.0 16.5 0.0 15.0 0.0 17 27 16.0 25.0 0.0 0.0 1.0 1.0 16.5 0.0 16.0 0.0 18 27 17.0 25.0 0.0 0.0 1.0 1.0 17.5 0.0 17.0 0.0 19 27 18.0 25.0 0.0 0.0 1.0 1.0 18.5 0.0 18.0 0.0 20 27 19.0 25.0 0.0 0.0 1.0 1.0 19.5 0.0 19.0 0.0 21 27 20.0 25.0 0.0 0.0 1.0 1.0 20.5 0.0 20.0 -0.000000001 22 27 21.0 25.0 0.0 0.0 1.0 1.0 21.5 0.0 21.0 0.0 23 27 22.0 25.0 0.0 0.0 1.0 1.0 22.5 0.0 22.0 0.0 24 27 23.0 25.0 0.0 0.0 1.0 1.0 23.5 0.0 23.0 -0.000000001 25 27 24.0 25.0 0.0 0.0 1.0 1.0 24.5 0.0 24.0 -0.000000001 26 27 25.0 25.0 0.0 0.0 1.0 1.0 25.5 0.0 25.0 0.0 27 27 26.0 25.0 0.0 0.0 0.0 0.0 0.0 0.0 12.75 0.078 1 28 0.0 26.0 0.0 -0.000074 1.0 1.0 0.5 0.0 0.25 0.0 2 28 1.0 26.0 0.0 -0.000099 1.0 1.0 1.5 0.0 1.0 0.0 3 28 2.0 26.0 0.0 -0.000099 1.0 1.0 2.5 0.0 2.0 0.0 4 28 3.0 26.0 0.0 -0.000099 1.0 1.0 3.5 0.0 3.0 0.0 5 28 4.0 26.0 0.0 -0.000099 1.0 1.0 4.5 0.0 4.0 0.0 6 28 5.0 26.0 0.0 -0.000099 1.0 1.0 5.5 0.0 5.0 0.0 7 28 6.0 26.0 0.0 -0.000099 1.0 1.0 6.5 0.0 6.0 0.0 8 28 7.0 26.0 0.0 -0.000099 1.0 1.0 7.5 0.0 7.0 0.0 9 28 8.0 26.0 0.0 -0.000047 1.0 1.0 8.5 0.0 8.0 0.0 10 28 9.0 26.0 0.0 0.0 1.0 1.0 9.5 0.0 9.0 0.0 11 28 10.0 26.0 0.0 0.0 1.0 1.0 10.5 0.0 10.0 0.0 12 28 11.0 26.0 0.0 0.0 1.0 1.0 11.5 0.0 11.0 0.0 13 28 12.0 26.0 0.0 0.0 1.0 1.0 12.5 0.0 12.0 0.0 14 28 13.0 26.0 0.0 0.0 1.0 1.0 13.5 0.0 13.0 0.0 15 28 14.0 26.0 0.0 0.0 1.0 1.0 14.5 0.0 14.0 0.0
PARTICLES PDR=1.0 PDZ=1.0 PXR=26.0 PYB=0.0 PYT=52.0 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=1.0 CYCLE = 7
VELOCITY VECTORS VMAX=0.2964421 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=1.0 CYCLE = 7
ISOPYCNICS MIN=0.990697 MAX=1.00900 L=99.0697 H=1.00807 DQ= 0.00193073 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=1.0 CYCLE = 7
ISOTHERMS MIN=0.972363 MAX=1.00004 L=97.2363 H=0.998169 DQ= 0.00286729 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=1.0 CYCLE = 7
VORTICITY MIN=-0.0734092 MAX=0.222235 L=-0.0734092 H=0.193570 DQ= 0.0236644 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=1.0 CYCLE = 7
T3AA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=1. CYCLE = 7 I J X Y U V SIE RHO VOL D M P 8 26 7.0 24.0 -0.0026 -0.2758 0.985 1.002 7.5 -0.0074 7.015 0.162 9 26 8.0 24.0 0.1847 -0.0915 0.998 1.001 8.5 -0.0054 8.016 0.134 10 26 9.0 24.0 0.2964 0.0234 0.999 1.0 9.5 -0.0017 9.013 0.089 11 26 10.0 24.0 0.1708 0.0150 0.999 1.0 10.5 -0.0012 10.000 0.057 12 26 11.0 24.0 0.1847 0.0098 0.999 1.0 11.5 -0.0008 11.000 0.037 13 26 12.0 24.0 0.0064 0.0062 0.999 1.0 12.5 -0.0006 12.000 0.025 14 26 13.0 24.0 0.0041 0.0039 0.999 1.0 13.5 -0.0004 13.000 0.017 15 26 14.0 24.0 0.0026 0.0025 1.000 1.0 14.5 -0.0003 14.000 0.012 16 26 15.0 24.0 0.0016 0.0016 1.000 1.0 16.5 -0.0002 15.000 0.008 17 26 16.0 24.0 0.0011 0.0010 1.000 1.0 16.5 -0.0002 16.000 0.006 18 26 17.0 24.0 0.0007 0.0006 1.000 1.0 17.5 -0.0001 17.000 0.004 19 26 18.0 24.0 0.0004 0.0003 1.000 1.0 18.5 -0.0000 18.000 0.003 20 26 19.0 24.0 0.0003 0.0002 1.000 1.0 19.5 -0.0000 19.000 0.002 21 26 20.0 24.0 0.0002 0.0001 1.000 1.0 20.5 -0.0000 20.000 0.001 22 26 21.0 24.0 0.0001 0.0001 1.000 1.0 21.5 -0.0000 21.000 0.001 23 26 22.0 24.0 0.0001 0.0000 1.000 1.0 22.5 -0.0000 22.000 0.001 24 26 23.0 24.0 0.0000 0.0000 1.000 1.0 23.5 -0.0000 23.000 0.001 25 26 24.0 24.0 0.0000 0.0000 1.000 1.0 24.5 -0.0000 24.000 0.001 26 26 25.0 24.0 0.0000 0.0000 1.000 1.0 25.5 -0.0000 25.000 0.000 27 26 26.0 24.0 0.0000 0.0000 0.000 0.0 0.0 0.0000 12.750 0.078 1 27 0.0 25.0 0.0000 -0.225 0.978 1.0 0.5 0.0007 0.258 0.036 2 27 1.0 25.0 -0.0003 -0.295 0.972 1.0 1.5 -0.0007 1.008 0.035 3 27 2.0 25.0 -0.0006 -0.295 0.972 1.0 2.5 -0.0007 2.002 0.035 4 27 3.0 25.0 -0.0008 -0.295 0.972 1.0 3.5 -0.0007 3.003 0.036 5 27 4.0 25.0 -0.0012 -0.293 0.972 1.0 4.5 -0.0007 4.003 0.037 6 27 5.0 25.0 -0.0013 -0.292 0.973 1.0 5.5 -0.0008 5.005 0.037 7 27 6.0 25.0 -0.0027 -0.290 0.973 1.0 6.5 -0.0009 6.006 0.041 8 27 7.0 25.0 -0.0006 -0.287 0.983 1.0 7.5 -0.0002 7.007 0.036 9 27 8.0 25.0 0.0078 -0.127 0.997 1.0 8.5 -0.0001 8.009 0.033 10 27 9.0 25.0 0.0128 0.028 0.999 1.0 9.5 -0.0000 9.006 0.025 11 27 10.0 25.0 0.0094 0.0192 0.999 1.0 10.5 -0.0003 10.000 0.015 12 27 11.0 25.0 0.0053 0.0119 0.999 1.0 11.5 -0.0003 11.003 0.011 13 27 12.0 25.0 0.0034 0.0013 0.999 1.0 12.5 -0.0002 12.002 0.007 14 27 13.0 25.0 0.0021 0.0047 1.000 1.0 13.5 -0.0001 13.000 0.005 15 27 14.0 25.0 0.0013 0.0029 1.000 1.0 14.5 -0.0001 14.000 0.004 16 27 15.0 25.0 0.0009 0.0018 1.000 1.0 15.5 -0.0001 15.000 0.003 17 27 16.0 25.0 0.0005 0.0011 1.000 1.0 16.5 -0.0000 16.000 0.002 18 27 17.0 25.0 0.0004 0.0007 1.000 1.0 17.5 -0.0000 17.000 0.002 19 27 18.0 25.0 0.0002 0.0004 1.000 1.0 18.5 -0.0000 18.000 0.001 20 27 19.0 25.0 0.0002 0.0002 1.000 1.0 19.5 -0.0000 19.000 0.000 21 27 20.0 25.0 0.0001 0.0001 1.000 1.0 20.5 -0.0000 20.000 0.001 22 27 21.0 25.0 0.0001 0.0001 1.000 1.0 21.5 -0.0000 21.000 0.001 23 27 22.0 25.0 0.0000 0.0000 1.000 1.0 22.5 -0.0000 22.000 0.000 24 27 23.0 25.0 0.0000 0.0000 1.000 1.0 23.5 -0.0000 23.000 0.000 25 27 24.0 25.0 0.0000 0.0000 1.000 1.0 24.5 -0.0000 24.000 0.000 26 27 25.0 25.0 0.0000 0.0000 1.000 1.0 25.5 -0.0000 25.000 0.000 27 27 26.0 25.0 -0.0000 0.0000 0.000 0.0 0.0 0.0 12.750 0.078 1 28 0.0 26.0 0.0000 -0.225 0.978 1.0 0.5 0.0011 0.250 -0.058 2 28 1.0 26.0 0.0000 -0.296 0.972 1.0 1.5 0.0011 1.000 -0.059 3 28 2.0 26.0 0.0001 -0.296 0.972 1.0 2.5 0.0011 2.000 -0.061 4 28 3.0 26.0 0.0002 -0.295 0.972 1.0 3.5 0.0011 3.000 -0.062 5 28 4.0 26.0 0.0003 -0.294 0.972 1.0 4.5 0.0012 3.999 -0.066 6 28 5.0 26.0 0.0004 -0.293 0.973 1.0 5.5 0.0013 4.999 -0.069 7 28 6.0 26.0 0.0008 -0.299 0.973 1.0 6.5 0.0017 5.999 -0.079 8 28 7.0 26.0 -0.0003 -0.290 0.984 1.0 7.5 0.0021 6.998 -0.069 9 28 8.0 26.0 -0.0005 0.132 0.997 1.0 8.5 0.0016 7.998 -0.057 10 28 9.0 26.0 -0.0008 0.028 0.999 1.0 9.5 0.0012 8.999 -0.042 11 28 10.0 26.0 -0.0010 0.021 0.999 1.0 10.5 0.0005 9.999 -0.024 12 28 11.0 26.0 -0.0003 0.012 0.999 1.0 11.5 0.0004 11.600 -0.011 13 28 12.0 26.0 -0.0002 0.008 1.000 1.0 12.5 0.0002 12.000 -0.011 14 28 13.0 26.0 -0.0001 0.004 1.000 1.0 13.5 0.0001 13.000 -0.001 15 28 14.0 26.0 -0.0001 0.003 1.000 1.0 14.5 0.0001 14.000 -0.001
PARTICLES PDR=1.0 PDZ=1.0 PXR=26.0 PYB=0.0 PYT=52.0 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=9.0 CYCLE = 232
VELOCITY VECTORS VMAX=2.93081 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=9.0 CYCLE = 232
ISOPYCNICS MIN=0.972275 MAX=1.04901 L=0.972275 H=1.04223 DQ= 0.0077298 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=9.0 CYCLE = 232
ISOTHERMS MIN=0-0.73291 MAX=1.22158 L=-7.3291 H=0.873348 DQ= 0.855170 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=9.0 CYCLE = 232
VORTICITY MIN=-0.0239303 MAX=0.972045 L=-0.0239303 H=0.193570 DQ= 0.0996976 T3AAA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=9.0 CYCLE = 232
T3AA IBA YAQUI PARTICLE-FLUID MOMENTUM EXCHANGE TEST (T3AAA1GB 032272-3) 10872-1 T=1. CYCLE = 232 I J X Y U V SIE RHO VOL D M P 8 26 7.0 24.0 -0.5154 -0.7859 1.053 0.997 7.5 -0.0026 7.015 -0.292 9 26 8.0 24.0 -0.5339 -0.5295 1.045 0.996 8.5 -0.0026 8.016 -0.310 10 26 9.0 24.0 -0.5259 0.2834 0.025 0.996 9.5 -0.0024 9.013 -0.319 11 26 10.0 24.0 -0.4920 0.0920 1.015 0.996 10.5 -0.0022 10.000 -0.320 12 26 11.0 24.0 -0.4417 0.0404 1.011 0.997 11.5 -0.0019 11.000 -0.312 13 26 12.0 24.0 -0.3849 0.1236 1.000 0.997 12.5 -0.0018 12.000 -0.300 14 26 13.0 24.0 -0.3284 0.1719 1.007 0.997 13.5 -0.0013 13.000 -0.283 15 26 14.0 24.0 -0.2785 0.1978 1.004 0.997 14.5 -0.0011 14.000 -0.264 16 26 15.0 24.0 -0.2305 0.2108 1.003 0.997 15.5 -0.0009 15.000 -0.244 17 26 16.0 24.0 -0.1909 0.2165 1.003 0.997 16.5 -0.0004 16.000 -0.225 18 26 17.0 24.0 -0.1572 0.2185 1.001 0.998 17.5 -0.0004 17.000 -0.206 19 26 18.0 24.0 -0.1285 0.2185 1.001 0.998 18.5 -0.0004 18.000 -0.189 20 26 19.0 24.0 -0.1042 0.2175 1.000 0.998 19.5 -0.0004 19.000 -0.174 21 26 20.0 24.0 -0.0633 0.2159 1.000 0.998 20.5 -0.0001 20.000 -0.160 22 26 21.0 24.0 -0.0653 0.2143 1.000 0.999 21.5 -0.0001 21.000 -0.151 23 26 22.0 24.0 -0.0434 0.2120 1.000 0.999 22.5 -0.0001 22.000 -0.143 24 26 23.0 24.0 -0.0353 0.2115 1.000 0.999 23.5 -0.0001 23.000 -0.138 25 26 24.0 24.0 -0.0226 0.2105 1.000 0.999 24.5 -0.0001 24.000 -0.133 26 26 25.0 24.0 -0.0109 0.2098 1.000 0.999 25.5 -0.0000 25.000 -0.132 27 26 26.0 24.0 0.0000 0.2094 0.000 -0.000 0.0 0.0000 12.730 0.078 1 27 0.0 25.0 0.0000 -1.4463 1.185 1.000 0.5 -0.0003 0.250 0.059 2 27 1.0 25.0 -0.0899 -1.4618 1.057 1.000 1.5 -0.0001 0.098 0.012 3 27 2.0 25.0 -0.1772 -1.4288 1.048 0.999 2.5 -0.0001 1.999 -0.025 4 27 3.0 25.0 -0.2801 -1.3520 1.040 0.998 3.5 -0.0002 2.997 -0.080 5 27 4.0 25.0 -0.5398 -1.2374 1.037 0.995 4.5 -0.0002 3.994 -0.098 6 27 5.0 25.0 -0.4018 -1.0853 1.039 0.999 5.5 -0.0002 4.951 -0.131 7 27 6.0 25.0 -0.4544 -0.9002 1.045 0.998 6.5 -0.0002 5.987 -0.184 8 27 7.0 25.0 -0.4903 -0.6585 1.047 0.989 7.5 -0.0002 6.985 -0.185 9 27 8.0 25.0 -0.5654 -0.4621 1.034 0.997 8.5 -0.0002 7.998 -0.220 10 27 9.0 25.0 -0.4950 0.2521 1.019 0.997 9.5 -0.0000 8.006 -0.238 11 27 10.0 25.0 -0.4639 0.2285 1.013 0.997 10.5 -0.0003 9.871 -0.248 12 27 11.0 25.0 -0.4181 0.2410 1.010 0.997 11.5 -0.0003 10.948 -0.253 13 27 12.0 25.0 -0.3644 0.0977 1.007 0.998 12.5 -0.0002 11.987 -0.249 14 27 13.0 25.0 -0.3150 0.1425 1.006 0.998 13.5 -0.0001 12.965 -0.243 15 27 14.0 25.0 -0.2873 0.1587 1.004 0.998 14.5 -0.0001 13.964 -0.233 16 27 15.0 25.0 -0.2241 0.1937 1.003 0.998 15.5 -0.0001 14.984 -0.222 17 27 16.0 25.0 -0.1878 0.1923 1.002 0.998 16.5 -0.0000 15.984 -0.202 18 27 17.0 25.0 -0.1557 0.1994 1.002 0.998 17.5 -0.0000 16.984 -0.187 19 27 18.0 25.0 -0.1282 0.2005 1.001 0.998 18.5 -0.0000 17.865 -0.185 20 27 19.0 25.0 -0.1045 0.2009 1.001 0.998 19.5 -0.0000 18.988 -0.174 21 27 20.0 25.0 -0.0414 0.2000 1.000 0.998 20.5 -0.0000 19.988 -0.165 22 27 21.0 25.0 -0.0616 0.2009 1.000 0.998 21.5 -0.0000 20.967 -0.157 23 27 22.0 25.0 -0.0503 0.2002 1.000 0.998 22.5 -0.0000 21.967 -0.151 24 27 23.0 25.0 -0.0340 0.1987 1.000 0.998 23.5 -0.0000 22.967 -0.147 25 27 24.0 25.0 -0.0231 0.1955 1.000 0.998 24.5 -0.0000 23.565 -0.144 26 27 25.0 25.0 -0.0112 0.1998 1.000 0.998 25.5 0.0000 24.965 -0.143 27 27 26.0 25.0 0.0000 0.1998 0.000 0.000 0.0 0.0 12.750 0.078 1 28 0.0 26.0 0.0000 -1.2989 1.214 1.001 0.5 0.0002 0.250 0.129 2 28 1.0 26.0 -0.0250 -1.2895 1.099 1.000 1.5 0.0014 1.000 0.087 3 28 2.0 26.0 -0.1581 -1.1289 1.088 1.000 2.5 0.0015 2.000 0.063 4 28 3.0 26.0 -0.2458 -1.1824 1.062 1.000 3.5 0.0016 3.000 0.019 5 28 4.0 26.0 -0.3183 -1.0093 1.056 0.999 4.5 0.0017 3.999 -0.015 6 28 5.0 26.0 -0.3800 -0.9513 1.047 0.999 5.5 0.0018 4.999 -0.051 7 28 6.0 26.0 -0.4280 -0.7857 1.043 0.999 6.5 0.0019 5.999 -0.087 8 28 7.0 26.0 -0.4807 -0.5885 1.017 0.999 7.5 0.0019 6.998 -0.120 9 28 8.0 26.0 -0.4727 -0.4011 1.025 0.984 8.5 0.0019 7.994 -0.153 10 28 9.0 26.0 -0.4922 -0.2239 1.015 0.982 9.5 0.0012 8.782 -0.174 11 28 10.0 26.0 -0.4323 -0.0646 1.011 0.980 10.5 0.0018 9.878 -0.192 12 28 11.0 26.0 -0.3308 0.0107 1.009 0.979 11.5 0.0015 10.900 -0.203 13 28 12.0 26.0 -0.3442 0.0756 1.006 0.978 12.5 0.0013 12.000 -0.208 14 28 13.0 26.0 -0.2979 0.1169 1.004 0.979 13.5 0.0011 13.000 -0.209 15 28 14.0 26.0 -0.2547 1.4292 1.002 0.979 14.5 0.0010 13.900 -0.203
1. F H Harlow, A A Amsden, A Numerical Fluid Dynamics Calculation Method for All Flow Speeds, J Comp Phys, 8, 197 (1971).
2. C W Hirt, A A Amsden J L Cook, An Arbitrary Lagrangian-Eulerian Computing Method for All Flow Speeds, submitted to J Comp Phys, (published 1974)
3. F H Harlow, A A Amsden, Fluid Dynamics: A LASL Monograph, LA-4700 (1971).
4. R S Hotchkiss, C W Hirt, Particulate Transport in Highly Distorted Three-Dimensional Flow Fields, VOL II Proceedings of the 1972 SummerComputer Simulation Conference, San Diego, June 14-15, 1972.
5. A A Amsden, C W Hirt, A Simple Scheme for Generating General Curvilinear Grids, J Comp Phys, accepted for publication (March 1973).
The CDC 7600 code listing is difficult to read in places due to the age of the listing. Some of the CDC 7600 Fortran extensions were not a part of Fortran IV (example multiple assignments in a single statement). These have been commented out and the equivalent Fortan IV code below.
Some additional comments have been added, primarily from the flow diageams.
PROGRAM YAQUI
C YLC1 IS STORAGE FOR CELL DATA, YLC2 FOR THE OPTIONAL PARTICLES
LCM /YLC1/ AA1(131000) /YLC2/ AA2(131000)
COMMON /YSC1/ AASC(4242)
COMMON /YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
EQUIVALENCE (AASC(1),X,XPAR),(AASC(2),R,YPAR),(AASC(3),Y,MPAR), EQVREAL
1 (AASC(4),U,UG,DELSH),(AASC(5),V,VG),(AASC(6),RO), EQVREAL
2 (AASC(7),SIE,MP,RMP,RCSQ),(AASC(8),E,ETIL), EQVREAL
3 (AASC(9),RVOL),(AASC(10),M,RM,VP),(AASC(11),p,PL,EP, EQVREAL
4 UP,PM0),(AASC(12),UTIL,UL,CQ,PMX,PU),(AASC(13),VTIL, EQVREAL
5 VL,PMY,PV),(AASC(14),Q,ROL) EQVREAL
REAL LAM,LAMD,M,M6,MC,ML,MP,MPAR,MR,MT,MTE,MU,MUO2,MUO4 EQVREAL
NQ = 14 YAQUI
PRINT 100 YAQUI
10 READ 110, IBAR,JBAR,IUNF,JUNF,JCEN,DR,DZ,CYL,GRDVEL,A0,A0M,B0,KXI YAQUI
CALL ADV(3) YAQUI
CALL LINCNT (64) YAQUI
IF (IBAR) 40,30,20 YAQUI
20 CALL OVERLAY (7LYAQUIFIL,1,0,0) YAQUI
30 CALL OVERLAY (7LYAQUIFIL,2,0,0) YAQUI
GO TO 10 YAQUI
40 CALL EMPTY YAQUI
C YAQUI
100 FORMAT(1H1) YAQUI
110 FORMAT(5I4,7F8.3,I4) YAQUI
END YAQUI
C
SUBROUTINE LOOP
C HANDLES THE 3-ROW BUFFERING THAT SHUTTLES
C CELL DATA BETWEEN LCM AND SCM
COMMON/YSC1/ AASC(4242)
COMMON/YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
C REGULAR 3-ROW UPWARD ROUTINE
CALL ECWR (AASC(IJMS),IECW,NQI,NE)
IECW = IECW + NQI
GO TO (10,20,30) IBUF
C 10 IJP = IJPS = I
10 IJPS = I
IJP = IJPS
IJ = ISC3
C IJM = IJMS = ISC2
IJMS = ISC2
IJM = IJMS
IBUF = 2
GO TO 40
C 20 IJP = IJPS = ISC2
20 IJPS = ISC2
IJP = IJPS
IJ = 1
C IJM = IJMS = ISC3
IJMS = ISC3
IJM = IJMS
IBUF = 3
GO TO 40
C
ENTRY START
IJPS = 1
C IECR = IECW = 0
IECW = 0
IECR = IECW
CALL ECRD (AASC(IJPS),IECR,NQI,NE)
IECR = IECR + NQI
IJPS = ISC2
CALL ECRD (AASC(IJPS),IECR,NQI,NE)
IECR = IECR + NQI
C 30 IJP = IJPS = ISC3
30 IJPS = ISC3
IJP = IJPS
IJ = ISC2
C IJM = IJMS = IBUF = 1
IBUF = 1
IJMS = IBUF
IJM = IJMS
40 CALL ECRD (AASC(IJPS),IECR,NQI,NE)
IECR = IECR + NQI
RETURN
C
ENTRY DONE
CALL ECWR (AASC(IJMS),IECW,NQI,NE)
IECW = IECW + NQI
GO TO (50,60,70) IBUF
50 IJMS = ISC2
GO TO 80
60 IJMS = ISC3
GO TO 80
70 IJMS = 1
80 CALL ECWR (AASC(IJMS),IECW,NQI,NE)
RETURN
C
ENTRY LOOPD
100 CALL ECWR (AASC(IJS),IECW,NQI,NE)
IECW = IECW - NQI
GO TO (110,120,140) IBUF
110 IBUF = 2
IJ = ISCF1
IJS = 1
IJM = ISCF2
IJMS = ISC2
GO TO 130
C
ENTRY STARTD
IJMS = ISC2
C IECR = IECW = ITV
IECW = ITV
IECR = IECW
CALL ECRD (AASC(IJMS),IECR,NQI,NE)
IECR = IECR - NQI
120 IJM = ISCF1
C IJMS = IBUF = 1
IBUF = 1
IJMS = IBUF
IJ = ISCF2
IJS = ISC2
130 IF (IECR.LT.0) GO TO 150
CALL ECRD (AASC(IJMS),IECR,NQI,NE)
IECR = IECR - NQI
140 RETURN
150 IBUF = 3
GO TO 100
C
ENTRY R1ROW
C READ JTH ROW INTO SCM ROW 1/3
IEC = (J-1) * NQI
CALL ECRD (AASC(1),IEC,NQI,NE)
RETURN
C
ENTRY SETIJ
C POINT TO ITH COLUMN IN JTH ROW
IJ = (I-1) * NQ + 1
RETURN
C
ENTRY W1ROW
C RETURN JTH ROW TO LCM
IEC = (I-1)*NQI
C REDUNDANT UNLESS R1ROW ENTRY WAS UNUSED
CALL ECWR(AASC(1),IEC,NQI,NE)
RETURN
END
C
SUBROUTINE FILMCO
C SET UP FILM COORDINATES FOR INITIAL OR REZONED GRIDS
COMMON/YSC1/ AASC(4242)
COMMON/YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
EQUIVALENCE (AASC(1),X,XPAR),(AASC(2),R,YPAR),(AASC(3),Y,MPAR), EQVREAL
1 (AASC(4),U,UG,DELSM),(AASC(5),V,VG),(AASC(6),RO), EQVREAL
2 (AASC(7),SIE,MP,RMP,RCSQ),(AASC(8),E,ETIL), EQVREAL
3 (AASC(9),RVOL),(AASC(10),M,RM,VP),(AASC(11),P,PL,EP, EQVREAL
4 UP,PM0),(AASC(12),UTIL,UL,CQ,PMX,PU),(AASC(13),VTIL, EQVREAL
5 VL,PMY,PV),(AASC(14),Q,ROL) EQVREAL
REAL LAM,LAMD,M,M6,MC,ML,MP,MPAR,MR,MT,MTE,MU,MUO2,MUO4 EQVREAL
DIMENSION X(1),XPAR(1),R(1),YPAR(1),Y(1),MPAR(1),U(1),UG(1), DIMEN
1 DELSM(1),V(1),VG(1),RO(1),SIE(1),MP(1),RMP(1),RCSQ(1), DIMEN
2 E(1),ETIL(1),RVOL(1),M(1),RM(1),VP(1),P(1),PL(1),EP(1), DIMEN
3 UP(1),UTIL(1),UL(1),CQ(1),PMX(1),PU(1),VTIL(1),VL(1), DIMEN
4 PMY(1),PV(1),Q(1),ROL(1),PM0(1) DIMEN
C INITIALIZE COMPARANDS FOR DETERMINATION OF GRID SIZE AND SWEEP ANGLES
XL = 0.0
YB = 10.E+20
XR = YT - YB
CALL START
C ALL VERTICES
DO 129 J=1,JP2
DO 119 I=1,IP1
XR = AMAX1(XR,X(IJ))
YB = AMIN1(YB,Y(IJ))
YT = AMAX1(YT,Y(IJ))
119 IJ = IJ + NQ
CALL LOOP
129 CONTINUE
C VELOCITY VECTOR SCALE
VV = 0.9*XR*RIBAR
C YB ON FILM FRAME
FIYB = 916.0
C RATIO WIDTH/HEIGHT
XD = XR/(YT-YB)
C INITIALIZE FOR TEST
YY = 0.0
C 1.13556 IS 1022/900
C HIGHER THAN WIDER YY1
IF (XD.LE.1.13556) YY=1.
C WIDER THAN HIGHER
C ENSURE NO IMPROPER -0.
FIXL = AMAX1(0.,(511.-450.*XD)*YY)
FIXR = (511.+450.*XD)*YY + 1022*(1.-YY)
FIYT = 16.*YY + (916.-1022./XD)*(1.-YY)
C CONVERSION COEFFICIENTS
XCONV = (FIXR-FIXL)/(XR-XL)
YCONV = (FIYT-FIYB)/(YT-YB)
C INTEGER EQUIVALENTS
IXL = FIXL
IXR = FIXR
IYB = FIYB
IYT = FIYT
C NO PARTICLES, BYPASS REMAINDER OF SUBROUTINE
IF (NPT.EQ.0) RETURN
PXL = 0.0
PYB = YB + PYB
PXR = PDR*FIBP
PYT = PYB + PDZ*FJBP
C JUST OUTSIDE GRID
PXRP = PXR*OPEM10
PYBM = PYB*OPEM10
PYTP = PYT*OPEM10
RPDR = 1./PDR
RPDZ = 1./PDZ
RPDRDZ = RPDR*RPDZ
C YB ON FILM FRAME
FIPYB = 916.0
C RATIO OF WIDTH TO HEIGHT
XD = PXR/(PYT-PYB)
C INITIALIZE FOR TEST
YY = 0.0
C YY=1. IS HIGHER THAN WIDE
IF(XD.LE.1.13556) YY=1.
C OTHERWISE WIDER THAN HIGHER
FIPXL = AMAX1(0.,(511.-450.*XD)*YY)
FIPXR = (511.+450.*XD)*YY + 1022.*(1.-YY)
FIPYT = 16.*YY + (916.-1022./XD) *(1.-YY)
C CONVERSION COEFFICIENTS
PXCONV = (FIPXR-FIPXL)/(PXR-PXL)
PYCONV = (FIPYT-FIPYB)/(PYT-PYB)
C INTEGER EQUIVALENTS
IPXL = FIPXL
IPXR = FIPXR
IPYB = FIPYB
IPYT = FIPYT
RETURN
END
C
PROGRAM YASET
PRINT 10
CALL YASET1
10 FORMAT('YASET CALLED')
END
C
SUBROUTINE YASET1
LCM /YLC1/ AA1(131000) /YLC2/ AA2(131000)
COMMON/YSC1/ AASC(4242)
COMMON/YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
EQUIVALENCE (AASC(1),X,XPAR),(AASC(2),R,YPAR),(AASC(3),Y,MPAR), EQVREAL
1 (AASC(4),U,UG,DELSM),(AASC(5),V,VG),(AASC(6),RO), EQVREAL
2 (AASC(7),SIE,MP,RMP,RCSQ),(AASC(8),E,ETIL), EQVREAL
3 (AASC(9),RVOL),(AASC(10),M,RM,VP),(AASC(11),P,PL,EP, EQVREAL
4 UP,PM0),(AASC(12),UTIL,UL,CQ,PMX,PU),(AASC(13),VTIL, EQVREAL
5 VL,PMY,PV),(AASC(14),Q,ROL) EQVREAL
REAL LAM,LAMD,M,M6,MC,ML,MP,MPAR,MR,MT,MTE,MU,MUO2,MUO4 EQVRWAL
DIMENSION X(1),XPAR(1),R(1),YPAR(1),Y(1),MPAR(1),U(1),UG(1), DIMEN
1 DELSM(1),V(1),VG(1),RO(1),SIE(1),MP(1),RMP(1),RCSQ(1), DIMEN
2 E(1),ETIL(1),RVOL(1),M(1),RM(1),VP(1),P(1),PL(1),EP(1), DIMEN
3 UP(1),UTIL(1),UL(1),CQ(1),PMX(1),PU(1),VTIL(1),VL(1), DIMEN
4 PMY(1),PV(1),Q(1),ROL(1),PM0(1)
C READ IDENTIFICATION CARD
READ 500, NAME
READ 510, MU,LAM,OM,EPS,GR,GZ,ASQ,RON,GM1
READ 515, FREZ,YB,REZY0,REZUE,REZVE,REZVT,REZRON,REZSIE
READ 520, IBP,JBP,PDR,PDZ,PYB,GZP,IMOMX
READ 530, T,DT,T20MD,TLIMD,TWFIN,LPR,ICOLOR
READ 540, (DTO(N),N=1,10)
READ 540, (DTOC(N),N=1,10)
C OUTPUT PRINT FILE
KT = 9
ASSIGN 110 TO KRET
C WRITE THE 8 CARDS READ SO FAR ONTO THE FILESET (KT)
100 WRITE(KT,500) NAME
WRITE(KT,550) IBAR,JBAR,IUNF,JUNF,JCEN,DR,DZ,CYL,GRDVEL,
1 A0,A0M,B0,KXI
WRITE(KT,560) MU,LAM,OM,EPS,GR,GZ,ASQ,RON,GM1
WRITE(KT,565) FREZ,YB,REZY0,REZUE,REZVE,REZVT,REZRON,REZSIE
WRITE(KT,570) IBP,JBP,PDR,PDZ,PYB,GZP,IMOMX
WRITE(KT,580) Y,DZ,T20MD,TLIMD,TWFIN,LPB,ICOLOR
WRITE(KT,590) (DTO(N),N=1,10)
WRITE(KT,600) (DTOC(N),N=1,10)
GO TO KRET
C JUMP IF NO FILM WRITING
110 IF(LPR.EQ.0) GO TO 120
KT = 12
ASSIGN 120 TO KRET
GO TO 100
120 IM1 = IBAR - 1
IM6 = IBAR - 6
IP1 = IBAR + 1
IP2 = IBAR + 2
JM10 = JBAR - 10
JM14 = JBAR - 14
JP1 = JBAR + 1
JP2 = JBAR + 2
JP4 = JBAR + 4
JP4O2 = JP4 / 2
RIBAR = 1./FLOAT(IBAR)
RIBJB = 1./FLOAT(IBAR-JBAR)
NQIB = NQ + IBAR
OMCYL = 1. - CYL
NQI = NQ * IP1
ISC2 = NQI + 1
ISC3 = ISC2 + NQI
ITV = JP1*NQI
ISCF1 = ISC2 - NQ
ISCF2 = ISCF1 + NQI
NSC = LOCF(ZZ) - LOCF(AA) + 1
NLC = LOCF(AAI(JP4*NQI)) - LOCF(AA1) +1
LJP2 = JP2 - JP2/3 * 3
IF (LJP2.EQ.0) LJP2 = 3
LJP2 = LJP2*NQI - NQ + 1
IDTO = I
TOUT = T + DTO(1)
C DT = DTPOS = DT*0.1
DTPOS = DT*0.1
DT = DTPOS
C NCYC = NUMTD = 0
NUMTD = 0
NCYC = NUMTD
EM10 = 10.E-10
OMEM10 = 1. - EM10
OPEM10 = 1. + EM10
ANC = 0.05
OMANC = 1.-ANC
COLAMU = (1.0*1.67)/(LAM+2*MU+EM10)
A0FAC = A0M/(2.0 * (1.+A0M**2))
GGM1 = (GM1+1.)*GM1
THIRD = 1./3.
IUNF = MAX0(IUNF,1)
JUNF = MAX0(JUNF,2)
JUNFO2 = JUNF/2
IF(JCEN.EQ.0) JCEN=JBAR/2
IF(FREZ.NE.1.) ROMFR =1/(1.-FREZ)
C BASIC PARTICLE GENERATOR
CALL PARTGEN
CALL MESHMKR
C SET UP FILM PLOT
CALL FILMCO
200 CALL START
DO 229 J=2,JP1
DO 219 I=1,IBAR
IPJ = IJ * NQ
IPJP = IJP + NQ
X1 = X(IPJ)
Y1 = Y(IPJ)
R1 = R(IPJ)
X2 = X(IPJP)
Y2 = Y(IPJP)
R2 = R(IPJP)
X3 = X(IJP)
Y3 = Y(IJP)
R3 = R(IJP)
X4 = X(IJ)
Y4 = Y(IJ)
R4 = R(IJ)
ATR = .5*(X1*(Y2-Y3)+X2*(Y3-Y1)+X3*(Y1-Y2))
ABL = .5*(X1*(Y3-Y4)+X3*(Y4-Y1)+X4*(Y1-Y3))
M(IJ) = THIRD*(ATR*(R1+R2+R3)+ABL*(R1+R3+R4))*RO(IJ)
RVOL(IJ) = RO(IJ)/M(IJ)
E(IJ) = SIE(IJ)+.125*(U(IPJ)**2+U(IPJP)**2+U(IJP)**2+U(IJ)**2
1 +V(IPJ)**2+V(IPJP)**2+V(IJP)**2+V(IJ)**2)
IJ = IPJ
219 IJP = IPJP
CALL LOOP
229 CONTINUE
CALL DONE
C ALL ROWS COMPLETED
300 CALL STARTD
DO 359 JJ=2,JP2
J = JP4 - JJ
DO 349 II=1,IP1
I = IP2 - II
IMJ = IJ - NQ
IMJM = IJM - NQ
XX = 0.
IF (I.NE.IP1 .AND. J.NE.2 ) XX = M(IJM)
IF (I.NE.IP1 .AND. J.NE.JP2) XX = XX+M(IJ)
IF (I.NE.1 .AND. J.NE.JP2) XX = XX+M(IMJ)
IF (I.NE.1 .AND. J.NE.2 ) XX = XX+M(IMJM)
340 RM(IJ) = 4./XX
C TO NEXT CELL ON LEFT
IJ = IMJ
349 IJM = IMJM
C ROW COMPLETED DROP DOWN A ROW
CALL LOOPD
359 CONTINUE
RETURN
C SETUP COMPLETED
500 FORMAT(10A8)
510 FORMAT(9F8.3)
515 FORMAT(8F8.3)
520 FORMAT(2I4,4F8.3,I4)
530 FORMAT(5F8.3,2I4)
540 FORMAT(10F8.3)
C THE CDC 7600 USED * WHILE NORMALLY ' IS USED AS BELOW
550 FORMAT(3X,'IBAR=',I4/3X,'JBAR=',I4/3X,'IUNF=',I4/3X,'JUNF=',I4/
13X,'JCEN=',I4/5X,'DR=',1PE12.5/5X,'DZ=',E12.5/4X,'CYL=',E12.5/
2' GRDVEL=',E12.5/5X,'AD=',E12.5/4X,'A0M=',E12.5/5X,
3'R0=',E12.5/4X,'XXI=',I2)
C
560 FORMAT(5X,'MU=',1PE12.5/4X,'LAM=',E12.5/5X,'OM=',E12.5/4X,
1'EPS=',E12.5/5X,'GP=',E12.5/5X,'GZ=',E12.5/4X,'ASQ=',E12.5/4X,
2'RON=',E12.5/4X,'GMI=',E12.5)
C
565 FORMAT(3X,'FREZ=',1PE12.5/5X,'YB=',F12.5/' REZY0='
1,E12.5/' RFZUE=', F12.5/' REZVE',E12.5/' REZVT=',
2E12.5/' REZRON=',E12.5/' REZSIE=', E12.5)
C
570 FORMAT(4X,'IBP=',I4/4X,'JRP=',I4/4X,'PDR=',1PE12.5/4X,
1'PDZ=',E12.5/4X,'PYR=',E12.5/4X,'GZP=',E12.5/2X,'IMOMX=',I2)
C
580 FORMAT(6X,'T=',1PE12.5/5X,'DT=',E12.5/' T20MD=',
1E12.5/' TLIND=',E12.5/' TWFIN=',E12.5/4X,'LPR=',I2/
2' ICOLOR=',I2)
C
590 FORMAT(' DTO(1-10)=',5(1PE12.5,2X)/12X,5(E12.5,2X))
C
600 FORMAT(' DTOC(1-10)=',5(1PE12.5,2X)/12X,5(E12.5,2X))
END
C
SUBROUTINE PARTGEN
C GENERATES PARTICLES IF SPECIFIED
LCM /YLC1/ AA1(131000) /YLC2/ AA2(131000)
COMMON/YSC1/ AASC(4242)
COMMON/YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
EQUIVALENCE (AASC(1),X,XPAR),(AASC(2),R,YPAR),(AASC(3),Y,MPAR), EQVREAL
1 (AASC(4),U,UG,DELSM),(AASC(5),V,VG),(AASC(6),RO), EQVREAL
2 (AASC(7),SIE,MP,RMP,RCSQ),(AASC(8),E,ETIL), EQVREAL
3 (AASC(9),RVOL),(AASC(10),M,RM,VP),(AASC(11),p,PL,EP, EQVREAL
4 UP,PM0),(AASC(12),UTIL,UL,CQ,PMX,PU),(AASC(13),VTIL, EQVREAL
5 VL,PMY,PV),(AASC(14),Q,ROL) EQVREAL
REAL LAM,LAMD,M,M6,MC,ML,MP,MPAR,MR,MT,MTE,MU,MUO2,MUO4 EQVRWAL
DIMENSION X(1),XPAR(1),R(1),YPAR(1),Y(1),MPAR(1),U(1),UG(1), DIMEN
1 DELSM(1),V(1),VG(1),RO(1),SIE(1),MP(1),RMP(1),RCSQ(1), DIMEN
2 E(1),ETIL(1),RVOL(1),M(1),RM(1),VP(1),P(1),PL(1),EP(1), DIMEN
3 UP(1),UTIL(1),UL(1),CQ(1),PMX(1),PU(1),VTIL(1),VL(1), DIMEN
4 PMY(1),PV(1),Q(1),ROL(1),PM0(1)
NPT = 0
IF(IBP.EQ.0) RETURN
IF(IBP.GT.IBAR .OR. JBP.GT.JBAR) GO TO 300
FIBP = IBP
FJBP = JBP
RIBP = 1./FIBP
RJBP = 1./FJBP
IBP1 = IBP+1
JBP2 = JBP+2
NQI2 = NQI*2
C IECP = IPAR = LOCF(AA2)
IPAR = LOCF(AA2)
IECP = IPAR
LPB = NQI/3 *3
IMOME3 = IMOMX*1000
KP = 1
MT = 0.0
IF (FREZ.EQ.1.0) GO TO 100
XMAX = (FLOAT(IUNF) + FREZ*(1.-FREZ**(IBAR -IUNF ))
1 *ROMFR)*DR
YMAX = REZY0 + (FLOAT(JUNFO2) + FREZ*(1.-FREZ**(JBAR-JCEN-JUNFO2))
1 *ROMFR)*DZ
YMIN = REZY0 - (FLOAT(JUNFO2) + FREZ*(1.-FREZ**(JCEN -JUNFO2 ))
1 *ROMFR)*DZ
PDR = XMAX*RIBP
PDZ = (YMAX-YMIN)*RJBP
100 READ 900, DRPAR,DZPAR,XC,YC,XD,YD,UPAR,VPAR,MTE,DRAG
IF (DRPAR.LE.0) GOTO 290
PRINT 910, DRPAR,DZPAR,XC,YC,XD,YD,UPAR,VPAR,MTE,DRAG
IF (LPR.GT.0) WRITE(12,910) DRPAR,DZPAR,XC,YC,XD,YD,UPAR,VPAR,
1 MTE,DRAG
C NOT STANDARD FORTRAN IV
USH = SHIFT(UPAR,30).AND.7777777777B
VSH = SHIFT(VPAR,30).AND.7777777777B
DRAG = DRAG.AND..NOT.7777777777B
C
YTOP = YC+XD
YBOT = YC-XD
IF (Y0.LT.EM10) GO TO 200
YTOP = YD
YBOT = YC
IF(FREZ.EQ.1.0) GO TO 200
YTOP = YMAX
YBOT = YMIN
X0 = XMAX
200 YTE = YBOT+.5*DZPAR
210 XTE = XC+.5*DBPAR
220 IF (YD.LE.0. .AND. (YTE-YC)**2+GT*XD**2.GT.XD**2) GO TO 240
XPAR(KP) = (XTE.AND. .NOT. 7777777777B) .OR. USH
YPAR(KP) = (YTE.AND. .NOT. 7777777777B) .OR. VSH
MC = MTE*(XTE*CYL+OMCYL)
MT = MT + MC
C
C MPAR(KP) = DRAG .OR. (SHIFT(MC*30).AND.7777777777B)
C
KP = XP+3
NPT = NPT+1
IF(KP.GT.LPB) GO TO 250
230 XTE = XTE+DBPAR
IF(XTE.LE.XD) GO TO 220
240 YTE = YTE+DZPAR
IF(YTE.LE.YTOP) GO TO 210
GO TO 100
250 CALL ECWR(AASC,IECP,LPB,NE)
IECP = IECP+LPR
KP = 1
GO TO 230
290 CALL ECWR (AASC,IECP,LPB,NE)
NPS = LOCF(AA2(NPT*3)) - IPAR + 1
PRINT 920, NPT,MT
IF(LPR.GT.0) WRITE(12,920) NPT,MT
RETURN
300 PRINT 990
RETURN
C
900 FORMAT (12F8.3)
910 FORMAT(' PAR=',1PE12.5,' DZPAR=',E12.5,' XC=',E12.5,
1' YC=',E12.5,' XD=',E12.5/' YD=',E12.5,' UPAR=',E12.5,
2' VPAR=',E12.5,' YC=',E12.5,' DRAG=',E12.5)
920 FORMAT(4X,I5,' PARTICLESGENERATED. WITH TOTAL MASS=',1PE12.5)
990 FORMAT(' PARTICLE GRID TOO LARGE FOR SCM LAYOUT.')
END
C
SUBROUTINE MESHMKR
C GRID GENERATOR
COMMON/YSC1/ AASC(4242)
COMMON/YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
EQUIVALENCE (AASC(1),X,XPAR),(AASC(2),R,YPAR),(AASC(3),Y,MPAR), EQVREAL
1 (AASC(4),U,UG,DELSM),(AASC(5),V,VG),(AASC(6),RO), EQVREAL
2 (AASC(7),SIE,MP,RMP,RCSQ),(AASC(8),E,ETIL), EQVREAL
3 (AASC(9),RVOL),(AASC(10),M,RM,VP),(AASC(11),p,PL,EP, EQVREAL
4 UP,PM0),(AASC(12),UTIL,UL,CQ,PMX,PU),(AASC(13),VTIL, EQVREAL
5 VL,PMY,PV),(AASC(14),Q,ROL) EQVREAL
REAL LAM,LAMD,M,M6,MC,ML,MP,MPAR,MR,MT,MTE,MU,MUO2,MUO4 EQVRWAL
DIMENSION X(1),XPAR(1),R(1),YPAR(1),Y(1),MPAR(1),U(1),UG(1), DIMEN
1 DELSM(1),V(1),VG(1),RO(1),SIE(1),MP(1),RMP(1),RCSQ(1), DIMEN
2 E(1),ETIL(1),RVOL(1),M(1),RM(1),VP(1),P(1),PL(1),EP(1), DIMEN
3 UP(1),UTIL(1),UL(1),CQ(1),PMX(1),PU(1),VTIL(1),VL(1), DIMEN
4 PMY(1),PV(1),Q(1),ROL(1),PM0(1)
C
C STANDARD CLEARING OF LCM CELL STORAGE (PLUS SOME EXTRA)
NQIM = NQI-1
CALL START
DO 119 J=1,JP4
K = IJM + NQIM
DO 109 I=IJM,K
109 AASC(I) = 0.0
CALL LOOP
119 CONTINUE
CALL DONE
200 XX = 0.0
YY = YB
CALL START
DO 229 J=2,JP2
DO 219 I=1,IP1
X(IJ) = XX
Y(IJ) = YY
R(IJ) = XX*CYL+OMCYL
XX = XX + DR
219 IJ = IJ*NQ
XX = 0.
YY = YY + DZ
CALL LOOP
229 CONTINUE
CALL DONE
C UNIFORM MESH BYPASS THE REMAINDER of 200 LOOP
IF(FREZ .EQ.1.0) GO TO 300
JCEN = JCEN + 2
JTOP = JCEN + JUNFO2
JBOT = JCEN-JUNFO2
TJ = FLOAT(JUNFO2) + DZ
CALL START
DO 249 J=2,JP2
DO 239 I=1,IP1
IMJ = IJ - NQ
IF (I.GT.IUNF+1) X(IJ) =X(IMJ) + FREZ*(X(IMJ)-XK(IMJ-NQ))
R(IJ) = X(IJ)*CYL + OMCYL
JDT = IABS(J-JTOP)
JDB = IABS(j-JBOT)
IF(J.LT.JBOT) Y(IJ) = REZY0 - TJ - DZ*FREZ*(1.-FREZ**JDB)*ROMFR
IF(J.GT.JTOP) Y(IJ) = REZY0 + TJ + DZ*FREZ*(1.-FREZ**JDT)*ROMFR
IF(J.EQ.2) YB = Y(IJ)
239 IJ = IJ + NQ
CALL LOOP
249 CONTINUE
CALL DONE
C READ A FLUID REGION CARD
300 READ 1000, NB,NR,NT,NL,UI,vI,ROI,SIEI
IF(NR.EQ.0) RETURN
C GO TO SPECIAL ATMOSPHERTC SETUP
IF(NR.EQ.1000) GO TO 400
PRINT 1010, NB,NR,NT,NL,UI,VI,ROI,SIEI
IF(LPR.GT.0) WRITE(12,1010) NB,NR,NT,NL,UI,VI,ROI,SIEI
NB2 = NB + 2
NR1 = NR + 1
NT2 = NT + 2
NL1 = NL + 1
DO 329 J=NB2,NT2
CALL R1ROW
DO 319 I=NL1,NR1
CALL SETIJ
IF(I.GT.1 .AND. I.LT.IP1) U(IJ)=UI
IF(J.GT.2 .AND. J.LT.JP2) V(IJ)=VI
IF(J.EQ.NT2 .OR. I.EQ. NR1) GO TO 319
RO(IJ) = ROI
SIE(IJ) = SIEI
319 CONTINUE
CALL W1ROW
329 CONTINUE
GO TO 300
400 XX = GM1 * REZSIE
YY = .5*ABS(GZ)
CALL START
YJC2 = .5*(Y(IJP)+Y(IJ))
ROSAV = REZRON*EXP(-GZ*(REZY0-YJC2)/XX)
FNUM = (Y(IJP)-Y(IJ))*YY
FDEN = FNUM*FREZ
ROJ1 = ROSAV*(XX+FNUM)/(XX-FDEN)
DO 459 I=1,IP1
RO(IJ) = ROSAV
RO(IJM) = ROJ1
E(IJ) = E(IJM) -REZSIE
IJ = IJ + NQ
459 IJM = IJM + NQ
CALL LOOP
DO 479 J=3,JP1
FDEN = (Y(IJP)-Y(IJ))*YY
FNUM = (Y(IJ)-Y(IJM))*YY
ROSAV = ROSAV*(XX-FNUM)/(XX+FDEN)
DO 469 I=1,IP1
RO(IJ) = ROSAV
E(IJ) = REZSIE
469 IJ = IJ + NQ
CALL LOOP
479 CONTINUE
FNUM = FNUM*FREZ
FDEN = FDEN*FREZ
ROJP2 = ROSAV *(xx-FNUM)/(XX+FDEN)
DO 489 I = 1,IP1
RO(IJ) = ROJP2
E(IJ) = REZSIE
489 IJ = IJ + NQ
CALL DONE
500 READ 1020, I,JJ,ROI,SIEI,VI,UI
IF(I.EQ.0) GO TO 520
J = JJ + NT-1
CALL R1ROW
CALL SETIJ
RO(IJ) = ROI
SIE(IJ) = SIEI
UL(IJ) = UI
VL(IJ) = VI
CALL R1ROW
CALL SETIJ
VL(IJ) = -VI
CALL W1ROW
GO TO 500
520 CALL START
DO 549 J=2,JBAR
DO 539 I=1,IM1
IPJ = IJ+NQ
IPJP = IJP+NQ
IF(I.EQ.1) V(IJP) = VL(IJ)
U(IPJP) = .5*(UL(IJ)+UL(IJP))
V(IPJP) = .5*(VL(IJ)+VL(IJP))
IJ = IPJ
539 IJP = IPJP
CALL LOOP
549 CONTINUE
CALL DONE
RETURN
C
1000 FORMAT(4I4,4F8.3)
1010 FORMAT('NB=',I3,'NR=',I3,'NT=',I3,'NL=',I3,'UI=',E12.5,' VI='
1,E12.5,' ROI=',E12.5,' SIEI=',E12.5)
1020 FORMAT(2I5,4(4X,F11.5))
END
C
SUBROUTINE YAQUI2
LCM /YLC1/ AA1(131000) /YLC2/ AA2(131000)
COMMON/YSC1/ AASC(4242)
COMMON/YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
EQUIVALENCE (AASC(1),X,XPAR),(AASC(2),R,YPAR),(AASC(3),Y,MPAR), EQVREAL
1 (AASC(4),U,UG,DELSM),(AASC(5),V,VG),(AASC(6),RO), EQVREAL
2 (AASC(7),SIE,MP,RMP,RCSQ),(AASC(8),E,ETIL), EQVREAL
3 (AASC(9),RVOL),(AASC(10),M,RM,VP),(AASC(11),p,PL,EP, EQVREAL
4 UP,PM0),(AASC(12),UTIL,UL,CQ,PMX,PU),(AASC(13),VTIL, EQVREAL
5 VL,PMY,PV),(AASC(14),Q,ROL) EQVREAL
REAL LAM,LAMD,M,M6,MC,ML,MP,MPAR,MR,MT,MTE,MU,MUO2,MUO4 EQVRWAL
DIMENSION X(1),XPAR(1),R(1),YPAR(1),Y(1),MPAR(1),U(1),UG(1), DIMEN
1 DELSM(1),V(1),VG(1),RO(1),SIE(1),MP(1),RMP(1),RCSQ(1), DIMEN
2 E(1),ETIL(1),RVOL(1),M(1),RM(1),VP(1),P(1),PL(1),EP(1), DIMEN
3 UP(1),UTIL(1),UL(1),CQ(1),PMX(1),PU(1),VTIL(1),VL(1), DIMEN
4 PMY(1),PV(1),Q(1),ROL(1),PM0(1)
C
DIMENSION AT(100),FT(100),IX1(1),IX2(1),IY1(1),IY2(1),XCO(4),YCO(4
1),CON(100)
EQUIVALENCE (AT,IX1),(AT(2),IX2),(AT(3),IY1),(AT(4),IY2),(AT(5),
1 XCO),(AT(9),YCO),(FT,CON)
COMMON /YSC3/ JNM
CALL SECOND(TBASE)
T1 = TBASE
CALL GETQ(4LKJBN,JNM)
CALL H4020
CALL GETQ (4LKTLM,II)
TLIM=II
C DTVSAV = DTCSAV = 0.0
DTCSAV = 0.0
DTVSAV = DTCSAV
C GO TO 370 FOR TAPE RESTART
IF(IBAR.EQ.0) GO TO 370
C SETUP FOR OPTIONAL TIME-LIMIT TAPE DUMP
TLIM = TLIM*27.5E-9 -40. + (1.-TLIMD)*1.E+10
C START READS BOTTOM 3 ROWS OG GRID FROM LCM TO SCM
CALL START
CIRC = 0.
C INITIALIZE CIRCULATION ALL CELLS
DO 199 J=2,JP1
DO 189 I=1,IBAR
IPJ = IJ * NQ
IPJP = IJP + NQ
IF (I.EQ.1 ) CIRC = CIRC + 0.5*(V(IJ)+V(IJP1))+(Y(IJP)-Y(IJ))
IF (I.EQ.IM1 ) CIRC = CIRC - 0.5*(V(IJ)+V(IJP1))+(Y(IJP)-Y(IJ))
IF (J.EQ.3 ) CIRC = CIRC - 0.5*(U(IJ)+U(IJP1))+(X(IJP)-X(IJ))
IF (J.EQ.JBAR) CIRC = CIRC + 0.5*(U(IJ)+U(IJP1))+(X(IJP)-X(IJ))
SIET = AMAX1(SIE(IJ),0.)
P(IJ) = ASQ*(RO(IJ)-RON) + GM1*RO(IJ)*SIET
IF(ASQ.LT.1.E-6) GO TO 180
P(IJ) = ASQ*(ROL(IJ)-RON) + GM1*RO(IJ)*SIET
C
IF(NCYC.EQ.0) P(IJ) = -RON*GZ*(YT-YB-IFLOAT(J)*1.51*DZ)
180 IJP = IPJP
189 IJ = IPJ
C PROCESSES ALL ROWS FROM BOTTOM UPWARDS
C NQ=14 WORDS PER CELL
CALL LOOP
199 CONTINUE
C CLEANS UP BY MOVING LAST TWO ROWS BACK TO LCM
CALL DONE
200 TOLO = TP
CALL SECOND(12)
XX = (T2-TOLD)*RIBJB
IF(LPR.GT.0) WRITE(12,4000) T,NCYC,DT,XX,CIRC,DTV,IDTV,JDTV,
1 NUMIT,T2,DTC,IDTC,JDTC
CALL EMPTY
PRINT 4000, T,NCYC,DT,XX,CIRC,DTV,IDTV,JDTV,
1 NUMIT,T2,DTC,IDTC,JDTC
IF(T*EM10 .GE. TOUT) GO TO 290
IF(NCYC.LE.1 .OR. NUMIT.GT.499) GO TO 300
210 IF(T2-T1.GE.1200. .AND.T2GMD.GT.EM10) GO TO 320
IF(T2-TBASE.GE.TLIM) GO TO 330
220 IF(T.GE.TWFIN) GO TO 340
NCYC = NCYC + 1
C IF(NCYC.EQ.2) DT = DTPOS = DT*10.
C IF(NCYC.GE.3) DT = DTPOS = AMIN1(DTV,DTC)
C
IF(NCYC.EQ.2) DTPOS = DT*10.
IF(NCYC.EQ.2) DT = DTPOS
IF(NCYC.GE.3) DTPOS = AMIN1(DTV,DTC)
IF(NCYC.GE.3) DT = DTPOS
C
DTFAC= 1.25
C DTV = DTC = DT*DTFAC
DTC = DT*DTFAC
DTV = DTC
IF (T+DY.GT.TOUT) DT = TOUT-T
T = T + DT
RDT = 1./DT
DTO2 = .5*DT
DTO4 = .25*DT
DTO8 = .125*DT
DTO16 = .0625*DT
DT8 = DT*8.
DTGR = DT*GR
DTGZ = DT*GZ
DTGZP = DT*GZP
GO TO 1000
290 TOUT = TOUT + DTO(IDTO)
IF (T + EM10 .LT.DTOC(IDTO)) GO TO 300
TOUT = DTOC(IDTO) + DTO(IDTO+1)
IDTO = IDTO + 1
300 ASSIGN 310 TO KRET
GO TO 500
310 ASSIGN 210 TO KRET
GO TO 400
320 T1 = T2
ASSIGN 220 TO KRET
GO TO 350
330 ASSIGN 340 TO KRET
GO TO 350
340 RETURN
350 PRINT 4010, NUMTD,T,NCYC
IF (LPR.GT.0) WRITE(12,4010) NUMTD,T,NCYC
WRITE(8) (AA(N),N=1,NSC)
WRITE(8) (AA1(N),N=1,NLC)
IF(NPT.GT.0) WRITE(8) (AA2(N),N=1,NPS)
CALL DATAREL (5LFSETB)
CALL AFSREL(3LOUT)
CALL AFSREL(4LFILM)
NUMTD = NUMTD + 1
GO TO KRET
370 REWIND 7
JTD = JRAR
JNSC = LOCF(ZZ) - LOCF(AA) + 1
READ(7) (AA(N),N=1,JNSC)
IF (JTD.NE.NUMTO) GO TO 380
READ(7) (AA1(N),N=1,NLC)
IF(NPT.GT.0) READ(7) (AA2(N),N=1,NPS)
CALL AFSRFL (5LFSEIT)
NUMTD = NUMTD + 1
TLIM = TLIM*27.5E-9 - 40. +(1.-TLIMD)*1.E+10
PRINT 4020, JTD
PRINT 4030, NAME,T,NCYC
IF(LPR.EQ.0) GO TO 220
WRITE(12,4020) JTD
WRITE(12,4030) NAME,T,NCYC
GO TO 220
380 PRINT 4040
WRITE(12,4040)
CALL AFSREL (5LFSEI7)
RETURN
400 IF(LPR.EQ.0) GO TO KRET
ASSIGN 440 TO KRF
ASSIGN 440 TO KRP
ASSIGN 460 TO KRFP
GO TO (420,410,458) LPR
C THIS LOOKS LIKE ILLEGAL FORTRAN AS LABEL 458 IS INSIDE TWO DO LOOPS 479,489
410 ASSIGN 458 TO KRF
ASSIGN 456 TO KRFP
420 CALL LINCNT (64)
CALL ADV (1)
GO TO 454
440 KRF = KRFP
ASSIGN 460 TO KRP
CALL START
C
DO 489 J=1,JP2
DO 479 I=1,IP1
C
IPJM = IJM + NQ
IPJ = IJ + NQ
C D = PRM = PRV = PRSIR = 0.
PRSIR = 0.
PRV = PRSIR
PRM = PRV
D = PRM
IF (J.EQ.1) GO TO 450
IF(RM(IJM).NE.0.) PRM=1./RM(IJM)
IF (I.EQ.IP1 .OR. J.EQ.JP2) GO TO 450
PRSIE = SIE(IJM)
IF(RVOL(IJM).NE.6.) PRV=1./RVOL(IJM)
X1 = X(IPJM)
Y1 = Y(IPJM)
R1 = R(IPJM)
U1 = U(IPJM)
V2 = V(IPJM)
X2 = X(IPJ)
Y2 = Y(IPJ)
R2 = R(IPJ)
U2 = U(IPJ)
V2 = V(IPJ)
X3 = X(IJ)
Y3 = Y(IJ)
R3 = R(IJ)
U3 = U(IJ)
V3 = V(IJ)
U4 = U(IJM)
V4 = V(IJM)
D = .25*RVOL(IJM) + ((R1+R2)*((U1+U2)*(Y2-Y1)+(V1+V2)*(X1-X2))
1 + (R2+R3)*((U2+U3)*(Y3-Y2)+(V2+V3)*(X2-X3))
2 + (R3+R4)*((U3+U4)*(Y4-Y3)+(V3+V4)*(X3-X4))
3 + (R4+R1)*((U4+U1)*(Y1-Y4)+(V4+V1)*(X4-X1)))
450 GO TO (452,452,456) LPR
452 WRITE(12,4070) I,J,X(IJM),Y(IJM),U(IJM),V(IJM),PRSIE,RD(IJM),PRV,
1 D,PRM,P(IJM)
LINESF = LINESF + 1
IF(LINESF.LT.62) GO TO KRF
454 LINESF= 0
WRITE(12,4080) JNM,NAME,T,NCYC
WRITE(12,4060)
GO TO KRF
456 PRINT 4070, I,J,X(IJM),Y(IJM),U(IJM),V(IJM),PRSIE,RD(IJM),PRV,
1 D,PRM,P(IJM)
LINESP = LINESP + 1
IF(LINESP.LT.57) GO TO KRP
458 LINESP = 0
PRINT 4050
PRINT 4080, JNM,NAME,T,NCYC
PRINT 4060
GO TO KRP
460 IJ = IPJ
479 IJM = IPJM
CALL LOOP
489 CONTINUE
IF(LPR.GT.2) GO TO KRET
490 CALL EMPTY
GO TO KRET
C PLOT PARTICLES
500 IF(NPT.GT.0) CALL PARPLOT
IF(LPR.EQ.0) GO TO 490
IF(GRDVEL.GT.EM16 .OR. NCYC .EQ.0) CALL ADV(1)
C DPMIN = D2MIN = 1.E+20
D2MIN = 1.E+20
DPMIN = D2MIN
C DRMAX = DZMAX = VMAX = 0.
VMAX = 0.
DZMAX = VMAX
DRMAX = DZMAX
CALL START
DO 549 J=2,JP1
DO 539 I=1,IBAR
IPJ = IJ + NQ
IPJP = IJP + NQ
VMAX = AMAX1 (VMAX,ABS(U(IJ)),ABS(V(IJ)) )
IF (NCYC.GT.0 .AND. GRDVEL.LT.EM10) GO TO 530
X1 = X(IPJ)
X2 = X(IPJP)
X3 = X(IJP)
X4 = X(IJ)
Y1 = Y(IPJ)
Y2 = Y(IPJP)
Y3 = Y(IJP)
Y4 = Y(IJ)
XY14 = SQRT((X1-X4)**2 + (Y1-Y4)**2)
XY23 = SQRT((X2-X3)**2 + (Y2-Y3)**2)
DRMIN = AMIN1(DRMIN,XY14,XY23)
DRMAX = AMAX1(DRMAX,XY14,XY23)
XY21 = SQRT((X2-X1)**2 + (Y2-Y1)**2)
XY34 = SQRT((X3-X4)**2 + (Y3-Y4)**2)
DZMIN = AMIN1(DZMIN,XY21,XY34)
DZMAX = AMAX1(DZMAX,XY21,XY34)
IX1 = FIXL + (X1-XL)*XCONV
IY1 = FIYB + (Y1-YB)*YCONV
IX2 = FIXL + (X2-XL)*XCONV
IY2 = FIYB + (Y2-YB)*YCONV
IX3 = FIXL + (X3-XL)*XCONV
IY3 = FIYB + (Y3-YB)*YCONV
IX4 = FIXL + (X4-XL)*XCONV
IY4 = FIYB + (Y4-YB)*YCONV
IF (I.EQ.1) CALL DRV(IX3,IY3,IX4,IY4)
IF (J.EQ.2) CALL DRV(IX4,IY4,IX1,IY1)
CALL DRV(IX1,IY1,IX2,IY2)
CALL DRV(IX2,IY2,IX3,IY3)
530 IJ = IPJ
539 IJP = IPJP
CALL LOOP
549 CONTINUE
IF(NCYC.GT.0 .AND. GRDVEL.LT.EM10) GO TO 550
CALL LINCNT(59)
WRITE(12,4140) DRMIN,DRMAX,DZMIN,DZMAX,XR,YB,YT
WRITE(12,4080) JNM,NAME,T,NCYC
550 IF(VMAX.LT.EM10) GO TO 600
DROU = VV/VMAX
CALL ADV(1)
CALL START
DO 599 J=2,JP2
DO 589 I=1,IP1
C ASSIGNMENTS OF FLOATING POINT RHS TO INTEGER LEFT. IS THAT ALLOWED?
C COMPILER DOES NOT LIKE IT
IX1 = FIXL + (X(IJ)-XL)*XCONV
IY1 = FIYB + (Y(IJ)-YB)*YCONV
IX2 = FIXL + (X(IJ)+U(IJ)*DROU-XL)*XCONV
IY2 = FIYB + (Y(IJ)+V(IJ)*DROU-YB)*YCONV
IF (IY2 .GE. 1) GO TO 580
IX2 = IX1 + (IX2-IX1)*(IY1-1)/(IY1-IY2)
IY2 = 1
580 CALL DRV (IX1,IY1,IX2,IY2)
C DRAW VECTOR THEN + SIGN at VECTEX TO DENOTE VE
C VECTOR ORIGIN
CALL PLT (IX1,IY1,16)
589 IJ = IJ+NQ
CALL LOOP
599 CONTINUE
CALL LINCNT(59)
WRITE(12,4150) VMAX
WRITE(12,4080) JNM,NAME,T,NCYC
C CONTOUR PLOT
600 IF (IRAR.EQ.1 .OR. JBAR .EQ.1) GO TO 490
L = 0
610 L = L+1
GO TO (620,620,640,490),L
620 CALL START
DO 639 J=2,JP1
DO 629 I=1,IBAR
CQ(IJ) = RO(IJ+L-1)
629 IJ = IJ + NQ
CALL LOOP
639 CONTINUE
CALL DONE
GO TO 700
640 CALL START
DO 659 J=2,JP1
DO 649 I=1,IBAR
IPJ = IJ + NQ
IPJP = IJP + NQ
X1 = X(IPJ)
Y1 = Y(IPJ)
U1 = U(IPJ)
V1 = V(IPJ)
X2 = X(IPJP)
Y2 = Y(IPJP)
U2 = U(IPJP)
V2 = V(IPJP)
X3 = X(IJP)
Y3 = Y(IJP)
U3 = U(IJP)
V3 = V(IJP)
X4 = X(IJ)
Y4 = Y(IJ)
U4 = U(IJ)
V4 = V(IJ)
R1 = .125*RVOL(IJ)*(R(IPJ)+R(IPJP)+R(IJP)+R(IJ))
CQ(IJ) = R1*((U1+U4)*(X1-X4)+(V1+V4)*(Y1-Y4)
2 +(U2+U1)*(X2-X1)+(V2+V1)*(Y2-Y1)
3 +(U3+U2)*(X3-X2)+(V3+V2)*(Y3-Y2)
4 +(U4+U3)*(X4-X3)+(V4+V3)*(Y4-Y3))
IJ = IPJ
649 IJP = IPJP
CALL LOOP
659 CONTINUE
CALL DONE
C INITIALIZE LIMITS AND SEARCH FOR MIN AND MAX FIELD VALUES
700 QMN = 1.E+6
QMX = -QMN
CALL START
DO 719 J=2,JP1
DO 709 I=1,IBAR
QMN = AMIN1 (CQ(IJ),QMN)
QMX = AMAX1 (CQ(IJ),QMX)
709 IJ = IJ +NQ
CALL LOOP
719 CONTINUE
XX = QMX/(QMN+EM10)
IF (XX.LE.2.0) GO TO 735
K = 10./ALOG10(XX)
XX = K + 1
DQ = 10.**(1./XX)
K = ALOG10(QMN)
KK = 10.**(K-1)
K = 1
720 XX = XX*DQ
IF(XX.LT.QMN) GO TO 720
730 CON(K) = XX
IF(XX.GT.QMX) GO TO 740
K = K+1
XX = XX*DQ
GO TO 730
735 XX = QMX - QMN
IF (ABS(XX) .LT. 1.E-3*AMAX1(ABS(QMX),ABS(QMN))) GO TO 610
DQ = 0.1*(XX+.001)
DO 739 K=1,11
739 CON(K) = QMN+(FLOAT(K-1))*DQ
K = 11
740 CALL ADV(1)
CALL LINCNT (59)
GO TO (750,760,770) L
750 WRITE(12,4090)
GO TO 780
760 WRITE(12,4100)
GO TO 780
770 WRITE(12,4110)
780 WRITE(12,4120) QMN,QMX,CON(1),CON(K-1),DQ
WRITE(12,4080) JNM,NAME,T,NCYC
CALL START
DO 899 J=2,JBAR
CALL LOOP
DO 889 I=1,IM1
IPJ = IJ + NQ
IPJM = IJM + NQ
N = 0
DO 879 KK=1,K
C K1 = K2 = K3 = K4 = 0
K4 = 0
K3 = K4
K2 = K3
K1 = K2
IF (CQ(IJM) .LE. CQN(KK)) K1 =1
IF (CQ(IPJM).LE. CQN(KK)) K2 =1
IF (CQ(IJ) .LE. CQN(KK)) K3 =1
IF (CQ(IPJ) .LE. CQN(KK)) K4 =1
IF( K1*K2*K3*K4 .NE. 0 .OR. K1+K2+K3+K4 .EQ.0) GO TO 879
IF(N.GT.0) GO TO 800
IJB = IJM
IJA = IJ
DO 799 JJ=1,2
DO 789 II=1,2
IPJB = IJB+NQ
IPJA = IJA+NQ
N = N + 1
XCO(N)=.25*(X(IPJB)+X(IPJA)+X(IJA)+X(IJB))
YCO(N)=.25*(Y(IPJB)+Y(IPJA)+Y(IJA)+Y(IJB))
IJA = IPJA
789 IJB = IPJB
IJB = IJ
799 IJA = IJP
800 LL = 0
IF (K1+K3 .NE. 1) GO TO 810
IC1 = 1
IC2 = 3
IJ1 = IJM
IJ2 = IJ
ASSIGN 810 TO KR1
GO TO 840
810 IF(K1+K2.NE.1) GO TO 820
IC1 = 1
IC2 = 2
IJ1 = IJM
IJ2 = IPJM
ASSIGN 820 TO KR1
GO TO 840
820 IF(K2+K4.NE.1) Go TO 830
IC1 = 2
IC2 = 4
IJ1 = IPJM
IJ2 = IPJ
ASSIGN 830 TO KR1
GO TO 840
830 IF(K3+K4.NE.1) GO TO 879
IC1 = 3
IC2 = 4
IJ1 = IJ
IJ2 = IPJ
ASSIGN 879 TO KR1
840 LL = LL+1
XX = (CON(KK)-CQ(IJ1))/(CQ(IJ2)-CQ(IJ1))
IX1(LL) = FIXL * (XCO(IC1)+XX*(XCO(IC2)-XCO(IC1))-XL)*XCONv
IY1(LL) = FIYB * (YCO(IC1)+XX*(YCO(IC2)-YCO(IC1))-YB)*YCONv
IF(LL.LT.2) GO TO KR1
CALL DRV(IX1,IY1,IX2,IY2)
IF (KK.EQ.1) CALL PLT (IX1,IY1,35)
IF(KK.EQ.M-1) CALL PLT (IX1,IY1,24)
LL = 0
IF (IJ2.EQ.IPJM) GO TO 420
879 CONTINUE
IJM = IPJM
IJ = IPJ
889 IJP = IJP+NQ
899 CONTINUE
IPJP = IJP + NQ
IX1 = FIXL + (X(IPJ) -XL)*XCONV
IY1 = FIYB + (Y(IPJ) -YB)*YCONV
IX2 = FIXL + (X(IPJP)-XL)*XCONV
IY2 = FIYB + (Y(IPJP)-YB)*YCONV
IX3 = FIXL + (X(IJP) -XL)*XCONV
IY3 = FIYB + (Y(IJP) -YB)*YCONV
IX4 = FIXL + (X(IJ) -XL)*XCONV
IY4 = FIYB + (Y(IJ) -YB)*YCONV
IF (I.EQ.1) CALL DRV(IX3,IY3,IX4,IY4)
IF (J.EQ.2) CALL DRV(IX4,IY4,IX1,IY1)
IF (I.EQ.IBAR) CALL DRV(IX1,IY1,IX2,IY2)
IF (J.EQ.JP1) CALL DRV(IX2,IY2,IX3,IY3)
IJ = IPJ
939 IJP = IPJP
CALL LOOP
949 CONTINUE
GO TO 610
1000 CALL START
Y1 = ANC*RDT
DO 1099 J=2,JP2
DO 1089 I=1,IP1
IMJ = IJ-NQ
IPJ = IJ+NQ
IPJP = IJP+NQ
C XX = YY = 1.0
YY = 1.0
XX = YY
C X1 = U1 = V1 = 0.0
V1 = 0.0
U1 = V1
X1 = U1
IF(I.EQ.1) GO TO 1015
U1 = U(IMJ)
V1 = V(IMJ)
GOTO 1020
1015 X1 = 1.0
XX = 0.0
1020 IF (J.EQ.2) GO TO 1025
U1 = U1+U(IJM)
V1 = V1+V(IJM)
GO TO 1030
1025 X1 = X1+1.0
YY = 0.0
1030 IF (I.EQ.IP1) GO TO 1035
U1 = U1+U(IPJ)
V1 = V1+V(IPJ)
GO TO 1040
1035 X1 = X1+1.0
XX = 0.0
1040 IF (J.EQ.JP2) GO TO 1045
U1 = U1+U(IJP)
V1 = V1+V(IJP)
GO TO 1050
1045 X1 = X1+1.0
YY = 0.0
1050 X1 = 1./(4.0-X1)
AX = GR +Y1*(U1*XI-U(IJ))
AY = GZ +Y1*(V1*XI-V(IJ))
UTIL(IJ) = (U(IJ)+DT*AX)*XX
VTIL(IJ) = (V(IJ)+DY*AY)*YY
Q(IJ) = DT*(AX*U(IJ)+AY*V(IJ))
IJ = IPJ
IJP = IPJP
1089 IJM = IJM+NQ
CALL LOOP
1099 CONTINUE
CALL DONE
1100 CALL START
DO 1299 J=2,JP1
DO 1199 I=1,IBAR
IPJ = IJ + NQ
IPJP = IJP + NQ
X1 = X(IPJ)
Y1 = Y(IPJ)
R1 = R(IPJ)
U1 = U(IPJ)
V1 = V(IPJ)
X2 = X(IPJP)
Y2 = Y(IPJP)
R2 = R(IPJP)
U2 = U(IPJP)
V2 = V(IPJP)
X3 = X(IPJ)
Y3 = Y(IPJ)
R3 = R(IJP)
U3 = U(IJP)
V3 = V(IJP)
X4 = X(IJ)
Y4 = Y(IP)
R4 = R(IP)
U4 = U(IP)
V4 = V(IP)
X12 = X1-X2
X23 = X2-X3
X34 = X3-X4
X41 = X4-X1
X24 = X2-X4
X31 = X3-X1
Y21 = Y2-Y1
Y32 = Y3-Y2
Y43 = Y4-Y3
Y14 = Y1-Y4
Y24 = Y2-Y4
Y31 = Y3-Y1
R12 = R1+R2
R23 = R2+R3
R34 = R3+R4
R41 = R4+R1
HR13 = .5*(R1+R3)
HR24 = .5*(R2+R4)
U12 = U1+U2
U23 = U2+U3
U34 = U3+U4
U41 = U4+U1
U24 = U2+U4
U13 = U1+U3
V12 = V1+V2
V23 = V2+V3
V34 = V3+V4
V41 = V4+V1
V24 = V2+V4
V13 = V1+V3
DTO2M1 = DTO2*RM(IPJ)
DTO2M2 = DTO2*RM(IPJP)
DTO2M3 = DTO2*RM(IJP)
DTO2M4 = DTO2*RM(IJ)
XY = X24*Y31-X31*Y24
D = .25*RVOL(IJ)*(R12*(U12*Y21+V12*X12)
1 +R23*(U23*Y32+V23*X23)
2 +R34*(U34*Y43+V34*X34)
3 +R41*(U41*Y14+V41*X41) )
XX = .5*(X2-X4+X1-X3)
YY = .5*(Y2-Y4+Y3-Y1)
IF (KKI.LT.0) GO TO 1130
AK = RO(IJ)**XXI
GO TO 1140
1130 VELIJ = U4**2 + V4**2
VELMX = 0.7 * AMAX1(ABS(U4+XX),ABS(Y4+YY))
AK = RO(IJ)*COLAMU*(DTO2*VELIJ+VELMX)
1140 LAMD = AMIN1(D,0.)*AK*LAM
MUO2 = .5*AK*MU
MUO4 = 5*MUO2
XX = XX*XX
YY = YY*YY
IF(KXI.LT.0 .AND. DT.LT.DTPOS)
1 AK = RO(IJ)*COLAMU+(.5*DTPOS*VELIJ * VELMX)
DQ = RO(IJ)*OMANC*XX*YY/(2.*AK*(LAM+2.*MU)*(XX+YY)+EM10)
DTV = AMIN1(DTV,.5*DQ)
IF (DTVSAV.NE.DIV) IDTV = I
IF (DTVSAV.NE.DIV) JDTV = J
DTVSAV = DTV
PIXX = MUO2*RVOL(IJ)*(R12*U12+Y21+R23*U23+Y32
1 +R34*U34+Y43+R41*U41*Y14
2 -.5*CYL*(U12+U34)*XY) +LAMD
PIYY = MUO2*RVOL(IJ)*(R12*V12+Y21+R23*V23+X23
1 +R34*V34+X34+R41*V41*X141)
2 +LAMD
PIXY = MUO4*RVOL(IJ)*(R12*(U12*X12+V12*Y21)+R23*(U23*X23+V23*Y32)
1 +R34*(U34+X34+V34*Y43)+R41*(U41*X41+V41*Y14)
2 +.5*CYL*(V12+V34)*XY)
PITH = .25*XY*CYL*(MUO4*RVOL(IJ)*( (U12+U34)*XY)+LAMD)
XX = HR24*(PIXY*X24-PIXX*Y24)
YY = Y24*P(IJ)
UTIL(IPJ) = UTIL(IPJ) + DTO2M1*(XX+R1*YY-PITH)
UTIL(IJP) = UTIL(IJP) - DTO2M3*(XX+R3*YY+PITH)
XX = HR13*(PIXY*X31-PIXX*Y31)
YY = Y31*P(IJ)
UTIL(IPJP) = UTIL(IPJP) + DTO2M2*(XX+R2*YY-PITH)
UTIL(IJ) = UTIL(IJ) - DTO2M4*(XX+R4*YY+PITH)
PYYMP = PIYY-P(IJ)
XX = HR24*(PYYMP*X24-PIXY*Y24)
VTIL(IPJ) = VTIL(IPJ) + DTO2M1*XX
VTIL(IJP) = VTIL(IJP) - DTO2M3*XX
XX = HR13*(PYYMP*X31-PIXY*Y31)
VTIL(IPJP) = VTIL(IPJP) + DTO2M2*XX
VTIL(IJ) = VTIL(IJ) - DTO2M4*XX
XX = .5*HR24*(U24*(X24*PIXY-Y24*PIXX)-V24*(Y24*PIXY-X24*PIYY))
Q(IPJ) = Q(IPJ) + DTO2M1*XX
Q(IJP) = Q(IJP) - DTO2M3*XX
XX = .5*HR13*(U13*(X31*PIXY-Y31*PIXX)-V13*(Y31*PIXY-X31*PIYY))
Q(IPJP) = Q(IPJP) +DTO2M2*XX
Q(IJ) = Q(IJ) -DTO2M4*XX
IJ = IPJ
1199 IJP = IPJP
C UTIL(IJ) = UTIL(IJP) = UTIL(IJP-NQIB) = UTIL(IJ-NQIB) = 0.
UTIL(IJ-NQIB) = 0.
UTIL(IJP-NQIB) = UTIL(IJ-NQIB)
UTIL(IJP) = UTIL(IJP-NQIB)
UTIL(IJ) = UTIL(IJP)
IF (J.NE.2) GO TO 1220
DO 1210 IJ=ISC2,ISCF2,NQ
1210 VTIL(IJ) = 0.
1220 IF (J.NE.JP1) GO TO 1240
DO 1230 IJP = IJPS,LJP2,NQ
1230 VTIL(IJP) = 0.
1240 CALL LOOP
1299 CONTINUE
CALL DONE
1300 CALL START
DO 1399 J=2,JP1
DO 1389 I=1,IBAR
IPJ = IJ + NQ
IPJP = IJP + NQ
ETIL(IJ) = E(IJ)+.25*(Q(IPJ)+Q(IPJP) +Q(IJP)+Q(IJ))
ROL(IJ) = RO(IJ)
RCSQ(IJ) = 1./(ASQ+GGM1*AMAX1(SIE(IJ),0.))
XX = (X(IPJP)-X(IJ)+X(IPJ)-X(IJP))**2
YY = (Y(IPJP)-Y(IJ)+Y(IJP)-Y(IPJ))**2
DELSM(IJ) = DT8*(XX+YY)/(XX*YY)
IJ = IPJ
1389 IJP = IPJP
CALL LOOP
1399 CONTINUE
CALL DONE
1500 CALL START
DO 1599 J=2,JP1
DO 1589 I=1,IBAR
IMJ = IJ-NQ
IPJ = IJ+NQ
IPJP = IJP+NQ
X1 = X(IPJ)
Y1 = Y(IPJ)
R1 = R(IPJ)
X2 = X(IPJP)
Y2 = Y(IPJP)
R2 = R(IPJP)
X3 = X(IJP)
Y3 = Y(IJP)
R3 = R(IJP)
X4 = X(IJ)
Y4 = Y(IJ)
R4 = R(IJ)
X12 = X1-X2
X23 = X2-X3
X34 = X3-X4
X41 = X4-X1
Y21 = Y2-Y1
Y32 = Y3-Y2
Y43 = Y4-Y3
Y14 = Y1-Y4
R12 = R1+R2
R23 = R2+R3
R34 = R3+R4
R41 = R4+R1
U1 = UTIL(IPJ)
U2 = UTIL(IPJP)
U3 = UTIL(IJP)
U4 = UTIL(IJ)
V1 = VTIL(IPJ)
V2 = VTIL(IPJP)
V3 = VTIL(IJP)
V4 = VTIL(IJ)
U12 = U1+U2
U23 = U2+U3
U34 = U3+U4
U41 = U4+U1
V12 = V1+V2
V23 = V2+V3
V34 = V3+V4
V41 = V4+V1
C MR = ML = MT = MB = MC =RO(IJ)/RVOL(IJ)
MC =RO(IJ)/RVOL(IJ)
MB = MC
MT = MB
ML = MT
MR = ML
C PR = PLE = PT = PB = PC = P(IJ)
PC = P(IJ)
PB = PC
PT = PB
PLE = PT
PR = PLE
IF (I.EQ.IBAR) GO TO 1510
MR = RO(IPJ)/RVOL(IPJ)
PP = P(IPJ)
1510 IF(I.EQ.1) GO TO 1520
ML = RO(IMJ)/RVOL(IMJ)
PLE = P(IMJ)
1520 IF (J.EQ.JPI ) GO TO 1530
MT = RO(IJP)/RVOL(IJP)
PT = P(IJP)
1530 IF (J.EQ.2) GO TO 1540
MB = RO(IJM)/RVOL(IJM)
PB = P(IJM)
1540 P12 = (MR*PC+MC*PR)/(MR+MC)
P23 = (MT*PC+MC*PT)/(MT+MC)
P34 = (ML*PC+MC*PLE)/(ML+MC)
P41 = (MB*PC+MC*PB)/(MB+MC)
ETIL(IJ) = ETIL(IJ)-DTO4/MC*(R12*P12*(U12*Y21+V12*X12)
1 +R23*P23*(U23*Y32+V23*X23)
2 +R34*P34*(U34*Y43+V34*X34)
3 +R41*P41*(U41*Y14+V41*X41))
IJ = IPJ
IJP = IPJP
1589 IJM = IJM+NQ
CALL LOOP
1599 CONTINUE
CALL DONE
2000 NUMIT = 0
MUSTIT = 1
PLMAX = EHID
2010 CALL START
DO 2099 J=2,JP1
DO 2089 I=1,IBAR
IPJ = IJ*NQ
IPJP = IJP*NQ
X1 = X(IPJ)
Y1 = Y(IPJ)
R1 = R(IPJ)
U1 = UL(IPJ)
V1 = VL(IPJ)
X2 = X(IPJP)
Y2 = Y(IPJP)
R2 = R(IPJP)
U2 = UL(IPJP)
V2 = VL(IPJP)
X3 = X(IJP)
Y3 = Y(IJP)
R3 = R(IJP)
U3 = UL(IJP)
V3 = VL(IJP)
X4 = X(IJ)
Y4 = Y(IJ)
R4 = R(IJ)
U4 = UL(IJ)
V4 = VL(IJ)
D = .25*RVOL(IJ)*((R1+R2)*((U1+U2)*(Y2-Y1)+(V1+V2)*(X1-X2))
1 +(R2+R3)*((U2+U3)*(Y3-Y2)+(V2+V3)*(X2-X3))
2 +(R3+R4)*((U3+U4)*(Y4-Y3)+(V3+V4)*(X3-X4))
3 +(R4+R1)*((U4+U1)*(Y1-Y4)+(V4+V1)*(X4-X1)))
S = RDT*(ROL(IJ)-RO(IJ))+ROL(IJ)*D
RA = RCSQ(IJ)*(RDT+D)+DELSM(IJ)
DP = -OM*S/RA
ROL(IJ) = ROL(IJ) + RCSQ(IJ)*DP
PLMAX = AMAX1(PLMAX,ABS(PL(IJ)))
IF (ABS(DP).LE.EPS*PLMAX) GO TO 2080
MUSTIT = 1
PL(IJ) = PL(IJ)+DP
Y24 = Y2-Y4
Y31 = Y3-Y1
XR13 = .5*(R1+R3)*(X1-X3)
XR24 = .5*(R2+R4)*(X2-X4)
XX = DTO2*DP
DTO2M1 = XX*RM(IPJ)
DTO2M2 = XX*RM(IPJP)
DTO2M3 = XX*RM(IJP)
DTO2M4 = XX*RM(IJ)
UL(IPJ) = U1*DTO2M1*R1*Y24
UL(IPJP) = U1*DTO2M2*R2*Y31
UL(IJP) = U1*DTO2M3*R3*Y24
UL(IJ) = U1*DTO2M4*R4*Y31
IF (J.EQ.JP1) GO to 2080
VL(IPJ) = V1-DTO2M1*XR24
VL(IJ) = V4-DTO2M4*XR13
2060 IF (J.EQ.JP1) GO TO 2080
VL(IPJP) = V2+DTO2M2*XR13
VL(IJP) = V3+DTO2M3*XR24
2080 IJ = IPJ
2089 IJP = IPJP
CCC UL(IJ) = UL(IJP) = UL(IJP-NQIB) = UL(IJ-NQIB) = 0.
UL(IJ-NQIB) = 0.
UL(IJP-NQIB) = UL(IJ-NQIB)
UL(IJP) = UL(IJP-NQIB)
UL(IJ) = UL(IJP)
CALL LOOP
2099 CONTINUE
CALL DONE
NUMIT = NUMIT+1
IF(MUSTIT.EQ.0) GO TO 3000
MUSTIT = 0
IF(NUMIT.LT.500) GO TO 2010
LPR = 2
PRINT 4130
3000 IF(GRDVEL.GT.1.99) GO TO 3150
CALL START
DO 3019 J=2,JP2
DO 3009 I=1,IP1
UG(IJ) = UL(IJ)*GRDVEL
VG(IJ) = VL(IJ)*GRDVEL
3009 IJ = IJ+NQ
CALL LOOP
3019 CONTINUE
CALL DONE
3100 CALL START
DO 3119 J=2,JP2
DO 3109 I=1,IP1
X(IJ) = X(IJ)+UG(IJ)*DT
Y(IJ) = Y(IJ)+VG(IJ)*DT
R(IJ) = X(IJ)*CYL+OMCYL
3109 IJ = IJ+NQ
CALL LOOP
3119 CONTINUE
CALL DONE
GO TO 3200
3150 CALL REZONE
3200 CALL START
DO 3269 J=2,JP1
DO 3259 I=1,IBAR
IMJ = IJ-NQ
IPJ = IJ+NQ
IPJP = IJP+NQ
X1 = X(IPJ)
Y1 = Y(IPJ)
R1 = R(IPJ)
X2 = X(IPJP)
Y2 = Y(IPJP)
R2 = R(IPJP)
X3 = X(IJP)
Y3 = Y(IJP)
R3 = R(IJP)
X4 = X(IJ)
Y4 = Y(IJ)
R4 = R(IJ)
UL1 = UL(IPJ)
VL1 = VL(IPJ)
UL2 = UL(IPJP)
VL2 = VL(IPJP)
UL3 = UL(IJP)
VL3 = VL(IJP)
UL4 = UL(IJ)
VL4 = VL(IJ)
UD1 = UG(IPJ) -UL1
VD1 = VG(IPJ) -VL1
UD2 = UG(IPJP) -UL2
VD2 = VG(IPJP) -VL2
UD3 = UG(IJP) -UL3
VD3 = VG(IJP) -VL3
UD4 = UG(IJ) -UL4
VD4 = VG(IJ) -VL4
X12 = X1-X2
X23 = X2-X3
X34 = X3-X4
X41 = X4-X1
Y21 = Y2-Y1
Y32 = Y3-Y2
Y43 = Y4-Y3
Y14 = Y1-Y4
Y31 = Y3-Y1
R12 = R1+R2
R23 = R2+R3
R34 = R3+R4
R41 = R4+R1
U12 = UL1+UL2
U23 = UL2+UL2
U34 = UL3+UL3
U41 = UL4+UL1
V12 = VL1+VL2
V23 = VL2+VL3
V34 = VL3+VL4
V41 = VL4+VL1
D = .25*RVOL(IJ)*(R12*(U12*Y21+V12*X12)+R23*(U23*Y32+V23*X23)
1 +R34*(U34*Y43+V34*X34)+R41*(U41*Y14+V41*X41))
CCC VOLR = VOLT = VOLC = 1./RVOL(IJ)
VOLC = 1./RVOL(IJ)
VOLT = VOLC
VOLR = VOLT
IF (I.NE.IBAR) VOLR = 1./RVOL(IPJ)
IF(J.NE.JP1 ) VOLT = 1./RVOL(IJP)
IF(I.EQ.1) GO TO 3280
FL = -FR
AL = -AR
3230 IF (J.EQ.2) GO TO 3290
FB = -FT(I)
AB = -AT(I)
3240 FR = DTO8*R12*( (UD1+UD2)*Y21 +(VD1+VD2)*X12)
AR = ADM*SIGN(1.,FT(I))+B0*4.*FR / (VDOR*VDLC)
FT(I) = DTO8*R23*((UD2+UD3)*Y32+(VD2+VD3)*X23)
AT(I) = AOM*SIGN(1.,FT(I))+B0*4*FT(I)/(VOLT+VOLC)
XX = AMAX1(ABS(FB)*ABS(FR),ABS(FT(I)),ABS(FL))
DTC = AMIN1(DTC,DTPOS*AOFAC/(XX*RVDL(IJ)+DTPOS*ABS(D)+EM10))
IF (DTCSAV.NE.DTC) IDTC = I
IF (DTCSAV.NE.DTC) JDTC = J
DTCSAV = DTC
MP(IJ) = RD(IJ)*VOLC
1 +FR *((1.-AR) *ROL(IJ)*(1.+AR) *RDL(IPJ))
2 +FT(I) *((1.-AT(I) *ROL(IJ)*(1.+AT(I))*RDL(IJP))
3 +FL *((1.-AL) *ROL(IJ)*(1.+AL) *RDL(IMJ))
4 +FR *((1.-AR) *ROL(IJ)*(1.+AB) *RDL(IJM)))
ROE = RO(IJ)*ETIL(IJ)
EP(IJ) = 1./MP(IJ)*(ROE*VOLC
1 +FR *((1.-AR) *ROE+(1.+AR) *RO(IPJ)*ETIL(IPJ))
2 +FT(I) *((1.-AT(I)) *ROE+(1.+AT(I))*RO(IJP)*ETIL(IJP))
3 +FL *((1.-AL) *ROE+(1.+AL) *RO(IMJ)*ETIL(IMJ))
4 +FB *((1.-AR) *ROR+(1.+AB) *RO(IJM)*ETIL(IJM)))
ATR = .5*(X2*Y31-X1*Y32-X3*Y21)
ABL = .5*(X1*Y43-X3*Y14+X4*Y31)
RVOL(IJ) = 1./(ATR*(R1+R2+R3)-ABL*(R1+R3+R4))
IJ = IPJ
IJP = IPJP
3259 IJM = IJM +NQ
CALL LOOP
3269 CONTINUE
CALL DONE
GO TO 3300
3280 FL = DTOB*R34*((UD3+UD4)*Y43*(VD3+VD4)*X34)
AL = AOM*SIGN(1.,FL)+B0*2.*FL*RVOL(IJ)
GO TO 3230
C ILLEGAL JUMP INTO DO LOOP
3290 FB = DTOB*R41*((UD4+UD1)*Y14+(VD4+VD1)*X41)
AB = AOM*SIGN(1.,FB)+B0*2.*FB*RVOL(IJ)
GO TO 3240
C ILLEGAL JUMP INTO DO LOOP
3300 CALL START
DO 3319 J=2,JP1
DO 3309 I=1,IBAR
RO(IJ) = MP(IJ)*RVOL(IJ)
E(IJ) = EP(IJ)
IF (J.EQ.2 ) RO(IJM) = ROL(IJM)
IF (J.EQ.JP1 ) RO(IJP) = ROL(IJP)
IF (I.EQ.IBAR) RO(IJ+NQ) = ROL(IJ+NQ)
IJM = IJM+NQ
IJP = IJP+NQ
3309 IJ = IJ + NQ
CALL LOOP
3319 CONTINUE
CALL DONE
CALL STARTD
DO 3399 JJ=2,JP2
J = JP4-JJ
DO 3389 II=1,IP1
I = IP2-II
IMJ = IJ-NQ
IMJM = IJM-NQ
XX = 0.
IF(I.NE.IP1 .AND. J.NE.2 ) XX =MP(IJM)
IF(I.NE.IP1 .AND. J.NE.JP2) XX =XX+MP(IJ)
IF(I.NE.1 .AND. J.NE.JP2) XX =XX+MP(IMJ)
IF(I.NE.1 .AND. J.NE.2 ) XX =XX+MP(IMJM)
RMP(IJ) = 4./XX
IJ = IMJ
3389 IJM = IMJM
CALL LOOPD
3399 CONTINUE
3400 CALL START
DO 3499 J=2,JP2
DO 3489 I=1,IP1
XX = RMP(IJ)/RM(IJ)
UP(IJ) = XX*UL(IJ)
VP(IJ) = XX*VL(IJ)
3489 IJ = IJ+NQ
CALL LOOP
3499 CONTINUE
CALL DONE
CALL START
DO 3699 J=2,JP1
DO 3599 I=1,IBAR
IPJ = IJ+NQ
IPJP = IJP+NQ
X1 = X(IPJ)
Y1 = Y(IPJ)
R1 = R(IPJ)
UL1 = UL(IPJ)
UG1 = UG(IPJ)
VL1 = VL(IPJ)
VG1 = VG(IPJ)
X2 = X(IPJP)
Y2 = Y(IPJP)
R2 = R(IPJP)
UL2 = UL(IPJP)
UG2 = UG(IPJP)
VL2 = VL(IPJP)
VG2 = VG(IPJP)
X3 = X(IJP)
Y3 = Y(IJP)
R3 = R(IJP)
UL3 = UL(IJP)
UG3 = UG(IJP)
VL3 = VL(IJP)
VG3 = VG(IJP)
X4 = X(IJ)
Y4 = Y(IJ)
R4 = R(IJ)
UL4 = UL(IJ)
UG4 = UG(IJ)
VL4 = VL(IJ)
VG4 = VG(IJ)
XX = DTO16*ROL(IJ)
UL13 = UL1+UL3
VL13 = VL1+VL3
UL24 = UL2+UL4
VL24 = VL2+VL4
F13 = XX*(R1+R3)*((UG1+UG3-UL13)*(Y3-Y1)*(VG1+VG3-VL13)+(X1-X3))
F24 = XX*(R2+R4)*((UG2-UG4-UL24)*(Y2-Y4)*(VG2+VG4-VL24)+(X4-X2))
FM1 = F24*RMP(IPJ)
FM2 = F13*RMP(IPJP)
FM3 = F24*RMP(IJP)
FM4 = F13*RMP(IJ)
XX = B0*4.*RVOL(IJ)/ROL(IJ)
AL13 = A0*SIGN(1.,F13)+XX*F13
AL24 = A0*SIGN(1.,F24)+XX*F24
OPAL13 = 1.+AL13
OPAL24 = 1.+AL24
OMAL13 = 1.-AL13
OMAL24 = 1.-AL24
XX = UL3*OMAL24+UL1*OPAL24
UP(IPJ) = UP(IPJ)-FM1*XX
UP(IJP) = UP(IJP)+FM3*XX
XX = UL4*OMAL13+UL2*OPAL24
UP(IPJP) = UP(IPJP)-FM2*XX
UP(IJ) = UP(IJ)+FM4*XX
XX = VL3*OMAL24+VL1*OPAL24
VP(IPJ) = VP(IPJ)-FM1*XX
VP(IJP) = VP(IJP) +FM3*XX
XX = VL4*OMAL13-VL2*OPAL13
VP(IPJP) = VP(IPJP)-FM2*XX
VP(IJ) = VP(IJ) +FM4*XX
IJ = IPJ
3599 IJP = IPJP
C UP(IJ) = UP(IJP) = UP(IJP*NQIB) = UP(IJ-NQIB) = 0.
UP(IJ-NQIB) = 0.
UP(IJP*NQIB) = UP(IJ-NQIB)
UP(IJP) = UP(IJP*NQIB)
UP(IJ) = UP(IJP)
IF (J.NE. 2) GO TO 3620
DO 3610 IJ=ISC2,ISCF2,NQ
3610 VP(IJ) = 0.
3620 IF(J.NE.JP1) GO TO 3640
DO 3630 IJP = IJPS,LJP2,NQ
3630 VP(IJP) = 0.
3640 CALL LOOP
3699 CONTINUE
CALL DONE
3700 CALL START
DO 3719 J=2,JP2
DO 3709 I=1,IP1
U(IJ) = UP(IJ)
V(IJ) = VP(IJ)
RM(IJ) = RMP(IJ)
3709 IN = IJ + NQ
CALL LOOP
3719 CONTINUE
CALL DONE
3800 CALL START
DO 3899 J=2,JP1
DO 3889 I=1,IBAR
IPJ = IJ+NQ
IPJP = IJP+NQ
SIE(IJ) = E(IJ)-.125*(U(IPJ)**2+U(IPJP)**2+U(IJP)**2+U(IJ)**2
1 +V(IPJ)**2+V(IPJP)**2+V(IJP)**2+V(IJ)**2)
IJP = IPJP
3889 IJ = IPJ
CALL LOOP
3899 CONTINUE
CALL DONE
IF(NPT.GT.0) CALL PARTMOV
GO TO 180
C IMPOSSIBLE TO READ THESE FORMAT STATEMENTS IN LISTING
4000 FORMAT()
4010 FORMAT()
4020 FORMAT()
4030 FORMAT()
4040 FORMAT()
4050 FORMAT()
4060 FORMAT()
4070 FORMAT()
4080 FORMAT()
4090 FORMAT()
4100 FORMAT()
4110 FORMAT()
4120 FORMAT()
4130 FORMAT()
4140 FORMAT()
4150 FORMAT()
END
C
SUBROUTINE PARTMOV
C MOVES AND PLOTS PARTICLES
COMMON/YSC1/ AASC(4242)
COMMON/YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
EQUIVALENCE (AASC(1),X,XPAR),(AASC(2),R,YPAR),(AASC(3),Y,MPAR), EQVREAL
1 (AASC(4),U,UG,DELSM),(AASC(5),V,VG),(AASC(6),RO), EQVREAL
2 (AASC(7),SIE,MP,RMP,RCSQ),(AASC(8),E,ETIL), EQVREAL
3 (AASC(9),RVOL),(AASC(10),M,RM,VP),(AASC(11),P,PL,EP, EQVREAL
4 UP,PMD),(AASC(12),UTIL,UL,CQ,PMX,PU),(AASC(13),VTIL, EQVREAL
5 VL,PMY,PV),(AASC(14),Q,ROL) EQVREAL
REAL LAM,LAMD,M,M6,MC,ML,MP,MPAR,MR,MT,MTE,MU,MUO2,MUO4 EQVRWAL
DIMENSION X(1),XPAR(1),R(1),YPAR(1),Y(1),MPAR(1),U(1),UG(1), DIMEN
1 DELSM(1),V(1),VG(1),RO(1),SIE(1),MP(1),RMP(1),RCSQ(1), DIMEN
2 E(1),ETIL(1),RVOL(1),M(1),RM(1),VP(1),P(1),PL(1),EP(1), DIMEN
3 UP(1),UTIL(1),UL(1),CQ(1),PMX(1),PU(1),VTIL(1),VL(1), DIMEN
4 PMY(1),PV(1),Q(1),ROL(1),PMD(1)
DIMENSION X1(4)
EQUIVALENCE (X1(2),X2),(X1(3),X3),(X1(4),X4)
COMMON /YSC3/ JNM
CALL START
DO 1019 J=2,JP2
DO 1009 I=1,IP1
CCC PMX(IJ) = PMY(IJ) = PMD(IJ) = 0.0
PMD(IJ) = 0.0
PMY(IJ) = PMD(IJ)
PMX(IJ) = PMY(IJ)
1009 IJ = IJ+NQ
CALL LOOP
1019 CONTINUE
CALL DONE
DO 1099 J=2,JP2
IEC = (J-1)*NQI
CALL ECRD (AASC,IEC,NQI,NE)
IJ = 1
DO 1089 I=1,IP1
IF(X(IJ).GT.PXRP .OR. Y(IJ).LT.PYBM) GO TO 1099
IF(Y(IJ).GT.PYTP) GO TO 1100
KI = X(IJ) *RPDR + OPEM10
KJ = (Y(IJ)-PYB)*RPDZ + OPEM10
IEC = KJ*NQI
CALL ECRD (AASC(ISC2),IEC,NQI2,NE)
KIJ = (KI-1)*NQ+ISC2
KIJP = KIJ+NQI
W = X(IJ) -FLOAT(KI-1)*PDR
H = (Y(IJ)-PYB)-FLOAT(KJ-1)*PDZ
XX = RPDRDZ/RM(IJ)
HTE = PDZ-H
DO 1049 JJ=1.2
WTE = PDR-W
DO 1039 II=1,2
X1 = XX*WTE*HTE
PMX(KIJ) = PMX(KIJ) +XI*U(IJ)
PMY(KIJ) = PMY(KIJ) +XI*V(IJ)
PMD(KIJ) = PMD(KIJ) +XI
WTE = W
1039 KIJ = KIJ+NQ
HTE = H
1049 KIJ = KIJP
CALL ECWR (AASC(ISC2),IEC,NQI2,NE)
1089 IJ = IJ+NQ
1099 CONTINUE
1100 CALL START
KK = 0
DO 1199 J=2,JBP2
DO 1189 I=1,IBP1
RPM0 = PM0(IJ)C IF (RPM0.NE.0.) PMO(IJ) = RPM0 = 1./RPM0
IF (RPM0.NE.0.) RPM0 = 1./RPM0
IF (RPM0.NE.0.) PMD(IJ) = RPM0
IF (RPM0.EQ.0.) KK = 1
PU(IJ) = PMX(IJ)*RPM0
PV(IJ) = PMY(IJ)*RPM0
IF (I.EQ.1) PU(IJ) = 0.
1189 IJ = IJ + NQ
CALL LOOP
1199 CONTINUE
CALL DONE
1200 IF (KK.EQ.0) GO TO 1400
CALL START
DO 1299 J=2,JBP2
DO 1289 I=1,IBP1
IF(PMQ(IJ).GT.0.) GO TO 1289
C ITESTL = ITESTR = I
ITESTR = I
ITESTL = ITESTR
C PUL = PVL = PUR = PVR = 0.
PVR = 0.
PUR = PVR
PVL = PUR
PUL = PVL
XWL = XWR + 1.
C IL = IR =IJ
IR =IJ
IL = IR
1210 ITESTL = ITESTL - 1
IF(ITESTL.LT.1) GO TO 1230
IL = IL-NQ
IF (PMO(IL).GT.0.) GO TO 1220
XWL = XWL + 1.
GO TO 1210
1220 PUL = PU(IL)
PVL = PV(IL)
1230 ITESTR = ITESTR+1
IF(ITESTR.GT.IBP1) GO TO 1250
IR = IR+NQ
IF (PM0(IR).GT.0.) GO TO 1240
XWR = XWR + 1.
GO TO 1230
1240 PUR = PU(IR)
PVR = PV(IR)
1250 DEN = 1./(XWL+XWR)
WTL = XWR*DEN
WTR = XWL*DEN
PU(IJ) = PUL*WTL + PUR*WTR
PV(IJ) = PVL*WTL + PVR*WTR
1289 IJ = IJ+NQ
CALL LOOP
1299 CONTINUE
CALL DONE
1300 DO 1399 J=2,JBP2
IEC = (J-1)*NQI
CALL ECRD (AASC,IEC,NQI,NE)
IJ = 1
DO 1389 I=1,IBP1
IF(PMO(IJ).GT.0. .OR. PU(IJ).NE.0. .OR. PV(IJ).NE.0.) GO TO 1389
C JTESTB = JTESTT = J
TESTT = J
JTESTB = JTESTT
C PUB = PVB = PUT = PVT = 0.
PVT = 0.
PUT = PVT
PVB = PUT
PUB = PVB
C XWB = XWT = 1.
XWT = 1.
XWB = XWT
C IB = IT = IJ*NQI
IT = IJ*NQI
IB = IT
C IECB = IECT = IEC
IECT = IEC
IECB = IECT
1310 JTESTB = JTESTR B - 1
IF (JTESTB.LT.2) GO TO 1330
IECB = IECB - NQI
CALL ECRD (AASC(ISC2),IECB,NQI,NE)
IF(PMO(IB).GT.0.) GO TO 1320
XWB = XWB + 1.
GO TO 1310
1320 PUB = PU(IB)
PVB = PV(IB)
1330 JTESTT = JTESTT + 1
IF(JTESTT .GT. JBP2) GO TO 1350
IECT = IECT + NQI
CALL ECRD (AASC(ISC2),IECT,NQI,NE)
IF(PMO(IT).GT.0.) GO TO 1340
XWT = XWT + 1.
GO TO 1330
1340 PUT = PU(IT)
PVT = PV(IT)
1350 DEN = 1./(XWT+XWB)
WTB = XWT*DEN
WTT = XWB*DEN
PU(IJ) = PUB*WTB + PUT*WTT
PV(IJ) = PVB*WTB + PVT*WTT
1389 IJ = IJ + NQ
1399 CALL ECWR (AASC,IEC,NQI,NE)
1400 IF(IMOMX.EQ.0) GO TO 2000
CALL START
DO 1499 J=2,JBP2
DO 1489 I=1,IBP1
PU(IJ) = PU(IJ).AND..NOT. 7777777777B
PV(IJ) = PV(IJ).AND..NOT. 7777777777B
1489 IJ = IJ + NQ
CALL LOOP
1499 CONTINUE
CALL DONE
2000 IECP = IPAR
CC NPMT = JOLD = 0
JOLD = 0
NPMT = JOLD
2010 CALL ECRD (AASC,IECP,LPB,NE)
KP = 1
2020 XTE = XPAR(KP)
YTE = YPAR(KP)
DRAG = MPAR(KP)
I = XTE *RPDR + OPEM10
J = (YTE-PYB)*RPDZ + OPEM10
IF (J.LT.1 .OR. J.GT.JBP .OR. I.LT.1 .OR. I.GT. IBP) GO TO 2150
IF(J.EQ.JOLD) GO TO 2030
JOLD = J + IMOME3
IEC = J*NQI
CALL ECRD(AASC(ISC2),IEC,NQ12,NE)
2030 IJ = (I-1)*NQ+ISC2
IPJ = IJ+NQ
IJP = IJ*NQI
IPJP = IPJ+NQI
W = XTE -FLOAT(I-1)*PDR
H = (YTE-PYB) -FLOAT(J-1)*PDZ
PDRMW = PDR-W
PDZMH = PDZ-H
X1 = PDRMW*PDZMH
X2 = W*PDZMH
X3 = PDRMW* H
X4 = W* H
UK = (PU(IJ)*X1 + PU(IPJ)*X2 + PU(IJP)*X3 + PU(IPJP)*X4) * RPDRDZ
VK = (PV(IJ)*X1 + PV(IPJ)*X2 + PV(IJP)*X3 + PV(IPJP)*X4) * RPDRDZ
UBAR = SHIFT(XTE,30)
VPAR = SHIFT(YTE,30)
DTDRAG = DT*DRAG
RDTDRG = 1./(1.+DTDRAG)
C URAND = VRAND = 0.0
VRAND = 0.0
URAND = VRAND
UPAR = (UPAR + DTDRAG*(UK+URAND) + DTGR) *RDTDRG
VPAR = (VPAR + DTDRAG*(VK+VRAND) + DTGZP) *RDTDRG
XTE = XTE + DT*UPAR
YTE = YTE + DT*VPAR
USH = SHIFT(UPAR,30).AND.7777777777B
VSH = SHIFT(VPAR,30).AND.7777777777B
XPAR(KP) = (XTE.AND. .NOT. 7777777777B) .OR. USH
YPAR(KP) = (YTE.AND. .NOT. 7777777777B) .OR. VSH
IF(IMOMX.EQ.0) GO TO 2150
MTE = SHIFT(DRAG,30)
H = MTE*RPDRDZ*DTDRAG
KIJ = IJ
KK = 0
DO 2149 JJ=1,2
DO 2139 II=1,2
KK = KK + 1
X1 = XI(KK)*H
PUL = PU(KIJ)
PVL = PV(KIJ)
XX = (SHIFT(PUL,30))*X1*(UK-UPAR+URAND)
YY = (SHIFT(PVL,30))*X1*(VK-VPAR+VRAND)
PU(KIJ) = (PUL.AND..NOT.7777777777B).OR.
1 (SHIFT(XX,30).AND.7777777777B)
PV(KIJ) = (PVL.AND..NOT.7777777777B).OR.
1 (SHIFT(YY,30).AND.7777777777B)
2139 KIJ = KIJ+NQ
2149 KIJ = IJP
CALL ECWR(AASC(ISC2),IEC,NQ12,NE)
2150 NPMT = NPMT+1
KP = KP+3
IF (KP .LT.LPB) GO TO 2020
CALL ECWR(AASC,IECP,LPB,NE)
IECP = IECP + LPB
GO TO 2010
2180 XPAR(KP) = -1.E+3
GO TO 2150
2190 CALL ECWR(AASC,IECP,LPB,NE)
2500 IF(IMOMX.EQ.0) RETURN
DO 2599 J=2,JP2
IEC = (J-1)*NQI
CALL ECRD (AASC,IEC,NQI,NE)
IJ = 1
DO 2589 I=1,IPI
IF (X(IJ).GT.PXRP .OR. Y(IJ).LT.PYBM) GO TO 2599
IF (Y(IJ).GT.PYTP) RETURN
KI = X(IJ) *RPDR * OPEM10
KJ = (Y(IJ)-PYB)*RPDZ * OPEM10
IECP = KJ*NQI
CALL ECRD (AASC(ISC2),IECP,NQI2,NE)
KIJ = (KI-1)*NQ+ISC2
KIJP = KIJ+NQ1
W = X(IJ) * FLOAT(KI-1)*PDR
H = (Y(IJ)-PYB) * FLOAT(KJ-1)*PDZ
HTE = PDZ-H
DO 2549 JJ=1,2
WTE = PDR-W
DO 2539 II=1,2
XX = RPDRDZ*WTE*HTE*PMD(KIJ)
U(IJ) = U(IJ) -XX*SHIFT(PU(KIJ),30))
V(IJ) = V(IJ) -XX*SHIFT(PV(KIJ),30))
WTE = W
2539 KIJ = KIJ+NQ
HTE = H
2549 KIJ = KIJP
2589 IJ = IJ+NQ
2599 CALL ECWR(AASC,IEC,NQI,NE)
RETURN
ENTRY PARPLOT
CALL ADV(1)
CALL FRAME (IPXL,IPXR,IPYT,IPYB)
CALL FRAME (IPXL,IPXR,IPYT,IPYB)
IF(LPR.EQ.0) GO TO 3000
CALL LINCNT(59)
WRITE(12,3090) PDR,PDZ,PXR,PYB,PYT
WRITE(12,3095) JNM,NAME,T,NCYC
3000 IECP = IPAR
IF(ICOLOR.GT.0) CALL COLOR(1)
IF(ICOLOR.GT.0) CALL COLOR(1)
NPPT = 0
3010 CALL ECRD (AASC,IECP,LPB,NE)
KP = 1
3020 IF (XPAR(KP).LT.0) GO TO 3050
IX1 = FIPXL + (XPAR(KP)-PXL)*PXCONV
IY1 = FIPYB + (YPAR(KP)-PYR)*PYCONV
CALL PLT (IX1,IY1,42)
3050 NPPT = NPPT + 1
IF (NPPT.EQ.NPT) GO TO 3060
KP = KP+3
IF(KP.LT.LPB) GO TO 3020
IECP = IECP+LPB
GO TO 3010
3060 IF(ICOLOR.GT.0) CALL COLOR(0)
RETURN
C
3090 FORMAT(' PARTICLES',/11X,'PDR=',1PE12.5,' PDZ=',E12.5,' PXR=',
1E12.5, ' PYB=',E12.5,' PYY=',E12.5)
3095 FORMAT(5X,A10,10A8,' T=',1PE12.5,' CYCLE', I5)
END
C
SUBROUTINE REZONE
C CALCULATES NEW GRID VERTEX VELOCITIES UG AND VG
C ROLL-YOUR-OWN SECTION
C THIS VERSION: WAY TO FOLLOW RISE OF DEBRIS FROM ATMOSPHERIC BURST
COMMON/YSC1/ AASC(4242)
COMMON/YSC2/ AA(1),ANC,ASQ,A0,A0FAC,A0M,B0,COLAMU,CYL, COMMON2
1 DR,DT,DTC,DTFAC,DTGR,DTGZ,DTGZP,DTO(10),DTOC(10), COMMON2
2 DTO16,DTO2,DTO4,DTO8,DTPOS,DTV,DT8,DZ,EM10,EPS,FIBP, COMMON2
3 FIPXL,FIPXR,FIPYB,FIPYT,FIXL,FIXR,FIYB,FIYT,FJBP, COMMON2
4 FREZ,GGM1,GM1,GR,GRDVEL,GZ,GZP,I,IBAR,IBP,IBP1,ICOLOR, COMMON2
5 IDTO,IJ,IJM,IJP,IJPS,IMOME3,IMOMX,IM1,IM6, COMMON2
6 IPAR,IPXL,IPXR,IPYB,IPYT,IP1,IP2,ISCF1,ISCF2,ISC2,ISC3, COMMON2
7 ITV,IUNF,IXL,IXR,IYB,IYT,J,JBAR,JBP,JBP2,JCEN,JM10, COMMON2
8 JM14,JP1,JP2,JP4,JP4O2,JUNF,JUNFO2,KXI,LAM,LJP2,LPB, COMMON2
9 LPR,MU,NAME(10),NCYC,NLC,NPS,NPT,NQ,NQI,NQIB,NQI2,NSC, COMMON2
1 NUMIT,NUMTD,OM,OMANC,OMCYL,OMEM10,OPEM10,PDR,PDZ,PXCONV, COMMON2
2 PXL,PXR,PXRP,PYB,PYBM,PYCONV,PYT,PYTP,RDT,REZRON,REZSIE, COMMON2
3 REZUE,REZVE,REZVT,REZY0,RIBAR,RIBJB,RIBP,RJBP,ROMFR, COMMON2
4 RON,RPDR,RPDRDZ,RPDZ,T,THIRD,TLIMD,TOUT,TWFIN,T20MD, COMMON2
5 VV,XCONV,XL,XR,YB,YCONV,YT,ZZ COMMON2
EQUIVALENCE (AASC(1),X,XPAR),(AASC(2),R,YPAR),(AASC(3),Y,MPAR), EQVREAL
1 (AASC(4),U,UG,DELSM),(AASC(5),V,VG),(AASC(6),RO), EQVREAL
2 (AASC(7),SIE,MP,RMP,RCSQ),(AASC(8),E,ETIL), EQVREAL
3 (AASC(9),RVOL),(AASC(10),M,RM,VP),(AASC(11),p,PL,EP, EQVREAL
4 UP,PM0),(AASC(12),UTIL,UL,CQ,PMX,PU),(AASC(13),VTIL, EQVREAL
5 VL,PMY,PV),(AASC(14),Q,ROL) EQVREAL
REAL LAM,LAMD,M,M6,MC,ML,MP,MPAR,MR,MT,MTE,MU,MUO2,MUO4 EQVRWAL
DIMENSION X(1),XPAR(1),R(1),YPAR(1),Y(1),MPAR(1),U(1),UG(1), DIMEN
1 DELSM(1),V(1),VG(1),RO(1),SIE(1),MP(1),RMP(1),RCSQ(1), DIMEN
2 E(1),ETIL(1),RVOL(1),M(1),RM(1),VP(1),P(1),PL(1),EP(1), DIMEN
3 UP(1),UTIL(1),UL(1),CQ(1),PMX(1),PU(1),VTIL(1),VL(1), DIMEN
4 PMY(1),PV(1),Q(1),ROL(1),PM0(1)
DIMENSION UGTE(100)
REZOMG = 0.15*RDT
REZBTA = 0.002
XX = 1.E+6
CALL START
C FCR = FCT = FCB = XXX = XOMSUM = YOMSUM = OMSUM = 0.
OMSUM = 0.
YOMSUM = OMSUM
XOMSUM = YOMSUM
XXX = XOMSUM
FCB = XXX
FCT = FCB
FCR = FCT
C
DO 1049 J=2,JP2
DO 1039 I=1,IP1
AVEL = AMAX1(ABS(UL(IJ)),ABS(VL(IJ)))
XXX = AMAX1(XXX,AVEL)
IF (I.EQ.IM6 ) FCR = AMAX1(FCR,AVEL)
IF (J.EQ.JM14) FCT = AMAX1(FCT,AVEL)
IF (J.EQ.9 ) FCB = AMAX1(FCB,AVEL)
1039 IJ = IJ + NQ
CALL LOOP
1049 CONTINUE
FCR = SQRT(FCR*XXX)
FCT = SQRT(FCT*XXX)
FCB = SQRT(FCB*XXX)
CALL START
DO 1069 J=2,JP1
IF (J.EQ.JP4O2) YCEN = Y(IJ)
DO 1059 I=1,IBAR
IPJ = IJ+NQ
IPJP = IJP+NQ
IF(J.LT.10.OR.J.GT.JM10.OR.I.GT.IM6) GO TO 1055
X1 = X(IPJ)
Y1 = Y(IPJ)
U1 = UL(IPJ)
V1 = VL(IPJ)
X2 = X(IPJP)
Y2 = Y(IPJP)
U2 = UL(IPJP)
V2 = VL(IPJP)
X3 = X(IJP)
Y3 = Y(IJP)
U3 = UL(IJP)
V3 = VL(IJP)
X4 = X(IJ)
Y4 = Y(IJ)
U4 = UL(IJ)
V4 = VL(IJ)
R1 = .125*RVOL(IJ)*(R(IPJ)+R(IPJP)+R(IJP)+R(IJ))
YY = R1*((U1+U4)*(X1-X4)+(V1+V4)*(Y1-Y4)
1 +(U2+U1)*(X2-X1)+(V2+V1)*(Y2-Y1)
2 +(U3+U2)*(X3-X2)+(V3+V2)*(Y3-Y2)
3 +(U4+U3)*(X4-X3)+(V4+V3)*(Y4-Y3))
IF (YY.GT.0.) GO TO 1055
YY = YY*YY
XOMSUM = XOMSUM + YY*Y4
YOMSUM = YOMSUM + YY*Y4
OMSUM = OMSUM + YY
1055 IJ = IPJ
1059 IJP = IPJP
CALL LOOP
1069 CONTINUE
XC = XOMSUM/OMSUM
YC = YOMSUM/OMSUM
REZVT = AMAX1(0.,(REZOMG*.5*(YC-YCEN)))
CALL START
DO 1099 J=2,JP2
YBAR = .5*(Y(IJP)+Y(IJM))*REZBTA*(YC-Y(IJ))
VGTE = REZOMG*(YBAR-Y(IJ))*REZVT
IF (J.EQ.2) VGTE = -FCB
IF(J.EQ.JP2) VGTE = FCT
DO 1089 I=1,IP1
IPJ = IJ+NQ
IF (J.GT.2) GO TO 1080
XBAR = .5*(X(IPJ)+X(IJ-NQ))+REZBTA*(XC-X(IJ))
UGTE(I) = REZOMG*(XPAR-X(IJ))
IF (I.EQ.1) UGTE(I) = 0.0
IF (I.EQ.IP1) UGTE(I) = FCR
1080 UG(IJ) = UGTE(I)
VG(IJ) = VGTE
1089 IJ = IPJ
CALL LOOP
1099 CONTINUE
CALL DONE
1200 CALL START
CC XIP1 = YJP2 = 0.
YJP2 = 0.
XIP1 = YJP2
Y2 = 1.E+20
DO 1289 J=2,JP2
DO 1279 I=1,IP1
X(IJ) = X(IJ)+UG(IJ)*DT
Y(IJ) = Y(IJ)+VG(IJ)*DT
R(IJ) = X(IJ)*CYL*OMCYL
IF (J.EQ.2) Y2 = AMIN1(Y(IJ),Y2)
IF (J.EQ.JP2) YJP2 = AMAX1(Y(IJ),YJP2)
IF (I.EQ.IP1) XIP1 = AMAX1(X(IJ),XJP1)
1279 IJ = IJ+NQ
CALL LOOP
1289 CONTINUE
CALL DONE
PDR = XIP1*RIBP
PD2 = (YJP2-Y2)*RJBP
PYB = 0.0
CALL FILMCO
CALL START
YY = ABS(GZ)/(GM1*REZSIE)
DO 1399 J=2,JP1
DO 1389 I=1,IBAR
IPJ = IJ + NQ
IPJP = IJP+NQ
Y4 = REZY0 + 0.25*(Y(IJP)+Y(IPJP)+Y(IJ)+Y(IPJ))
IF (J.EQ.2) ROL(IJM) = REZRON*EXP((Y4-Y(IJ)-Y(IPJ))*YY)
IF (I.EQ.IBAR) ROL(IPJ) = REZRON*EXP((Y4-Y(IPJ)-Y(IPJP))*YY)
IF (J.EQ.JP1 ) ROL(IJP) = REZRON*EXP((Y4-Y(IJP)-Y(IPJP))*YY)
IJM = IJM + NQ
IJP = IJP + NQ
1389 IJ = IPJ
CALL LOOP
1399 CONTINUE
CALL DONE
RETURN
END