True interpolation

Although the primary variables are determined at the element centers the governing equations require their values at the element faces, as well as their gradient at the centers and faces. Let us have a look at a typical element P and one of its neighbors F (Figure [*]; P and F will be used to denote the elements as well as their centers).

Figure: Mesh parameters
\begin{figure}\epsfig{file=figF1.eps,width=12cm}\end{figure}

Based on the values of a variable \phi at the element centers P and F the value in the center f of face f is looked for. To that end the following vectors are defined: \boldsymbol{i}_{\xi} is a unit vector connecting P with f, while \boldsymbol{n} is the normal vector on face f at its center f (pointing way from the element). The angle between both is \alpha. Furthermore, \boldsymbol{j}_{\xi} is a vector connecting P with F, \beta is the angle between \boldsymbol{n} and \boldsymbol{j}_{\xi}, Q is the orthogonal projection of f on the straight line connecting P with F and \boldsymbol{r}_f connects Q to f. The value \phi_f can now be approximated by

\phi_f \approx \phi _Q \approx \frac{[\phi_P d(Q,F) + \phi_F d(Q,P)]}{d(P,F)}, (553)

where d(Q,F) is the Euclidean distance between Q and F. The next better approximation takes \boldsymbol{r}_f into account:

\phi_f \approx \phi_Q + (\nabla \phi)_Q \cdot \boldsymbol{r}_f. (554)

To this end the gradient is approximated by

(\nabla \phi)_Q \approx \frac{(\nabla \phi)_P d(Q,F)+(\nabla \phi)_F d(Q,P)}{d(P,F)}. (555)

The value of the gradient at the element centers can be obtained from the facial values. Indeed, by Gauss the following applies:

\int_{V}^{} \nabla \cdot \boldsymbol{a} d V = \int_{A}^{} \boldsymbol{a} \cdot \boldsymbol{n} d A (556)

for an arbitrary vector \boldsymbol{a}. Setting \boldsymbol{a}=\phi \cdot \boldsymbol{i}_k, where \boldsymbol{i}_k is a unit vector in carthesian direction k, one obtains

\int_{V}^{} \nabla \cdot (\phi \boldsymbol{i}_k) d V = \int_{A}^{} \phi \boldsymbol{i}_k \cdot \boldsymbol{n} dA \approx \sum_{f}^{} \phi_f \boldsymbol{i}_k \cdot \boldsymbol{n}_f A_f = \sum_{f}^{} \phi_f (n_k)_f A_f. (557)

and

\int_{V}^{} \nabla \cdot (\phi \boldsymbol{i}_k) d V = \int_{V}^{} (\nabla \phi) \cdot \boldsymbol{i}_k d V = \int_{V}^{} \phi,_k d V \approx (\phi,_k)_P V_P. (558)

From which one obtains

(\phi,_k)_P \approx \frac{1}{V_P} \left [ \sum_{f}^{} \phi_f (n_k)_f A_f \right], (559)

where the sum is over all faces f belonging to element P, \phi_f is the value of \phi at the center of face f, A_f is the area of face f and V_P is the volume of element P.

The calculation of \phi_f is usually done iteratively using the following steps:

  1. use Equation ([*]) to determine \phi_f except if \phi_f is known by a boundary condition (BC). In the latter case take the value of the BC.
  2. use Equation ([*]) to determine (\nabla \phi)_P and (\nabla \phi)_F.
  3. by Equation ([*]) determine (\nabla \phi)_Q.
  4. Use Equation ([*]) to improve \phi_f except if \phi_f is know from a BC (in that case take the value of the BC).

Steps 2 up to 4 can be repeated to improve the value of \phi_f further. Usually one such extra iteration is sufficient. Notice that (\nabla \phi)_Q can be used as an approximation of (\nabla \phi)_f.

Figure: Mesh parameters for a boundary element
\begin{figure}\begin{center}
\epsfig{file=figF2.eps,width=10cm}\end{center}\end{figure}

For a boundary element (Figure [*]) the situation is slightly different. The value at Q is now extrapolated from the center value at the neighboring element W (if any) bordering the face opposite to face f:

\displaystyle \phi_Q \displaystyle = \phi_P + \frac{\phi_P - \phi_W}{d(P,W)} d(P,Q) (560)
  \displaystyle = \phi_P \left ( 1 + \frac{d(P,Q)}{d(P,W)} \right) - \phi_W \frac{d(P,Q)}{d(P,W)}. (561)

Apart from this, the situation is completely analogous. If no element on the opposite side exists the center value is taken, i.e. \phi_f=\phi_P. This also amounts to (\nabla \phi)_f \cdot \boldsymbol{i}_\xi =0. Notice also that for a boundary element \boldsymbol{i}_\xi = \boldsymbol{j}_\xi.

The preceding discourse is generic and, although explained here for 2-dimensional elements applies in 3 dimensions equally well. In fact, CalculiX only knows three-dimensional hexahedral, pentahedral and tetrahedral fluid elements. In the following the values of \phi_f and (\nabla \phi)_f obtained in the above way will be denoted by \overline{\phi }_f and (\overline{\nabla \phi })_f, respectively.

Figure: Contributing elements if performing steps 1-2-3
\begin{figure}\begin{center}
\epsfig{file=figF3.eps,width=10cm}\end{center}\end{figure}

(\overline{\nabla \phi })_f, however, can be further improved. The value of (\overline{\nabla \phi })_e (Figure [*]), as calculated in step 3 above, depends on (\nabla \phi)_P and (\nabla \phi)_E. These depend on the values \overline{\phi }_f of all faces belonging to elements P and E, respectively (step 2). Each of these depends on the neighboring element values, e.g. \overline{\phi }_n depends on \phi_N and \phi_P (step 1). Overall, all dashed elements in addition of element P and E in the figure contribute to (\overline{\nabla \phi })_e. If, instead of just performing steps 1 through 3, steps 1-2-3-4-2-3 are taken, the values in the vertically dashed elements in Figure [*] will also contribute.

Figure: Contributing elements if performing steps 1-2-3-4-2-3
\begin{figure}\begin{center}
\epsfig{file=figF4.eps,width=10cm}\end{center}\end{figure}

However, for consistency it has proven vital that the simple neighbor relationship

(\nabla \phi)_e \cdot \boldsymbol{j}_{PE} = \frac{\phi_E-\phi_P}{l_{PE}} (562)

be satisfied ( l_{PE}=d(P,E)). This can be obtained by setting

(\nabla \phi)_e = (\overline{\nabla \phi })_e + \left ( \frac{\phi_E-\phi_P}{l_{PE}}- (\overline{\nabla \phi })_e \cdot \boldsymbol{j}_{PE} \right ) \boldsymbol{j}_{PE}. (563)

(\nabla \phi)_e is an improved value of the facial gradient. It will be denoted without bar on top.