Morphing

The boundaries of a region can change position and shape over time—for example, due to the motion of a contacting body. In Simcenter STAR-CCM+, the morphing algorithm redistributes the mesh vertices in response to the displacement of the boundaries, allowing for non-rigid deformations of the mesh.

The morpher algorithm redistributes the mesh vertices—it does not act on cell volumes, faces, or edges. The morpher does not optimize the grid, since it cannot compute mesh-quality metrics.

The initial movement that the morpher imposes on the mesh is defined through a set of control points. These control points originate from the mesh vertices on the morphing boundary, or from a control point table. Each control point is associated with a known displacement vector. The morpher uses these displacements to construct an interpolation field, which is used to calculate the displacements of all mesh vertices.

The morpher uses the control point displacements to generate an interpolation field for the region, by which the morpher motion is defined. The morpher then uses the interpolation field to translate the mesh vertices to their new positions.

Interpolation Algorithms

Simcenter STAR-CCM+ provides two interpolation methods—BSpline and Radial Basis Functions (RBFs). The locality of the BSpline coefficients and the parallelization schemes makes the BSpline morpher more scalable for parallel computations than RBF—potentially reducing computational expense.

BSpline

Following [18], Simcenter STAR-CCM+ uses a hierarchical, multilevel BSpline approach based on progressively adaptive grids.

The interpolation method is C²-continuous—producing a smoother surface than either C⁰ or C¹-continuous curves.

Consider a simple, 1D example. The diagram below shows a line segment, with four BSpline coefficients at the 'knots' (labeled 0-3), spaced equidistantly along it. The parameter $s$ is the normalized distance in the interval wherein one or more control points may lie:

In this example, the BSpline curves are fitted through control points (not shown) along a line that is superimposed onto the geometry domain. These points are not, by definition, coincident with the model geometry, but instead surround the mesh vertices. A cubic BSpline curve across the central interval (1)-(2) requires a minimum of four coefficients:

Figure 1. EQUATION_DISPLAY

f (x) = \sum_{k = 0}^{3} B_{k} (s) ϕ_{k}

(4845)

where:

$f (x)$ is the fitting curve.
$s$ is the normalized distance along the intervals
$ϕ_{k}$ is the value of BSpline coefficient at knot $k$ , with $k$ =0, 1, 2, and 3.

The weight functions

B_{k} (s)

are Uniform Basis Functions, given as:

Figure 2. EQUATION_DISPLAY

\begin{matrix} B_{0} (s) = {(1 - s)}^{3} / 6 \\ B_{1} (s) = (3 s^{3} - 6 s^{2} + 4) / 6 \\ B_{2} (s) = (- 3 s^{3} + 3 s^{2} + 3 s + 1) / 6 \\ B_{3} (s) = s^{3} / 6 \end{matrix}

(4846)

These act as weighting factors. Each has a proportionate influence on $ϕ_{k}$ , such that $\sum_{k = 0}^{3} B_{k} = 1$ .

The approach is easily extended to 2D and 3D cases:

Figure 3. EQUATION_DISPLAY

\begin{matrix} f (x, y) = \sum_{k = 0}^{3} \sum_{l = 0}^{3} B_{k} (s) B_{l} (t) ϕ_{k l} \\ f (x, y, z) = \sum_{k = 0}^{3} \sum_{l = 0}^{3} \sum_{m = 0}^{3} B_{k} (s) B_{l} (t) B_{m} (v) ϕ_{k l m} \end{matrix}

(4847)

where $B_{k} (s)$ , $B_{l} (t)$ , and $B_{m} (v)$ are the basis functions in the $x, y$ , and $z$ coordinate directions, respectively; these describe knots on a 1D line, 2D square, or a 3D cube, as appropriate.

The BSpline approach implemented in Simcenter STAR-CCM+ uses a hierarchical process—progressively refining the control-point grid until the residual error is within the limits specified by the morpher tolerance. A sequence of corrections is applied at different grid-scales. For simplicity, let the coordinates be named $x_{i}$ , where $i = 1, 2, ... n$ denotes the grid refinement step, and $t$ the target coordinates, with the residual being the difference between these, $r_{i} = t - x_{i}$ . The residual is, in fact, the displacement, which is known at the control points, but unknown elsewhere.

The residual can be represented as a sum of interpolation $I_{i}$ and the remaining residual $r_{i - 1} = I_{i} + r_{i}$ until the desired displacement is reached.

An equation for $t$ is formed as:

Figure 4. EQUATION_DISPLAY

t = x_{0} + \sum_{i = 1}^{n} I_{i} + r_{n}

(4848)

where $x_{0}$ is the coordinates of the point at the start of the process, $I_{i}$ is the interpolation at step $i$ , and $n$ is the last refinement step at which the desired displacement is reached.

The refinement is performed until the morphing error is converged to within the precision of the morpher tolerance. The convergence of such an algorithm is guaranteed, as the proximity data sets are getting smaller and, at a certain refinement level, each of them contains just one control point and exact interpolation is realized.

The hierarchical approach is inefficient, as the number of calculation points grows exponentially with refinement, while many of the calculations are for regions where there is no data. A more efficient approach is to remove these points from the calculation. As follows from Eqn. (4848), the coefficient $ϕ$ is only equal to zero when either:

all contributions from points in the proximity data are zero, or
the proximity data set is empty.

In both cases, it can be assumed that the spline correction at these points is already converged, and that no further corrections are required. Thus, at subsequent refinement levels, the algorithm collects just non-zero coefficients—a process named adaptive interpolation. After refinement, the number of points in the control lattice grows, but, conversely, there are less points where any correction is required.

The picture below shows the control lattice at refinement level 0 (left image), and refined control lattice at level n (right image):

If the control data points are defined on the boundary surface (black line), after a certain refinement level, no inner-region vertices are corrected. The number of operations is of the order $O (M + N) \times n$ , where $M$ is the number of control points (that is, vertices with displacements), $N$ is the number of evaluation points (that is, all vertices in the geometry), and $n$ is the number of refinement levels.

Radial Basis Functions (RBF)

The RBF morpher generates an interpolation field via a system of equations; these are created using the control points, and their specified displacements, where the known displacement ${d'}_{i}$ of each control vertex $i$ is expanded as

Figure 5. EQUATION_DISPLAY

{d'}_{i} = \sum_{j = 1}^{N} f_{b, j} (r_{i j}) λ_{j} + α

(4849)

in which $f_{b, j} (r_{i j})$ is a radial basis function of the form

Figure 6. EQUATION_DISPLAY

f_{b, j} (r_{i j}) = \sqrt{{r_{i j}}^{2} + {c_{j}}^{2}}

(4850)

and $r_{i j}$ is the magnitude of the distance between two vertices:

Figure 7. EQUATION_DISPLAY

r_{i j} = | x_{i} - x_{j} |

(4851)

In these equations, $λ_{j}$ is the expansion coefficient, $x_{i}$ is the position of vertex $i$ , N is the number of control vertices, and $c_{j}$ is the basis constant. The basis constant is set to zero in Simcenter STAR-CCM+. The constant vector $α$ is used to satisfy the additional constraint that

Figure 8. EQUATION_DISPLAY

\sum_{j = 1}^{N} λ_{j} = 0

(4852)

This additional constraint bounds the expansion for large $x$ . The system that is described in Eqn. (4849) is non-singular if $x_{i} \neq x_{j}$ for all $i \neq j$ .

Equations Eqn. (4849) and Eqn. (4852) are solved to give the Cartesian components of all $λ_{j}$ and the components of constant vector $α$ . This results in the desired interpolation field:

Figure 9. EQUATION_DISPLAY

d' (x) = \sum_{j = 1}^{N} f_{b, j} (r) λ_{j} + α

(4853)

where $f_{b, j} (r)$ applies at all vertices in the mesh, that is,

\begin{matrix} f_{b, j} (r) = \sqrt{r^{2} + {c_{j}}^{2}} \\ r = | x- x_{j} | \end{matrix}

This equation can now be used to move mesh vertices by the calculated displacement $d'$ .

The matrices that are involved in equations Eqn. (4849) and Eqn. (4853) are dense, and a straightforward matrix operation would require $O (N^{2})$ operations. An $O (N \log (N))$ method has been developed to solve the linear equation Eqn. (4849) and evaluate equation Eqn. (4853) at many vertex locations. The RBF morpher accelerates the matrix-vector multiplication using a conjugate gradient method, with a special pre-conditioner ([15], [16]), along with a fast multipole method [17].

Boundary Plane Conditions (RBF Morpher)

Simcenter STAR-CCM+ provides two constraints for the vertices that lie on the morphing boundary—the Boundary Plane and the Fixed Boundary Plane constraints. For the fixed boundary plane condition, the vertices which lie on the specified plane are fixed, while for the boundary plane condition they are permitted to move within the plane itself. An effect similar to the fixed boundary plane condition can be realized by using control vertices, but the fixed boundary plane condition itself can be more efficient. The movement of vertices near the constraint planes is such that they join smoothly with the movement of vertices on the plane.

For the case of the RBF morpher only, the movement near the constraint planes can be controlled using the damping distance $L_{p}$ that is described below.

These constraints are imposed before solving Eqn. (4849). In preparation, user-specified displacements $d_{i}$ are transformed to ${d'}_{i}$ by solving

Figure 10. EQUATION_DISPLAY

d_{i} = \prod_{p = 1}^{N_{p}} [s_{p} n_{p} n_{p}^{T} + t_{p} (1 - n_{p} n_{p}^{T})] {d'}_{i}

(4854)

where $n_{p}$ is the normal of the p-th plane and $N_{p}$ is the number of planes. The displacements ${d'}_{i}$ are transformed back to $d_{i}$ together with other mesh vertices at the end of the morphing process. $s_{p}$ and $t_{p}$ are functions that depend on the plane type and the continuity condition you choose. These functions are given as follows:


	$C^{0}$ -continuous		$C^{2}$ -continuous
	$s_{p}$	$t_{p}$	$s_{p}$	$t_{p}$
Fixed Boundary Plane	$f_{0} (l)$	$f_{0} (l)$	$f_{2} (l)$	$f_{2} (l)$
Boundary plane	$f_{0} (l)$	1	$f_{2} (l)$	1

Here:

Figure 11. EQUATION_DISPLAY

f_{0} (l) = {\begin{matrix} 1 \\ 1 - {(1 - l)}^{3} \\ 0 \end{matrix} \begin{matrix} ; l \geq 1 \\ ; 1 > l \geq 0 \\ ; l < 0 \end{matrix}

(4855)

Figure 12. EQUATION_DISPLAY

f_{2} (l) = {\begin{matrix} 1 \\ l^{3} & (10 - 15 l + 6 l^{2}) \\ 0 \end{matrix} \begin{matrix} ; l \geq 1 \\ ; 1 > l \geq 0 \\ ; l < 0 \end{matrix}

(4856)

in which

Figure 13. EQUATION_DISPLAY

l = \frac{n_{p} (x - x_{0, p})}{L_{p}}

(4857)

is the ratio of the signed distance from the p-th plane to $x_{i}$ and the damping distance $L_{p}$ . $x_{0, p}$ is a position vector on the p-th plane and $L_{p}$ is either entered as a finite distance or calculated automatically as $(L_{m i n} + L_{m a x}) / 2$ . $L_{m i n}$ and $L_{m a x}$ are respectively the distances from the plane to the nearest and farthest control vertices on the active side of the plane.

Functions $f_{0}$ and $f_{2}$ in the range from 0 to 1 are plotted below for comparison. A continuity level of $C^{2}$ results in a smooth transition.

The steps that are involved in the RBF morpher stage within a time-step are:

The control vertex displacements $d_{i}$ are transformed to ${d'}_{i}$ according to Eqn. (4854).
Eqn. (4849) is solved with Eqn. (4852) to provide an interpolation field.
The displacement vectors are transformed from ${d'}_{i}$ to $d_{i}$ using Eqn. (4854).

Linear Fitter

The linear fitter can be used with both the BSpline and RBF morphers.

Each displacement component $(d x, d y, d z)$ is represented as a function of coordinates. Using least squares, a hyperplane equation is found for each displacement component.

Consider a problem of fitting the x-displacement, $d x$ :

$d x (x, y, z) = A x + B y + C z + D$

An over-determined system of linear equations can be constructed

$(\begin{matrix} ⋮ & ⋮ & ⋮ & ⋮ \\ x_{i} & y_{i} & z_{i} & 1 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋮ \end{matrix}) (\begin{matrix} A \\ B \\ C \\ D \end{matrix}) = (\begin{matrix} ⋮ \\ {d x}_{i} \\ ⋮ \\ ⋮ \end{matrix})$

which, written in matrix form, is

$P c = d x$

where $P$ is a matrix with point coordinates, $c$ contains the hyper-plane coefficients, and $d x$ is a displacement column.

To solve the problem using the least-squares approach, the normal equations are constructed:

$c = P^{+} d x$ , where $P^{+}$ is the pseudo-inverse matrix.

The RBF morpher uses the Cholesky decomposition for solving the linear equations. For the BSpline Linear Fitter implementation, the matrix is decomposed with singular-value decomposition (SVD).