chapter. Boundary Conditions 第8章 边界条件[cfd-0009]

Before we proceed with the discussion of the boundary conditions, we should mention the concept of dummy cells or dummy points. This approach is very popular on structured grids. However, dummy cells offer some advantages also on unstructured grids. The dummy cells are additional layers of grid cells or points outside the physical domain. This is sketched in Fig. 8.1 for the case of a 2-D structured grid. As we can see, the whole computational domain is surrounded by two layers of dummy cells (denoted by dashed line). The dummy cells (points) are usually not generated as the grid inside the domain (except on interfaces between grid blocks). Rather, the cells are only virtual, although geometrical quantities like volumes or face vectors are associated with them.

The purpose of the dummy cells is to simplify the computation of the fluxes, gradients, dissipation, etc. along the boundaries. This is achieved by the possibility to extend the stencil of the spatial discretisation scheme beyond the physical boundaries. As we can see in Fig. 8.1, the same discretisation scheme can be employed at the boundaries as inside the physical domain. Thus, we can solve the governing equations in the same way for all "physical" grid points. This makes the discretisation schemes much easier to implement. Furthermore, all grid points of a structured grid can be accessed in a single loop, which is of significant advantage particularly on vector computers. The condition is of course that the dummy cells (points) contain appropriate values of the conservative variables as well as of the geometrical quantities. Clearly, the number of dummy cell layers must be such that the part of the stencil outside the physical domain is completely covered. The conservative variables in the dummy cells (points) are obtained from boundary conditions. The geometrical quantities are usually taken from the corresponding control volume at the boundary. In the case of boundaries between multiple grid blocks (like in Fig. 3.4), all flow variables and the geometry are transferred from the neighbouring block.

The grey-shaded dummy cells (and the associated points) in Fig. 8.1 represent a certain problem, since it is not quite clear how to set their values if there is no adjacent grid block. Their values are not required by the standard cross-type discretisation stencil. However, the cells (points) are important for the computation of gradients (viscous fluxes - see Section 4.4), or for the transfer operators within multigrid (Section 9.4). The simplest way to solve the problem is to average the values from the adjacent "regular" dummy cells, as indicated in Fig. 8.1 by arrows. However, this does not work satisfactorily for wall or symmetry boundaries. In these cases, it is better to extend the physical boundary into the dummy cell layers (see the structured 2-D code on the CD-ROM).

The numerical simulation of external flows past airfoils, wings, cars and other configurations has to be conducted within a bounded domain. For this reason, artificial farfield boundary conditions become necessary. The numerical implementation of the farfield boundary conditions has to fulfil two basic requirements. First, the truncation of the domain should have no notable effects on the flow solution as compared to the infinite domain. Second, any outgoing disturbances must not be reflected back into the flow field [9]. Due to their elliptic nature, sub- and transonic flow problems are particularly sensitive to the farfield boundary conditions. An inadequate implementation can lead to a significant slow down of convergence to the steady state. Furthermore, the accuracy of the solution is likely to be negatively influenced. Various methodologies were developed, which are capable of absorbing the outgoing waves at the artificial boundaries [10]-[15]. An review of different non-reflecting boundary conditions can be found in [16].

In the following two subsections, we shall discuss the concept of characteristic variables as it was described by Whitfield and Janus [12]. We shall also present an extension of the farfield boundary conditions for lifting bodies.

Various approaches were devised for the implementation of numerical inlet, and in particular, of outlet (also named open) boundary conditions for the Navier-Stokes equations [22]-[26]. Here, we will concentrate on methodologies, which were developed for turbomachinery applications. Suitable non-reflecting inlet and outlet boundary conditions were described, e.g., in [27]-[30]. Giles [31], and Hirsch and Verhoff [32] suggested non-reflecting boundary conditions for the Euler equations, which are intended for domains with a short distance between the body and the inlet or the outlet plane.

In certain cases, the inlet and outlet boundary are additionally periodic with respect to the velocity as well as the pressure and temperature gradient. This type of flow is encountered, for example, in the simulation of heat exchangers [33]. The implementation of periodic inlet and outlet boundary conditions was presented in [34]-[36] for LES in channels.

A common procedure consists of the specification of the total pressure, total temperature, and of two flow angles. One characteristic variable has to be interpolated from the interior of the flow domain. One possibility is to employ the outgoing Riemann invariant [30], which is defined as

where the index \(d\) denotes the state inside the domain (cf. Fig. 8.6a). The Riemann invariant is used to determine either the absolute velocity or the the speed of sound at the boundary. In practice, it was found that selecting the speed of sound leads to a more stable scheme, particularly for low Mach-number flows. Therefore, we set

with \(\theta\) being the flow angle relative to the boundary, and \(c_0\) denoting the stagnation speed of sound. Hence,

and

where \(\|\vec{v}_d\|_2\) denotes the total velocity at the interior point \(d\) (Fig. 8.6a). The unit normal vector \(\vec{n}\) in Eq. (8.34) was assumed to point out of the domain.

Quantities like the static temperature, pressure, density, or the absolute velocity at the boundary are evaluated as follows

where \(T_0\) and \(p_0\) are the given values of total temperature and pressure, \(R\) and \(c_p\) represent the specific gas constant and the heat coefficient at constant pressure, respectively. The velocity components at the inlet are obtained by decomposing \(\|\vec{v}_b\|_2\) according to the two (one in 2D) prescribed flow angles.

In turbomachinery, the static pressure is usually prescribed at the outlet. The subsonic outlet boundary can be treated in a way quite similar to the outflow condition in Eq. (8.23). Only the ambient pressure \(p_a\) is replaced here by the given static exit pressure.

Flow variables in the dummy cells (points) can be obtained by linearly extrapolating the states at the boundary and at the interior point \(d\).

The injection velocity and other flow variables are computed from the given mass flow rate \(\dot{m}\) and the injection temperature \(T_{inj}\). We will assume here that the mass is injected perpendicular to the boundary. The convective fluxes are evaluated using Eq. (2.21) with the contravariant velocity \(V\) set equal to the injection velocity, i.e.,

Velocity components are evaluated using the face normal vector as

(see Fig. 8.8). The density \(\rho_{inj}\) in Eq. (8.37) is a function of \(T_{inj}\) and of the pressure at the boundary \(p_b\). Hence, in the case of perfect gas we have

For the cell-centred scheme sketched in Fig. 8.8a, \(p_b\) can be set identical to the pressure in the boundary cell, i.e., \(p_b=p_2\). The median-dual cell-vertex scheme requires an interpolation of \(p_b\) from the points defining the particular face area (i.e., from \((i,2)\) and \((i+1,2)\) in Fig. 8.8b). Formulae like (8.8), (8.11)-(8.13) can be utilised for this purpose. Values in the dummy cells (points) are obtained by 0-th order or by linear extrapolation from the interior.

If the flow is to be symmetrical with respect to a line or a plane, the first condition which must be met is that there is no flux across the boundary. This is equivalent to the requirement that the velocity normal to the symmetry boundary is zero. Furthermore, the following gradients have to vanish:

• gradient of a scalar quantity normal to the boundary,
• gradient of the tangential velocity normal to the boundary,
• gradient of the normal velocity along the boundary (since \(\vec{v}\cdot\vec{n}=0\)).

We can write these conditions as

where \(U\) stands for a scalar variable and \(\vec{t}\) denotes a vector tangential to the symmetry boundary.

The implementation of the symmetry boundary condition can be largely simplified by employing dummy cells. The flow variables in the dummy cells are obtained using the concept of reflected cells. This means that scalar quantities like density or pressure in the dummy cells are set equal to the values in the opposite interior cells, i.e.,

The notation corresponds to that in Fig. 8.2. The velocity components are reflected with respect to the boundary as indicated in Eq. (8.10). The normal gradient of the normal velocity in the dummy cell equals to that in the opposite interior cell, but it has a reversed sign.

Two different approaches can be followed. One possibility is to construct the missing half of the control volume by mirroring the grid on the boundary. The fluxes and the gradients are then evaluated like in the interior using reflected flow variables (see above). The second methodology computes the fluxes for the halved control volume (but not across the boundary). The components of the residual normal to the symmetry plane are then zeroed out. It is also necessary to correct normal vectors of those faces of the control volume, which touch the boundary (like at point \(2^{*}\) in Fig. 8.4). The modification consists of removing all components of the face vector, which are normal to the symmetry plane. The gradients have also to be corrected according to Eq. (8.40).

This type of boundary condition is encountered only in the case of structured grids. The coordinate cut represents an artificial, not a physical, boundary. It is a line (plane in 3D) composed of grid points with different computational coordinate(s) but the same physical location. This means that the grid is folded such that it touches itself. As we shall see in Subsection 11.1.1, the coordinate cut appears for the so-called C- (Fig. 11.5) or O-grid topology (Fig. 11.9). The flow variables and their gradients have to stay continuous across the cut.

The best way to implement the cut boundary condition is to employ dummy cells (points). The situation is sketched in Fig. 8.9. As we can see, the dummy layers here are not virtual, but they coincide with the grid on the opposite side of the cut. Hence, the values of physical quantities in the dummy cells (cell-centred scheme), or in the dummy points (cell-vertex scheme), are obtained directly from the opposite cells (points). In the case of the cell-centred scheme, the fluxes across the faces of the boundary cell (shaded in Fig. 8.9a) are evaluated exactly like in the interior field.

The cut boundary can be treated in two different ways for the cell-vertex scheme. One possibility is to generate a complete control volume at the cut (the second part is denoted by a dashed line in Fig. 8.9b). Using the dummy points, the fluxes can be calculated in the same way as inside the domain. If the implementation is done correctly, the flow quantities at the points 2 (upper grid part) and 5 (lower part) will be equal. The second approach is to integrate the fluxes separately for each half of the control volume. The residuals at the points 2 and 5 in Fig. 8.9b are then added. It is important that the partial control volumes at the points 2 and 5 are summed up as well.

There are certain practical applications where the flow field is periodic with respect to one or multiple coordinate directions. In such a case, it is sufficient to simulate the flow only within one of the repeating regions. The correct interaction with the remaining physical domain is enforced via periodic boundary conditions.

We can distinguish between two basic types of periodic boundaries. The first one covers translational periodicity. This means that one periodic boundary can be transformed into the other boundary by pure coordinate translation. The second type represents periodic boundaries, which were generated by coordinate rotation. Thus, we speak of rotational periodicity.

In the following, we shall describe the implementation of the periodic boundary conditions for the cell-centred and the cell-vertex scheme. We shall also consider the case of rotational periodicity. Further details of the treatment of periodic boundaries can be found in Refs. [37], [38].

The utilisation of the dummy-cells concept enables a simple implementation of the periodic boundary condition. Let us consider the example from turbomachinery in Fig. 8.10. The configuration is periodic in the vertical direction. The shaded cells 1 and 2 are located on the lower and the upper periodic boundary, respectively. Due to the periodicity condition, the first dummy-cell layer corresponds to the boundary cells at the opposite periodic boundary. The second dummy-cell layer communicates with the second layer of the physical cells and so on. Hence, all scalar quantities (density, pressure, etc.) in the dummy cells are obtained directly from the corresponding physical cells, i.e,

The same relations hold also for the vector quantities (velocity, gradients) in the case of translational periodicity. Rotational-periodic boundaries require a correction of the vector variables. This will be discussed further below.

This situation is sketched in Fig. 8.11. One approach for the treatment of periodic boundaries consists of the integration of the fluxes around the faces of the shaded control volumes. The residuals at the points 1 and 2 in Fig. 8.11 are then summed in order to obtain the complete net flux. Thus,

The partial control volumes at the points 1 and 2 (shaded in Fig. 8.11) have to be added up as well. In the case of translational periodicity, the residuals from the opposite boundary remain unchanged, i.e., \(\vec{R}_{1'}=\vec{R}_1\) and \(\vec{R}_{2'}=\vec{R}_2\). This results in \(\vec{R}_{1,sum}=\vec{R}_{2,sum}\). Rotationally periodic boundaries require a transformation of the momentum equations before Eq. (8.43) can be applied.

The rotational periodicity condition is based on a rotation of the coordinate system. Therefore, all vector quantities like velocity or gradients of scalars have to be transformed accordingly. Scalar quantities like pressure or density, which are invariant with respect to coordinate rotation, remain unchanged. If we assume the rotational axis is parallel to the \(x\)-axis (see Fig. 8.12), the rotation matrix becomes

where the angle \(\phi\) between the periodic boundaries \(A\) and \(B\) is positive in the clockwise direction. Hence, for example, the velocity vector transformed from boundary \(A\) to \(B\) (cells \(1'\), \(2'\) in Fig. 8.10 and points \(1'\)-\(4'\) in Fig. 8.11) reads

It is easy to show that the \(x\)-component of \(\vec{v}_A\) (i.e., \(u_A\)) is not changed by the rotation. Thus, \(u_B=u_A\). The gradients of all flow quantities are transformed in similar way.

As stated above, in the case of the cell-vertex scheme the residuals of the momentum equations must be corrected before the summation in Eq. (8.43) can take place. The application of the rotation matrix Eq. (8.44) leads to

The superscript \(u,v,w\) in Eq. (8.46) denotes the three momentum equations.

During the discussion of the spatial discretisation with structured grids in Section 3.1, it became evident that it is usually not possible to generate a single grid inside a geometrically complex domain (see Fig. 3.4). We mentioned two possible methodologies how to solve the problem. The first one was the multiblock approach and the second one was the Chimera technique. In the following, we shall describe the basic implementation issues of the multiblock approach. For more throughout discussion, the reader is referred to Refs. [39]-[45]. A very helpful introduction to the multiblock methodology was presented in [46]. Details of the Chimera technique, which is not treated here, can be found in Refs. [47]-[51].

Within the multiblock technique, the physical domain is split into a certain number of virtual parts. Consequently, the computational domain becomes also divided into the same number of blocks. In a general case, the physical solution in a particular block will depend on the flow in one or multiple neighbouring blocks. Therefore, we have to provide a data structure which allows for an efficient exchange of information between the blocks. The structure is also required for communication, if different processors are used to solve the governing equations in the grid blocks.

The first part of the data structure consists of the numbering of the block boundaries. One particular numbering scheme is displayed in Fig. 8.13. The numbering strategy in Fig. 8.13 can be summarised as follows:

It is important that all blocks employ the same numbering scheme. The indices \(i\), \(j\), \(k\) of the grid points in the computational space are defined in the ranges

The cell indices \(I\), \(J\), \(K\), which are required by the cell-centred scheme are defined in a similar way. Since the multiblock approach is usually implemented using dummy cells/points, the physical cells/points will have a certain offset from the start or the end of each range (see Fig. 8.1).

The boundary of each block is divided into a number of non-overlapping patches. This allows the specification of different boundary conditions on the same block boundary. The situation is depicted in Fig. 8.14. For a unique identification of each patch it is necessary to store the number of the corresponding block and the number of the block boundary. Furthermore, the origin, the height and the width of the patch must be stored. For this purpose, the coordinates \(L1BEG\), \(L1END\), \(L2BEG\) and \(L2END\) are used in Fig. 8.14. It is suggested to orient the coordinate system of the patch according to the cyclic directions. This means, that if we consider the \(i\)-coordinate, \(j\) and \(k\) will be the first and the second cyclic direction. In the case of the \(j\)-coordinate, the cyclic directions will become \(k\) and \(i\), respectively. Therefore, since the patch in Fig. 8.14 is on the \(j=JBEG\) boundary, the \(l_1\)-coordinate is oriented in the \(k\)-direction and \(l_2\) in the \(i\)-direction. The application of the cyclic directions allows for a unique definition of the orientation of each patch.

The remaining part of the data structure makes sure that data can be exchanged between those patches, which represent interfaces between the blocks (we assume here that the blocks communicate only across their faces). For this purpose, it is required to extend the above patch data structure by the numbers of the adjacent block and patch. It is furthermore necessary to code the orientation of the communicating patches with respect to each other.

The exchange of flow quantities between two blocks is sketched in Fig. 8.15. The procedure consists of two steps. In the first step, variables from the part of the domain, which is overlapped by the dummy layers of the adjacent patch are written to the own dummy cells/points or to a temporary storage (\(A'\) and \(B'\) in Fig. 8.15). This is done for all blocks. In the second step, the data in \(A'\) and \(B'\) is exchanged between both blocks. This means that \(A'\) is written to the dummy layers of block B and \(B'\) to the dummy layers of block A. If the two patches have a different orientation, the data must be transformed accordingly. In cases where the grid lines do not match at the block interface, further operations are required as described, e.g., in [52], [53].

We already stated in Subsection 5.3.4 that the evaluation of the flow gradients requires some care in the case of the median-dual scheme. If the gradients are calculated on triangular or tetrahedral grids using the Green-Gauss approach in Eq. (5.50), the contributions from the boundaries of the domain (except at symmetry or periodic boundaries) must be evaluated similar to Eq. (8.12) or Eq. (8.13), instead of the arithmetic average. Otherwise, the gradient will not be accurate. Considering the notation in Fig. 8.4, the contribution to the boundary node 1 reads

where \(\Delta S_{12}\) is the length of the boundary face between node 1 and 2 (therefore halved). Corresponding to Eq. (8.13), the contribution of the triangular face 1-2-3 to node 1 in Fig. 8.5 becomes

with \(\Delta S_{123}/3\) being the grey area in the triangle 1-2-3. On mixed grids, it is more appropriate to employ the least squares approach with virtual edges [8] (see Fig. 5.15).

The cell-centred scheme requires no special provisions at symmetry or periodic boundaries. The implementation is identical to that discussed for the fluxes in Section 8.6 or 8.8. This holds also for the median-dual scheme, if the gradients are evaluated using the least-squares approach. The only additional work required is to set certain gradients to zero as described previously in Section 8.6 (cf. Eq. (8.40)).

If the Green-Gauss approach is employed within the median-dual scheme (i.e., if Eq. (5.50) is applied), it is necessary to correct normal vectors of those faces of the control volume, which touch the boundary (like at point \(2^{*}\) in Fig. 8.4). This is done by setting all components of the face vector to zero, which are normal to the symmetry plane. Finally, the gradients are corrected as discussed in Section 8.6. At periodic boundaries, the gradients and the volumes from both sides of the boundary have to be summed up as presented in Section 8.8 for the fluxes (Eq. (8.43)). In the case of rotational periodicity, the gradients needs to be transformed by applying the rotation matrix in Eq. (8.44).