The LVEL turbulence model for conjugate heat transfer at low Reynolds numbers

The LVEL turbulence model for conjugate heat transfer at low Reynolds numbers
by
Dereje Agonafer (IBM), Liao Gan-Li (CHAM) and Brian Spalding (CHAM)

Abstract

A simple method is described for calculating distances from walls in circumstances of complex geometry.

This is then used as the basis of a comparative study of conjugate heat transfer by means of three low-Reynolds-number turbulence models.

One of these models is new.

It is argued that the new model, called LVEL, performs as well as the older Lam-Bremhorst-Yap and 2-layer-k-epsilon models.

Since it is computationally less expensive, and easily applied to three-dimensional problems, it is recommended for practical use, especially for electronics-cooling applications.

Contents list

The problem considered
Method of solution
Tests for a two-dimensional case
Exemplification for a three-dimensional case
Conclusions
References
Appendices

1. The problem considered

Turbulence models exist in the literature for computing steady flow and heat transfer at low and transitional Reynolds numbers; but two difficulties attend their application to circumstances arising in engineering practice, namely:-

several of them require the distance from the nearest solid wall to be computed, for every point in the field of flow; and
it is rarely possible to employ the fine computational grid necessary for the terms in the equations to be accurately evaluated.

Difficulty (1) has both practical and conceptual aspects. The first arises from the fact that the domain of interest often contains many solids, of various shapes and sizes, so that calculating the distance to the wall involves a search procedure.

The second arises from doubt as to how the differing distances, which are found when the search is carried out along lines of differing direction, should be weighted in order to assign a single distance is to a selected point.

Difficulty (2) also results from the very large number of objects which may be present in the domain, rendering it often necessary to represent the gaps between them with only two or three grid intervals.

In these circumstances, the employment of methods which require the accurate computation of (say) velocity gradients cannot be justified.

The present paper reports a method which circumvents these difficulties and which has the additional advantage of being computationally inexpensive.

It utilises what has been called the LVEL turbulence model (because it requires knowledge only of the wall distances and the local velocities).

The LVEL model is to be regarded as providing a practical solution to heat-transfer-engineers' problems, and not as one providing new scientific insight.

2. Method of solution

2.1 The calculation of the distance from the wall

The LVEL model requires, like the Two-Layer-k-epsilon and Lam-Bremhorst-Yap models to be referred to later, knowledge of the distance from the wall.

It computes this by means of a method reported at the 1994 International Heat Transfer Conference (Spalding, 1994).

This method deduces the wall distance from the solution of the differential equation:

div grad W = - 1

with the boundary condition: W = 0

within the whole of the space occupied by fluid.

W is not itself the wall distance; indeed its dimensions are those of distance squared.

However the distance to the wall, and also the distance between walls, are deducible from the local value and the local gradient of W, by way of a formula which is exact when two extensive parallel walls are present (the only case where the concepts have precise significances).

For all other circumstances investigated so far, this formula gives results which are plausible.

The solution of the above differential equation is of course a very simple task for any computer code which is capable of solving the much more complex Navier-Stokes equations.

As implemented in PHOENICS (Spalding, 1996) therefore, the solution of the W equation is a swiftly-conducted preliminary to the main flow-simulating calculation which uses its results.

2.2 The effective-viscosity calculation

The LVEL model is of the "effective-viscosity" variety.

This means that it provides the values of the "effective" transport properties of momentum and energy which are necessary for the solution of the time-averaged Navier-Stokes and heat-transport equations.

The calculation of the effective viscosity, here given the symbol V, proceeds as follows.

(a) The main idea

For channel-like flows it is presumed that a universal velocity profile exists which connects:

the dimensionless distance from the wall,

y+ (defined as distance * square root(shear stress/density)/ laminar viscosity) with:

the dimensionless velocity parallel to the wall,

u+ (defined as velocity / square root(shear stress/density)).

By differentiation of this u+_versus_y+ expression, a universal relationship can be derived between:

(1) the dimensionless effective viscosity,

v+ (defined as effective viscosity/laminar viscosity),

and

(2) the local Reynolds number, which equals: u+ * y+ .

The use of this effective viscosity in the momentum equations amounts to the employment of a 'zero-equation turbulence model'.

This is bound to yield accurate results for channel and pipe flows; and it can be expected to be approximately correct for more complex interactions between solids and fluids.

This is indeed all that the LVEL model aims to do.

(b) The law of the wall

For the u+ versus y+ relation, a differentiable formula is used which covers the entire laminar and turbulent regimes (Spalding, 1961). It is as follows:

the dimensionless quantities y+ and u+ are linked by:


**************************************************
*                                                *
*  y+ = u+ + (1/E) * [ exp(K*u+) - 1 - K*u+ -    *
*      (K*u+)**2/2 - (Ku+)**3/6 - (K*u+)**4/24 ] *
*                                                *
**************************************************

Here K is the von Karman constant (0.417) and E is another constant (8.6) needed to fit the well-known logarithmic law.

The expression inside the square bracket is the exponential function less the first five terms of its Taylor-series expansion.

This law is known to fit experimental data very closely in: - the laminar sub-layer close to the wall, - the logarithmic region far away from the wall, and - the intermediate transition layer.

It entails that the dimensionless effective viscosity v+, which must equal d(y+)/(dv+) by definition, can be deduced by differentiation as:

v+ = 1 + (K/E) * [ exp(K*u+) - 1 - K*u+ - (K*u+)**2/2 - (K*u+)**3/6 ]

This implies that v+ = 1 close to the wall, where u= is very small; and, where u+ is large, far away from the wall, it reduces to the well-established result:

v+ = K * y+

v+ - 1 is, of course, the turbulent contribution to the (non-dimensional) effective viscosity.

(c) Solving the equations

Let the local Reynolds number R be defined as:

R = the local value of the absolute velocity of the fluid * the distance from the nearest wall / the laminar viscosity.

It follows that R = u+ * y+ . Therefore R can be computed for every point in the flow from the formula:

R = u+ * {u+ + (1/E) * [ exp(K*u+) - 1 - K*u+ - (K*u+)**2/2 - (K*u+)**3/6 - (K*u+)**4/24 ]} .

Further, u+ can be computed for every point in the flow, albeit by an iterative Newton-Raphson procedure, such as:

u+ = guessed_u+ + (R_actual - R(guessed_u+))/(dR/du+) .

As a consequence, v+, and so the effective viscosity V, can also be computed for every point in the flow, which is what the turbulence model is supposed to permit.

3. Tests for a two-dimensional case

3.1 The physical problem considered

In order to test the utility of the method, calculations have been made of the flow and heat transfer in the situation illustrated in Fig.1, through which air flows steadily.

The two steel blocks are uniformly heated. The problem is to compute the maximum temperature within the solid at various rates of flow of air.


********************************************************************
*                                                                  *
*       ///////////////////////////////////////////                *
*       --------------- adiabatic wall ------------                *
*   --->                   duct                    ----->          *
*       -------------                 -------------                *
*       // steel ///|     cavity      |/// steel //                *
*       -------------------------------------------                *
*       ////////////// aluminium //////////////////                *
*       -------------------------------------------                *
*                                                                  *
* Fig.1  Horizontal dimension 3 cm, vertical dimension 1.2 cm      *
*        The two steel blocks and the cavity are each 1 cm wide.   *
*        The blocks, the aluminium base and the air-flow duct      *
*        are each 0.4 cm thick.                                    *
*                                                                  *
********************************************************************

3.2 The models employed

Predictions were made, with the above-indicated input data, by activating three of the low-Reynolds-Number models which are supplied as standard in the PHOENICS computer code, namely:

(1) the LVEL model, as described in the present paper;

(2) the Lam-Bremhorst-Yap model (Lam & Bremhorst, 1981; Yap, 1987);

(3) the two-layer k-epsilon model (ElHadidy, 1980; Rodi, 1991).

Since all three of these models require, as an input, the distance from the wall, this was computed by solving the W equation as described above.

Four different uniform computational grids were used, namely:

(a) 6 * 6 ; (b) 15 * 15 ; (c) 30 * 30 ; and (d) 60 * 60 .

Only the last of these can be regarded as likely to provide grid- independent results.

However, the coarser grids are still interesting because, as has been explained above, they are much more likely to represent the numbers of cells which can be devoted to flow elements such as are shown in Fig.1, when these form parts of much larger assemblies.

3.3 The results for a single run. LVEL is in use

The grid is 60 * 60 . The Reynolds Number is 2000. Flow is from left to right.

Velocity vectors
Wall-distance contours
Wall-gap contours
The temperature field
The effective-viscosity distribution
Contours of left-to-right velocity component
Contours of pressure

Table 1 contains the computed values of the maximum temperatures in the metal, organized according to Reynolds number, grid and model, the latter being distinguished by the abbreviations, with obvious meanings: lvel, 2-layr and lam-br.



------------------------------------------------------------------
| Table 1. Maximum temperature rises in the metal, in deg Celsius|
------------------------------------------------------------------
|  Re  grid lvel  2-layr lam-br  |   Re  grid lvel 2-layr lam-br |
|--------------------------------|-------------------------------|
|  100  a  12.96  12.98  13.01   |  1000   a  8.66  7.90   8.76  |
|       b  17.15  17.32  17.40   |         b  5.03  5.03   5.70  |
|       c  17.43  17.61  17.32   |         c  4.04  4.32   5.12  |
|       d  17.57  17.76  17.84   |         d  4.62  4.97   6.01  |
|                                |                               |
 COMPARE NUMBERS UNDER LVEL, 2-LAYR AND LAM-BR ON ANY HORIZONTAL

|  200  a  12.40  12.43  12.32   |  2000   a  6.58  5.72   6.13  |
|       b   8.81   9.13   9.28   |         b  3.61  3.57   4.32  |
|       c  12.58  12.93  13.18   |         c  2.89  2.62   3.80  |
|       d  12.75  13.12  13.36   |         d  2.70  2.82   3.61  |
|                                |         ----------------------|
|  500  a  10.84  10.65  11.11   |  Computer times for attainment|
|       b   6.39   6.67   7.03   |  of steady state to within 1%,|
|       c   6.26   6.82   7.36   |  in minutes on a Pentium 200  |
|       d   7.90   8.27   9.01   |   2000  d   <1     10    >10  |
------------------------------------------------------------------

The computer times in the bottom right-hand corner of the table were deduced by visual observation of the PHOENICS "graphical monitor", and should therefore be taken as approximate.

The >10 entry for the Lam-Bremhorst model appears because oscillations of maximum temperature of amplitude in excess of 1% were still continuing at the 10-minute mark, and seemed likely to continue.

Generally, convergence was more difficult to procure for the Lam-Bremhorst model than for the other two.

3.4 Discussion

Four features of the results are relevant to the argument of the present paper, namely:

The grid-fineness effect is indeed appreciable, the temperature rise varying by a factor of more than 2 between the coarsest and the finest in the worst cases.
Moreover, the effect was surprisingly irregular, in that grid b was sometimes better and sometimes worse than grid a.
The Reynolds-number effect is appreciable.
The differences of the LVEL-based results from those of the other two models are small compared with the differences which arise from the coarseness of the grid.
The differences of the LVEL results from those of the 2-layer model are, in most cases, smaller than those between the 2-layer and Lam-Bremhorst-Yap values.

The computer times for attainment of convergence are much smaller for LVEL than for the other models.
Feature (1) needs to be carefully considered by those users of (say) electronics-cooling packages who find that the passages between solids are represented by very few grid intervals.
The second feature is of course expected, and requires no other comment than that a double effect is present, in that, at the lower Reynolds numbers, the temperature-rise in the air is appreciable.
Features (3) and (4) are the ones of most relevance to the aim of the present paper; for they lead to the obvious question:
Why use either of the computationally-expensive models when the LVEL gives results which are equally reliable (or, of course, unreliable)?

4. Exemplification for a three-dimensional case

The geometry of the two-dimensional study just reported is immensely simpler than the three-dimensional geometries encountered in practice. It is therefore important to emphasise that the LVEL method is very well-suited to the solving of conjugate-heat-transfer problems in complex three-dimensional geometries.

Space limitations preclude extensive demonstrations of this assertion; but the point is perhaps sufficiently made by the presentation of Figs 2 and 3, which show the geometry, and an associated surface of constant temperature, of a three-dimensional conjugate heat-transfer study, in which both forced and free convection are active.

5. Conclusions

It has been shown that the LVEL model of turbulence is capable of providing predictions of conjugate-heat-transfer phenomena which are as suitable for engineering purposes as are those of more elaborate low-Reynolds-Number models.

Its use for both two- and three-dimensional studies has been demonstrated; and it has proved to be both robust and economical.

It deserves consideration in all those circumstances in which complexities of geometry and limited computer memory prevent the use of the extremely fine grids which numerical accuracy requires, especially when the Reynolds number are transitional.

6. References

Spalding DB, 1961. "A single formula for the law of the wall" Trans ASME Series A, J Appl Mech vol28, no 3, pp 444-458
Spalding DB, 1994, Poster session. International Heat Transfer Conference, Inst. Chem. Eng. London
Spalding DB, 1996, The PHOENICS Encyclopaedia, Concentration, Heat and Momentum Ltd, London
C.K.G.Lam and K.Bremhorst, 1981. ' A modified form of the k-e model for predicting wall turbulence', ASME J. Fluids Engng., Vol 103, p 456.
C.Yap, 1987, 'Turbulent heat and momentum transfer in recirculating and impinging flows', PhD Thesis, Faculty of Technology, University of Manchester.
M.A.Elhadidy, 1980,'Applications of a low-Reynolds-number turbulence model and wall functions for steady and unsteady heat-transfer computations', PhD Thesis, University of London.
W.Rodi, 1991, 'Experience with two-layer models combining the k-e model with a one-equation model near the wall', AIAA-91-0216, 29th Aerospaces Sciences Meeting, January 7-10, Reno, Nevada, USA

7. Appendices

The parametric study was performed in a single multi-run, in which the new PIL command incl( ) was employed, the included files being data0 and data.

Appendix 1 The data0 file

This defines variables used for distinguishing one run from another

REAL(REYNO,UIN,TKEIN,EPSIN,QIN) integer(unit) char(mod)

Appendix 2 The data file

!!!! This is mainly a conventional Q1; but some novel features are utilised.
GROUP 1. Run title and other preliminaries
TITLE DISPLAY The problem concerns 2d incompressible, turbulent flow and heat transfer in the following geometry
-----------adiabatic----------------------- ---> -----> ----------- ----------- //heated // //heated // ------------------------------------------- /////////////////////////////////////////// -------------------------------------------
ENDDIS
!!!! Note the use of # in place of load( ) and of unigrid, which saves writing grdpwr commands UIN=1.0; qin=1.e5 nx=6*unit; xulast=0.03 ; ny= 3*unit ; yvlast = 0.006 yvlast=yvlast*2.0 #unigrid
GROUP 7. Variables stored, solved & named SOLVE(P1,U1,V1,TEM1);SOLUTN(P1,Y,Y,Y,N,N,N) STORE(ENUT,PRPS,WDIS,WGAP)
!!!! LVEL is here the default model turmod(lvel)

IF(:MOD:.EQ.LB) THEN + STORE(FMU,FTWO) + TURMOD(KEMODL-LOWRE-YAP);STORE(FONE,REYT) + KELIN=3 + RELAX(KE,LINRLX,0.25);RELAX(EP,LINRLX,0.25) ELSE + IF(:MOD:.EQ.2L) THEN + TURMOD(KEMODL-2L) + KELIN=3;RELAX(KE,LINRLX,0.25);RELAX(EP,LINRLX,0.25) + ELSE + IF(:MOD:.EQ.LV) THEN + TURMOD(LVEL) + ENDIF + ENDIF ENDIF MOD SOLVED STORED
GROUP 8. Terms (in differential equations) & devices TERMS(TEM1,N,P,P,P,P,P) egwf=t GROUP 9. Properties of the medium (or media) !!!! fluidmat and solidmat are always-declared and -set character variables, enabling materials to be introduced by name
# could have been used in place of l( ,which does not now need a closing bracket
l(fluidmat
GROUP 11. Initialization of variable or porosity fields
FIINIT(PRPS) = air20
l(solidmat
iniadd=f PATCH(BLOCK1,INIVAL,1,2*UNIT,UNIT+1,UNIT*2,1,1,1,1) COVAL(BLOCK1,PRPS,0.0,copper)
PATCH(BLOCK2,INIVAL,4*UNIT+1,6*UNIT,UNIT+1,UNIT*2,1,1,1,1) COVAL(BLOCK2,PRPS,0.0,copper)
PATCH(BASE,INIVAL,1,UNIT*6,1,UNIT,1,1,1,1) COVAL(BASE,PRPS,0.0,alumin)

GROUP 13. Boundary conditions and special sources
UIN=reyno * 1.e-5 /( 0.3333 * YVLAST ) TKEIN=0.001*UIN*UIN; EPSIN=0.1643*TKEIN**1.5/(0.1*YVLAST) FIINIT(U1)=UIN; FIINIT(V1)=0.0

PATCH(INLET,WEST,1,1,2*UNIT+1,NY,1,1,1,1) COVAL(INLET,P1,FIXFLU,1.189*UIN); COVAL(INLET,U1,ONLYMS,UIN) COVAL(INLET,TEM1,ONLYMS,0.0)
IF(:MOD:.EQ.LB.OR.:MOD:.EQ.2L) THEN FIINIT(EP)=EPSIN; FIINIT(KE)=0.04*UIN*UIN COVAL(INLET,KE,ONLYMS,TKEIN); COVAL(INLET,EP,ONLYMS,EPSIN) ENDIF
PATCH(OUTLET,EAST,NX,NX,1,NY,1,1,1,1);COVAL(OUTLET,P1,1.0E+05,0.0) COVAL(OUTLET,U1,ONLYMS,0.0);COVAL(OUTLET,V1,ONLYMS,0.0) COVAL(OUTLET,KE,ONLYMS,0.0);COVAL(OUTLET,EP,ONLYMS,0.0) COVAL(OUTLET,tem1,ONLYMS,SAME)
PATCH(TOP,NWALL,1,NX,NY,NY,1,1,1,1) COVAL(TOP,U1,1.0,0.0)
PATCH(HEAT1,VOLUME,1,UNIT*2,UNIT+1,UNIT*2,1,1,1,1) COVAL(HEAT1,TEM1,FIXFLU,QIN)
PATCH(HEAT2,VOLUME,4*UNIT+1,6*UNIT,UNIT+1,UNIT*2,1,1,1,1) COVAL(HEAT2,TEM1,FIXFLU,QIN)
GROUP 15. Termination of sweeps LSWEEP=1000;TSTSWP=-1 GROUP 16. Termination of iterations LITER(U1)=2;LITER(V1)=2 GROUP 17. Under-relaxation devices REAL(DTF);DTF=0.1*xulast/UIN relax(p1,linrlx,0.5) RELAX(U1,FALSDT,DTF); RELAX(V1,FALSDT,0.1*DTF) relax(tem1,linrlx,0.25); relax(enut,linrlx,0.1) varmax(enut)=0.1*uin*yvlast varmax(ke)=0.01*uin*uin
!!!! adddif, in its new role (see encyclopaedia) promotes convergence in low-Re duct flows adddif=t resfac=0.1
GROUP 21. Print-out of variables ixprf=nxs+1;ixprl=nxs+5;NXPRIN=1;NYPRIN=1;iyprl=10 patch(prof,profil,nxs+3,nxs+3,1,nys,1,1,1,1) coval(prof,tem1,0.0,0.0);coval(prof,u1,0.0,0.0) coval(prof,wdis,0.0,0.0);coval(prof,wgap,0.0,0.0) GROUP 22. Monitor print-out IXMON=NX-1;IPLTL=2000;IYMON=1;NPRMON=10000;YPLS=T LIBREF=208 echo=f
!!!! minimal print-out is obligatory, when 60 runs are to be performed, lest RESULT becomes enormous TEXT(mod = :mod: reyno= :reyno: unit=:unit: output(p1,n,n,n,n,n,n) output(u1,n,n,n,n,n,n) output(v1,n,n,n,n,n,n) output(prps,n,n,n,n,n,n) output(ltls,n,n,n,n,n,n) output(wdis,n,n,n,n,n,n) output(wgap,n,n,n,n,n,n) output(enut,n,n,n,n,n,n) if(.not.:mod:.eq.lv) then output(ke,n,n,n,n,n,n) output(fone,n,n,n,n,n,n) output(ftwo,n,n,n,n,n,n) output(ep,n,n,n,n,n,n) output(epke,n,n,n,n,n,n) output(reyt,n,n,n,n,n,n) output(fmu,n,n,n,n,n,n) endif

Appendix 3

The multi-run q1, with reyno (ie Reynolds number), mod (ie turbulence model) and unit (ie grid fineness) as parameters
Note the repeated use of incl(data0) and incl(data) talk=f;run(1,60) Re = 100 incl(data0) mod=lv; reyno=100; unit=2 incl(data) stop
incl(data0) mod=lb; reyno=100; unit=2 incl(data) stop incl(data0) mod=2l; reyno=100; unit=2 incl(data) stop
incl(data0) mod=lv; reyno=100; unit=5 incl(data) stop incl(data0) mod=lb; reyno=100; unit=5 incl(data) stop incl(data0) mod=2l; reyno=100; unit=5 incl(data) stop
incl(data0) mod=lv; reyno=100; unit=10 incl(data) stop incl(data0) mod=lb; reyno=100; unit=10 incl(data) stop incl(data0) mod=2l; reyno=100; unit=10 incl(data) stop
incl(data0) mod=lv; reyno=100; unit=20 incl(data) stop incl(data0) mod=lb; reyno=100; unit=20 incl(data) stop
incl(data0) mod=2l; reyno=100; unit=20 incl(data) stop
Re = 200
incl(data0) mod=lv; reyno=200; unit=2 incl(data) stop incl(data0) mod=lb; reyno=200; unit=2 incl(data) stop incl(data0) mod=2l; reyno=200; unit=2 incl(data) stop
incl(data0) mod=lv; reyno=200; unit=5 incl(data) stop incl(data0) mod=lb; reyno=200; unit=5 incl(data) stop incl(data0) mod=2l; reyno=200; unit=5 incl(data) stop
incl(data0) mod=lv; reyno=200; unit=10 incl(data) stop incl(data0) mod=lb; reyno=200; unit=10 incl(data) stop incl(data0) mod=2l; reyno=200; unit=10 incl(data) stop
incl(data0) mod=lv; reyno=200; unit=20 incl(data) stop incl(data0) mod=lb; reyno=200; unit=20 incl(data) stop incl(data0) mod=2l; reyno=200; unit=20 incl(data) stop
Re = 500
incl(data0) mod=lv; reyno=500; unit=2 incl(data) stop incl(data0) mod=lb; reyno=500; unit=2 incl(data) stop incl(data0) mod=2l; reyno=500; unit=2 incl(data) stop
incl(data0) mod=lv; reyno=500; unit=5 incl(data) stop incl(data0) mod=lb; reyno=500; unit=5 incl(data) stop incl(data0) mod=2l; reyno=500; unit=5 incl(data) stop
incl(data0) mod=lv; reyno=500; unit=10 incl(data) stop incl(data0) mod=lb; reyno=500; unit=10 incl(data) stop incl(data0) mod=2l; reyno=500; unit=10 incl(data) stop
incl(data0) mod=lv; reyno=500; unit=20 incl(data) stop incl(data0) mod=lb; reyno=500; unit=20 incl(data) stop incl(data0) mod=2l; reyno=500; unit=20 incl(data) stop
Re = 1000
incl(data0) mod=lv; reyno=1000; unit=2 incl(data) stop incl(data0) mod=lb; reyno=1000; unit=2 incl(data) stop incl(data0) mod=2l; reyno=1000; unit=2 incl(data) stop
incl(data0) mod=lv; reyno=1000; unit=5 incl(data) stop incl(data0) mod=lb; reyno=1000; unit=5 incl(data) stop incl(data0) mod=2l; reyno=1000; unit=5 incl(data) stop
incl(data0) mod=lv; reyno=1000; unit=10 incl(data) stop incl(data0) mod=lb; reyno=1000; unit=10 incl(data) stop incl(data0) mod=2l; reyno=1000; unit=10 incl(data) stop

incl(data0) mod=lv; reyno=1000; unit=20 incl(data) stop incl(data0) mod=lb; reyno=1000; unit=20 incl(data) stop incl(data0) mod=2l; reyno=1000; unit=20 incl(data) stop

Re = 2000
incl(data0) mod=lv; reyno=2000; unit=2 incl(data) stop incl(data0) mod=lb; reyno=2000; unit=2 incl(data) stop incl(data0) mod=2l; reyno=2000; unit=2 incl(data) stop
incl(data0) mod=lv; reyno=2000; unit=5 incl(data) stop incl(data0) mod=lb; reyno=2000; unit=5 incl(data) stop incl(data0) mod=2l; reyno=2000; unit=5 incl(data) stop
incl(data0) mod=lv; reyno=2000; unit=10 incl(data) stop incl(data0) mod=lb; reyno=2000; unit=10 incl(data) stop incl(data0) mod=2l; reyno=2000; unit=10 incl(data) stop
incl(data0) mod=lv; reyno=2000; unit=20 incl(data) varmax(wdis)=1.e-5 lsweep=1 output(prps,y,y,y,y,y,y) stop incl(data0) mod=lb; reyno=2000; unit=20 incl(data) stop incl(data0) mod=2l; reyno=2000; unit=20 incl(data) stop