9. Supported Boundary Conditions
9.1. Inflow Boundary Condition
9.1.1. Continuity
Continuity uses a flux boundary condition with the incoming mass flow rate based on the user specified values for velocity,
As this is a vertex-based code, at inflow and Dirichlet wall boundary locations, the continuity equation uses the specified velocity within the inflow boundary condition block.
9.1.2. Momentum, Mixture Fraction, Enthalpy, Species, \(k_{sgs}\), k and \(\omega\)
These degree-of-freedoms (DOFs) each use a Dirichlet value with the specified user value. For all Dirichlet values, the row is zeroed with a unity placed on the diagonal. The residual is zeroed and set to the difference between the current value and user specified value.
9.2. Wall Boundary Conditions
9.2.1. Continuity
Continuity uses a no-op.
9.2.2. Momentum
When resolving the boundary layer, Momentum again uses a no-slip Dirichlet condition., e.g., \(u_i = 0\).
In the case of a wall model, a classic wall function is applied. The wall shear stress enters the discretization of the momentum equations by the term
Wall functions are used to prescribe the value of the wall shear stress rather than resolving the boundary layer within the near-wall domain. The fundamental momentum law of the wall formulation, assuming fully-developed turbulent flow near a no-slip wall, can be written as,
where \(u^+\) is defined by the the near-wall parallel velocity, \(u_{\|}\), normalized by the wall friction velocity, \(u_{\tau}\). The wall friction velocity is related to the turbulent kinetic energy by,
by assuming that the production and dissipation of turbulence is in local equilibrium. The wall friction velocity is also computed given the density and wall shear stress,
The normalized perpendicular distance from the point in question to the wall, \(y^+\), is defined as the following:
The classical law of the wall is as follows:
where \(\kappa\) is the von Karman constant and \(C\) is the dimensionless integration constant that varies based on authorship and surface roughness. The above expression can be re-written as,
or simplified to the following expression:
In the above equation, \(E\), is referred to in the text as the dimensionless wall roughness parameter and is described by,
In Nalu-Wind, \(\kappa\) is set to the value of 0.42 while the value of \(E\) is set to 9.8 for smooth walls (White suggests values of \(\kappa=0.41\) and \(E=7.768.\)). The viscous sublayer is assumed to extend to a value of \(y^+\) = 11.63.
The wall shear stress, \(\tau_w\), can be expressed as,
where \(\lambda_w\) is simply the grouping of the factors from the law of the wall. For values of \(y^+\) less than 11.63, the wall shear stress is given by,
The force imparted by the wall, for the \(i_{th}\) component of velocity, can be written as,
where \(A_w\) is the total area over which the shear stress acts.
The use of a general, non-orthogonal mesh adds a slight complexity to specifying the force imparted on the fluid by the wall. As shown in Equation (9.11), the velocity component parallel to the wall must be determined. Use of the unit normal vector, \(n_j\), provides an easy way to determine the parallel velocity component by the following standard vector projection:
Carrying out the projection of a general velocity, which is not necessarily parallel to the wall, yields the velocity vector parallel to the wall,
Note that the component that acts on the particular \(i^{th}\) component of velocity,
provides a form that can be potentially treated implicitly; i.e., in a way to augment the diagonal dominance of the central coefficient of the \(i^{th}\) component of velocity. The use of residual form adds a slight complexity to this implicit formulation only in that appropriate right-hand-side source terms must be added.
9.2.3. Mixture Fraction
If a value is specified for each quantity within the wall boundary condition block, a Dirichlet condition is applied. If no values are specified, a zero flux condition is applied.
9.2.4. Enthalpy
If the temperature is specified within the wall boundary condition block, a Dirichlet condition is always specified. Wall functions for enthalpy transport have not yet been implemented.
The simulation tool supports multi-physics coupling via conjugate heat transfer and radiative heat transfer. Coupling parameters required for the thermal boundary condition are post processed by the fluids or PMR Realm. For conjugate and radiative coupling, the thermal solve provides the surface temperature. From the surface temperature, a wall enthalpy is computed and used.
9.2.5. Thermal Heat Conduction
If a temperature is specified in the wall block, and the surface is not an interface condition, then a Dirichlet approach is used. If conjugate heat transfer is included, then the boundary condition applied is as follows,
where \(h\) is the heat transfer coefficient and \(T^o\) is the reference temperature. The details of how these quantities are computed are currently omitted in this manual. In general, the quantities are post processed from the fluids temperature field. A surface-based gradient is computed on the boundary face. Nodes on the face augment a heat transfer coefficient field while nodes off the face contribute to a reference temperature.
For radiative heat transfer, the boundary condition applied is as follows:
where \(H\) is again the irradiation provided by the RTE solve.
If no temperature is specified or an adiabatic line command is used, a zero flux condition is applied.
9.2.6. Species
If a value is specified for each quantity within the wall boundary condition block, a Dirichlet condition is applied. If no values are specified, a zero flux condition is applied.
9.3. Atmospheric Boundary Layer Surface Conditions
9.3.1. Monin-Obukhov Theory
Consider atmospheric flow over a flat but non-smooth surface; the coordinate system convention is that flow is along the \(x\)-axis, while the \(z\)-axis is oriented normal to the surface. The surface layer is the relatively thin layer near the surface where strong wind and temperature gradients exist. Turbulence within this layer can be generated through mechanisms of both shear and thermal convection; the relative contributions of these two mechanisms is determined by the stability state of the atmosphere. The stability state is characterized by the Monin-Obukhov length:
\(u_\tau\) is the friction velocity, defined as the square root of the magnitude of the Reynolds shear stress at the surface, or
\(\theta_{ref}\) is a reference (virtual potential) temperature associated with the air within the surface layer; for example, the average temperature within the surface layer. \(\kappa \approx 0.41\) is the von Karman constant, and \(g\) is the acceleration of gravity. \(\overline{w^\prime \theta^\prime}_s\) is the surface turbulent temperature flux. Both the turbulent shear stress and turbulent temperature flux are approximately constant within the surface layer.
Applying a gradient diffusion model for the turbulent temperature flux leads to:
The sign of \(L\) is then connected to the sign of the temperature gradient within the surface layer. Three regimes are delineated:
\(\frac{1}{L} > 0, \quad \frac{\partial \theta}{\partial z} > 0\), stable stratification
\(\frac{1}{L} = 0, \quad \frac{\partial \theta}{\partial z} = 0\), neutral stratification
\(\frac{1}{L} < 0, \quad \frac{\partial \theta}{\partial z} < 0\), unstable stratification
Monin-Obukhov theory postulates the following similarity laws for mean velocity parallel to the surface and temperature,
where the forms of the non-dimensional functions \(\phi_m\) and \(\phi_h\) are determined from empirical observations. Analytical functions have been fit to the data; these are not given here, rather, we present the integrated form of ((9.15)) and ((9.16)), since these are the forms required by the code implementation.
For neutral stratification, \(\phi_m = 1\) and we recover the logarithmic profile for a “fully rough” surface,
where \(z_0\) is the characteristic roughness height. Note that viscous scaling involving surface viscosity and density properties is not required with this form of the logarithmic profile, since the roughness height is large enough to eliminate the presence of a laminar sublayer and buffer layer.
For stable stratification, the surface layer profiles take the form
\(\theta_*\) is calculated from the temperature flux and friction velocity as \(\theta_* = -\frac{\overline{w^\prime \theta^\prime}_s}{u_\tau}\), and \(\gamma_m\), \(\alpha_h\), and \(\gamma_h\) are constants specified below.
For unstable stratification, the surface layer profiles take the form
where
The constants used in ((9.18)) – ((9.23)) are [Dye74]
9.3.2. ABL Wall Function
The equations from the preceeding section can be used to formulate a wall function boundary condition for simulation of atmospheric boundary layers. The user-specified inputs to this boundary condition are: roughness length, \(z_0\), and surface heat flux, \(q_s = \rho C_p (\overline{w^\prime \theta^\prime})_s\). The surface layer profile model is evaluated for each surface boundary flux integration point; the wall-normal distance of the “first point off the wall” is taken to be one fourth of the length of the nearest edge intersecting the boundary face. The boundary condition is specified weakly through the imposition of a surface shear stress and surface heat flux.
The procedure for applying the boundary condition is as follows:
Determine the stratification state of the boundary layer by calculating the sign of the Monin-Obukhov length scale.
Solve the appropriate profile equation, either ((9.17)), ((9.18)), or ((9.20)), for the friction velocity \(u_\tau\). For the neutral case, \(u_\tau\) can be solved for directly. For the stable and unstable cases, \(u_\tau\) must be solved for iteratively because \(L\) appears in these equations and \(L\) depends on \(u_\tau\).
The surface shear stress is calculated as \(\tau_s = \rho_s u_\tau^2\). For calculating left-hand-side Jacobian entries, the form (9.24) is used, where \(\psi^\prime\) is zero for a neutral profile, \(-\gamma_m z/L\) for a stable profile, and \(\psi_h(z/L)\) for an unstable profile. The Jacobian entries follow directly from this form.
The user specified surface heat flux is applied to the enthalpy equation. Evaluation of surface temperature is not required for the boundary condition specification. However, if surface temperature is required for evaluation of other quantities within the code, the appropriate surface layer temperature profile should be used, either ((9.19)) or ((9.21)).
9.4. Atmospheric Boundary Layer Top Conditions
The abltop option is intended for the upper boundary of atmospheric boundary layer simulations. Currently it has two functions: (1) to provide a fully-automated open boundary condition that allows for inflows and outflows generated by terrain features or obstacles placed within the domain interior, and (2) to allow for mean temperature gradients at the upper boundary. These two features will be described separately below.
9.4.1. Inflow-Outflow Capability
The inflow-outflow capability is facilitated by a potential flow solution on a sub-domain consisting of a thin slab extending from the upper boundary a short distance into the domain. The vertical velocity at the slab bottom face is measured directly from the solution at each time step and this information is used to form a potential flow solution for the entire slab, which effectively provides all three velocity components at the main domain upper boundary. The potential flow solution is achieved at very little computational cost via Fast Fourier Transforms. The mathematical aspects of the boundary condition solution procedure are sketched below.
We start by assuming that the flow near the upper boundary is either laminar or weakly turbulent, as is the norm in atmospheric boundary layer simulations due to the presence of capping inversions (a region of increased stability that effectively terminates the boundary layer). If any turbulence within the slab is weak then the flow is approximately irrotational and any density fluctuations are weak. Furthermore, if the slab is thin (compared with the density scale height \(\sim 8\) km) then the mean density variation is also negligible. Thus we idealize the flow as being irrotational and of constant density, which leads to the following potential flow framework
where \(\phi\) is the disturbance potential and \(\langle\vec{u}\rangle\) is the mean velocity.
A Fourier transform of Laplace’s equation in \(x\) and \(y\) leads to
where \(k\) and \(l\) are the wavenumber components in the \(x\) and \(y\) directions respectively, and where \(\hat{\phi}\) is the Fourier transform of \(\phi\).
The bounded solution to the above equation is
where \(\hat{\phi}_0\) is the Fourier transform of \(\phi\) on the horizontal plane \(z=z_0\). The vertical velocity component (in Fourier \(x-y\) space) is
When evaluated at the position \(z=z_0\) the above relation yields the following specification for \(\hat{\phi}_0\)
Making use of this result, the solution for \(\hat{\phi}\) can be written as
If we identify \(z_0\) with the position of the potential flow slab lower face (the sampling plane), then the above prescription for \(\hat{\phi}\) yields the solution anywhere in the slab using just knowledge of \(\hat{w}\) at its lower boundary. The three velocity components are
Notice that all three velocity components contain a common wavenumber-dependent exponential damping term. This feature indicates that high wavenumber (small-scale) information on the sampling plane is strongly filtered in constructing the velocity field at the upper boundary. Thus the boundary condition velocity field will always be smoother that the velocity on the sampling plane. The filtering also increases with larger separations between the sampling plane and the upper boundary. This fact poses the need to compromise between the additional computational cost of adding a thicker potential flow region and the reduced filtering resulting from a thin slab. We have found that a potential flow region thickness of 10% of the overall computational domain height provides a good compromise between these two competing aspects of the solution.
The solution procedure is rather simple:
Measure \(w(x,y,z=z_0)\), at the position of the sampling plane.
Take the 2D Fourier transform of \(w\) to get \(\hat{w}_0(k,l)\).
Solve for the three velocity components given above using the computational domain upper boundary position for \(z\) and the sampling plane position for \(z_0\).
Inverse Fourier transform all three solution components.
While this solution procedure is simple and can be computed at very low cost using Fast Fourier Transforms (FFTs), it does have some obvious restrictions, namely
Data must be sampled on a uniform Cartesian grid at a fixed elevation.
Periodic or symmetry boundary conditions must exist in the lateral directions.
The simplest way to achieve the first constraint is to add a Cartesian mesh block to the upper portion of the main computational domain. The thickness of the addition can be rather small, say about 10% of the original domain height. The current implementation assumes that a plane at constant elevation with a uniform Cartesian mesh exists near the upper boundary. However, if required in the future, interpolation can be used to sample the solution on a horizontal plane. The BdyLayerVelocitySampler class can do this.
While the second constraint may seem restrictive, it allows for any combination of periodic or inflow/outflow boundary conditions in either of the two lateral directions. Half-wave instead of periodic transforms are used for an inflow/outflow direction. For example, for inflow/outflow in the \(x\) direction, we take \(u\sim\cos(\pi k x/xL)\), \(v\sim\sin(\pi k x/xL)\), \(w\sim\sin(\pi k x/xL)\). These prescriptions dictate that \(\partial u/\partial x=0\) at \(x=0\) (the inlet) and at \(x=xL\) (the outlet) and that \(v=0\), \(w=0\) at the inlet and outlet. In order to meet these constraints, the terrain should be horizontal for a short distance downstream of the inlet and a short distance upstream of the outlet. The appropriate combination of half-wave and full-wave (periodic) transforms in the two coordinate directions allow for various combinations of inflow/outflow or periodic conditions in these directions.
9.4.2. Temperature Gradient Capability
The temperature almost always varies with height within the atmosphere. The temperature gradient option allows the simulation to be consistent with this fact specifying the rate of change in temperature at the upper boundary.
9.4.3. Implementation Details
The ABL top boundary condition is activated by the keyword abltop_boundary_condition in the input file. The following is an example block of an ABL top boundary condition specification:
- abltop_boundary_condition: bc_uppertarget_name: topabltop_user_data:potential_flow_bc: truegrid_dimensions: [121, 2, 61]horizontal_bcs: [1, -1, 0, 0]z_sample: 0.85normal_temperature_gradient: 0.01
The potential flow and temperature gradient options can be used independently or together. The Boolean input potential_flow_bc: true activates the potential flow feature, whereas the presence of the keyword normal_temperature_gradient: value activates the temperature gradient feature. If only the temperature gradient feature is called for, then symmetry conditions (i.e. a slip wall) are used continuity and momentum boundary conditions at the upper boundary. If only the potential flow boundary condition is called for then the condition \(d\langle T\rangle/dz=0\) is used at the upper boundary (where \(\langle\rangle\) indicates a horizontal average).
If the potential flow option is selected, the user must then specify three additional user inputs: grid_dimensions, horizontal_bcs and z_sample. The grid_dimensions input specifies the number of mesh points in the three coordinate directions for the structured Cartesian portion of the domain containing the potential flow slab. This may be the entire domain, or just a sub-portion of it. The current implementation assumes that the grid points within the structured Cartesian region are tagged with a mesh index flag that indicates the relative position of each point within the structured mesh system. For example, if (imax,jmax,kmax) are the grid dimensions in the three coordinate directions, then the mesh index for the point (i,j,k) is k*imax*jmax + j*imax + i. The abl_mesh program will include the mesh index tag in the exodus grid file.
The horizontal_bcs input specifies the lateral boundary conditions in use in the i and j grid directions. Inflow is specified with +1, outlfow -1, and periodic with zero. Thus in the example above, inflow is used at the i=0 boundary, outflow at the i=imax boundary, and periodic conditions are used in the j direction.
The input z_sample specifies the elevation of the sampling plane. If this input is missing, the default position of 90% of the distance between the lower and upper boundary will be used.
9.5. Turbulent Kinetic Energy, \(k_{sgs}\) LES model
When the boundary layer is assumed to be resolved, the natural boundary condition is a Dirichlet value of zero, \(k_{sgs} = 0\).
When the wall model is used, a standard wall function approach is used with the assumption of equal production and dissipation.
The turbulent kinetic energy production term is consistent with the law of the wall formulation and can be expressed as,
The parallel velocity, \(u_{\|}\), can be related to the wall shear stress by,
Taking the derivative of both sides of Equation (9.26), and substituting this relationship into Equation (9.25) yields,
Applying the derivative of the law of the wall formulation, Equation (9.2), provides the functional form of \({\partial u^+ / \partial y^+}\),
Substituting Equation (9.2) within Equation (9.27) yields a commonly used form of the near wall production term,
Assuming local equilibrium, \(P_k = \rho\epsilon\), and using Equation (9.29) and Equation (9.3) provides the form of wall shear stress is given by,
Under the above assumptions, the near wall value for turbulent kinetic energy, in the absence of convection, diffusion, or accumulation is given by,
This expression for turbulent kinetic energy is evaluated at the boundary faces of the exposed wall boundaries and is area-assembled to the nodal value for use in a Dirichlet condition.
9.5.1. Turbulent Kinetic Energy and Specific Dissipation SST Low Reynolds Number Boundary conditions
For the turbulent kinetic energy equation, the wall boundary conditions follow that described for the \(k_{sgs}\) model, i.e., \(k=0\).
A Dirichlet condition is also used on \(\omega\). For this boundary condition, the \(\omega\) equation depends only on the near-wall grid spacing. The boundary condition is given by,
which is valid for \(y^{+} < 3\).
9.5.2. Turbulent Kinetic Energy and Specific Dissipation SST High Reynolds Number Boundary conditions
The high Reynolds approach uses the law of the wall assumption and also follows the description provided in the wall modeling section with only a slight modification in constant syntax,
In the case of \(\omega\), an analytic expression is known in the log layer:
which is independent of \(k\). Because all these expressions require \(y\) to be in the log layer, they should absolutely not be used unless it can be guaranteed that \(y^{+} > 10\), and \(y^{+} > 25\) is preferable. Automatic blending is not currently supported.
9.5.3. Solid Stress
The boundary conditions applied are either force provided by a static pressure,
or a Dirichlet condition, i.e., \(u_i = u^{spec}_i\), on the displacement field. Above, \(F^n_i\) is the force for component \(i\) due to a prescribed [static] pressure.
9.5.4. Intensity
The boundary condition for each intensity assumes a grey, diffuse surface as,
9.6. SST of the Atmospheric Boundary Layer
The following boundary conditions simulate the Atmospheric Boundary Layer, as described in Bautista, [Bau11] and [BDM15]. The Nalu-Wind SST RANS implementation matches the Monin-Obukhov profile when used with the model constants from Table-A I-1 (Boundreault, 2011) in [Bau11] and the meshing method described in [BDM15]. The mesh described in [BDM15] gives the Monin-Obukhov profile for roughness height 0.1. When the roughness height is decreased, the mesh must be refined near the wall. For example, for the [BDM15] ABL test case using roughness height 0.001 instead of 0.1, the mesh size needs to be halved near the wall.
The \(k\) and \(\omega\) boundary conditions are the same as in the Turbulent Kinetic Energy and Specific Dissipation SST High Reynolds Number Boundary Conditions:
and
The friction velocity, \(u_{\tau}\) is calculated from a reference velocity, \(u_h\), at height, \(h\). \(h\) and \(u_h\) could, for example, be the hub height of a wind turbine and the velocity measured at that height.
where \(\kappa\) is the Von K{'a}rm{'a}n constant and \(z_0\) is the roughness height.
The momentum boundary condition is a no-slip Dirichlet condition, \(u_i=0\), as described in the momentum wall boundary conditions.
The streamwise and spanwise boundary conditions are periodic, as described in periodic boundary conditions.
The \(k\), \(\omega\), and \(u\) wall boundary conditions are set in the input file by specifying a wall boundary condition with RANS_abl_bc. The input file must also specify a height and the velocity at that height with reference_height and reference_velocity.
realms:
- name: fluidRealm
boundary_conditions:
- wall_boundary_condition: bc_lower
wall_user_data:
RANS_abl_bc: yes
reference_velocity: 6.6
reference_height: 90.0
Rather than specifying a momentum source term, \(dp/dx\), the velocity is set to the geostrophic (freestream) value near the top of the domain using the ABL forcing method described in ABL Forcing Source Terms. The ABL forcing term is turned on and the desired wind velocity and height must be specified in the input file:
realms:
- name: fluidRealm
solution_options:
options:
- source_terms:
momentum:
- abl_forcing
abl_forcing:
output_format: "abl_%s_sources.dat"
momentum:
type: computed
relaxation_factor: 1.0
heights: [2500.0]
velocity_x:
- [0.0, 17.5]
- [1000000.0, 17.5]
velocity_y:
- [0.0, 0.0]
- [1000000.0, 0.0]
velocity_z:
- [0.0, 0.0]
- [1000000.0, 0.0]
9.7. AMS of the Atmospheric Boundary Layer
The following boundary conditions are meant for Active Model Split (AMS) of the Atmospheric Boundary Layer (ABL) when using SST as the mean (RANS) contribution to AMS. The boundary conditions on \(u_i\), \(k\), and \(\omega\) are the same as for SST of the ABL, as described in SST of the Atmospheric Boundary Layer. When using this boundary condition for AMS, the energy transfer (LES) components of AMS should be turned off below the sand grain roughness height, \(k_s = 30 z_0\), where \(z_0\) is the roughness height. Below this height the solution from the boundary condition is not necessarily accurate. Therefore including energy transfer components of AMS, such as the active forcing, in this region could introduce errors.
Specifically, the energy transfer components of AMS below \(k_s\) are turned off by setting the scaling coefficient on the energy transfer contribution to AMS to one, resolution adequacy indicator to one, and active forcing to zero. To turn off the energy transfer components the input file should include:
realms:
- name: fluidRealm
solution_options:
options:
- user_constants:
roughness_height: 0.3
rans_below_ks: true
9.8. Open Boundary Condition
Open boundary conditions require far more care. In general, open bcs are assembled by iterating faces and the boundary integration points on the exposed face. The parent element is also required since oftentimes gradients are used (for momentum). For an open boundary condition the flow can either leave or enter the domain depending on what the computed mass flow rate at the exposed boundary integration point is. Two options are available computing the velocity of the entrained flow—either the normal velocity at the integration point is used or a specified normal velocity is used. The tangential components are always specified.
9.8.1. Continuity
For continuity, the boundary mass flow rate must also be computed. This value is stored and used for the other equations that require advection. The same formula is used for the pressure-stabilized mass flow rate. However, the local pressure gradient for each boundary contribution is based on the difference between the interior integration point and the user-specified pressure which takes on the boundary value. This can optionally be modified to be a “total pressure”—removing the kinetic energy associated with entrainment at the open. The interior integration point is determined by linear interpolation. For CVFEM, full elemental averaging is used while in EBVC discretization, the midpoint value between the nearest node and opposing node to the boundary integration point is used. In both discretization approaches, non-orthogonal corrections are required. This procedure has been very important for stability for CVFEM tet-based meshes where a natural non-orthogonality exists between the boundary and interior integration point.
For wind energy applications, the usage of the standard open boundary mass flow rate expression, which includes pressure contributions, is not appropriate due to complex temperature/buoyancy specifications. In these cases, a global correction algorithm is supported. Specifically, pressure terms are dropped at the open boundary mass flow rate expression in favor or a pre-processing algorithm that uniformly distributes the continuity mass flow rate (and possible density accumulation) “error” over the entire set of open boundary conditions. The global correction scheme may perform well with single open boundary condition specification, e.g., multiple inflows with a single open location, however, it is to be avoided if the flow leaving the domain is complex in that a simulation includes multiple open boundary conditions. A complex situation might be an open jet with entrainment from the side (open boundary that allows for inflow) and a top open that allows for outflow. However, a routine case might be a backward facing step with a single inflow, side periodic, top wall and open boundary. Not that the ability for the continuity solve to be well conditioned may require an interior Dirichlet on pressure as the open pressure specification for the global correction algorithm is lacking. In most cases, a Dirichlet condition is not actually required as the nullspace of the continuity system may not be found in the solve.
9.8.2. Momentum
For momentum, the normal component of the stress is subtracted out we subtract out the normal component of the stress. The normal stress component for component i can be written as \(F_k n_k n_i\). The tangential component for component i is simply, \(F_i - F_k n_k n_i\). As an example, the tangential viscous stress for component x is,
which can be written in general component form as,
Finally, the normal stress contribution is applied based on the user specified pressure,
For CVFEM, the face gradient operators are used for the thermal stress terms. For EBVC discretization, from the boundary integration point, the nearest node (the “Right” state) is used as well as the opposing node (the “Left” state). The nearest node and opposing node are used to compute gradients required for any derivatives. This equation follows the standard gradient description in the diffusion section with non-orthogonal corrections used. In this formulation, the area vector is taken to be the exposed area vector. Non-orthogonal terms are noted when the area vector and edge vector are not aligned.
For advection, if the flow is leaving the domain, we simply advect the nearest nodal value to the boundary integration point. If the flow is coming into the domain, we simply confine the flow to be normal to the open boundary integration point area vector. The value entrained can be the nearest node or an upstream velocity value defined by the edge midpoint value or by a specified value.
9.8.3. Mixture Fraction, Enthalpy, Species, \(k_{sgs}\), k and \(\omega\)
Open boundary conditions assume a zero normal gradient. When flow is entering the domain, the far-field user supplied value is used. Far field values are used for property evaluations. When flow is leaving the domain, the flow is advected out consistent with the choice of interior advection operator.
9.9. Strong Symmetry Boundary Condition
There are two implementations of the symmetry boundary condition: strong and weak. In the strong symmetry implementation, the normal velocity is set to zero at the boundary. Strong symmetry has only been implemented for a cartesian mesh, meaning it can be used for flat surfaces that are aligned with the principle cartesian directions. It cannot be used for curved surfaces or flat surfaces that are not aligned with the principle cartesian directions. Both the strong and weak symmetry boundary conditions have an associated error. In the strong form the associated error lies not on the boundary but in the domain.
9.10. Weak Symmetry Boundary Condition
9.10.1. Continuity, Mixture Fraction, Enthalpy, Species, \(k_{sgs}\), k and \(\omega\)
Weak symmetry applies zero diffusion at the boundary for scalar quantities, which effectively sets the boundary-normal gradients of these quantities to zero. This means that, unlike for strong symmetry, in the weak symmetry implementation, normal velocity can be non-zero. This is possible because the cell averaged quantities in Nalu-Wind’s discretization are stored at the nodes and therefore reside on the boundaries. In general, a non-zero normal velocity can cause net inflow or outflow. The Poisson solve prevents this by enforcing mass conservation globally to the order of linear solver convergence.
9.10.2. Momentum
A symmetry boundary is one that is described by removal of the tangential stress. For weak symmetry this is done in the momentum equation by applying only the normal component of stress:
which can be written in general component form as,
The momentum equation also penalizes non-zero normal velocity. The strength of this penality depends on the penalty factor. One can enforce a stronger representation of the boundary by amplifying the penalty factor through the symmetry_bc_penalty_factor variable in solution_options. Its default value is 2.0 which is the minimum required for stability. Amplifying the penalty factor will enforce the boundary in a stronger sense, but can also lead to a more difficult matrix solve as with any large penalty term.
Both strong and weak symmetry boundary conditions have associated error. In the weak form the error manifests as non-zero local velocities at the boundary. As the mesh is refined, the weak symmetry boundary condition converges to the strong symmetry boundary condition (zero normal velocity) with a first order rate. Note that switching from weak to strong symmetry boundary condition does not make the error go away; the error just moves off the boundary and into the domain. The errors for both weak and strong symmetry boundary conditions do go away as the mesh is refined, decreasing with first order rate.
Weak symmetry has not been implemented for Active Model Split (AMS).
If the symmetry type is not specified in the input file then the code defaults to weak symmetry. If the ABL top boundary condition is used and it defaults to symmetry, as described in ref:Atmospheric Boundary Layer Top Conditions <theory_abltop_bc>:, then strong or weak symmetry can be specified explicitly in the input file:
realms:
- name: fluidRealm
boundary_conditions:
- abltop_boundary_condition: bc_upper
target_name: upper
symmetry_user_data:
symmetry_type: "z_direction_strong"
9.10.3. Specified Boundary-Normal Temperature Gradient Option
The standard symmetry boundary condition applies zero diffusion at the boundary for scalar quantities, which effectively results in those scalars having a zero boundary-normal gradient. There are situations, especially for atmospheric flows in which the user may desire a finite boundary-normal gradient of temperature. For example, the atmospheric boundary layer is often simulated with a stably stratified capping inversion in which the temperature linearly increases with height all the way to the upper domain boundary. We apply symmetry conditions to this upper boundary for momentum, but we specify the boundary-normal temperature gradient on this boundary to match the capping inversion’s gradient.
This is an option in the symmetry boundary condition specification, which appears in the input file as:
- symmetry_boundary_condition: bc_upper
target_name: upper
symmetry_user_data:
normal_temperature_gradient: -0.003
In this example, the temperature gradient normal to the symmetry boundary is set to -0.003 K/m, where the boundary-normal direction is pointed into the domain.
Nalu-Wind does not solve a transport equation for temperature directly, but rather it solves one for enthalpy. Therfore, the boundary-normal temperature gradient condition is applied internally in the code through application of a compatible heat flux,
where \(q_n\) is the heat flux at the boundary, \(\kappa_{eff}\) is the effective thermal diffusivity (the molecular and turbulent parts), \(c_p\) is the specific heat, and \(\partial T / \partial n\) is the boundary-normal temperature gradient.
9.11. Periodic Boundary Condition
A parallel multiple-periodic boundary condition is supported. Mappings are created between master/slave surface node pairs. The node pairs are obtained from a parallel search and are expected to be unique. The node pairs are used to map the slave global id to that of the master. This allows the linear system to include matrix rows for only a subset of the overall set of nodes. Moreover, a periodic assembly for assembled quantities is managed via: \(m+=s\) and \(s=m\), where \(m\) and \(s\) are master/slave nodes, respectively. For each parallel assembled quantity, e.g., dual volume, turbulence quantities, etc., this procedure is used. Periodic boxes and periodic couette and channel flow have been simulated in this code base. Tow forms of parallel searches exist and are supported (one through the Boost TPL and another through the STK Search module).
9.12. Non-conformal Boundary Condition
A surface-based approach based on a DG method has been discussed in the 2010 CTR summer proceedings by Domino, [Dom10]. Both the edge- and element-based formulation currently exists in the code base using the CVFEM and EBVC approaches.
9.1 Two-block example with one common surface, \(\Gamma_{AB}\).
Consider two domains, \(A\) and \(B\), which have a common interface, \(\Gamma_{AB}\), and a set of interfaces not in common, \(\Gamma \backslash \Gamma_{AB}\) (see Figure 9.1), and assume that the solution of the time-dependent advection/diffusion equation is to be solved in both domains. Each domain has a set of outwardly pointing normals. In this cartoon, the interface is well resolved, although in practice this may not be the case.
An interior penalty approach is constructed at each integration point at the exposed surface set. The numerical flux for a general scalar \(\phi\) is constructed at the current integration point which is based on the current (\(A\)) and opposing (\(B\)) elemental contributions,
where \(q_j^A\) and \(q_j^B\) are the diffusive fluxes computed using the current and opposing elements and normals are outward facing. The penalty coefficient \(\lambda^A\) contains the diffusive contributions averaged over the two elements,
Above, \(\Gamma^k\) is the diffusive flux coefficient evaluated at current and opposing element location, respectively, and \(L^k\) is an elemental length scale normal to the surface (again for current and opposing locations, \(A\) and :math`B`). When upwinding is activated, the value of \(\eta\) is unity.
As written in Equation (9.35), the default convection and diffusion term is a Galerkin approach, i.e., equally averaged between the current and opposing face. The standard advection term is given by,
For surface A, the form is as follows:
with the nonconformal mass flow rate given by,
In the above set of expressions, the consistent definition of \(\hat{u}_j\), i.e., the convecting velocity including possible pressure stabilization terms, is retained.
As with the interior advection scheme, the mass flow rate involves pressure stabilization terms. The value of \(\gamma\) defines whether or not the full pressure stabilization terms are included in the mass flow rate expression. Equation (9.39) also forms the continuity nonconformal boundary contribution.
With the substitution of \(\eta\) to be unity, the effective convective term is as follows:
Note that this form reduces to a standard upwind operator.
Since this algorithm is a dual pass approach, a numerical flux can be written for the integration point on block \(B\),
As with Equation (9.41), \(\dot{m}^B\) (see Equation (9.42)) is of similar form to \(\dot{m}^A\),
For low-order meshes with curved surface, faceting will occur. In this case, the outward facing normals may not be (sign)-unity factors of each other. In this case, it may be adventageous to define the opposing outward normal as, \(n_j^B = -n_j^A\).
Domino, [Dom10] provided an overview of a FEM fluids implementation. In such a formulation, the interior penalty term appears, i.e.,
and
Although the sign of this term is often debated in the literature, the above set of expressions acts to increase penalty term stencil to include the full element contribution. As the CVFEM uses a piecewise-constant test function, this term is currently neglected.
Average fluxes are computed based on the current and opposing integration point locations. The appropriate DG terms are assembled as boundary conditions first with block \(A\) integration points as \(current\) (integrations points for block B are \(opposing\)) and then with block \(B\) integration points as \(current\) (surfaces for block A are, therefore, \(opposing\)). Figure 9.1 graphically demonstrates the procedure in which integration point values of the flux and penalty term are computed on the block \(A\) surface and at the projected location of block \(B\).
A parallel search is conducted to project the current integration point location to the opposing element exposed face. The search, therefore, provides the isoparametric coordinates on the opposing element. Elemental shape functions and shape function derivatives are used to construct the numerical flux for both the edge- and element-based scheme. The location of the Gauss points on the current element are either the Gauss Labatto or Gauss Legendre locations (input file specification). For each equation (momentum, continuity, enthalpy, etc.) the numerical flux is computed at each exposed non-conformal surface.
As noted, for most equations other than continuity and heat condition, the numerical flux includes advection and diffusion contributions. The diffusive contribution is easily provided using elemental shape function derivatives at the current and opposing surface.
9.2 Description of the numerical flux calculation for the DG algorithm. The value of fluxes and penalty values on the current block (\(A\)) and the opposing block (\(B\)) are used for the calculation of numerical fluxes. \(\tilde \varphi\) represents the projected value.
Above, special care is taken for the value of the mass flow rate at the non-conformal interface. Also, note that the above written form does not upwind the advective flux, although the code allows for an upwinded approach. In general, the advective term contains contributions from both elements identified at the interface, specifically.
The penalty coefficient for the mass flow rate at the non-conformal boundary points is again a function of the blended inverse length scale at the current and opposing element surface location. The form of the mass flow rate above provides the continuity contribution and the form of the mass flow rate used in the scalar non-conformal flux contribution.
The full connectivity for element integration and opposing elements is within the linear system. As such, for sliding mesh configurations, the linear system connectivity graph changes each time step. Recent prototyping of the dG-based and the overset scheme has allowed this method to be used across both disparate low-order topologies (see Figure 9.3).
9.3 Discontinuous Galerkin non-conformal interface mixed topology (hex8/tet4).
Footnotes