This report presents a detailed description of the Marker-And-Cell (MAC) technique developed at Los Alamos for solving fluid flow problems. The method is appropriate for use with a high-speed digital computer.
The fluid is incompressible, viscous, and moves through large-amplitude contortions in several space dimensions. There may be a free surface upon which waves can form and break, or the flow may be entirely confined by walls. Development of the fluid configuration can be investigated through as much elapsed time as desired. The motions are calculated by using the complete Navier-Stokes equations, including all nonlinear terms. The only approximations arise from the finite-difference representation. Pressure and velocity are used directly as the dependent variables; neither stream function nor vorticity enters except as results derived from the velocity field.
Part 1 describes the basic method, its theoretical properties, and some extensions. Part 2 presents a detailed description of the logic of the present code. In Part 3, there are numerous calculational examples to show the wide scope of applicability. The extension of the method to computations involving two fluids is described in Part 4.
We are grateful to Wolfgang Rasmussen for his analysis of the bore calculations described in Part 1, and to Jacob Fromm for details of an iterative technique to solve the finite-difference Poisson's equation.
Rewriting the text below to HTML presents problems due to the lack of support for mathematical symbols having both subscripts and superscripts immediately to the right of the symbol. Also MATHML is poorly supported in modern browsers. One solution is to make changes to the text via Javascript. A second is to use SVG for the mathematical text, which is what has been used in this document. It does mean small adjustments to the text to ensure that all the mathematics is in SVG and none embedded in the HTML.
The numerical solution of problems in fluid dynamics has been greatly facilitated in the last decade by the development of high-speed, large-memory electronic computers. To take advantage of these facilities, this hardware development has been accompanied by the invention of numerical analysis techniques.
Anyone experienced in the field will confirm that the difficulties involved in formulating useable methods for complicated problems are formidable. It is suspected that, at the present time, our technique capabilities for the existing machines are far less sophisticated than those that eventually will be developed. Even so, each new approach is broadening the scope of allowable investigations to such an extent that the term computer revolution is beginning to be a realistic description of the present situation.
Through necessity, the greatest advances in fluid dynamics computing have come in the field of compressible (high-speed) flows, in which shocks and rarefactions are prominent features. Particularly spectacular have been the results for large-distortion, time-dependent processes in several space dimensions. Many investigators have contributed to the development of techniques and to their successful applications.
The subject of incompressible (low-speed) flows has only recently received the same intensive effort towards the development of techniques for transient-flow solutions. The method developed by Fromm has been particularly broad in its applicability to confined, incompressible flows [1].
Recently, we proposed a different approach, by which it is possible to solve for incompressible flows with a free surface, as well as for flows that are confined [2][3]. The original description of this method was somewhat sketchy in its presentation of details, and a subsequent paper emphasized only the applications [4]. Therefore, this report is intended to present a detailed description of our method for computing time-dependent, viscous, incompressible fluid flows in several space dimensions, a method we have termed the Marker-And-Cell (MAC) technique.
This part of the report discusses the theoretical background, together with a description of the basic technique, its properties, and some of its allowable modifications. The method of solution and the details of computer logic are given in Part 2. In Part 3 we have gathered a potpourri of results from actual calculations in order to illustrate the applicability of the method. Part 4 shows the technique extended to calculations of two different fluids in contact with each other.
Before discussing details of the MAC technique, it is useful to record the form of the differential equations from which the difference equations are derived, and to discuss some properties of the equations that are pertinent to the numerical manipulations employed in the computer.
For the basic field variables we use the symbols
The properties of the material are represented by the kinematic viscosity coefficient, which is taken to be constant at this stage of the discussion. Finally, it is convenient to define a discrepancy term:
Equation (1.2) expresses the conservation of mass for an incompressible fluid. If D fails to vanish anywhere, there is a discrepancy in this conservation requirement - hence the name for this term.
Equation (1.3) describes the local production of momentum in terms of the following sources:
For the present, we neglect the effects of density layering, such as would occur in the atmosphere of the earth. Also neglected are the effects of temperature variations that can produce buoyancy (free-convection problems) as well as variations in the kinematic viscosity. Eddy viscosity effects, whose importance can be great, are likewise neglected; ν is considered to represent molecular viscosity only, and the flows are imagined to be non-turbulent.
Numerous methods have been used for the analytical solution of these equations. In many cases, the approach has been tailored to take advantage of any helpful features of the particular initial or boundary conditions. For example, it has often been useful to linearize the equations (slow-motion and low-amplitude approximations), or to neglect the effects of viscosity, or to impose restrictions of periodicity in space and/or time. Our solution technique uses none of these assumptions. It is based upon the full Navier-Stokes equations. The only approximations are related to:
All three of these matters are considered in detail in this report.
One of the most common solution techniques, both analytical and numerical, involves the introduction of a stream function, through which Eq (1.2) becomes identically satisfied. In addition, it is often useful to introduce vorticity, thereby eliminating pressure as a variable. These techniques are particularly useful for confined flows in two space dimensions. For flows with a free surface, boundary condition difficulties arise with these substitute variables; likewise, if extension to cylindrical coordinates is desirable, stream function and vorticity become somewhat awkward to handle. The MAC technique, therefore, uses pressure and velocity directly as the primary variables; and no difficulties arise with respect to satisfying Eq (1.2).
Equation (1.3) can be written in two modified forms, both useful for the numerical method, and both derived with the help of Eq. (1.2):
The reason for the change in the transport term is that the finite-difference form of this equation retains rigorous momentum conservation, while the finite-difference form of the transport term in Eq (1.3) does not. This conservative property can be exhibited directly in the differential equation. Thus, integrating Eq (1.4) over a fixed volume V, whose surface is S, we get:
This equation shows that the only contributions to the momentum within the volume come from the body force in the interior and/or from fluxes through the surface. Equation (1.2) would not have allowed the transformation of the transport term to a surface integral. The situation in finite difference form is directly analogous. The form of Eq (1.4) is perfectly appropriate for problems in Cartesian coordinates; it has, in fact, the advantage over Eq (1.5) of relating shear stresses directly to the velocity differences producing them. The derivations of technique and the applications described in this report for Cartesian coordinates are all based on Eq (1.4).
For other coordinate systems, the proper starting point is Eq (1.5). We shall show in the discussion on cylindrical coordinate applications that there actually is an unexpected advantage in using Eq (1.5) that would be realized even in Cartesian-coordinate applications. This has been discovered only recently, so that all the tests reported in this discussion were based upon Eq (1.4). Further discussion of this advantage is also included in the section on cylindrical coordinates.
It should be mentioned, incidentally, that conservation of energy can be demonstrated from the differential equations. In the absence of viscous or body forces, kinetic energy is separately conserved. (The body force contributes a coupling to the potential energy, while the viscous force couples to the internal, or heat, energy.) We have not been able to extend these energy conservation properties rigorously to the finite-difference approximation; but, fortunately, the discrepancies can easily be made negligible.
The first term on the right is a combination of convective and work terms. The second term exhibits the coupling with the potential energy of the fluid, through work done by or against the body force. The third term gives the dissipation from kinetic to heat energy; it can be shown that the internal contribution from this term is negative definite. To see this, it is necessary to evoke the identity
The internal contribution is, thus, obviously negative definite; the other term contributes to viscous work through the surface of the volume.
For any specific problem, it is necessary to supply an appropriate set of initial and boundary conditions. We shall be concerned particularly with a prescribed set of rigid walls that may be no-slip or free-slip, and with inflow and outflow boundaries. The rigid walls may partially confine the fluid, or they may define an obstacle about which the fluid flows. Inflow boundaries have prescribed conditions of fluid influx through them, while outflow boundaries are arranged in such a way that fluid outflux through them will occur with minimal disturbance to the fluid remaining in the calculation region.
In addition to the prescribed boundary specifications, there will be boundary conditions to apply at the free surface, whose position varies with time in a manner not a priori known. One of the major features of the MAC technique is the method by which the free-surface boundary conditions are applied.
The rigid-wall boundary conditions are the simplest to derive; they follow directly from the momentum equations. For a free-slip wall, the normal velocity component must vanish; for a no-slip wall, the tangential components must, in addition, vanish. Corresponding boundary conditions on pressure, obtained through Eq (1.3), relate the normal derivative of φ to the body and viscous forces in a straightforward manner. These differential boundary conditions need not be written in detail here, as, for the numerical calculations, it is necessary to derive the finite-difference analogies to the boundary conditions directly from the finite-difference momentum equations.
Conditions along an inflow boundary are similarly derived; the only difference is that the velocity components are prescribed in some arbitrary manner, rather than forced to vanish. Pressure boundary conditions then follow from Eq (1.3) in such a way as to again insure consistency with the momentum balance.
An outflow boundary condition, in contrast, is very difficult to derive because there are no unique criteria to aid in its formulation. There is some consolation in the fact that if the outflow velocity exceeds the wave speed, then it should not matter how the boundary is treated. (If, however, there is a viscous boundary layer along the approach to the outflow line, then even this consolation may disappear.) Any technique that one applies at the outflow boundary must be amenable to a physical interpretation. Thus, for example, if the tangential velocity components are forced to vanish, and the normal derivative of the normal velocity component is likewise set to zero, then the model presumably represents an outflow into numerous frictionless pipes normal to the boundary. If the normal component is, instead, prescribed as a function of position and time, then there is a representation of an outflow pump at each of the tubes whose flow rate is controlled.
Boundary conditions at the free surface remain the most interesting of them all. For the numerical method, they are the most difficult conditions to apply accurately. It is not difficult to state the principles that form the basis for the free-surface boundary conditions:
The second principle implies a slight generalization from a strictly free boundary. It allows for pressure-driven surface motions to be produced by an atmosphere whose inertial contribution to the dynamics is negligible. The first principle implies, however, that such an atmosphere be incapable of exerting a shear stress. If the latter is actually of importance, or if Bernoulli pressure variations in the atmosphere are a significant effect of its dynamics, then the one-fluid, free-surface treatment breaks down; and a two-fluid calculation is required.
To exhibit the complications arising from the free-surface boundary conditions, it is useful to record them for the special case of plane motions in two space dimensions. Let nx and ny be the components of a unit outward vector, normal to the surface. Then, if z(x,t) is the equation for the surface height,
Thus, it is clear that accurate application of these conditions requires accurate knowledge of the free-surface orientation. The difficulties of this requirement will become apparent when the numerical method is presented.
It should be mentioned, incidentally, that inclusion of surface tension into the free-surface boundary conditions could be useful for some types of problems. If T the surface tension (stress per unit length) and R1, R2 are the principal radii of curvature, then Eq (1.10) must be modified by a term
At present, we have not found an effective way to incorporate this effect into the numerical calculations.
To illustrate the MAC method for solving free-surface problems, we confine our attention in this section to two-dimensional motions in a plane. Extensions to cylindrical and other coordinates are discussed in a subsequent section.
In the development of any computing method for fluid dynamics problems, there are two interacting considerations that must be taken into account:
Many representations can be visualized for calculating the flow of an incompressible fluid with a free surface. Several approaches were investigated before the specific MAC arrangement was discovered. In most cases, ideas were discarded because of difficulties encountered in the incorporation of boundary conditions. One particularly promising approach proved unworkable because of the impossibility of achieving rigorous mass and momentum conservation. Thus, the representation we have chosen is not as arbitrary as might at first be supposed; it is the result of considerable trial-and-error experimentation to assure the inclusion of many properties that seemed necessary for successful computing.
There are, in effect, two coordinate systems used in MAC-method calculations: The primary one covers the entire domain of interest with a rectangular grid of cells, each of dimensions δx by δy. The cells are numbered by indices i and j, with i counting the columns in the x direction and j counting the rows in the y direction. The field-variable values describing the flow field are directly associated with these cells. Their points of definition, relative to a cell, are shown in Fig 1.1.
Actually, the true fluid would, in general, have a different set of field-variable values for every infinitesimal point in the fluid. The representation used for computing, however, must be restricted to a finite number of values, each approximating an average through the immediately adjacent region. It follows that the accuracy of the representation depends strongly upon the fineness of the mesh compared to the macroscopic structure of the flow.
The placement of the field-variable quantities relative to the mesh is of considerable importance to the matter of conservation. No arrangement other than the one shown in the figure appears to be workable.
For example, if the field variables are placed at the cell centers, the following complication is introduced: To attain rigorous finite-difference mass conservation, the finite-difference equation for pressure would require the involvement of the next layer of cells beyond that which immediately surrounds any central cell. Such an involvement adds enormously to the complexity of the solution technique. Even more important is the expectation from previous experience that the use of far-distant quantities would introduce inaccuracies, at best, or even instabilities that could reduce the results to nonsense.
In addition to the primary-coordinate system of finite-difference cells, there is a coordinate system of particles whose motions describe the trajectories of fluid elements. These particles serve two purposes: First, they show which cells are surface cells, into which the surface boundary conditions should be applied. Second, they show the motion of the fluid and all its distortions as it passes through the computing region. For the first purpose, it would be necessary to have particles only near the fluid surface. The second purpose, however, is also of considerable value, particularly if facilities are available for easily plotting particle positions at various times through the progress of the calculation. We have found that complete configuration plots carry in concise form much of the information required for the analysis of results.
It should be emphasized that these particles serve only as mass-less markers of the centers of mass of the elements of fluid. They contribute nothing to the dynamics, and enter into the calculations only insofar as they designate which are the surface cells.
The cell-and-particle system enables an instantaneous representation of the fluid for any particular time during the evolution of the dynamics. In addition, it is necessary to have a means of actually calculating the changes with time of the fluid representation. What is needed is a computing technique whereby the prescribed initial conditions can develop, within the limitations imposed by the boundary conditions, into that subsequent set of configurations that most nearly represent the behavior of a true fluid.
As in most other fluid dynamics computing methods for transient problems, the MAC technique works with a time cycle, or movie frame, point of view. This means that the calculation proceeds through a sequence of cycles, each advancing the entire fluid configuration through a small, but finite, increment of time, δt. The results of each cycle act as initial conditions for the next one, and the calculation continues for as many cycles as the investigator wishes. Each cycle is itself subdivided into phases:
By the end of the cycle, the results have been arranged in the computer memory in such a way that the next cycle can immediately begin.
The finite difference analogy to D is
which we require for every cell at every time step.
Using Eq (1.4) as the basic form for this two-dimensional Cartesian example of the method, we obtain the difference equations
The superscript n or n+1 refers to a value at time nδt or (n+1)δt, so that n counts the number of time cycles. Where the superscript is omitted, n is implied, i.e., the value of the quantity at the beginning of the cycle.
A number of undefined quantities appear in Eq (1.15) and Eq (1.14). These are velocity components at localities for which values have not been stored in the machine memory. In each case, a simple average is to be used. For example,
Where a product of such quantities appears, it is to be understood that each is to be averaged first, and then the product is to be formed.
It may be noticed that as soon as the pressures are known for all the cells, then Eq (13) and Eq (1.14) are immediately appropriate for the calculation of new velocities, a process accomplished by simple algebraic substitution. To find an equation for the pressures, it is only necessary to manipulate Eq (13) and Eq (14) into an expression for the rate of change of Di j. First, define:
Note that the resulting values of φ lack superscripts; that is, they are time-centered at the beginning of the cycle, appropriate for their use in Eq (13) and Eq (14).
The solution of Eq (18) requires much of the computing time, because some type of iterative process is usually necessary. The process we use is described in Part 2 of this report; the details have no bearing on the basic MAC methodology, except insofar as the matter of efficiency is concerned. The solution of Eq (18) will increase the required computer time quite rapidly as the accuracy criterion is tightened. It follows that any process that decreases the accuracy requirement will speed the computing time. It also must be noticed that errors in the pressure solution, resulting from iterating to a coarse convergence criterion, will result in non-vanishing values of Di j. Thus, the reason for keeping the D terms in Eq (19) becomes clear. We have, in fact, compared the results of calculations run in the following three ways:
The conclusion is that the roughness in alternative (2) can be surprisingly great before the results of the comparisons show appreciable differences from those of alternative (1), but that number (3) differs strongly from the other two even if the roughness of the solution is not nearly as great as in number (2). This means that alternative (2) introduces considerable efficiency in the computer usage and is recommended for all MAC-method calculations.
Thus, Equations Eq (13), Eq (14) and Eq (18) are the three equations required for calculating the new quantities in each cycle. Once the pressures have been obtained from Eq (1.18), the velocities can be found directly from Eqs (1.13) and (1.14) without further iteration.
The marker particles enter into the solutions of the cell quantity equations only insofar as they serve to show where the moving, free surface is located, and accordingly which cells should have free-surface boundary conditions imposed in them. In order to keep this information on free-surface position current, it is necessary to move the marker particles each cycle in such a way that they accurately represent the fluid motion.
The technique we have been using is to find a velocity for the movement of each particle by using a simple area-weighted interpolation method among the nearby cell velocities. Details are given in Part 2 of this report. The principal accuracy criterion for the method of moving particles is based upon the requirement for conserving the volume of the region they fill. This is discussed in detail in the section on accuracy.
At rigid walls, the basic boundary conditions are simply that the normal velocity component vanishes and, in addition, that the tangential component vanishes if no slippage is to be allowed. With regard to the latter, the determination as to whether or not slippage is to be allowed depends upon the thickness of the boundary layer that would be expected to develop in the true fluid. If this is much less than the dimensions of a finite-difference cell, then a free-slip condition is appropriate; if it is larger than a cell, then a no-slip condition is required. For intermediate cases, the proper condition to use depends upon the exact circumstances, and in some cases it is appropriate to try both ways and compare the results.
In present applications of the MAC method, we restrict rigid walls, as well as influx and outflux walls, to follow cell boundaries. This means that all side or obstacle walls must be formed of horizontal or vertical segments only. (It is difficult to relax this restriction in Eulerian fluid dynamics calculations. We have found in other codes that diagonal walls or even circular ones can be built into a rectangular mesh, but the complexities are appreciable. The use of a curvilinear coordinate system might be helpful, but applicability would be restricted in each case to a relatively small class of well-fitting examples.)
The natural arrangement is to have vertical walls passing through horizontal velocity (u) points, and horizontal walls passing through vertical velocity (v) points. This allows for direct incorporation of the boundary condition that normal velocity vanishes.
If the wall is to allow for free slip, then whenever an equation calls for use of an exterior tangential velocity, the calculation simply uses the value of the tangential velocity at the image point back in the computation region. Likewise, if an exterior normal velocity is required, the negative of the image value is used. (Thus, the average of the exterior and interior normal components is zero, consistent with the vanishing of the normal component at the wall.)
If the wall is to have a no-slip condition, then the exterior tangential velocity component must be the negative of the interior image value. This, then, leads us to the following strange procedure for exterior normal velocity components: In order that D vanish for the exterior cell, we require that exterior normal velocity components have the same value as they have at their image interior points. The values at the wall still vanish, however, even if the interpolated values do not.
The reason that D must vanish for exterior cells adjacent to the boundary can be seen in reference to Eq (17), which shows that D will otherwise diffuse into interior cells. It should be noted that a set of difference equations based on Eq (1.5) instead of on Eq (1.4) would lack the viscous diffusion terms of Eq (17), and the requirement of a vanishing exterior D would no longer be necessary. We have not yet tried this approach, but will do so in the expectation that it will have several useful properties. (For cylindrical coordinates, the alternative of starting from Eq (1.5) appears to be the only acceptable one.)
In addition, we require boundary conditions for the solution of the φ equation. At rigid valls, the momentum equations supply precisely the needed guide for finding these conditions. They show how the normal difference of the pressure is balanced by the normal component of the body force and by the viscous diffusion of normal momentum. These balance expressions give exactly the information from which any equation needing an exterior φ value can find such a value in terms of available data. Examples of this and related boundary condition matters are given in detail in Part 2, which, in particular, shows how to handle the special problems arising at corners.
Boundary conditions along a line of prescribed constant-inflow rate differ only slightly from those at a rigid wall; the modification is straightforward. At an outflow wall, there is no unique prescription to apply. We have had good results in several examples by requiring the normal derivatives of the field variables to vanish. Further discussion of the details is given in Part 2.
It is the handling of free surfaces which particularly distinguishes the MAC method from all others for fluid dynamics computing. The general differential conditions have been stated in Eq (1.9) and Eq (1.10); they express the vanishing of tangential stress and the balance of normal stress with that which is externally applied.
If the viscous effects are negligible, then the normal stress condition is all that is required; and we simply put the appropriate local value of the applied pressure into each surface cell. Often, this value would simply be a fixed constant that we choose to be zero, representing the truly free surface adjacent to an atmosphere with negligible inertial or dynamic influence.
If viscous effects are important, then it is, in principle, necessary to construct finite-difference analogies to the stress balance conditions stated in Eqs (1.9) and (1.10). It is easy to see, however, that the construction of these analogies is troublesome because of the difficulty of sensing accurately the local orientation of the surface. As a result, we proceed somewhat differently. It is, of course, of primary importance that Di j = 0 for every surface cell, this inspite of the apparent contradiction that each surface cell is in the process of changing its total amount of mass. The vanishing of D refers to the interior of the fluid near the surface, with the resulting velocity components extrapolated to the exterior edges of the cell. Such a condition, however, is unique only if one side of a surface cell faces the exterior. If two sides face the exterior, then the vanishing of tangential stress ought to be invoked in order to determine the two boundary velocity values. Because of the difficulties described above, however, we require instead simply that
separately. Surprisingly, we find that even for very viscous fluid calculations, there is good evidence that this procedure results in negligible free-surface tangential stress. We have found no way, however, to derive rigorously a measure of the error to be expected from this treatment.
The balance of normal stress in a strongly viscous fluid is likewise approximated in our present MAC calculations. It can be seen that the viscous correction in Eq (1.10) is proportional to the normal derivative of the normal component of velocity. Experience has shown that in many circumstances this gives a very small stress correction, in comparison with the other stresses present near the surface. As a result, we simply neglect that term and use the externally applied stress only. Again, we have only experimental, but not analytical, justification for the approximation. For both of these stress-balance boundary conditions, we hope to experiment with more accurate representations of the true conditions. Meanwhile, we have been pleased to note that in many examples of interest our approximation model has been adequate.
If the flow is completely confined by rigid walls, including inflow and outflow walls, so that there is no free surface, then the MAC method is still applicable. The procedure is actually simplified in several respects:
In the absence of the free-surface boundary conditions, there is no reference pressure in the system to guide the convergence of the φ equation iterations. Convergence does not seem to be affected adversely by this lack of reference. In order to compare pressure results from cycle to cycle, it is useful to normalize the pressures after convergence. This can be accomplished by adding any desired constant to all of the φ values.
For confined flows, the MAC technique still appears to exhibit advantages over the stream-function-and-vorticity methods. Boundary conditions at rigid walls do not involve derived variables and, thus, can be applied easily and accurately without any special techniques or half-cell ambiguities in the definition of wall position. Examples of some confined flows are presented in Part 3.
The solution of initial-value problems by finite-difference approximations is almost always plagued by potential difficulties with numerical instability. The MAC method, however, is fortunate in that the stability restrictions are not very stringent for most applications. In addition, we have found that viscosity is not required as a stabilizer, an unexpected contrast to many of the Eulerian computing techniques in which stagnant regions become hashy in the absence of real or artificial viscosity.
It appears that there are two principal stability requirements, but we have not been able to demonstrate the sufficiency of either. One is directly analogous to the Courant condition that occurs in compressible flow calculations. In place of the sound speed, however, the wave speed appears:
In both cases, the result is effectively a restriction on the size of δt, the time increment per cycle. Because both conditions should be satisfied, it is necessary for any calculation to pay attention to whichever is the more restrictive.
For calculations with ν = 0, it is noticed that the particle arrangements usually become slightly irregular, in comparison to the neat line-up arrangement that persists when ν is not zero. Apparently the non-viscous calculations experience slight cell-to-cell velocity fluctuations whose principal manifestation is in the slight particle disarrangements. This effect has not been serious in any of our calculations.
It does not behave like a disastrous instability inasmuch as its effects are bounded to a very small magnitude.
Obvious requirements for accuracy include the necessity for cells fine enough to resolve the features of interest and for time steps small enough to prevent instability. A precise statement concerning cellwise resolution appears impossible to give. We generally say that the cells must be so small that no field variable changes by much across any cell. We find, however, that the usual interpretation of this is too stringent; very useful results have been obtained when the resolution was much coarser than might have been thought desirable.
One aspect of cell resolution concerns the boundary layer that can be expected along a rigid wall. It is not at all necessary to resolve the details of such a layer if the effect on the external flow would really be negligible. In such cases it is, in fact, appropriate to ignore an unresolvable boundary layer by incorporating a free-slip condition at the wall.
If the time step per cycle is just small enough to prevent instability, then the calculation may be quite accurate, in the sense that any smaller value of δt produces a negligible change in the results. This is not always the case, but is more likely to be so if the viscous stability condition is the more restrictive of the two.
The number of particles per cell used to define the fluid configuration has little effect on accuracy. However, for this statement to be true, it is necessary that cells interior to the fluid be considered as full, even if particle-position fluctuations should momentarily remove all particles from such a cell. (This, incidentally, can easily happen, especially in circumstances where compression in one direction and the corresponding expansion in the other result in the spreading of lines of particles.) For example, in defining the surface position, an average of four particles per cell is not much worse than sixteen.
Experience with many different types of computing techniques has shown that one of the most important contributors to accuracy is the degree of rigorous conservation of the finite-difference equations. For incompressible-flow calculations of the type discussed in this report, the primary quantities to conserve are mass (or volume) and momentum.
Mass or volume conservation is assured if, at every cycle, D= 0 for every cell, and if the motions of the marker particles well represent the velocities that contribute to the vanishing of D. It has been shown that a test of this mass conservation can be accomplished by keeping account of the total number of cells that have any particles in them [3]. If, for example, the computing region is bounded by rigid walls, so that the total number of particles is constant, then statistically one should expect the number of cells with any particles to be
in which A is the true area that the fluid should occupy, P is the perimeter of the area, and λ is the ratio of mean interparticle spacing to cell size.
In a typical example of the use of Eq (1.22), we suppose fluid to be laid out initially to fill an area ten cells by ten cells. Then A=100δxδy and initially, with a neat layout, one sees that the number of cells with any particles is 100. As the fluid moves, however, and the particles become less regularly placed, the number of cells increases. By assuming the glob of fluid to be removed from any wall, the initial perimeter of the glob is P=20(δx + δy); and its subsequent perimeter will continue to be approximately the same. With four particles per cell, λ= 0.5; and we see that ultimately
If the cells were square, then the number of cells with any particles would be 112.7. Examples of this sort have often been checked against the results of actual computer calculations, with excellent agreement.
Momentum conservation is also strongly required for accuracy of calculations. It was for this reason that terms in the transport part of the momentum equations of the form
were transformed (using the incompressibility condition) to
The difference form of such an expression can then be written as pure differences, so that the flux of momentum out of one side of a cell exactly equals the flux into that same side of the adjacent cell.
To see the importance of these conservation requirements, consider the example of the hydraulic jump. It is well known that the basic properties of the jump are completely determined by the overall conservation conditions from one side to the other, regardless of the detailed nature of the jump itself. In similar fashion, the overall properties of a flow can be given correctly as the result of rigorous conservation of the crucial physical quantities, even if the detailed structure is not accurately resolved.
In the absence of viscous dissipation or body forces, the differential equations of the MAC technique also conserve kinetic energy. This is not the case, however, with the difference equations; it does not appear to be possible to conserve mass (or volume), momentum, and kinetic energy in the finite-difference approximation. As in the analogy to the true hydraulic jump, which does not conserve kinetic energy, the energy seems to take care of itself. The basic source of non-conservation in the finite difference form comes from smoothing the velocity fluctuations caused by the finite cell size. The process resembles the production of entropy in compressible fluid calculations. We speculate, but have not proved, that the kinetic energy non-conservation rate is negative definite. If a heat energy equation is carried along in the calculations (for such purposes as determining heat flow rates or temperatures for buoyancy terms) then it should be possible to assure rigorous conservation of total energy.
There is one more aspect of accuracy that is worth mentioning. If the boundary conditions are incorporated into all phases of the calculation with precise consistency, and if the computer executes all instructions accurately, then the size of D for every cell will be bounded by a number that decreases as the convergence criterion for the φ equation is made more stringent. Therefore, the value of D should, at no time or place, exceed a specified upper bound. To check this, it is useful to build a comparison routine into the calculation that will print diagnostic information in case of a violation.
With u and v the velocity components in the r and z directions, respectively, we may write the appropriate differential equations for cylindrical coordinates in the form
Letting i and j count cells in the r and z directions, respectively, we may then define
This, then, completes the essential part of the derivation. The φ equation follows from Eq (1.30) with Dijn+1 =0. Note the absence of viscous diffusion of D in Eq (1.30), in contrast to the analogous Eq (1.17) for the Cartesian-coordinate description. The reason is not related to the difference in coordinate systems, but rather to the form of the differential equation from which the finite-difference momentum equations were derived. If the Cartesian-coordinate derivation were based upon Eq (1.5), the same type of modification to Eq (1.17) would have resulted. We have not yet actually tried this in calculations, but anticipate doing so in the hopes of simplifying the boundary conditions and improving the accuracy of the results.
Boundary conditions for the cylindrical-coordinate calculations would be very similar to those in Cartesian coordinates. Care must be taken, as before, in applying the conditions in precisely the same way throughout the calculation steps. Factors of r will be liberally sprinkled throughout the boundary conditions in order to accurately account for momentum balances.
In situations involving large motions of a fluid, it would be advantageous to have the coordinate mesh follow the mean fluid motion. If there is a constant translation rate, then the only procedural change involves the boundary conditions at rigid walls. Translation parallel to a no-slip wall means that the wall velocity relative to the mesh must be considered. Translation normal to a wall can be computed only if the translation speed equals the wall speed.
If the translation is to a non-inertial (accelerating) coordinate system, then the acceleration enters the momentum equations as an addition to the body acceleration term (gravity). For translation parallel to a wall, the same considerations arise as if the translation rate were constant. For translation normal to a rigid wall, in which the translation speed must always equal the wall speed, care must be taken in the pressure boundary conditions to include a balance against the accelerative term.
The only admonition to be offered in this report is that the investigator must be careful to assure rigorous mass (or volume) and momentum conservation in the difference equations. He must balance with care the appropriate momentum fluxes to derive boundary conditions that are physically meaningful; and he must take care to use the boundary conditions consistently in the same way through all steps of the calculation, in order to assure the vanishing of D. The apparent accelerations related to the constraints of curved mesh lines (non-inertial paths) must receive particular attention.
Problems in three-dimensional Cartesian coordinates appear to require only simple and obvious modifications to the present technique for two space dimensions. An iterative solution is required for but a single variable, the pressure, in contrast to stream-function techniques where it appears that three iterative solutions per cycle would be required.
The principal difficulties in three-dimensional calculations relate to computer limitations and to effective display of the results. At present, too much computing time is required if enough cells are used to give reasonable resolution in all three dimensions. The display problems are somewhat more mundane, related to human inefficiency for three-dimensional visualization - a problem beyond the scope of this report.
In many examples of interest, the full energy equation is not required; the dynamics of the incompressible flow are essentially unaffected by the distribution of heat, so that the latter simply accommodates itself to the changes in fluid configuration. This simple uncoupling is by no means a universal property of incompressible flows, however. The most common exception arises from variations of viscosity with temperature. Buoyancy effects may also be an important manifestation of temperature variations. Even without coupling back onto the dynamics, the study of heat flow can be of considerable practical interest in itself [5].
The incorporation of heat transport into the MAC calculations is quite a straightforward process. It is necessary only to add a step to the sequence of processes that advance all field variables each cycle. Even if the heat flow effects couple back to the dynamics, the required modifications are relatively simple to include.
The basic equation for heat transport can be specialized directly from the general energy equation of fluid dynamics
Here I is the internal (heat) energy per unit mass, often expressed as a product of specific heat and temperature; T is the temperature; k is the heat conduction coefficient; and φ is proportional to the rate of heat production caused by viscous dissipation. The other symbols retain their previous meanings.
The first term on the right side of Eq (1.31) represents the convection of heat produced by fluid motion. (Thus, even in problems for which the heat flow does not couple back into the dynamics, we see that the dynamics always couples strongly into the heat flow.) The second term is the heat production rate from compressive effects, which we ignore for incompressible flows. The third term is of considerable importance: it gives the heat flow from conduction processes. Finally, the last term in Eq (1.31) is negligible for most flows that can be called incompressible. To see that this is so, it is necessary to recall that the definition of incompressible flow is more accurately related to a comparison of fluid speed with sound speed than it is to the simple statement, ρ = constant. If the Mach number, M = u/cs, in which cs is local sound speed, is everywhere much less than unity, then the flow is incompressible. We can estimate the magnitude of φ by using only the fact that it is formed by a sum of the squares of velocity derivatives. Thus, if s is the distance over which any flow quantities vary appreciably, and
This is to be compared with either the convection term or the conduction term. The former can be estimated as follows:
that is, provided that the Reynolds number is very large in comparison to the square of the Mach number. This, of course, is almost invariably the case for the types of flow that the MAC technique can handle.
Comparison of the dissipative term with the heat conduction term shows the former to be negligible, provided that M22 Pr ≪ ≪1 (in which Pr is the Prandtl number, for most substances of order unity). Again, this criterion is certainly required anyway if the MAC technique is to be applicable.
Thus, for most purposes, it is sufficient to write for heat transport
Many techniques have been proposed for the solution of this equation by finite differences. In the MAC technique, the equation can be adapted readily to the mesh, as already defined. Temperatures are associated with cell centers; and the finite difference equation can be written, for example, as
The stability requirement for this expression is only slightly more stringent than that for the dynamics if Pr<1; otherwise, it is less stringent.
Appropriate boundary conditions for the heat equation can be derived directly from the physical requirements of insulation or constant temperature at a wall or free surface; there even could be a radiation condition at the free surface.
Incorporation of the effects into the dynamics are relatively simple. Variations of the viscosity coefficient with temperature involve calculation of the coefficient wherever required, using an appropriate function of the local temperature value. Buoyancy is incorporated, generally with the Boussinesq approximation, by allowing density variations in the bodyforce term, but not in the inertial terms. Thus, instead of ρ0g, one has (ρ0+δρ)g in the momentum equation. The density variation, δ, is calculated using a simple bulk expansion coefficient,
δρ/ρ = β(T-T0)
in which T0, is the temperature at which ρ= ρ0.
The groundwork for including chemical or nuclear energy production in the form of internal heat sources has already been laid in the above discussion. The right side of Eq (1.31) is simply modified by adding an appropriate source term, whose strength may be a function of temperature or of mixture concentration. The diffusion and convection of materials dissolved in the fluid can be followed by equations much like that for heat transport. Slight density variations produced by dissolved substances can affect the dynamics through buoyancy terms exactly analogous to those produced by heat. The importance of such computations for oceanographic or chemical engineering purposes would seem to be great.
The MAC technique was devised to handle the calculation of multi-dimensional problems in the time-dependent dynamics of viscous, incompressible fluids. It is, therefore, ideally suited for the study of the development of laminar instabilities for as far into the nonlinear phases of motion as may be desired. Free surfaces can be included, so that surface waves may form, break, recede, and form again, as many times as desired. Through all of this, the results are as accurate as the equations and their finite-difference approximations. This means that for a wide variety of phenomena whose important features can be resolved by the mesh, the MAC calculations are pertinent, useful, and representative of true physical processes.
The matter of resolution can, however, be a seriously limiting feature in at least one important type of phenomenon: turbulence. In almost every conceivable circumstance, the scale of turbulence is very small in comparison to the size of the flow-field structures in which the turbulence is created. Thus, to examine both the detailed structure of the turbulence itself, together with the corresponding effects on the macroscopic flow, it would be necessary to have a very fine resolution in the finite-difference mesh of cells. Add to this the fact that true physical turbulence is necessarily a three-dimensional phenomenon, and one must conclude that the actual representation of true, fully-developed turbulence cannot now, with present computer limitations, be included as a part of macroscopic fluid flow problems.
These considerations do not require that the important problems of turbulence be completely ignored. There are at least two ways in which meaningful studies can be carried on, even at the present time.
One of these is the study of the onset of turbulence. The earliest stages of turbulent growth from a laminar shear layer must be closely related to the problems of laminar instability. The question of infinitesimal versus finite perturbation is a significent one; and the effects of nonlinearity are crucial. The calculations are not easy and will require much careful thought both in formulation and interpretation.
The calculations will require much computer time, because the instability of interest is typically associated with large Reynolds numbers, of the order of 104, To see this, suppose that the conditions are achieved by means of unit velocity, unit flow region size and, thus, a kinematic viscosity of 10-4 . (This can be accomplished by the use of appropriate scaling.) In this set of units, laminar instabilities grow roughly as exp(10-2t). (This estimate is taken from the theory for plane Poiseuille flow.) This means that t=100 is required for significant amplification, during which time the fluid will have moved about 100 distance units. For resolution, each distance unit must have about 10 cells or more, meaning that the fluid will have traveled past at least 1000 cells through the development of the instability. Because the time step per cycle must be distinctly less than the time required for a cell-width of motion, it is seen that many computer steps will be required.
There are, of course, several ways to relax this pessimistic estimate. Commencing the flow with a large initial perturbation is the most obvious one. Other approaches will undoubtedly come to mind and be attempted; and, in the end, considerable information on the early stages of transition to turbulence may be forthcoming. It must always be remembered, however, that if the calculation is conducted in two-dimensional space, the results must not be given physical significance beyond the time when the line vortices would actually experience significant three-dimensional instability.
A second way in which turbulence problems can be studied with meaning is probably more applicable to practical problems. It involves the simulation of turbulence effects on the mean flow-field, through the use of transport coefficients: eddy viscosity and eddy conductivity. Prandtl mixing-length theory has formed the basis for a number of semi-empirical techniques for incorporating momentum and heat transport resulting from turbulent fluctuations. For many simple, steady-state situations, it is possible to derive expressions for an eddy viscosity coefficient which, when incorporated into the Navier-Stokes equations for the mean flow quantities, gives results in excellent agreement with experiments. Such eddy viscosity coefficients, however, are no longer simple fluid properties (like molecular viscosity); instead they are functionals of the mean velocity field itself (i.e., formed from various spatial derivatives of the velocity components).
For non-steady problems, the simulation of turbulent transport effects becomes slightly more complicated. The approrpiate wvay to proceed seems to involve two steps:
In steady-state flows, the first step can be accomplished through semi-empirical relationships between the turbulent flow properties and those of the local mean flow. In non-steady flows, however, the transport of the turbulent properties themselves will have to be calculated. Just as in the case of heat transport, an equation can probably be written for the transport of ψ, a measure of the turbulent properties,
The second step probably would not differ appreciably from the analogous step in the steady-state interpretation of turbulent transport.
The computational region is composed of a rectangular Eulerian mesh of cells in two-dimensional Cartesian coordinates, and a set of marker particles that define the fluid configuration. The cells are numbered by the indices i and j which refer to the cell center, where 1 ≤ i ≤ I, and 1 ≤ j ≤ J. Cell boundaries are designated by i+½, i−½, j+½, and j−½. For example, ui+½ j would be defined at the center of the right-hand boundary of cell i j, whereas ui j would be defined at the center of the cell i j. (See Fig 2.1).
The complete mesh of cells consists of IJ cells. In addition to cells falling inside the boundaries that define the system, there are cells falling outside the boundaries. This arrangement permits easier handling of boundary conditions. Because a free surface may be present, there is a necessity for different cell types within the system. Each cell is, therefore, flagged,as the different cells must be handled differently (see Fig 2.2). Cells are flagged as follows (where flagging means any selected method of identifying different cell types):
A given cell may have several flags. For example, a cell could be FULL, OB, and URON. Another cell might be BIND, FRSLP, COR, and URON.
Because the fluid is not stationary, it is obvious that a cell may change from EMP to SUR and then to FULL. It is, therefore, necessary to check at the end of each time cycle to see if cells need to be re-flagged.
Marker particles are numbered by the index k, where 1 ≤ k ≤ K. These particles are used to define the position of the surface and to give a visual representation of the fluid. Each marker particle is moved at the end of each time cycle, with a weighted average of the four nearest cell velocities. It is important to understand that these particles do not enter directly into the calculation but are used merely to define the position of the fluid and, in particular, the position of the free surface. Whether a cell is flagged as EMP, SUR, or FULL depends on whether or not the cell contains any marker particles.
In order to have fluid input and output, we must have some way to create and destroy marker particles. Therefore, the particles are flagged in the following manner:
The cell variables for the problem are positioned within the cells, as shom in Fig 2.3. Di j, Ri j, and φi j are positioned at the cell center. Velocities in the x direction are positioned at the left- and right-hand boundaries of the cell, and velocities in the y direction are at the top and bottom boundaries. In the formulas there are velocities needed that do not fall at these points; in this case an average is used.
For example:
The fluid flow is advanced through a series of time cycles, each of finite length δt. One time cycle consists of the following series of steps:
This completes one time cycle. The process can be repeated for as long as the problem is of interest, usually for several hundred cycles.
The type of boundary conditions applied depends on the type of boundary under consideration. The boundary conditions to be used for a boundary at the left wall will be discussed for each boundary type. The conditions at other walls are analogous. The indices i j will refer to the cell inside the system, and i-1 j will refer to the BND cell outside the system (see Fig 2.4).
a) IN: An input wall allows fluid to move into the system at a constant velocity; this velocity never changes throughout the run.
b) INCOR: Same as IN cell except:
c) OUT: An output wall allows for fluid to leave the system. Velocities are calculated for each time cycle.
d) FRSLP: Free-slip boundaries represent a line of symmetry or a non-adhering (greased) surface.
e) NOSLP: No-slip boundaries represent a viscous boundary.
The free surface is defined by a set of SUR cells. The treatment of the free surface is relatively simple and is as follows (ij = SUR cell):
The only exception is in the case where we have an applied pressure at the surface. In this case φi j =φa, a prescribed function of position and time.
Velocities at the free surface can be handled in a variety of ways, the main consideration being the requirement that Di j=0 for each SUR cell. For a SUR cell that is open on one side only, we calculate the other three velocities in the usual manner, and calculate the fourth by using Eq (2.1) and setting Di j=0. For a cell that has two open sides, we merely set the velocity for each of the open sides equal to the velocity opposite. For a cell with three open sides, the open side opposite the fluid side has the velocity of the fluid side; and the other two remain unchanged, except for the effect of body forces (gravity). A SUR cell with four open sides merely follows a free-fall trajectory.
When Eq (2.4) and Eq (2.5) are used to calculate velocities in SUR cells, velocities from EMP cells are sometimes referred to. Because these velocities are not defined, it is necessary to apply appropriate boundary conditions. This is done by using the velocities for the SUR cell each time an EMP cell is referred to.
For FULL cells only, Ri j is calculated with Eq (2.2). The only problems that might arise occur in cells next to a boundary (OB cells). In this case, the boundary conditions discussed earlier apply wherever needed. It is important to assure that the boundary conditions used in the Ri j calculations are rigorously consistent with those used in the velocity calculation. This is particularly true at corners of obstacles.
The pressure field is represented by a Poisson's equation, solved by an iterative procedure. Each iteration consists of solving Eq (2.3) for every FULL cell in the system, starting at the lower left-hand corner and working across and upward, respectively. The boundary condition at the free surface is merely φi j=0, or φi j=φa for each SUR cell. The boundary conditions for pressure at a wall are applied for each iteration. Iteration continues until the following conditions are met:
Velocities are calculated for all FULL and SUR cells in the system, using Eq (2.4) and Eq (2.5). When calculating for a SUR cell, care must be taken to assure that for each EMP cell velocity, some other appropriate velocity is substituted, in order to avoid having the EMP cells exert a false drag on the fluid. The treatment of COR cells, which must be taken into consideration, will be discussed under the heading Corners. Having found the new velocities for all FULL and SUR cells, we now calculate the velocities for BND cells, using appropriate boundary conditions and assuring that for all SUR cells Di j = 0 (as discussed earlier under Free Surface).
Velocities needed for calculating ui+½ j and vi j+½ are shown in Fig 2.5.
Marker particles are moved with a velocity that is a weighted average of the nearest cell velocities. velocity uk is calculated as an interpolated value of the four nearest horizontal cell velocities, and vk is calculated as an interpolated value of the four nearest vertical-cell velocities (see Fig 2.6).
The location of particle k with respect to the cells can be found as follows:
i = the integer part of (xk/δxi + 2)
j = the integer part of (yk/δyj + 2)
This locates particle k as being in cell i j. When calculating uk it is also necessary to know if the particle is in the upper half of the cell or in the lower half of the cell, inasmuch as the four nearest cell values for u are different in the two cases. If we consider
fy = the fractional part of (yk/δy + 2)
then, if fy < ½, we are in the lower half of the cell; and, if fy > ½ we are in the upper half of the cell. Likewise, when calculating vk it is necessary to know if we are in the left or right half of the cell.
One other thing needs to be taken into consideration: the movement of a particle near a COR cell. Again, reference is made only to the velocity in the x direction (uk). We have found from experience that a particle below LINE A (Fig 2.6, lower drawing) must be moved with an x-direction velocity which assumes that the wall extends upward (i.e., u2=0), if we are to avoid having a particle move into the COR cell. As soon as the particle is above LINE A, the particle should be moved in the x direction, as if the obstacle did not exist (i.e., u4=u2). In the case of a NOSLP wall, u4=-u2 as soon as the particle has moved into the cell above the COR cell. These same arguments can be applied to all COR cell orientations, and to the vk calculation.
After uk and vk have been found, the particles are moved as follows:
Corners of obstacles present numerous problems that must be considered throughout the program. The most obvious difficulty is that the cell pressure (φi j) is not uniquely defined for a COR cell, but depends upon which adjacent cell is being considered. When applying boundary conditions for cell i-1 j the value of φi j is different from that in the case of cell i j+1 (see Fig 2.7). This presents no real problem, but must be kept in mind when calculating pressures.
Another problem arises with regard to velocity calculations. In calculating vi j+1 (see Fig 2.8) for a cell next to a corner, the stored value of vi+1 j+½ (zero) should not be used in the equation. The appropriate values to be used are:
vi+1 j+½ = vi j+½ for free-slip condition,
vi+1 j+½ = -vi j+½ for no-slip condition,
Analogous conditions apply for both u and v, and for the four different orientations of the corners.
Problems arising during particle movement near a COR cell were discussed in the preceding section.
The most elusive problem arising near a corner is making the Ri j calculation consistent with the ui+½ j and vi j+½ calculation.
If this is not carefully handled, the Di j's for the three cells bordering a COR cell will not be nearly as small as they should be; and the system will not conserve mass properly. Let us consider the velocity calculation for cell i j in Fig 2.8. Conditions for the calculation of vi j+1/2 have already been discussed. As part of the boundary condition, ui+1/2 jn+1 is set to zero. Then, the important question arises: In setting ui+1/2 jn+1 equal to zero, what have we assumed? We have assumed that if we were to actually calculate ui+1/2 j, we would, indeed, get zero. Therefore, we must look at the equation for ui=1/2 jn+1, and see what set of conditions would give us zero if we performed the calculations. A careful analysis of the equation shows that the only part of the equation that does not vanish naturally is a term involving ui+1/2 j+1. We have, therefore, been assuming that:
ui+1/2 j+1 = 0
when we set ui+1/2 jn+1 equal to zeo.
Therefore, we must go to the Ri j calculation and make this same assumption. The only place in the Ri j calculation where ui+1/2 j+1 occurs is in the Di j+1 term. Thus, we must subtract ui+1/2 j+1/δx from the value of Di j+1 before we use it. Analogous arguments pertain to any cell lying above, below, or to the right of a COR cell.
Cells that are diagonal to a corner also require careful treatment. Referring to cell i j+1 in Fig 2.8, what must we assume when calculating R1 j+1? In the velocity calculation for ui+1/2 j+1n+1 we have assumed that
ui+1/2 j = ui+1/2 j+1 or ui+1/2 j = −ui+1/2 j+1
We must, therefore, make the same assumption in the calculation of Ri j+1, this time by adding or subtracting ui+1/2 j/δx from Di j before using it in the Ri j+1 calculation, depending on whether it is a free-slip or a no-slip wall, respectively.
Another approach is to assume that the φi j−φi+1 j/δx term must balance the term incorporating ui+1/2 j+1, when calculating ui+1/2 jn+1. This leads to an additional term that must be added to the boundary conditions in the pressure calculation. When calculating φi j, the condition for φi+1 j must have the additional term
+ (δxν ui+1/2 j+1) /δy2
This method is easier to apply than the other, but increases the required computer time, as additional tests are necessary within the pressure-iteration loops.
The setup for the problem comprises three main regions. The regions are gone through once at the beginning of a run, and are not included in the time cycle.
The first part of the setup is termed a GENERAL SETUP. This region reads input cards containing general information, such as the number of cells, viscosity, time step, and other parameters pertaining to the problem. Then, index words and constant parameters needed for the problem are set up and stored for later use. The x and y coordinates for the cells are calculated and stored, and finally, information needed by the output routines is calculated and stored.
The second part pertains to the setup of the CELLS. Cards are read that contain information concerning the system boundaries, including the shape of the system of cells and the types of boundaries represented. This information is then used to flag all cells in the system, as discussed earlier. Because at this time there are no particles in the system, all interior cells are flagged as EMP cells. The coordinates for the boundary points are stored in such a way that the plot routines can use them to draw the boundaries of the system whenever necessary.
Once all the cells are flagged, we can proceed to the third portion of the setup, which has to do with MARKER PARTICLES. Data concerning the marker particles are read into the machine; then the array of particles that will represent the fluid is created, and the coordinates of the particles are stored for future reference. The velocity field of the fluid is calculated and stored into the appropriate cells in the system, and the cells that contain particles are flagged as FULL or SUR cells. The setup is then essentially finished. Only two things remain to be done: A print is made of the initial conditions, and a plot is drawn to show how the cells are flagged. We are now ready to go into the first time cycle.
The program output is in the form of on-line and off-line listings, and Stromberg Carlson SC-4020 plots of particle configurations, velocity vectors, and pressure contours. A one-line, on-line print is taken for every cycle, showing the time, time step, number of iterations, number of cells containing particles, and other information pertaining to the iterative procedure.
There are two types of off-line prints, one presenting information pertaining to each marker particle, and the other containing cell data, such as velocities and pressures. The prints are made only occasionally, according to a predetermined print interval incorporated in the program or upon demand by the operator through the use of sense switches.
The best output from the program is in the form of the SC-4020 plots. Three types of plots are taken at given time intervals, which are input to the program.
In the programs employing marker particles, perhaps the most descriptive pictures and the easiest to obtain are those of particle configurations. Each marker particle has an x and y coordinate stored in the computer memory. By plotting these coordinates and drawing the boundaries that define the system, we get a picture showing the shape of the fluid and its relationship to the confining walls of the system.
Figure 3.2,in the next part of the report, shows a set of particle plots representing the flow of water from a broken dam. The first frame in the picture shows the fluid configuration an instant after the dam has broken. (The dam pieces have been removed.) Subsequent frames show the water as it flows downstream and collides with an immovable object. The mesh of cells is not shown in these plots.
Another type of picture produced by using particles is a plot of smoke (or streak) lines. In Fig 3.14 there is a continuous input of fluid from the left, and a continuous output of fluid on the right. As the fluid comes in from the left, it flows past an obstruction. The smoke lines shown in the top picture are formed by feeding in particles from the left, and allowing these particles to move downstream with the motion of the fluid. The effect is the same as that of injecting lines of smoke into a wind tunnel, or jets of dye into water flowing down a channel. As the fluid flows past the obstruction and into the channel, eddies are formed. The smoke particles visually demonstrate this quite well.
Particle plots, then, have the advantages of giving a nice visual effect, and of being relatively easy to obtain. They do not, however, convey complete information with respect to details of the flow. A single particle plot, for example, does not show the direction of flow, or any information about pressures or velocities.
For showing the direction of flow and the velocity field of the fluid, we use velocity vector plots (see Fig 3.14). For each cell in the system, we draw one velocity vector starting at the cell center, with a length proportional to the cell velocity and in the direction of the local flow. For each cell there is an x component of velocity (u), and a y component of velocity (v). We create a velocity vector by plotting two points, and connecting them with a straight line in the following manner:
x1 = xcell
y1 = ycell
x2 = x1 + kucell
y2 = y1 + kvcell
where k is chosen in such a way as to give the vectors a reasonable length for display.
A useful method for showing the pressure field is to make contour plots (again, see Fig 3.14). This is done by plotting lines of constant φ.
That is, each contour line represents a given value of φ; and a given contour plot will show lines for several values of φ, separated by a prescribed contour interval. The effect is the same as that of geographical contour maps, where each contour represents a certain altitude.
As long as the data are relatively smooth, the contour plots provide very useful and informative pictures. However, a knowledge of the problem under consideration is necessary in order to decide whether the lines show increasing values or decreasing values of φ, as the values of the individual lines are not printed on the plots. We do, however, know the value of the lowest and highest contours, and the contour interval. From these values it is possible to determine the values of the other lines.
Contour plots are not easy to produce. It is necessary to first find the maximum and minimum values for the pressure. From these values the contour interval, δφ, is calculated in such a way as to give a desired number of lines. It is also desirable that δφ be a rounded number, so that the contour plots can be more easily read and interpreted. One solution is to allow the number of contour lines for a given plot to fall somewhere between N and 2N, where N is an integer to be chosen. If we calculate
δφ = (φmin − φmax)/N
and then change δφ to the next lower power of 2, we have accomplished two things. First, δφ has been rounded off, and second, the number of lines will fall somewhere between N and 2N. Because δφ was originally calculated to give exactly N lines, use of the next lowest power of 2 will never give more than 2N lines. An appropriate value of N depends largely on the type of problem under consideration.
Having found the contour interval, we next find the locations of the contour lines within our system of cells. Inasmuch as φ is defined at the cell center for each cell in the system, the positions of these values form a rectangular array of points. The points can be thought of as forming a series of triangles. It can be shown that if a contour line passes between the two points that form one side of a triangle, it must also pass between the two points forming one of the other sides of the triangle (but never both of the other sides). The point where the contour line enters and the point where the line exits can be found by a simple linear interpolation; the two points can then be connected by a straight line. Therefore, if we consider each triangle individually and draw the short segment for each contour line that passes through the triangle, we will have a completed contour plot.
In a given computer run, there are usually several hundred time-cycles. If plots are made for every time-cycle, we have a motion picture that runs for several seconds. The motion pictures are useful demonstration techniques and also give additional information concerning the nature of the flow.
The most graphic movies are those made from marker particle plots. Marker particle movies have been made for the flow of water from a broken dam (see Fig. 3.2) the flow of water under a sluice gate (Fig. 3.5), and the von Karman Vortex Street (Fig. 3.14).
The movies have been compiled into a short film, with appropriate titles, available on request from the LASL Report Library.
Flow diagrams of the SPLASH code are on the following pages of this section.
A few results of MAC-method calculations have been published, but these fail to show anything like the full scope of applicability [2] [3] [4]. This part of the report presents an album of calculational results, designed to illustrate some additional types of problems for which the technique is suited. None of the examples is analyzed in detail. Previous publications have shown, by comparison, that the MAC-method results are accurate; and only a few additional comments in this regard are given here.
All of the calculations vere performed on the IBM 7030 (Stretch) Computer. The plots were processed directly from computer output through the Stromberg Carlson SC-4020 Microfilm Recorder, and are not retouched or otherwise altered.
The figures are tilted to give downward direction to gravity; in the actual calculation, the bottom of the tank was level, and there was a negative horizontal component of the body acceleration. The wave was generated by dropping the blob of fluid shown to the left at t=0. By the time t=6.0, the resulting wave reached the tip of water; subsequently, it ran up on the shore with decreasing amplitude. The bottom allows free slip, and the viscosity is negligible.
g = -1.0 Height of mesh = 2.1 ν = 0.01 No. cells in vertical direction = 23
The dam holding the reservoir is removed at t=0. Subsequently, the water falls away toward the obstacle in its path. The collision and splash over the top of the obstacle are shown in the last two frames. The calculated pressure history on the obstacle can then be used to predict damage. The bottom allows free slip, and the viscosity is negligible. Comparison of similar calculations (but lacking an obstacle) have been made with the results of experiments, and the agreement is well within experimental error.
g = -1.0 ν = 0.01 Length of mesh = 4.8 Height of mesh = 4.0 No. cells across bottom = 50
The wave is generated, as it was in Fig 3.1, by dropping a blob of fluid at the left, outside the computing region shown. The relatively empty region just under the wave crest is treated, calculationally, as being full of water. As in Fig 3.1 and Fig 3.2, all rigid boundaries allow free slip; and the viscosity is small enough to have a negligible effect on all results, except for a slight smoothing of the particle arrangements.
g = -1.0 Length of mesh = 9.0 Height of mesh = 2.1 ν = 0.01 No. cells = 92 × 23 (2116)
The calculation is similar to that in Fig 3.3, but the breakwater has been widened into a shelf or reef. Note that the wave breaks at late times.
g = -1.0 Length of mesh = 9.0 Height of mesh = 2.1 ν = 0.01 No. cells = 92 × 23 (2116)
Water in a reservoir behind a sluice gate is forced out under pressure into a shallow, quiescent reservoir. The resulting wave breaks backward, toward the gate. The walls allow free slip, and the viscosity is negligible.
g = -1.0 Length of mesh = 4.8 Applied φ in deep reservoir = 2.5 ν = 0.01 Height of mesh = 3.0 No. cells = 50 × 32 (1600)
The calculation resembles that of Fig 3.5. Note that the shallower downstream reservoir permits the wave to break forward, at first.
g = -1.0 Length of mesh = 4.8 Applied φ in deep reservoir = 2.5 ν = 0.01 Height of mesh = 3.0 No. cells = 50 × 32 (1600)
From an opening at the upper left, water pours down onto the bottom and splashes off to the right. The resulting wave becomes highly irregular. The left wall allows free slip; but the bottom has a no-slip condition, and viscous drag is important.
g = -1.0 Length of mesh = 7.5 Input velocity = -1.0 ν = 0.01 Height of mesh = 2.5 No. cells = 77 × 27 (2079)
This is the same calculation as shown in Fig 3.7. The line segments are computer-plotted velocity vectors. Note the effect of the no-slip bottom boundary condition.
The run shown in Fig 3.7 and Fig 3.8 was continued to very late times, giving these marker-particle and velocity vector configurations.
The box of water has rigid ends, and the water is initially moving to the right. It piles up on the end, and a jump progresses to the left, at first into water of uniform depth, then into shallowing water. The walls and bottom allow free-slip, and the viscosity was chosen to be just great enough to prevent breaking of the wave. (For the effects of varying viscosity, see Fig 3.12.) In calculations of this type, the height and speed of the jump agree very well with analytical predictions. Tests of this type were made for a variety of situations, as a useful way to demonstrate the validity of the calculations.
g = -1.0 Length of mesh = 9.6 Initial water velocity = 1.0 ν = 0.10 Height of mesh = 2.2 No. cells = 98 × 25 (2450)
These are velocity vectors (one for each computational cell) and isobars for the same calculation illustrated in Fig 3.10. The interval between lines of equal φ is δ φ= 0.10 with the top line at φ=0.05.
Effects of viscosity are shown in these comparative frames, all of which are at a time t=4.0 after the jump was input. In each calculation, quiescent water was present at t=0; and the jump, with a vertical face, was fed in thereafter from the left. The input conditions were chosen to match the analytical jump predictions. The comparison between no-slip and free-slip bottom boundary is designed to show that viscosity can affect the appearance of the jump in two different ways: both internally in the fluid and as a result of drag. Many aspects of these results have been compared qualitatively with experimental data, and no disagreements were observed. Unfortunately, the available experimental data are too restricted for such comparisons to be really crucial tests of computing accuracy.
Another set of tests of these particular runs compared the number of cells containing any fluid with the Santalo prediction given in Eq (1.21). In every case, the comparison was accurate to well within the fluctuational uncertainty of the statistical equation.
g = -1.0 Length of mesh = 9.8 Height of mesh = 2.3 No. cells = 100 × 25 (2500) Initial depth = 0.5 * Input depth = 1.5 * Input velocity = 1.0 *
* Note: These conditions satisfy only the small ampitude bore equations! Actually all of the three values should be slightly different to be strictly correct.
A very viscous fluid is input from the left onto a flat plate. A comparison between no-slip and free-slip bottom boundary conditions shows a very large difference in results.
g = 1.0 Length of mesh = 9.8 Height of mesh = 2.2 ν = 0.10 No. cells = 100 × 24 (2400) Input depth = 1.5 Input velocity = 1.0
Water pumped in from the central region rises against gravity and falls sideways in both directions. The particle plot is shown above; velocity vectors, below.
g = -1.0 ν = 0.01 Input velocity = 1.0 Length of mesh = 3.0 No. cells across bottom = 30
The applicability of the MAC technique to confined flows is illustrated by this calculation of the von Karman vortex street formed behind a rectangular cylinder. The results have been compared with those obtained previously and the agreement is excellent [5]. The example shown here is for a downstream Reynolds number of 100. Both the obstacle and the confining walls had free-slip boundary conditions, so that vorticity had to be created at the obstacle corners. Production of the proper amount confirms the validity of the present treatment at corners (see discussion in Part 2). The figures show different ways of plotting the fluid state, all for the same time after the street was well developed.
The top frame illustrates the streakline configuration created by lines of particles fed in from the left. The middle frame is a plot of velocity vectors; the streamlines can be visualized easily from this plot. The bottom frame gives lines of constant pressure, with equal intervals between lines. Lift and wall stresses can be obtained accurately by integrating such pressure details.
g = 0.0 Distance between walls = 3.0 Height of obstacle = 1.0 ν = 0.01 Input velocity = 1.5 Output velocity = 1.0 Interval between isobars = 0.03125
The calculation resembles that of Fig 3.15, but there is a free surface above the obstacle. The pictures show late-time configurations of particles (above), velocities (middle), and isobars (below).
g = -1.0 Length of mesh = 5.625 Height of mesh = 4.0 ν = 0.01 No. cells = 45 × 32 (1440) Input velocity = 1.5 Obstacle height = 2.0
A thick, viscous fluid is pumped in from the upper right. It falls onto a rigid plane surface, some splashing backward and some forward, where it is lost from the computing region. A late-time particle configuration is shown above. Corresponding velocity vector and isobar plots are also given (middle and below, respectively).
g = -1.0 Length of mesh = 4.5 Height of mesh = 3.8 ν = 0.10 No. cells = 1710 Input velocity = 1.0
A bubble, initially circular in shape (bottom), deforms as it rises through a confined fluid. The particle configurations are shown for t=0, 1.6, 3.5, and 4.5 (from bottom to top, respectively).
g = -1.0 Height of mesh = 4.7 ν = 0.1 No. cells = 47 x 12 Input bubble radius = 10.4 Bubble pressure = 0.0 (constant in time)
Velocity vectors (at t=1.6, 3.2 and 4.5) and the pressure at t=4.5 are shown for the calculation described in Fig 3.18.
The MAC method has had its greatest application in the investigation of the effects of inhomogeneity on the dynamics of an incompressible fluid. The inhomogeneity discussed previously was of the extreme form represented by a free surface. In this part of the report we describe an extension of the MAC technique which studies the effects of more moderate density discontinuities. The change is one more of emphasis than of degree and, therefore, requires a somewhat different procedure for obtaining a numerical solution. For this reason, the calculation of the flow of mildly inhomogeneous fluids is being considered as a separate part of this report.
The crucial difference in the two techniques concerns their treatment of the fluid interface. As implied by its name, a free-surface calculation considers the interface to be a boundary line between a fluid region and an empty region. Its position is determined through the use of marker particles, which have no other function in the calculation. Computations are made only for the occupied parts of the mesh, and the fluid in these cells is considered to be uniform.
In the two-material technique, however, provision must be made for gradations in density in the cells that mark the interface. Furthermore, this treatment must be such that the position of the interface remains well defined throughout the course of the calculation. This dual requirement has been satisfied by enlarging the role which the particles play in the calculation. They are now used not only to mark the density discontinuity, but also to determine the values of cellular densities and viscous coefficients in all cells in the mesh.
The list of field variables for calculating the dynamics of a heterogeneous fluid must be somewhat expanded over that of Part 1 to include:
The independent variables are the time and the Eulerian coordinates.
The density and viscosity coefficient are considered to be constant in homogeneous regions of the fluid. At the interface these constant values change abruptly; but we assume that there is always a functional relationship between ρ and μ so that the above list of dependent variables can be reduced by one. In general, however, this relationship is not a direct proportionality, so that we cannot assume a constant kinematic viscosity, ν=μ/ρ, throughout the system.
The equations describing the flow of a heterogeneous, incompressible, viscid fluid are most often written in the forms
Equation (4.1) relates local density changes to fluid transport. Equations (4.2) and (4.3), expressing volume and momentum conservation, respectively, are identical to their one-material counterparts, Eqs (1.2) and (1.3), except that the two-material momentum equation contains a more complicated viscous diffusion term. This complication is necessary to account for spacial variations in the coefficient of viscosity. In the case where μ is constant throughout the fluid, it can be shown - with the aid of Eq (4.2) - that Eq (4.3) reduces to Eq (1.3).
The similarity between these two equations extends to their conservation properties. Just as Eq (1.3) when written in finite difference form, fails to rigorously conserve momentum, so also does Eq (4.3). Furthermore, the latter equation is inappropriate to a two-material calculation because it does not take into account the changes in momentum that result from density variations. Both of these difficulties can be remedied by combining Eqs (4.1) and (4.3) and using Eq (4.2) to obtain
It is this form of the momentum equation that is used throughout Part 4.
The conservative nature of Eq (4.4) can be demonstrated in the same way as was that of Eq (1.4). Integrating over a fixed volume and applying the divergence theorem to the appropriate volume integrals, we can show that, neglecting gravitational accelerations, momentum changes result entirely from external forces.
Although Eulerian finite difference techniques have previously been applied to multifluid flow calculations [6], they have been limited in application by their tendency to smear density discontinuities. This characteristic difficulty results from the fact that the method makes no provision for resolving interfaces. An Eulerian cell containing two fluids is usually assumed to have the mixture distributed uniformly throughout. Therefore, when mass fluxes are calculated, it is this mixture that is transported to the neighboring cells. This results in a calculational mass diffusion that quickly eradicates sharp density discontinuities.
Therefore, a technique was needed which would maintain sharp fluid interfaces without becoming so complicated and time-consuming in operation as to preclude the calculation of meaningful physical problems. It was discovered that this could be accomplished by extending the marker-particle-and-cell concept of the MAC method.
The approach is similar to that used in Part 1 in that the calculations are performed relative to the Eulerian mesh of cells, while the fluid is represented within these cells by marker particles. But, whereas in a typical one-fluid problem the particles differentiate fluid cells from empty cells, in a two-fluid problem they must distinguish one fluid from the other. This is accomplished by using different types of particles to represent each fluid. We differentiate between them by flagging the data associated with each type of particle in a characteristic way.
In addition to the fact that the particles are no longer uniform, there is another important manner in which they differ from those used in a single-fluid calculation. While the latter particles merely mark the position of the fluid but do not directly affect the calculations, the particles with which we deal in Part 4 are used to determine the values of ρ and μ needed for the cellwise calculations. Specifically, if a given Eulerian cell contains nx particles of x fluid and ny particles of y fluid, then ρ and μ in that cell are given by
At the fluid-fluid interface one could expect to find cells with various values of ρ and μ intermediate between those of the homogeneous cells.
In many respects, this type of interface treatment is far easier to incorporate into a numerical technique than is the free-surface treatment of Part 1. For example, there is no ambiguity concerning the volume of fluid contained in an interface cell; it is entirely filled with fluid just like any other cell in the mesh. Therefore, the criterion for volume conservation holds without modification, even though the mass of the cell may be changing.
Theoretically, the normal and tangential stress conditions are also automatically satisfied at the interface (at least to the extent that the finite-difference resolution and our knowledge of the viscosity of inhomogeneous fluids will permit). Experience has shown, however, that there may be difficulty in achieving continuity of pressure at an interface when the fluids have markedly different densities, and that drag will occur along strong slip lines even if the fluids are inviscid. The latter problem results from the particle movement procedure described below. Each of these difficulties is calculational in nature, and neither is considered insurmountable. They do, however, require further study.
The points of definition of the velocity components relative to the cellular mesh are the same as those indicated by Fig 1.1. The other field variables - pressure, density, and viscous coefficient - are cell centered.
With these definitions, we can write the finite-difference analogies of the two components of Eq (4.4):
As in Part 1, the absence of a superscript indicates that the quantity should be evaluated at time nδt.
When values of the velocity components, the density, or the viscosity coefficient are required at localities other than their point of definition, a simple average is generally used. The only exception to this rule involves the terms that express momentum transport normal to the calculated direction. In such cases, the momentum flux is taken in the direction indicated by the flow field. We illustrate the procedure with the term (ρuv)i+1/2 j-1/2, from Eq (4.6). This term expresses the flux of x-direction momentum in the y-direction. It could assume either of two possible forms:
A similar prescription applies to the other momentum flux terms.
The incompressibility condition, Eq (1.12), is the third and final difference equation that is required. It is used, in the manner described below, to determine the pressure and density fields necessary for the solution of Eq (4.6) and Eq (4.7).
Before explaining how these difference equations are applied, let us briefly review the computing procedure in Part 1. The incompressibility condition is applied in each cell in such a way as to be effective at time (n+1)δt. This process yields a finite difference Poisson equation for the pressure with a known source term in each cell. The Poisson equation is solved by means of an iterative technique. Finally, the resulting pressure field is used in the momentum equations to evaluate the advanced values of the velocity components. This velocity field is automatically conservative.
The same technique is used here, except that we now must solve for the pressure field and the (advanced time) density field simultaneously. This is accomplished by means of a double relaxation method: By guessing at an advanced time density field, the associated pressure field is determined by an iterative procedure. These pressures and densities are then used to find the velocities. Particle trajectories are calculated from the velocities in order to determine the new density field. This dual iteration process is repeated until, eventually, the densities remain static and the pressure field is sufficiently relaxed. Essential to this technique is the marker particle method of determining densities, for it insures that cell densities will change only by discrete increments.
Let us now go through the procedure in detail. For convenience we abbreviate Eq (4.6) and Eq (4.7) as
For a given density field, the incompressibility condition, Eq (1.12), can be written
The density field at time not is used as the starting guess for the advanced densities. With these density values, Eq (4.8) is solved (to a rather crude degree of accuracy) for the pressures. This provides us with sufficient data to compute the advanced velocity values from Eq (4.6) and Eq (4.7) and to determine particle trajectories. The particles are not actually moved at this time. We merely note any trajectories that would cause a particle to cross a cell boundary and contribute to a density change. The appropriate changes are made by means of Eq (4.5), and the resulting density field is used as the second guess in solving for the pressures in Eq (4.8). The process is repeated until there are no changes in density. Then, the pressures are relaxed very accurately; the new velocities are calculated; and the particles are moved. The final values of the densities and the viscous coefficients are calculated from Eq (4.5). A flow diagram showing the order of these steps is presented in Fig 4.1.
In general, the final density values will be identical to those used in the last pressure iteration; a discrepancy would be an indication that the crude pressure iteration was not sufficiently accurate. Occasionally, however, a cell density will fluctuate from one value to another on each succeeding iteration as a result of the fact that a particle has landed almost exactly on the cell boundary. In such a situation the density field would, of course, never converge. Therefore, it is necessary to sense this condition when it occurs and terminate the iterations. In this case a density change after the final iteration is permitted.
At first reading, the finite difference scheme described above may appear so complicated and time-consuming as to preclude its application to meaningful physical problems. In fact, however, such is not the case. Calculation times appear to be roughly equivalent to the time required for single-fluid MAC problems. The reason is that there are ways to optimize the efficiency of the two-part relaxation technique, which consumes the largest percentage of the calculation time. The iteration scheme and its important timesaving features are described in the following paragraphs.
Prior to the start of the iterations, the ξ and ζ terms are computed; thereafter, these terms remain fixed throughout the entire calculation cycle. Also during this preliminary period, we combine the source terms and the coefficients of the pressures in Eq (4.8) in order to write that equation in the simplified form,
These coefficients and the source term, Ai j, will also remain fixed for the entire cycle, unless the cell or one of its neighbors undergoes a density change in the course of the iteration. Notice that Eq (4.9) is a Poisson equation only in those (homogeneous) regions of the fluid where the B's are all equal to unity.
To obtain the solution of Eq (4.9), we employ the relaxation technique described in Part 2 of this report. Successive sweeps of the mesh (from left to right and from bottom to top) are made, computing the pressure in each cell by solving Eq (4.9) in terms of the most recent pressures from the neighboring cells. Convergence is tested periodically to determine whether additional iterations are required.
Two types of pressure iterations have been discussed in Part 4:
The crude iterations are simply a timesaving device, because absolute accuracy is not required at that stage of the calculation. They do, however, provide a headstart toward the solution and thereby shorten the fine iterations. Both types make use of the method described above, but they differ in respect to their convergence criteria and the frequency with which these criteria are checked. The convergence of the crude iteration is tested after four sweeps of the mesh against the criterion
where g is the magnitude of the body acceleration and h is the depth of the fluid. The fine iteration is checked after every 13 sweeps of the mesh against the number 0.0002.
The process of calculating density changes is greatly facilitated by a constraint placed on the time increment, δt. We require, for the accuracy of the numerical method, that the time step be sufficiently small so that no particle could traverse more than a single cell in one time cycle. As a result, we may be sure that the interface will not move more than a cell width in a cycle. This permits us to restrict our attention to a relatively few cells in the vicinity of the interface for the purpose of calculating density changes. All other cells will remain homogeneous during this time cycle.
In order to take advantage of this restriction, we define two particular types of cells:
Before any iterations begin, a search is made through the cells to determine the interface and contributing cells for that time cycle. These cells are then suitably marked by flagging the data associated with them.
In calculating particle trajectories after the crude pressure iterations, we may limit our attention to the contributing cells, inasmuch as these are the only particle movements that could cause density changes. The velocities used to calculate these trajectories use the latest pressure and density data by applying the following simple formulas:
We take note only of those trajectories that cross cell boundaries, and register their effect by tallying the gains and losses for each type of particle for the cells involved.
When all the required particle movements have been accounted for, it is time to determine the resulting density changes, if any. For this purpose we consider only the interface cells. The number and type of particles gained or lost by each interface cell are combined with the number of particles of each type that were contained in the cell at the beginning of the cycle, in order to determine the cell density by means of Eq (4.5). If the new cell density is different from its previous value, then new values of the B's and of Ai j must be calculated for that cell and its four neighbors for use in Eq (4.9). If the density is unchanged, we proceed to the next interface cell.
If none of the interface cells experience a density change, then we consider the density iterations to be complete and proceed with the final pressure iterations. If, however, there were density changes, then the crude pressure iteration and density change calculation must be repeated. After three such repetitions, we review all interface cells to determine which have undergone density changes. The interface cell flags are turned off for those cells that have experienced no changes during the three iteration cycles. The contributing cells are similarly re-evaluated in order to progressively reduce the number of cells for which these calculations are required.
The particle movement technique is exactly the same as that used for the single-fluid calculations. The components of the particle velocity are computed as a weighted average of the four nearest cellular velocity components. The purpose of the weighting is to give the greatest emphasis to the cell velocities nearest the particle.
Figure 2.6 illustrates the manner in which the weights are determined for each velocity component. Consider the x-component. A cell-size rectangle is imagined to be centered about the coordinates of the four nearest u-velocity components, and a similar rectangle is placed about the particle in question. The proportion of area overlap of this latter box on the other four determines the weighting to be applied to each of the four velocities.
For the purposes of numerical calculation, the whole process may be reduced to a single algebraic expression (see the flow diagrams for Region 60 in Part 2 for a detailed description of this numerical treatment; a similar procedure for the y-component of the particle velocity is also described there).
The calculations described here use rather simple boundaries. In all cases, the calculational regions are rectangles, bounded by rigid walls (coincident with mesh lines) with no interior obstacles. At the rigid walls, the normal component of velocity vanishes. The tangential component may also vanish; or it may be unaffected by the wall, depending on whether no-slip or free-slip boundary conditions are used. As mentioned in Part 1, the choice of the boundary condition depends upon the width of the boundary layer anticipated in the actual problem. If the boundary layer were on the order of a cell width or greater, a no-slip boundary condition would probably be required. If it were smaller, a free-slip boundary could be used.
The finite-difference implications of the no-slip and free-slip boundary conditions have been described in some detail in Part 1 and Part 2; and in this part of the report, we merely summarize the points of addition to it.
The Eulerian calculation method requires the consideration of an extra row (or column) of boundary cells outside the region of calculation (see Fig 4.2). The following summary discusses the treatment of the flow variables in those cells for no-slip and free-slip boundary conditions.
Density and viscosity coefficient: In order to make the boundary treatment of velocities and pressures consistent with that used in the single-fluid calculations, it is necessary to choose the density and the viscosity coefficient in the external boundary cell equal to the values of those quantities in the image cell inside the boundary. This is done at both no-slip and free-slip walls.
Velocity components: Rigid walls are always chosen to correspond to mesh lines, along which normal velocity components are defined. The vanishing of the normal component of velocity at rigid walls is therefore easily specified. The tangential components and the normal component outside the wall are chosen in such a way as to simultaneously satisfy the boundary condition and the incompressibility condition in the exterior cell.
In the following examples, consider a rectangular mesh with horizontal cell indices varying from i = 1 to I and vertical cell indices from j = 1 to J. The external boundary cells have a horizontal index of 0 on the left or I+1 on the right and a vertical index of 0 on the bottom or J+1 on the top (see Fig 4.2).
No-slip walls: Tangential velocities reverse, while the external normal velocity is equal to that at its image point.
At the left-hand wall,
At the right-hand wall,
At the bottom wall,
At the top wall:
Free-slip walls: The tangential velocity is equal to that at its image point while the external normal velocity reverses.
At the left wall,
At the right wall,
At the bottom wall,
At the top wall:
Pressures: With densities, viscosity coefficients, and velocities determined as above, the pressures are defined in such a way as to satisfy Eq (4.6) or Eq (4.7) identically (depending upon whether the boundary is vertical or horizontal).
No-slip walls:
At the left wall,
At the right wall,
At the bottom wall,
At the top wall,
Free-slip walls:
At the left wall,
At the right wall,
At the bottom wall,
At the top wall,
The stability requirements for one-fluid calculations, Eq (1.20) and Eq (1.21), apply equally well in homogeneous regions of multifluid problems. At the interface, however, the analytic determination of stability criteria is complicated by variations in density and viscosity coefficient. Here one must rely on experience, which to date indicates that no difficulties are encountered at the interface when the homogeneous fluid requirements are satisfied.
These stability criteria are, in effect, limitations on the size of the time increment, δt. Another limitation, required for accuracy, restricts fluid motion through the cellular mesh by means of the requirement
where Umax is the maximum velocity component in the mesh and δX is the smaller of δx and δy.
The numerical examples described below make use of the maximum value of δt, consistent with these stability and accuracy requirements, in order to reduce calculation times. The time step is tested at each cycle of calculation to determine whether an increase or a decrease is in order. In the event that one or more of the above conditions is violated, δt is halved; and then the tests are repeated as many times as necessary to avoid violations. The time increment is doubled if the left-hand expression in the inequality (4.10) is less than 0.125 and if the stability requirements are satisfied using the new value of δt. Only a single doubling is permitted in each time cycle.
As mentioned in Part 1, an important measure of accuracy is the degree of rigorous mass and momentum conservation of the finite-difference equations. The discussion of momentum conservation given there is equally true here and will not be repeated. The subject of mass conservation in a multifluid closed system is, however, sufficiently different from that in single-fluid, free-surface calculation to warrant additional comment.
Volume conservation for the system as a whole follows from the fact that it is closed, plus the requirement that all cells in the system be full of fluid at all times. The latter requirement holds even in those cells temporarily void of particles. In such a case, the cell retains the density and viscosity coefficient associated with it at the time it became empty of particles. Such instances are not uncommon, especially where the fluid as a whole is undergoing large distortions, but they are generally shortlived.
Volume conservation on a cell basis results from the stipulation that Di j vanish in every cell at each time cycle, so that the net flow into a cell is exactly balanced by the outward flow. However, the mass in a particular cell may vary as a result of a change in the proportionate number of particles of each type that the cell contains.
Overall mass conservation, therefore, depends upon maintaining the initial areal distribution of the marker particles of each fluid. Consequently, mass is not rigorously conserved at every time cycle; but we do expect only minor fluctuations, because the motion of the marker particles is governed by a conservative velocity field. As a test of conservation, we sum the mass in all cells at every time cycle and compare it with the initial mass. Typically, the percentage absolute mean deviation from the initial mass is about 0.1%. The maximum deviations for particular examples are given in the following section.
A reduction in the size of the mass deviation could no doubt be achieved by increasing the number of marker particles. However, the small discrepancies noted in the calculations performed to date do not seem to warrant the added calculation time that would result from such an increase. All of these calculations used an average of four particles per cell.
Examples of calculations performed with the multifluid extension of the MAC technique are presented in the form of particle plots printed directly from the computer data, with the aid of a Stromberg-Carlson SC 4020 microfilm recorder. The particles of the two fluids are differentiated by being printed with different degrees of darkness. For clarity, the interface between the fluids has been sketched in by hand; but otherwise the photos are not retouched.
Two fluids rest side by side in a box, separated by a thin diaphragm. The lighter fluid, marked by the heavy black dots, is on the left; and the heavier fluid is on the right. At t=0 the diaphragm is removed, the heavy fluid falls under the force of gravity, which is directed vertically downward. The two fluids revolve until the interface is nearly a horizontal line in the final figure of the sequence. Further running of the problem indicates that the interface overshoots this equilibrium line, reaching to the top of the box.
The walls of the box are rigid boundaries at which free-slip boundary conditions apply. Both fluids are inviscid. The apparent drag along the interface is the result of using area-weighted velocities for particle movement.
Number of calculation cells = 30 × 20 (600) Number of marker particles = 2400 Density radius = 2 Kinematic viscosity (both fluids) = 0 Maximum deviation from initial mass = 0.5%
A heavy fluid is superposed over a lighter fluid in a gravitational field. Initially, the fluid is perturbed by a low-amplitude, conservative-velocity field whose vertical component varies in the form of a cosine wave along the interface. This small perturbation grows exponentially for a time. When nonlinear effects become important, the interface develops a spike and bubble shape characteristic of Taylor instability. At late times, the effects of Helmholtz instability are visible along the edge of the spike.
Number of calculation cells = 20 × 60 (1200) Number of marker particles = 4800 Density ratio = 2 Kinematic viscosity (both fluids) = 0 Maximum deviation from initial mass = 0.1%
This problem is the same as that illustrated in Fig. 4.4, except for the added complication of viscosity. The primary effect of viscosity is to retard the growth of both Taylor and Helmholtz instabilities:
Number of calculation cells = 20 × 60 (1200) Number of marker particles = 4800 Density ratio = 2 Kinematic viscosity (both fluids) = 0.001 Maximum deviation from initial mass = 0.2%
1. J E Fromm, A Method for Computing Nonsteady, Incompressible, Viscous Fluid Flows, Los Alamos Scientific Laboratory Report IA-2910 (September 1963).
2. F H Harlow, J P Shannon, J E Welch, Liquid Waves of Computer, Science 149, 1092 (1965).
3. F H Harlow, J E Welch, Numerical Calculation of Time Dependent Viscous Incompressible Flow of Fluid with Free Surface, Phys. Fluids, , 2182-2189 (1965).
4. F. H. Harlow, J E Welch, Numerical Study of Large Amplitude Free Surface Motions, Phys. Fluids, submitted for publication.
5. F. H. Harlow, J E Fromm, Dynamics and Heat Transfer in the von Karman Wake of a Rectangular Cylinder, Phys. Fluids 7, 11 4 7 (1964).
6. A. Blair, et al, A Study of a Numerical Solution to a Two-Dimensional Hydrodynamical Problem, Los Alamos Scientific Laboratory Report IA-2165 (October 1958).