- Introduction
- Definition
- Mathematical aspects
- Finite-volume equations
- Integration procedure
- Storage implications

- Implementation in PHOENICS
- Graphical display
- Examples

When computational fluid dynamics first engaged the attention of engineers, during the 1960s, parabolic flows were a natural focus, because of their relatively small demands on computing power. For example, the Patankar-Spalding program of 1967, which was later developed into GENMIX, concerned two-dimensional parabolic flows, such as boundary layers, jets and wakes.

A later Patankar-Spalding contribution to CFD, of 1973, namely the SIMPLE algorithm, was also applied to three-dimensional parabolic flows in the first instance.

It is therefore not surprising that PHOENICS is equipped with a "parabolic option", the nature and utilisation of which will now be described.

- steady (
*i.e.*independent of time), - predominantly in one direction, defined as that in which the velocity vector nowhere has a negative component; and
- with only negligible diffusion fluxes in that direction.

These conditions often prevail in boundary layers, shear layers, jets, wakes, pipes and other conduits, when the Reynolds and Peclet numbers, based upon the cross-stream dimension, are large.

Smoke plumes, flows in not-too-winding rivers, and jet-engine exhausts are practical examples.

Such flows are called parabolic because mathematicians have applied that adjective to the differential equations which describe them, contrasting them with the related but different "elliptic" and "hyperbolic" ones.

The main reason for distinquishing parabolic flows from elliptic ones is that their special features may be exploited so as to reduce greatly the time and expense of the calculations.

This exploitation is effected, in PHOENICS, by setting PARAB = T in the Q1 file. The same applies, incidentally, to hyperbolic flows, which also allow differences to flow downstream only.

**Hyperbolic** flows are those in which the velocity in the
main direction of flow is supersonic, which circumstance prevents
the upstream transmission of information in a steady flow. Then the
same economies as are associated with parabolic flows, namely
unlimited grid-refinement in the z-direction, can be enjoyed.

It is sometimes then said that the "parabolized Navier-Stokes equations are being solved".

Whatever it is called, PHOENICS provides for it.

Click here for the "mathematics of PHOENICS lecture"

The finite-volume equations for parabolic problems differ from the general ones in three respects.

The **first** is that the term

phiH*aH

In parabolic calculations, the terms expressing variations with time are also absent; so the general equation reduces to the following form when PARAB=T:

phiP = aE.phiE + aW.phiW + aN.phiN + aS.phiS + aL.phiL + S ---------------------------------------------------- aE + aW + aN + aS + aL + aP

A **further**, less obvious, difference is the omission of the diffusive component
from the coefficient in the term

phiL*aL

The **third** difference concerns the representation of the pressure
force on the z-directed velocity resolute w. The z-wise variations
of pressure level, p, are **decoupled** from the lateral (*i.e.* x and y)
variation of the pressures (which continue to be determined by
reference to cell-wise continuity).

The force on the w velocities is taken as

-Dp/Dz .

Because of the absence of the A_{H} term, a single sweep through the integration domain suffices. However, in order that the non-linearity and inter-connectedness of the equations can be adequately
represented, a sufficient number of **iterations** must be conducted at each slab, in
order that imbalances in the equations have been adequately reduced;
for single-sweep parabolic calculations allow no 'second chance'.

Therefore, instead of the LSWEEP = 50 (say) which is set for an elliptic-flow problem, the corresponding setting in a parabolic-flow calculation would be LITHYD = 50 .

To make one forward step in the integration sweep, it is necessary to hold in computer memory the variables relating to only two slabs, namely (1) the local one, and (2) its immediately-upstream neighbour.

This means
that an **unlimited** number of slabs can be employed without any any
increase of computer storage.

This has the advantage that a very-fine-grid solution can be obtained; for all available storage can be used for the cross-stream directions.

In PHOENICS, the predominant direction of a parabolic flow is always the z-direction. Storage is therefore provided only for the current and the low slabs.

Attempts to access "high" values by way of In-Form statements or GROUND coding will therefore fail. So, probably will attempts to access "low" values; for, as soon as a new forward step starts, these are moved into current-step storage as initial guesses.

To instruct PHOENICS to simulate a flow in the parabolic manner, it is first necessary to set PARAB = T in the Q1 file.

LITHYD should be set to a sufficiently high value to ensure convergence. LSWEEP, which is used for that purpose in elliptic problems, has no significance for parabolic flows.

The integer IPARAB also needs to be set in order to indicate how the downstream pressure is to be computed. This differs as between confined and unconfined flows. For the former (IPARAB=0), EARTH calculates automatically what pressure gradient will allow mass-continuity to be satisfied.

IPARAB.GT.1 also allows certain hyperbolic flows to be simulated economically (See PHENC entry: IPARAB). The important case of IPARAB=1 is discussed at length in section 4.4 below.

It is necessary that the lateral dimensions of the domain, *i.e.* XULAST or YVLAST, vary with with downstream distance so as to:

- contain the growing extent of the diffusion-influenced region; and
- expand or contract according to whether the main-flow velocity is decreasing or accelerating.

Users are therefore advised to use instead the In-Form grid-specification statements described here.

Caution is needed when selecting the width-enlargement formulae. Over-rapid enlargement places some cells too far away from the physically-interesting region to be economical; and too slow a rate leads to implausibly-steep profiles near the free-stream boundary. Inspection of the computed profiles, and adjustment of the enlargement-rate formula so as to avoid the extremes, is always advisable until a user has found formulations which always work satisfactorily in his or her circumstances.

It is best to employ **solution-adaptive** In-Form formulae such as:

This increases YVLAST in proportion to the step size, DZ; and also to a user-chosen factor yrfac times the solution-generated expression :

When this velocity aproaches the free-stream velocity, the enlargement of the width of the grid is much reduced, as can be seen from the following table, which shows the total enlargement of the grid,when this formula is used, for core-library case 190.

yrfac | w ratio | YVLAST/YINLET |

0.1 | 9.419406E-01 | 1.132596E+00 |

1.0 | 9.897953E-01 | 1.450421E+00 |

10.0 | 9.985036E-01 | 1.944726E+00 |

100.0 | 9.998467E-01 | 2.314851E+00 |

The first of the next two contour diagrams. corresponding to yrfac=0.1, reveals that indeed that factor is too small; for the contours intersect the upper edge. The second, for yrfac=100.0, shows that the enlargement is rather excessive; for a significant region of the upper region has uniform velocity. nevertheless even these extremely different grids produce not-very-different solutions in the important region close to the wall.

Allowance for changes in main-stream velocity are best accomplished by including an additional multiplier in the YRAT (or XRAT) formula, namely wfree_last/wfree_current. This will be illustrated below.

The former may be represented by way of In-Form thus:

What are dwfz and wfre? The following in-Q1 lines answer:MAKE DPDZ is 0.0) (STORE1 DPDZ is -:rho1:*(:dwfz:)*wfre) patch(prgrd,volume,1,nx,1,ny,1,nz,1,1) (source of w1 at prgrd is -dpdz)

wherein dwfz must have been declared as a PIL variable to represent the (uniform) increase of w1 with z.(MAKE wfre is :wfree:) (STORE1 wfre is :WFREE:+(:dwfz:)*ZZW)

As mentioned in section 4.3, the width of the grid should vary with the mainstream velocity. Specifically the width*velocity product should remain constant. The following yrat formula effects this:

Human error being always to be expected, it is wise to check the compatibility of the dpdz and wfre settings. The following image provides such a check by supplying a uniform upstream profile and de-activating the solid-wall boundary condition. The free-stream velocity is caused, by the widening of the grid without diffusional effects to further enlarge the boundary layer, to fall linearly from 33 to 27 m/s.

The fact that the contour lines are precisely vertical shows that the *wfre *formula, which controls the upper edge, is entirely consistent with the *dpdz* formula, which controls the values below it. This satisfactory check lends credibility to the quite-different contours which are generated when the **non**-uniform upstream profile **and** frictional-wall boundary condition are reactivated. It is shown below; and it reveals that the adverse pressure gradient has the effect of causing what is called 'boundary-layer separation', *i.e.* a zero-velocity region close to the wall.

In the case in question, which was derived from input-library case 190, the grid-expansion factor accounted for **both** diffusional and velocity-reduction effects; and was:

The imposed deceleration value,(YGRID of YRAT is 1.0-(:dwfz:*dz/wfre)+:yrfac:*(1.-w1[,ny-1,]/$ :wfree:)*DZ/YVLAST)

In general, because the cells of the parabolic and elliptic grids are likely to be of different sizes, as shown below, interpolation is needed.

In such cases the values of *dpdz* and *dwfz* extracted by way of interpolation will vary with the distance *z* along the surface and will appear as tables of numbers. In-Form allows these to be input into the calculation by way of the *PWL* (piece-wise-linear functions); but care is necessary in formulating them to ensure that the two functions are wholly consistent.

The above sketch containers a reminder that a **two-way** interchange of information may take place between the elliptic and parabolic calculations. Thus the elliptic calculation may take place at first, with the assumption that friction at the solid surfaces is absent. Its predicted pressure distribution is for that reason not quite correct. However the ensuing parabolic calculation takes detailed account of friction and can report the so-called 'displacement thickness' of the boundary layer. If this is
transmitted back to the elliptic solver, that can repeat its calculation on the assumption that the effective size ot the solid object is larger than it first supposed. Its second flow prediction will be correspondingly more accurate.

Iterative exchanges between the two solvers can ensue; and, because of the relatively weak interactions between them, convergence is likely to be quickly achieved.

These two simple lines have the effect of placing in the 3D store of the variable called(stored var dis is dyv*(1.-w1/:wfree:)) (stored var idis is dis+idis[,-1,]!if(iy.gt.1))

Evidently the displacement thickness can be taken as approx. 1.58E-03 m in this case. [By contrast, the displacement thickness for the aboveIY IDIS DIS 22 1.583E-03 0.0 20 1.589E-03 1.090E-05 18 1.573E-03 2.331E-05 16 1.529E-03 3.747E-05 14 1.456E-03 5.342E-05 12 1.349E-03 7.103E-05 10 1.205E-03 9.013E-05 8 1.020E-03 1.105E-04 6 7.927E-04 1.320E-04 4 5.208E-04 1.551E-04 2 1.946E-04 1.946E-04

No doubt, if the thus-computed value is fed back to the elliptic solver, which thereupon transmits a different pressure distribution to the parabolic solver, the latter will calculate a slightly different displacement thickness next time. Further refinement will however not be justifiable. We are, after all, employing a turbulence model, not to mention the finite number of grid cells; and these are practices are notoriously unreliable.

It is convenient to include the quantity *idis* among those of which contour diagrams are plotted. Those for the cases cited above are as follows:

The different behaviours of the separating (left) and non-separating (right) boundary layers are visible at a glance.

Since, by definition, downstream events can have no upstream influence in a parabolic situation, no downstream boundary conditions are needed. Moreover, attempts to provide them may have undesirable effects, because PHOENICS will attempt to apply them at the last slab that it knows about, which is the current one.

Consequently, simply to set PARAB = T, at the end of a Q1 file that had been earlier set up for an elliptic-mode solution, would give an undesired result. For example, the IZ=NZ boundary condition corresponding to a fixed-pressure outlet would prevent any pressure variation at all.

Users who have hitherto used only elliptic-mode (*i.e.* PARAB=F)
settings, are therefore advised to study the parabolic examples in
the Input File Libraries.

One consquence of there being no outlet patch at the last high-z surface is that no values of the outflows through it are printed in the RESULT file. Therefore, for parabolic calculations, the significant differences between the "neg sum" and "pos sum" values in the "nett flux" print out might disquiet users familiar only with elliptic flows. The differences do not however indicate poor convergence; they point only to an insufficiency in the print-out provisions.

If PARAB is T, EARTH storage is provided only for the current z-slab
and the one preceding it (*i.e.* the lower slab); as a result, the data
file used for PHOTON and EARTH restarts (PHI or PHIDA) contains only
data from the last slab.

It can therefore be used for restarts but is not very useful for viewing results with PHOTON. However, a GX... subroutine, GXPARA, has been provided for the creation of PHOTON- readable output files for the graphical representation of the results of parabolic runs. GXPARA is called only when the following setting are made in the Q1 file:

PARAB = T and IDISPA = 1 or greater.

Setting IDISPA to a value greater than one will cause GXPARA to group the slabs in groups of IDISPA, and dumping only the values at the central slab of each group to the output file, thus allowing for economy of space and processing time for problems with a large number of slabs.

Users can also set IDISPB and IDISPC to the first and last slabs, respectively, for which results are to be dumped. The defaults are: IDISPB=1; IDISPC=NZ.

Formatted or direct-access files are created by GXPARA according to whether PHIDA=T or F in has been set in the PREFIX file by the user. Formatted files are named PARPHI, direct-access files PARADA.

If either XULAST or YVLAST is a function of IZ, then a geometry file PARXYZ (formatted) or PARZDA (direct-access) file is also written. PHOTON treats these files as it would the files generated by EARTH for a body-fitted-coordinate problem.

If there are any variables of which the values are not required to be written to PARPHI, the variable name should be set to a 4-character string, the last 2 characters of which are '**'; eg. NAME(R2)=R2**

(NOTE that only PHIDA or PHI should be used for restart runs, NOT PARADA or PARPHI)

An early example: a buoyant smoke plume in a tunnel

The turbulent Bunsen burner with a 14-fluid turbulence model

The plane turbulent mixing layer, with a multi-fluid turbulence model

More recent **compressible-flow** developments are described below.

Because the flow is partly sonic and partly subsonic, the PIL variable IPARAB has to be set equal to 5 in the Q1 file; and STORE(MACZ) should appear, also, in order that the appropriate testing for Mach Number can be carried out.

The typical succession of higher- and lower-pressure regions is clearly seen in the following pictures, in which the flow is from right to left.

Mach number contours

pressure contours

axial-velocity contours

velocity vectors

The flow is three-dimensional. The IPARAB and MACZ settings are as before.

The typical succession of higher- and lower-pressure regions is clearly seen in the following pictures.

Mach number

contours

pressure

contours

axial-velocity

contours

velocity vectors

When the flow is wholly supersonic, IPARAB may be set to 4 in the Q1 file; and STORE(MACZ) is not needed.

The following pictures relate to such a flow. The gas enters from the left, in a direction parallel to the lower wall.

A shock wave starts from the upper wall, in order to change the flow direction to accord with the inclination of that wall. This wave crosses the duct and is then reflected from the lower wall.

Mach number contours

pressure contours

density contours

velocity vectors