Convection/diffusion

Thermal equilibrium equation

The thermal equilibrium equation for a continuum in which a fluid is flowing with velocity $v$ , is

\int δ θ [ρ c {\frac{\partial θ}{\partial t} + v \cdot \frac{\partial θ}{\partial x}} - \frac{\partial}{\partial x} \cdot (k \cdot \frac{\partial θ}{\partial x}) - q] d V + \int_{S_{q}} δ θ [n \cdot k \cdot \frac{\partial θ}{\partial x} - q_{s}] d S = 0,

where $θ (x, t)$ is the temperature at a point, $δ θ (x, t)$ is an arbitrary variational field, $ρ (θ)$ is the fluid density, $c (θ)$ is the specific heat of the fluid, $k (θ)$ is the conductivity of the fluid, q is the heat added per unit volume from external sources, $q_{s}$ is the heat flowing into the volume across the surface on which temperature is not prescribed ( $S_{q}$ ), $n$ is the outward normal to the surface, $x$ is spatial position, and t is time. Although most fluids will have isotropic conductivity, so that $k = k I$ (where $k (θ)$ is a scalar and $I$ is a unit matrix), we provide for anisotropic conductivity to cover such cases as that of fluid flowing through a set of baffle plates whose conductivity is smeared into that of the fluid.

The boundary conditions are that $θ (x)$ is prescribed over some part of the surface, $S_{θ}$ , and that the heat flux per unit area entering the domain across the rest of the surface, $q_{s} (x)$ , is prescribed or is defined by convection and/or radiation conditions. For example, the boundary layer between fluid convection elements and solid elements might be modeled by DINTERx-type elements. The boundary term in the thermal equilibrium equation defines

q_{s} = - n \cdot k \cdot \frac{\partial θ}{\partial x} .

This implies that $q_{s}$ is the flux associated with conduction across the surface only—any convection of energy across the surface is not included in $q_{s}$ . This makes no difference if the surface is part of a solid body (where $q_{s}$ would be defined by heat transfer into the adjacent body), since then the normal velocity into that body, $v \cdot n$ , is zero. But it does make a difference when there is fluid crossing the surface, as—for example—on the upstream and downstream boundaries of the mesh. In this case the choice of $q_{s}$ for the natural boundary condition (instead of using the total flux crossing the surface) is desirable because it avoids spurious reflections of energy back into the mesh as the fluid flows through the surface.

These equations are discretized with respect to position by using first-order isoparametric elements. The fluid velocity, $v$ , is assumed to be known. (Abaqus actually requires that the mass flow rate of the fluid per unit area be defined, because this is generally more convenient for the user. The velocity is computed from the mass flow rate and the density of the fluid.)

The time discretization generates the solution at time $t + Δ t$ from the known solution at time t.

The interpolation for the temperature, $θ$ , is defined over an element and over a time increment as

θ (x, t) = N^{N} (x) A^{n} (t) θ^{(N, n)}, for N = 1, 2, \dots, n = t, t + Δ t,

where the $N^{N}$ are standard isoparametric functions and the time interpolation, $A^{n}$ , is linear:

A^{t} = 1 - \frac{τ}{Δ t}, A^{t + Δ t} = \frac{τ}{Δ t},

where $Δ t$ is the time increment and $0 \leq τ \leq Δ t$ .

The Petrov-Galerkin discretization proposed by Yu and Heinrich couples this linear interpolation with the weighting functions

\begin{matrix} δ θ & ​ = W^{N} δ θ^{N} \\ ​ = [N^{N} \bar{A} + \frac{h}{2} (α \bar{A} + β \frac{Δ t}{2} \frac{d \bar{A}}{d t}) \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x}] δ θ^{N}, \end{matrix}

where

\bar{A} = 6 \frac{τ}{Δ t^{2}} (1 - \frac{τ}{Δ t});

$v$ is the average fluid velocity over the element; $| v |$ is its magnitude; and h is a characteristic element length measure, defined below. $α$ and $β$ are control parameters. The $α$ term in the weighting is introduced to eliminate artificial diffusion of the solution, while the $β$ term is introduced to avoid numerical dispersion. Yu and Heinrich show that the optimal choices are

α = \coth \frac{γ}{2} - \frac{2}{γ} and β = \frac{C}{3} - \frac{2}{γ} \frac{α}{C},

where $γ$ is the local Péclet number in an element and C is the local Courant number, defined as

γ = | v | h \frac{ρ c}{k} and C = | v | \frac{Δ t}{h} .

The above expression for $β$ yields negative values for very small fluid velocities, which may destabilize the solution; hence, for low velocities dispersion control is switched off.

The characteristic element length measure, h, is defined by Yu and Heinrich as follows.

Let $h_{α}$ be the $α$ isoparametric line across the element passing through its centroid. The projection of $h_{α}$ in the direction of the fluid velocity vector at the element's centroid is

h_{α} = h_{α} \cdot \frac{v}{| v |} .

Then we define h as

h = \sum_{α} | h_{α} | .

When $β$ is nonzero, these elements require that $C \leq 1$ for numerical stability.

Since the weighting functions are biased (“upwinding”), they are discontinuous from one element to the next. Some care is, therefore, required in manipulating the weak form of the thermal equilibrium equation (see Hughes and Brooks, 1982). In particular, the usual integration by parts of the conduction term

δ θ \frac{\partial}{\partial x} \cdot (k \cdot \frac{\partial θ}{\partial x})

can be performed only for the continuous part of the weighting functions used to discretize $δ θ$ : otherwise, continuity of heat flux between elements is not assured. For convenience we write the discontinuous part of the weighting as

P^{N} = \frac{h}{2} (α \bar{A} + β \frac{Δ t}{2} \frac{d \bar{A}}{d t}) \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} .

The weak form of thermal equilibrium is

\begin{matrix} \int [W^{N} ρ c {N^{M} \frac{d A^{n}}{d t} + v \cdot \frac{\partial N^{M}}{\partial x} A^{n}} + \bar{A} \frac{\partial N^{N}}{\partial x} \cdot k \cdot \frac{\partial N^{M}}{\partial x} A^{n} - P^{N} k : \frac{\partial^{2} N^{M}}{\partial x \partial x} A^{n}] d V \\ - \int W^{N} q d V - \bar{A} \int_{S} N^{N} q_{s} d S & ​ = 0. \end{matrix}

This can be rewritten as

\begin{matrix} \bar{A} \frac{d A^{n}}{d t} \int ρ c (N^{N} + α \frac{h}{2} \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x}) N^{M} d V \\ + \bar{A} A^{n} \int ρ c (N^{N} + α \frac{h}{2} \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x}) v \cdot \frac{\partial N^{M}}{\partial x} d V \\ + β \frac{h Δ t}{4} \frac{d \bar{A}}{d t} (\frac{d A^{n}}{d t} \int ρ c \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} N^{M} d V + A^{n} \int ρ c \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} v \cdot \frac{\partial N^{M}}{\partial x} d V) \end{matrix}

\begin{matrix} + \bar{A} A^{n} \int (\frac{\partial N^{N}}{\partial x} \cdot k \cdot \frac{\partial N^{M}}{\partial x} - α \frac{h}{2} \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} k : \frac{\partial^{2} N^{M}}{\partial x \partial x}) d V \\ - β \frac{h Δ t}{4} \frac{d \bar{A}}{d t} A^{n} \int \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} k : \frac{\partial^{2} N^{M}}{\partial x \partial x} d V \\ - \bar{A} \int_{V} [N^{N} + α \frac{h}{2} \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x}] q d V - \frac{d \bar{A}}{d t} \int_{V} β \frac{Δ t}{2} \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} q d V - \bar{A} \int_{S} N^{N} q_{s} d S = 0. \end{matrix}

We now integrate this equation from time t to $t + Δ t$ to provide an average equilibrium statement for the increment. We use the results

\int_{Δ t} \bar{A} d t = 1, \int_{Δ t} \frac{d \bar{A}}{d t} d t = 0,

\int_{Δ t} \bar{A} \frac{d A^{t + Δ t}}{d t} d t = - \int_{Δ t} \bar{A} \frac{d A^{t}}{d t} d t = \frac{1}{Δ t}, \int_{Δ t} \bar{A} A^{t + Δ t} d t = \int_{Δ t} \bar{A} A^{t} d t = \frac{1}{2},

and

\int_{Δ t} \frac{d \bar{A}}{d t} \frac{d A^{t + Δ t}}{d t} d t = \int_{Δ t} \frac{d \bar{A}}{d t} \frac{d A^{t}}{d t} d t = 0, \int_{Δ t} \frac{d \bar{A}}{d t} A^{t + Δ t} d t = - \int_{Δ t} \frac{d \bar{A}}{d t} A^{t} d t = - \frac{1}{Δ t},

to give

\begin{matrix} \frac{1}{Δ t} \int ρ c (N^{N} + α \frac{h}{2} \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x}) N^{M} d V (θ^{M, t + Δ t} - θ^{M, t}) \\ + \frac{1}{2} \int ρ c (N^{N} + α \frac{h}{2} \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x}) v \cdot \frac{\partial N^{M}}{\partial x} d V (θ^{M, t + Δ t} + θ^{M, t}) \\ - β \frac{h}{4} \int ρ c \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} v \cdot \frac{\partial N^{M}}{\partial x} d V (θ^{M, t + Δ t} - θ^{M, t}) \\ + (\frac{1}{2} \int \frac{\partial N^{N}}{\partial x} \cdot k \cdot \frac{\partial N^{M}}{\partial x} d V - α \frac{h}{4} \int \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} k : \frac{\partial^{2} N^{M}}{\partial x \partial x} d V) (θ^{M, t + Δ t} + θ^{M, t}) \\ + β \frac{h}{4} \int \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x} k : \frac{\partial^{2} N^{M}}{\partial x \partial x} d V (θ^{M, t + Δ t} - θ^{M, t}) \\ - \int (N^{N} + α \frac{h}{2} \frac{v}{| v |} \cdot \frac{\partial N^{N}}{\partial x}) q d V - \int_{S} N^{N} q_{s} d S = 0. \end{matrix}

For the steady-state case the third term in this equation is omitted. In both transient and steady-state forms the contribution of such a convective element to the system of equations for the heat transfer model is not symmetric, requiring the use of the nonsymmetric matrix storage and solution scheme.