Main Subroutines of FASTEST/Boundary Conditions

Subroutine bcin_exp

Description/Structure:

The subroutine bcin_exp sets velocity, density and viscosity at the inlet, outlet and on the walls in case of explicit time discretization. At the inlets, bcin2_exp is called. This subroutine sets $u$, $v$, $w$, density and viscosity values taking from the user defined boundary conditions and properties subroutines funcbcs and funcprp.

At the outlets, density and viscosity values are copied from the neighbour control volume in vertical direction to the outlet face.

In case of walls, the velocities on the boundaries are taken from the "wall velocities" section of funcbcs. Viscosity is taken from funcprp.

Subroutine bpres

Description/Structure:

The subroutine bpres is called from the subroutines fmg3d, calcdp, inter and restr, and it calls setind which sets indices of actual block.

This subroutine sets boundary pressures by extrapolating pressure from interior points.

In case of "static pressure exit boundary condition", outlet pressure is taken from fpstat psrt of subroutine funcbcs.

In "fluid-structure interaction" problems, subroutine bpresedge is called. bpresedge sets pressures on edges and corners by extrapolating pressures from interior wall points.

Subroutine btau

Description/Structure:

ONLY USED IN CASE OF FLUID-STRUCTUR INTERATCION!

The subroutine btau sets wall shear stresses. For each wall, tau_face is called. It calculates wall shear stress and $y$+ values for a given node. For calculating these on edges and corners, btauedge is called. This subroutine extrapolates the values from interior wall points.

In tau_face, first the angles between global $x$-, $y$-, $z$-directions and the normal vector to the face are calculated. Using these angles, paralel component of the velocity to the wall is calculated ($up$, $vp$ and $wp$). Distance between the first cell in vertical direction and the wall ($deln$) is found. The coefficient $tcoef$ which is used for the calculation of shear stresses is calculated in terms of viscosity and $deln$. At the end, $y^^++$ is calculated and the shear stress on the walls are found.

Subroutine funcbcs

Description/Structure:

funcbcs is the user defined boundary conditions subroutine.

This subroutine includes several functions defining some variables on the boundaries.

• functions fu, fv, fw:

  These functions defines the inlet velocity at x-, y- and z-directions respectively.

• function fexbc:

  This function defines the convective velocity on the outlet. It is only needed in case of "convective outlet condition".

• functions fte, fed:

  Turbulent kinetic energy and dissipation of turbulent kinetic energy at the inlet boundary.

• functions frsc[xx,yy,zz,xy,yz,zx]:

  These functions are used to define reynolds stresses in xx, yy, zz, xy, yz and zx planes at the inlet boundary.

• function ft:

  Defines temperature at the inlet.

• function fconc:

  Defines mass concentration of the species (icphi) at the inlet

• function fphf:

  Used to define the value of the phase function at the inlet.

• function fpstat:

  In case of static pressure boundary condition, static pressure at the outlet is defined here.

• functions f[u,v,w]w:

  Wall velocities in x-, y- and z-directions.

• functions f[t,q,htc]w, ftsht:

  .Temperature (K) at the wall in case of temperature bc SOL;

  .Heat flux through the wall in case of temperature bc FLX; (+ if flux from outside, - if flux to outside the computation area)

  .Heat transfer coefficient on the outside of the computation area for temperature bc HTC;

  .Surrounding temperature for heat transfer on the outside of the computation area for temperature bc HTC.

• functions femisk, ftransk, ftsrad:

  In case of radiation,

  .Emisivity coefficient of the radiating wall;

  .Transmisivity coefficient of the radiating wall;

  .Surrounding temperature for heat transfer on the outside of the computation area for radiation

• functions flbc, fcwa:

  -Concentration boundary condition in case of:

  .walls with temperature

  .walls with heat flux

  .wall with reactions

  -Dependence of the lbc to the:

  .concentration on the wall

  .flux through the wall

  .reaction speed on the wall

• function ftsurr:

  It is needed for the radiation bc on free surfaces. ftsurr is the temperature of the illuminated surrounding media. (free surface is radiating to the surrounding media)

• function fmixt, fmixtvar:

  Variables for combustion progress.

• functions fsrs, fsrscorx, fsrscory, fsrscorz:

  Single reactive scalar, and its -u, -v and -w correlations at the inlet boundary.

Subroutine modc

Description/Structure:

The subroutine modc is called from the subroutines calcc and prepstep, and sets the boundary conditions for concentration for the inlet, symmetry and for wall boundaries.

a. Inlet Boundary Calculations for Concentration: The subroutine inlbcc3 serves as follows:

1. Computation of start/end indices of patch in local $i$, $j$ direction
2. Retrieve block and face number
3. Get orientation of face
4. Loop over all CVs on current patch
5. Compute index of CV on the inlet plane
6. Compute index of CV adjacent to the boundary
7. Compute indices of surrounding CVs
8. Calculates area of the boundary face
9. Computation of inlet diffusion coefficient
10. Sum convective and diffusive fluxes
11. Projection of vectors 'ar2xyz' and 'ar3xyz' onto 'ar1xyz'
12. Interpolate concentrations to faces
13. Maintain diagonal dominance
14. Add contributions to source terms and central coefficients and set $A$ coefficient

b. Symmetry Boundary Calculations for Concentration: The subroutine symt3 serves as follows:

1. Computation of start/end indices of patch in local $i$, $j$ direction
2. Retrieve block and face number
3. Get orientation of face
4. Loop over all CVs on current patch
5. Compute index of CV on the inlet plane
6. calls subroutine symt

c. Wall Boundary Calculations for Concentration: The subroutine wallcbc3 expand patch of wall type for concentration and serves as in the above case.

Subroutine modke

Description/Structure:

This subroutine is called from the calcke.

The treatment of boundary conditions for the $k - \varepsilon$ equations differs from the usual approach for scalars. The differences are detailed and explained below. In the following we describe the steps to build up the coefficients ($A$ matrix) ("source term"):

• Compute coefficients of $A$ matrix and contributions to the source term from the convective and diffusive fluxes (inside the whole computational domain, except for boundaries)
  1. the coefficients of the A matrix contain the convective fluxes (UDS) and the diffusive fluxes (CDS).
  2. the source term contains contributions from the convective fluxes due to the deferred correction and from the diffusive fluxes in case of non-orthogonality of the grid.
  3. the matrix coefficients are stored in $aw$, $as$, $ab$, $at$, $an$, and $ae$. The source term is stored in $bt$. done in subroutines celke and celnbke.
• Compute generation and destruction terms for $k - \varepsilon$ equation (not on boundaries)
  1. the destruction term is always added to 'sp' which in turn is later added up (together with the other off-diagonal elements of matrix A) to give the diagonal element 'ap'.
  2. the generation term for the 'k' equation consists of two parts: one due to the mean strain rates (stored in 'apu') and the other due to buoyancy effects (stored in 'apv'). These two parts together with the source term of the convective/diffusive fluxes is summed up to give the source term $su$
     $$su(inp)=bt(inp)+apu(inp)+apv(inp)$$

  3. the generation term of the $\varepsilon$ equation is added to $bt$ to give the source term $su$.
  4. done in subroutine calcke

• Boundary conditions are often implemented by adding an extra source term to the equations of the control volumes next to the boundary. This extra source term is added to the value of $su$ already computed. The inflow boundary condition (with mass flux prescribed) is treated this way whereas the wall boundary condition in the $k - \varepsilon$ equation differs from this scheme.

  1. Calculate the source term for the wall boundary condition. Here the 'law of the wall' is used to calculate the source term which is stored in $suu$. Contrary to other boundary conditions this source term is not an extra source term but it replaces the generation term in the $k$ equation. Therefore, you should not add the source term 'suu' to the already computed value of the source term $su$ (like it is done for the inflow boundary condition) since it contains already the generation term due to the main strain rates and bouyancy effects (stored in $apu$ and $apv$). For control volumes adjacent to the wall the source term 'su' has to be computed in the following manner wallkbc3

     $$su(inp)=bt(inp)+suu$$

     It can be easily verified that this procedure will overwrite the source terms 'suu' from all other boundary conditions of the same control volume if they happen to have a lower patch number $ipa$ than the wall boundary condition. It is therefore essential that there is first a loop over all patches to work on all wall boundary conditions and then a second loop to do the remaining boundary conditions. Thus it is ensured that the standard generation terms do not enter the source term in control volumes adjacent to a wall but other boundary conditions can still contribute their share to the source term. This procedure is implemented in this subroutine modke.

     However, it should be mentioned that there remains one unresolved issue: If a control volume is adjacent to more than one wall, the different source terms $suu$ will not add up but the wall with the highest patch number will determine the value of the source term. It is correct that the source terms do not add up since they are 'volumetric' source terms (like the $k$ generation term). But since the source term is computed from the parallel velocity $up$ and the wall shear stress $tcoef$ it is not clear which of the different values would be the best one (or which mean one would have to choose).

  2. Calculate the source term for inflow boundary condition. Since all values are known at the inflow plane (except the pressure) the fluxes can be computed.
  3. the source term contains contributions from convective plus diffusive fluxes times the value of the scalar at the inflow plane. An extra term is added from the diffusive fluxes in case of non-orthogonality of the grid
  4. the source term is stored in $uu$ and added to $su$
     $$su(inp)=su(inp)+suu$$

  5. done in subroutine inlbcke and inlbcke3

  6. The symmetry boundary condition contains just an extrapolation of the value of the scalar in the cell to the boundary point. The source term remains unchanged.

  7. Solve system of equation using SIP solver.

Subroutine modpr

Description/Structure:

The subroutine modpr sets the outlet boundary conditions for the static pressure. This subroutine is called from the iprops which provides modifications of fluid properties if these are changing during the solution process (i.e changes due to the temperature dependence of several fluid properties) and calls outbcpstat3 which serves for the calculation of the followings:

1. Computation of start/end indices of patch in local $i$, $j$ direction
2. Retrieve block and face number
3. Get orientation of face
4. Loop over all CVs on current patch
5. Compute index of CV on the outlet plane
6. Compute index of CV adjacent to the outlet plane
7. Extrapolated reverse coefficient of momentum equations
8. Interpolation factor used for extrapolation
9. Extrapolation of velocities
10. Mass flux

Subroutine modrsm

Description/Structure:

The subroutine modrsm serves for the calculation of boundary conditions for reynolds stress components and called from the subroutines caledrsm and calsrs. It calls the subroutines walledrsm and inlbcrsm

modrsm does the followings :

1. Computation of start/end indices of patch in local $i$, $j$ direction
2. Retrieve block and face number
3. Loop over all CVs on current patch
4. Compute index of CV on the boundary
5. Get orientation of face

walledrsm calculates the followings additional to the above:

1. Vector connecting the cell center in the face and the cell center of the CV adjacent to the boundary face
2. Area of the boundary face
3. Components of the unit normal vector
4. Compute normal distance by scalar product
5. Turbulent energy dissipation
6. Set source term, contribution to central coefficient, and $A$ coefficient

inlbcrsm calculates the inflow boundary turbulence:

1. Computation of start/end indices of patch in local $i$, $j$ direction
2. Retrieve block and face number
3. Loop over all CVs on current patch
4. Compute index of CV on the inlet plane
5. Compute index of CV adjacent to the boundary
6. Compute indices of surrounding CVs
7. Area of the boundary face
8. Inlet laminar viscosity
9. Sum convective and diffusive fluxes
10. Projection of vectors 'ar2xyz' and 'ar3xyz' onto 'ar1xyz'
11. Interpolate variable $fi$ and $tij$ to faces
12. Maintain diagonal dominance
13. Add contributions to source terms and central coefficients and set A coefficient to zero

Subroutine modspec

Description/Structure:

The subroutine modspec calculates the heat conductivity at boundaries between fluid regions and called from the calmef (makes the effective heat conductivity available to the boundary array for w,e,s,n,b,t)

Subroutine modsrs

Description/Structure:

The subroutine modsrs calculates the inlet boundary conditions for the single reactive scalar and called from calsrscor and calsrs, and the solution procedure is like in modrsm

Subroutine modsrscor

Description/Structure:

The subroutine modsrscor calculates the boundary conditions for the single reactive scalar/velocity correlation for the inlet (inlbcsrs), wall (wallsrscor and symmetry (symsrscor) cases.

Subroutine modt

Description/Structure:

The subroutine modt is for the temperature boundary conditions. This subroutine is called from the calct and restr. The boundary condition variables are introduced here as lb-t-[w,e,n,s,t,b]. They include:

• 1- inlet
• 2- outlet (defined again to provide confusion in the grid-input-file)
• 3- symmetry or adiabatic wall
• 4- wall with given temperature
• 5- wall with given heatfluxdensity.
• 6-
• 7-
• 9- wall with heat transfer coeff. and surr. temp.

• For -5 and -9: the symmetry routine is called to interpolate the temperature onto the normal to the wall. The temperature difference is then provided by adding $qw/lam*dx/area$ to the wall temperature.

a. Inlets: ifac (= 1 for w, s, b and = -1 for e, n, t) is transfered from calling routine. It is used as sign for the flux and to calculate the "displacement" in array when the flux is calculated; eg. for w the flux in i=1 is used, whereas in e, the flux in i=nim is used.

CAUTION: In all subroutines for the boundary conditions ar1x, ar1y and ar1z contain the metric normal to the boundary plane. In the inner field these variables contain the derivatives $d( ksi)/d(xi)$ (with xi=x,y,z).

b.Symmetry: look from modc

c. Wall boundary conditions: look from the walltbc3

Subroutine moduvw_exp

Description/Structure:

The subroutine moduvw_exp is called from the runge-kutta and sets the inlet inlbc3exp, symmetry symuvw3exp and wall wallbc3exp boundary conditions with explicit time discretization.

Subroutine outbc_exp

Description/Structure:

In case of explicit time discretization, in exptimedis, runge_kutta subroutine is called. runge_kutta uses outbc_exp for calculating outlet boundary condition.

For each face, om3 is called. It calculates the outlet velocities and outlet flux. First, start/end indices of patch in local i-j plane is calculated. Then block number, face number, face orientation and auxilliary variables for patch orientation are found. Necessary cv indices on the outlet, adjacent to the boundary node are calculated.

If there is "convective exit boundary condition" then indices for the surrounding cvs are calculated. Velocities are interpolated to the faces and the convective velocities for the boundary condition is set ($u_{conv}$). After that the velocity gradients for $u_{conv}$ are calculated and using it, the new velocities are found. At the end the entire mass flux for the outlet plane is calculated.

If there is zero-gradient exit boundary condition, then the velocities for zero-gradient bounadry is set ($u$(adjacent)=$u$(outlet)). After that the entire mass flux for the outlet plane is calculated.

Subroutine outbcmix

Description/Structure:

The subroutine outbcmix is called from the fmg3d and it calculates the outlet boundary for combustion progress variables.

Subroutine setbv

Description/Structure:

The subroutine setbv initializes the variables on the boundaries.

a. Initialization of the Inlets:

At the inlet faces, bvin* (* : vel, turb, temp, rsm, srs, srscor, conc, mix, phf) is/are called depending on the problem and model type. This subroutine takes the initial variable values from the user defined boundary condition subroutine funcbcs. For example, in case of a temperature dependent $k - \varepsilon$ turbulent model, bvinvel, bvinturb and bvintemp are called. bvinvel initializes $u$, $v$ and $w$ values as $fu$, $fv$ and $fw$ respectively. bvintemp takes the inlet temperature as $ft$, and bvinturb initializes the turbulent kinetic energy ($te$) and energy dissipation ($ed$) at the inlet as $fte$ and $fed$. Variables $fu$, $fv$, $fw$, $ft$, $fte$ and $fed$ are defined in funcbcs.

b. Initialization of the Outlets:

At the outlets, bvout* (* : vel, turb, temp, rsm, srs, srscor, conc, mix, phf) is/are called depending on the problem and model type. bvout* type of subroutines are used to correct the variables on the outlet. Control volume indices on the outlet plane (inh) and adjacent to the boundary (inp) are calculated. For the above example, bvoutvel, bvoutturb and bvouttemp are called. These subroutines copy the corresponding variables from the cv adjacent to the boundary, to the cv on the outlet boundary.

c. Initialization of the Symmetric Boundary Planes:

In order to correct the variables on the symmetic boundaries, bvsym* (* : vel, turb, temp, rsm, srs, srscor, conc, mix, phf) is/are called depending on the problem and model type. For correcting temperature, turbulent kinetic energy, energy dissipation and concentration subroutine symscal (symmetry boundary conditions for scalar) is used.

By using the assumption from Peric, the scalar values on the face are changed during the iteration process. This influence the convergence when multigrid is used, because the values are changed during each restriction cycle. In order to overcome this, the best solution would be to find the point where the normal from symmerty plane intersects the area between local NW, NE, SW and SE and then interpolate the right value.

In subroutine symscal, first start and end indices of the boundary patch is found. After setting the necessary neighbor control volume indices, variable value on the local east, west, north and south points (according to the cv adjacent to the inlet plane cv) are found. Grid points are calculated on the boundary. Necessary lenghts are calculated to find the discretization function. Then the variable value on the boundary cv $\phi_B$ is found in terms of variable values in adjacent cv $\phi_P$ and its neighbors ($\phi_E$, $\phi_W$, $\phi_N$ and $\phi_S$).

d. Initialization of the Walls:

On the walls, bvwal* (* : vel, turb, rsm, srscor, phf) is/are called depending on the problem and model type. For the example of a temperature dependent $k - \varepsilon$ turbulent model, bvwalvel assigns wall velocities read from funcbcs on the walls, bvwalturb sets turbulent kinetic energy ($te$) and energy dissipation ($ed$) to zero, and bvwaltemp sets the wall temperature defined in funcbcs.

For initialization of the single reactive scalar and combustion variables on the wall, bvsym* subroutines are used.

Subroutine setradbcv

Description/Structure:

Sets radiation parameters (emisivity, transmissivity and surrounding temperature) on boundary planes. setradw is called. This subroutine assigns values of radiation parameters according to the boundary type of the block.

• Inlet, Outlet:

  Inner emisivity = 0.

  Outer emisivity = 0.

  Transmisivity = 1.

  Surrounding temperature = ftsrad (defined surr. temp. in funcbcs)

• Block : Solid; Neighbour : Fluid:

  Inner emisivity = 0.

  Outer emisivity = femisk (defined emisivity coeff. of the wall in funcbcs)

  Transmisivity = ftransk (defined transmisivity coeff. of the wall in funcbcs)

  Surrounding temperature = ftsrad (defined surr. temp. in funcbcs)

• Block : Fluid; Neighbour : Solid:

  Inner emisivity = femisk (defined emisivity coeff. of the wall in funcbcs)

  Outer emisivity = 0.

  Transmisivity = ftransk (defined transmisivity coeff. of the wall in funcbcs)

  Surrounding temperature = ftsrad (defined surr. temp. in funcbcs)

• Both blocks are fluid or wall with temperature:

  Inner emisivity = femisk (defined emisivity coeff. of the wall in funcbcs)

  Outer emisivity = femisk (defined emisivity coeff. of the wall in funcbcs)

  Transmisivity = ftransk (defined transmisivity coeff. of the wall in funcbcs)

  Surrounding temperature = ftsrad (defined surr. temp. in funcbcs)

Subroutine symvel

Description/Structure:

In case of symmetric boundary condition, symvec is called. It calculates the cartesian components of the parallel velocity on the boundary.

Subroutine wallfluxbc

Description/Structure:

This subroutine is called in case of "fluid-structure interaction" problems with "moving grids".

After calwallvel is called (Calculation of wall velocities on a moving or distorted structure), wallfluxbc takes place.

Wall flux at each cell face is calculated in terms of wall velocity (coming from calwallvel) and the fluid density. After that these fluxes are summed up to total flux of the control volume.