# Towards industrial use of anisotropic hexahedral mesh adaptation for hydrodynamic flows

Jeroen Wackers, Ganbo Deng, Emmanuel Guilmineau, Alban Leroyer, Patrick

Queutey, Jules Richeux, Michel Visonneau

► To cite this version:

Jeroen Wackers, Ganbo Deng, Emmanuel Guilmineau, Alban Leroyer, Patrick Queutey, et al. Towards industrial use of anisotropic hexahedral mesh adaptation for hydrodynamic flows. SIAM International Meshing Roundtable 2023, Mar 2023, Amsterdam, Netherlands. [hal-03978422](https://hal.science/hal-03978422)

HAL Id: hal-03978422

https://hal.science/hal-03978422

Submitted on 8 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
TOWARDS INDUSTRIAL USE OF ANISOTROPIC
HEXAHEDRAL MESH ADAPTATION FOR
HYDRODYNAMIC FLOWS

Jeroen Wackers

Ganbo Deng

Emmanuel Guilmineau

Alban Leroyer

Patrick Queutey

Jules Richeux

Michel Visonneau

LHEEA Lab, Ecole Centrale de Nantes / CNRS UMR 6598, Nantes, France.
jeroen.wackers@ec-nantes.fr

ABSTRACT

This paper discusses the preparation of the flow solver ISIS-CFD for industrial simulation of free-surface hydrody-
namic flows with adaptive mesh refinement. To accurately capture the free surface and boundary layers, anisotropic
hexahedral mesh refinement and derefinement is performed by local cell division. The adaptation is based on metric
criteria which combine Hessians of the flow variables and free-surface detection. Analytical mesh deformation is used
to account for body movement, while a dedicated local deformation algorithm handles the projection of the adapted
mesh onto the geometry. It is shown that computational setup protocols can be defined to choose the simulation
parameters, without user intervention, for large classes of simulations. Tests on different ships and lifting hydrofoils
show the efficiency and versatility of the method.

Keywords: mesh adaptation, anisotropy, hexahedral, mesh deformation, automatic setup, ship hydrodynamics

1. INTRODUCTION

Industrial application of mesh adaptation is gaining
momentum, thanks to the increasing power and relia-
bility of mesh adaptation algorithms on one side and
the growing industry awareness of their benefits on the
other side. Mesh adaptation was identified as a future
key technology in NASA's CFD Vision 2030 [1] and
is currently used in aerospace by for example NASA
and Boeing [2, 3], and Safran[4, 5]. For hydrodynamic
simulations, mesh adaptation is used with flow solvers
such as OpenFOAM [6, 7], ReFresco [8], StarCCM+
[9], and FINE/Marine [10, 11, 12].

The ISIS-CFD Navier-Stokes solver, developed at Cen-
trale Nantes / CNRS and commercialised by Cadence
Design Systems as part of FINE/Marine, includes an
adaptive grid capability since 2009. Industrial use of
this feature rose sharply around 2018 when the soft-
ware had matured and accumulated experience al-

lowed the user interface to be simplified; it is now
used routinely by about half of the industrial users.
The principal application is mesh refinement to cap-
ture water surfaces [13].

Despite research progress [12, 14, 15], Hessian-based
mesh adaptation with ISIS-CFD has so far seen little
industrial application. Routine simulation of realistic
flows, where the entire mesh is created using adaptive
refinement, presents several mesh-related challenges:

* Most hydrodynamic flow solvers use hexahedral meshes with highly stretched cells to capture boundary layers, wakes, and the free surface. Thus, for mesh adaptation, anisotropic refinement of hexahedral meshes is essential.

* Even in steady simulations, the position of a ship is rarely fixed; at the very least, the freedom to move to an equilibrium position is often needed.
Furthermore, refined meshes must be adjusted to
conform exactly to the CAD topology. Both is-
sues require mesh deformation procedures.

* For industrial simulation to be successful, the simulation setup should be simple and user guidelines should be reliable and able to handle a large variety of different cases, without requiring expert knowledge.

This article discusses the work done to prepare Hessian-based adaptation for industrial application. The mesh refinement and deformation capabilities are available today, while on-going work intends to show that straightforward computational setup protocols can be established for all the typical classes of hydrodynamic simulations. As examples, protocols have been presented for displacement hulls like cargo ships and military vessels [12] and for lifting hydrofoils used on modern high-speed sailing craft [16].

The paper is organised as follows. Section 2 presents the flow solver and its specific features for hydrodynamic problems. Section 3 then describes the grid adaptation method, concentrating on those aspects that are particular for hexahedral meshes and hydrodynamic flows. Mesh deformation algorithms are discussed in section 4, while section 5 studies automatic computation setup and parameter choices for large classes of computations. Finally, in section 6 the method is demonstrated for ships and hydrofoils. The conclusion (section 7) discusses the readiness and the usefulness of the procedure for industrial application.

## 2. THE ISIS-CFD FLOW SOLVER

ISIS-CFD is an incompressible unstructured finite-
volume solver for multifluid flow [17]. The veloc-
ity field is obtained from the momentum conserva-
tion equations and the pressure field is extracted from
the mass conservation constraint transformed into a
pressure equation. These equations are similar to the
Rhie and Chow SIMPLE method [18], but have been
adapted for flows with discontinuous density fields.
The method features turbulence models like Menter's
$k-\omega$ SST [19] but also an anisotropic EASM model
[20] and hybrid RANS/LES models of the DES fam-
ily. Imposed and resolved motions of rigid or flexible
bodies can be computed as will be shown in section
4, while large-amplitude multibody motion is possible
with sliding or overset meshes. The code is fully paral-
lel using the message passing interface (MPI) protocol.

Free-surface flow is simulated with a water-air mixture
model [17]. Instead of explicitly reconstructing the
water surface, its position is obtained implicitly from a
conservation equation for the volume fraction of water.
This equation is discretised with specific compressive

discretisation schemes to keep the transition between
air and water close to a discontinuity (see figure 1).
At the surface, regular and fine meshes are required to
accurately resolve the volume fraction.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Water volume fraction on the hull of a ship and its centre symmetry plane. Red: volume fraction 1, blue: volume fraction 0.</figcaption></figure>

The unstructured discretisation of all the flow equa-
tions is face-based. While the unknown state variables
are cell-centered, the fluxes which form the systems of
equations are constructed face by face, so cells with
an arbitrary number of arbitrarily-shaped constitutive
faces are accepted. This enables for example adaptive
mesh refinement by local cell division: a cell with a
divided face and hanging nodes is treated naturally as
several independent faces.

Since free-surface flow simulations need regular meshes for good accuracy, ISIS-CFD simulations are mostly performed on unstructured hexahedral meshes like the ones generated by Hexpress from Cadence Design Systems. These meshes feature local refinement with hanging nodes and body-fitted boundary layer grids. An example of such a mesh is shown in figure 3.

## 3. MESH ADAPTATION

### 3.1 Anisotropic hexahedral refinement

The natural mesh adaptation approach for hexahedral
meshes is adaptive refinement by local cell division,
starting from a coarse initial mesh (figure 3). Hy-
drodynamic flows contain many flat or stretched flow
features, such as vortical wakes or the water surface.
To efficiently capture these, anisotropic mesh adapta-
tion is a necessity [21]. Therefore, the adaptation in
ISIS-CFD [14, 15, 22] is performed by anisotropic di-
vision: cells can either be split in one or in several
directions. For unsteady flows or to accommodate the
convergence of a steady flow, existing refinement can
be undone. This derefinement is also anisotropic and
refinement can be created and removed in a different
order [23].

During a simulation, the refinement procedure is called
repeatedly, until the flow has converged and the mesh
is no longer changed by the adaptation steps. Since a
cell can only be refined or derefined once during each
<figure><img src="image_3.png" /><figcaption>Figure 2: An example of an adapted mesh. Top to bottom: free surface, details of the Y-symmetry plane, X-cut at the stern. The case is the DTMB 5415 at Fr = 0.41 of section 6.2, with $T_{r,H}$ = 0.07.</figcaption></figure>

step, the adaptation is an iterative process which requires many calls (40–50) to converge. Therefore, the flow is not fully solved on each adapted mesh. Instead, the adaptation procedure is called every 10 – 25 time steps and the flow converges together with the mesh.

Division-based mesh adaptation has its limits. Specifically, the initial grid has a large influence on the adapted mesh, since (1) the division procedure can only reduce the initial cell sizes by a power of two, so not all cell sizes can be obtained, and (2) the adaptation cannot modify the orientation of the cells. This means that the initial grid with its body-aligned boundary layer grid has to be created with care; coarse but orthogonal grids and thick boundary layer meshes produce optimal results. An illustration of an initial and adapted ship mesh is provided in figures 3 and 2, while figure 4 shows a mesh for a hydrofoil.

<figure><img src="image_4.png" /><figcaption>Figure 3: Initial mesh corresponding to the bottom image of figure 2.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: An adapted mesh for a lifting hydrofoil geometry, coloured with the water volume fraction.</figcaption></figure>

## 3.2 Refinement criterion

The decision where to refine is based on a metric refinement criterion [24, 25, 26], which is applied to division-based refinement of hexahedrals as follows. Starting from a symmetric tensor field $C(x, y, z)$ computed from the flow, the grid is adapted until the three dimensions $\mathbf{d}_{i,j}$ of each hexahedral cell $i$, which are the vectors between the face centres for the three pairs $(j = 1, 2, 3)$ of opposing faces, satisfy:

$$ |C_i \mathbf{d}_{i,j}| = T_r, \quad \forall i, \forall j \in (1, 2, 3), \quad (1) $$

as closely as possible. Here $|\cdot|$ indicates a vector length. The constant threshold $T_r$ functions as a global specification of the mesh fineness. Cells are refined in a direction $j$ whenever the modified distance $|C_i \mathbf{d}_{i,j}|$ exceeds $T_r$ and derefined whenever it is below $T_r/2.5$, to make sure the derefined cells are not immediately refined again. This definition leads to anisotropic refinement whenever the cells are aligned with a non-uniform metric, which is typically the case at the water surface. Although the formulation (1) differs slightly from the approach of George et al. [24], who use a metric $M_P$ to directly compute a modified vector length $\mathbf{u}^T M_P \mathbf{u}$, this difference is cosmetic. If $C$ is chosen such that
$CC = M_P$, the resulting modified distances are the same between the approaches.

To create fine meshes entirely using adaptive refinement, the refinement criterion $C$ must react to all the flow features which are relevant for ship resistance. For the orbital flow fields in the waves, the Hessian of the pressure is a suitable indicator [14]. The accurate resolution of the boundary layer and wake is needed to correctly predict both viscous and pressure forces, which suggests refinement based on the velocity Hessians as well. Thus, the criterion is based on the Hessians of both the pressure and the velocity, weighted in the way in which they appear in the flux [15]. A common weight $\rho V$ is assigned to all the velocity Hessians, where $V = \sqrt{u^2 + v^2 + w^2}$:

$$C_H = \left( \max(\|\mathcal{H}(p)\|, \rho V \|\mathcal{H}(u)\|, \rho V \|\mathcal{H}(v)\|, \rho V \|\mathcal{H}(w)\|) \right)^\varphi. \quad (2)$$

$\mathcal{H}$ is the Hessian operator. $\| \cdot \|$ indicates a matrix having the same eigenvectors as the original one and the absolute values of its eigenvalues; the power $\varphi$ is again applied to the eigenvalues, where $\varphi = 0.5$ is used in practice. Finally, the maximum of two tensors is computed following the procedure in [22], which is based on [25]. The convection-derived weight $\rho V$ for the velocity Hessians is open to discussion, since these mainly indicate diffusion-based errors. This is a subject of on-going study.

The power $\varphi$ can be used to give more or less weight to the strongest flow features. For $\varphi = \frac{1}{2}$, the equivalent metric of [24] is $M_P = C_H^2 \sim \|\mathcal{H}(p)\|, \|\mathcal{H}(u)\|$, which is what [24] themselves use. Hessian metrics are linked with interpolation error control and, following [26], the use of a Hessian without weighting by its determinant implies a minimisation of the interpolation error in the $L^\infty$ norm, i.e. an equidistribution of the error. However, for our hexahedral finite-volume discretisation, the interpolation error may not be proportional to the discretisation error, so the latter is not equidistributed.

Finally, a smooth anisotropic mesh at the free surface is required for the volume fraction equation. Since the pressure gradient is discontinuous at the free surface due to hydrostatic forces, the pressure Hessian is undefined there. To solve this issue, the Hessian criterion is extrapolated from below through the surface region and a second criterion [22] creates the refined mesh at the free surface: this criterion is non-zero whenever the water volume fraction $\alpha$ is neither 0 nor 1. Directional refinement normal to the surface is obtained from normal vectors $\mathbf{v} = \nabla \alpha_A / |\nabla \alpha_A|$, where $\alpha_A$ is a smeared volume fraction field created with Laplacian

smoothing. The criterion then becomes:

$$C_S = \begin{cases} \mathbf{v} \otimes \mathbf{v} & \text{if } 0.1 \le \alpha_A \le 0.9, \\ \mathbf{0} & \text{otherwise.} \end{cases} \quad (3)$$

The choice is based on $\alpha_A$ to create a buffer layer of a few cells around the surface position for safety.

The two tensor criteria are combined into one [14] by taking a weighted maximum of the tensors. Since the free-surface criterion always has a unit eigenvalue, a weighting constant $c$ is applied only to the Hessian criterion:

$$C_C = \max(C_S, cC_H). \quad (4)$$

In practice, choosing $c$ directly is inconvenient, so the user specifies separate thresholds $T_{rS}$ and $T_{rH}$ for the free-surface and Hessian criterion. These are then used to set $T_r = T_{rS}$ and $c = \frac{T_{rS}}{T_{rH}}$. The choice of $T_{rS}$ and $T_{rH}$ is discussed in section 5.2.

In figure 2, combined refinement at and below the surface can be seen in waves, while the mesh around the hull is refined at the bow and in the boundary layer. Figure 4 shows refinement in the pressure peaks at the leading and trailing edge, as well as anisotropic refinement in the wake.

## 3.3 Quality measures

Certain measures are applied in the mesh adaptation process, to preserve the quality of the adapted mesh [22]. The two main aspects for this article concern post-treatment of the refinement criterion. First, the eigenvalues of the criteria $C$ can be limited to a user-specified value. This imposes a minimum cell size, below which the cells are never refined. This option is useful, for example, to prevent spurious refinement near flow discontinuities. Furthermore, the metric can be set to zero in the $x-$, $y-$ and $z$-direction on the outside of limiting boxes, to prevent refinement close to the border of the domain. The specific use of these features is explained in section 5.

Other quality measures act directly on the refined mesh, like the difference of one refinement level allowed between neighbour cells, or the obligation to keep the same refinement in a column of boundary layer cells. For brevity, these are not further described here.

## 4. MESH DEFORMATION

For realistic flow simulation in hydrodynamics, mesh deformation is essential. The first reason for this is, that the forces on a ship are highly dependent on its position with respect to the water surface: a ship that sits slightly deeper in the water has a much higher drag. Thus, even for steady simulations, it is more realistic to let a ship move to its dynamic equilibrium
position, than to fix the position. For fast ships and hydrofoils, this can lead to significant motion. Furthermore, for hydrofoils, naval architects are observing that the deformation of the geometry needs to be taken into account to get accurate forces. This requires mesh deformation to follow the motion of rigid and flexible bodies. Finally, division-based mesh refinement creates meshes where the newly inserted nodes are not body-fitted. To correct this, the mesh is locally deformed to project it back onto the geometry.

Since these mesh deformations do not have the same requirements, different deformation methods are employed in ISIS-CFD for rigid-body motion, flexible-body deformation, and for the projection of adapted meshes. This section briefly describes these methods. For the coupling with the flow equations, ISIS-CFD uses an Arbitrary Lagrangian Eulerian formulation [17] which allows to solve the flow on moving meshes.

## 4.1 Rigid-body motion

For rigid-body ship motions [27], the challenge is to preserve a good mesh at the free surface. Therefore, rigid-body motions may be divided in those that do not displace the free surface (x- and y-translation, z-rotation, with z the vertical axis), and those that do (z-translation, x- and y-rotation). For single-body flows, the former can be treated by rigid-body displacement of the entire mesh. The latter require mesh deformation where the outer boundaries of the mesh are not deformed, to make sure the mesh at the free surface moves as little as possible. If the free-surface grid is created with adaptive refinement, z-translation can also be treated with rigid-body movement, but x- and y-rotation always require deformation.

The deformation is performed with analytical weighted regridding: the required motion is applied to all the nodes of the mesh, multiplied by a weighting coefficient $w(x)$ which goes to zero on the outer boundaries so that only the inner part of the mesh moves (figure 5). An initial coefficient $w_0(x)$ is computed as a solution of Laplace's equation with Dirichlet boundary conditions 1 on the body surface and 0 on the outer boundaries. This coefficient is then modified to produce the actual weighting:

$$w(x) = \max(w_0(x)^p/c_{\max}, 1). \quad (5)$$

The power $p$ reduces the zone where the mesh deforms significantly, to leave the free surface far away from the ship undisturbed. The factor $c_{\max}$ on the other hand, imposes a zone of rigid-body motion $w = 1$ close to the hull, to preserve the shape of the boundary layer mesh. For ship simulations with mesh refinement, the default values are $p = 3$ and $c_{\max} = 0.8$.

<figure><img src="image_6.png" /><figcaption>Figure 5: Deformed mesh for rotation of the central body (left) and the weighting coefficient $w(x)$ (right).</figcaption></figure>

## 4.2 Flexible-body deformation

The deformation of flexible bodies cannot be performed with the analytical weighting presented above. To perform fluid-structure interaction (FSI), two deformation approaches are included in ISIS-CFD. The first, presented by [28, 29], is to divide the mesh into layers around the hull. The displacement of each face in a layer is copied to the layer around it, defined using the cell-face connectivity, with a diffusion effect to homogenise the deformation and a weighting similar to the one above, to obtain zero deformation on the outer boundaries. The second approach, used with the modal FSI mentioned in section 6.4, is to deform the mesh with radial basis functions distributed over the geometry, to distribute the motion of each of these parts [30]. These techniques go beyond the scope of this article and are not described in detail.

## 4.3 Post-refinement projection

In division-based refinement, new nodes are inserted in the centres of existing edges and faces. For curved geometries, these nodes are not placed on the geometry, so the refined mesh has to be deformed to fit the geometry. These deformations are highly local, since each misplaced boundary node is surrounded by existing nodes which are placed correctly on the geometry. The layer-based flexible-body mesh deformation mentioned above can be (and has been) used for this projection, but it is not ideal. First, it is inefficient to use such an expensive global deformation to solve essentially local problems. And second, the method does not propagate the local deformations very far from the boundary nodes, which leads to strongly deformed and potentially twisted cells.

Therefore, a dedicated local mesh deformation has been introduced for post-refinement projection. This method starts with the projection of the boundary
<figure><img src="image_7.png" /><figcaption>Figure 6: Ellipsoid around a boundary node, with wall-normal vector, projection displacement $dx_1$, and orthogonality-preserving displacement $dx_2$.</figcaption></figure>

nodes onto the CAD geometry, represented by a triangulation. Each node is projected orthogonally onto the triangle which gives the smallest displacement from the current node position. In the second step, the interior mesh is deformed. For this, ellipsoids (figure 6) are defined around each boundary node to project, with their long axes oriented along the normal direction to the geometry and with two short axes, of equal length, parallel to the wall. The ellipsoids are dimensioned such that they contain no other boundary nodes. This is obtained by starting from an infinite ellipsoid for each point and then searching the cells nearby. Each time a boundary node is found, either the long or the short axis length is reduced to fit the ellipsoid surface through this node. The axis to adjust is chosen as the one which maintains the largest volume for the ellipsoid.

Once all the ellipsoids are defined, the boundary node deformations ($dx_1$ in figure 6) are applied to all the nodes in each ellipsoid, with a linear weighting so that the outer boundary of the ellipsoid does not deform. Furthermore, a secondary shear translation $dx_2$ is applied to the interior nodes with a quadratic weighting, such that the interior node closest to the boundary node is projected towards the normal vector. Without this secondary lateral deformation, when faces are not aligned with the wall-normal direction, large normal displacements may reduce the orthogonality of the cells to the point that faces can become inverted. Let $r(\mathbf{x})$ be the non-dimensional radius of a point in the ellipsoid, such that $r = 0$ in the centre and $r = 1$ on the boundary. Then the total displacement of $\mathbf{x}$ is:

$$dx(\mathbf{x}) = (dx_1 + dx_2r(\mathbf{x})/r_2)(1 - r(\mathbf{x})) \quad (6)$$

with $r_2$ the radius of the closest node.

Figure 7 gives an example of the process. It shows the undeformed mesh with the ellipsoids and the final deformed mesh. Around the leading edge of the airfoil, secondary displacement is visible.

<figure><img src="image_8.png" /><figcaption>Figure 7: Example of a refined grid with ellipsoids (left) and corresponding projected grid (right).</figcaption></figure>

## 5. REFINEMENT PROTOCOLS AND COMPUTATIONAL SETUP

To use these algorithms for routine simulations, default values for the simulation parameters are needed, which must be straightforward, reliable, and valid over a large range of ship lengths, velocities, hull shapes, etc. These guidelines are based on physical arguments, notably dimensional analysis, to reduce the number of input variables that have to be considered, combined with experience to find sensible default settings for the remaining parameters [23].

Simulation protocols have been defined in [12] for displacement ships and in [16, 31] for lifting hydrofoils. Without repeating these protocols fully, we provide examples here of how the guidelines are defined, focusing on the mesh adaptation. Figures 2 and 4 show meshes created with these protocols.

### 5.1 Reference length

For non-dimensional guidelines, a key issue is the choice of the reference length, since many parameters such as the threshold $T_r$ of equation (1), but also the flow domain size, the boundary layer mesh, certain motion laws, etc. depend on a length scale. To obtain similar meshes for similar flow features, it is important that the length scale represents the characteristic size of the flow, which is typically the length over which the boundary layer develops.
<figure><img src="image_9.png" /><figcaption>**Figure 8:** KCS, Free-surface elevation with limiting box $0.3L$ behind the stern (left) and $2L$ behind the stern (right). Hessian threshold $T_{rH} = 0.05L$. The dashed line indicates the aft edge of the limiting box.</figcaption></figure>

For ships, this is the hull length $L$ which is also the main dimension of the ship. For hydrofoils however, the largest dimension is the span but the characteristic boundary layer length is the chord. Therefore, contrary to existing practices which take the span, [16] chooses the maximum chord $c$ as the reference.

Below, all lengths are expressed in terms of $L$ for ships and of $c$ for hydrofoils, to underline that the physics of these cases is different. However, both parameters denote the reference length for their respective cases.

## 5.2 Thresholds

The mesh density is specified here by the threshold $T_r$. Many authors adjust such mesh density parameters to obtain a specified number of cells, for example by normalising the metric in such a way that the threshold directly specifies the target number of cells [26].

To capture free-surface waves, we prefer a different approach. For non-adapted ship meshes, established ISIS-CFD guidelines specify target cell sizes at the free surface in terms of the ship length $L$. To obtain the same behaviour for adaptive refinement, the free-surface refinement criterion is defined with unit vectors in equation (3), so its threshold $T_{rS}$ equals the desired cell size normal to the free surface, which is chosen as $L/1000$ like for non-adapted meshes. For hydrofoils, [16] suggests $c/64$ which gives accurate forces, without reproducing all the details of the free-surface flow. The cell sizes of the initial mesh (section 3.1) should be chosen such that this target cell size can actually be obtained by cell division.

To keep the criteria compatible, a similar behaviour is sought for the Hessian refinement criterion: the mesh

<figure><img src="image_10.png" /><figcaption>**Figure 9:** KCS, effect of the minimum cell size (a) and limiting box (b) on the resistance convergence with $T_{rH}$. The blue lines show the uncertainty interval of the experimental results [32].</figcaption></figure>

density on the geometry should be proportional to the threshold value, independent of the case parameters. Therefore, the criterion is non-dimensionalised using the reference length and velocity. As a result, the choice for $T_{rH}$ no longer depends on the hull length and velocity; only a (weak) dependence on $Fr$, $Re$, and the hull shape remains, which can be ignored for a specific class of geometries (such as displacement ships, fast ships, or hydrofoils). For ships, the tests in [12] suggest the range $T_{rH} \in [0.2L, 0.025L]$ for coarse to fine grids. For hydrofoils, [16] recommends $T_{rH} \in [0.5c, 0.0625c]$.

## 5.3 Cutoff filter and limiting box

Adaptive refinement based on the combined refinement criterion captures flow details in the entire domain. However, for the accurate evaluation of resistance, not all these details may be required. Thus, the criterion is modified in two ways.

The resolution can be controlled by specifying a minimum cell size, below which cells are no longer refined.
A large minimum cell size works like a cutoff filter: while the adaptive refinement still captures the main flow features, very fine cells to resolve small details are not created. For the KCS test case of section 6.1, figure 9a shows the dependence of the resistance on the minimum cell size. This parameter (chosen small in earlier tests, around $10^{-4}L$) has a pronounced influence: for a larger minimum cell size, the number of cells is strongly reduced, for only a modest loss of precision. Also, the solutions for all minimum cell sizes converge to the same value when $T_{rH}$ is reduced. Thus, for resistance simulations, a large minimum cell size, i.e. $L/1000$ like the free-surface threshold, is preferable. For hydrofoils, the free-surface size was found to be too large so a minimum size of $c/256$ is recommended.

Furthermore, resolving the wake and the wave field behind the ship may not be needed for resistance computations. Therefore, limiting boxes were tested for the KCS which forbid horizontal refinement from a given distance behind the stern. Preventing refinement aft of $0.3L$ behind the stern removes the entire far wake field (figure 8). However, compared to a $2L$ box which preserves most of the waves, the box at $0.3L$ does not significantly alter the forces while reducing the number of cells by up to 40% (figure 9b). An even shorter box ($0.1L$ behind the stern) changes the converged resistance, so it is not adopted. The box at $0.3L$ is the best compromise. Similar tests suggest a box at $1c$ behind the trailing edge for hydrofoils.

# 6. TEST CASES

Without attempting a full test of the simulation approach, this section provides four examples of the variety of test cases which it can handle. Where available, the results are compared with experiments and simulations on non-adapted meshes. More complete tests, including further details of the cases presented here, are provided in [12] for ships and in [16] for hydrofoils.

<figure><img src="image_11.png" /><figcaption>**Figure 10:** Hull shapes of the KCS (a) and DTMB 5415 (b) models.</figcaption></figure>

## 6.1 KRISO Container Ship

A modern cargo ship, the KRISO Container Ship (KCS) [32], is simulated in model-scale towed condition at $Re = 1.257 \cdot 10^7$ and $Fr = 0.260$. The $k-\omega$

SST turbulence model with a wall law is used. Figure 11 shows the convergence of the resistance with the Hessian threshold $T_{rH}$, as well as the numerical uncertainty evaluated with the grid-convergence based method of Eça and Hoekstra [33] (where $T_{rH}$ is used as the measure for the grid size, instead of an averaged cell size). A power-law like convergence is obtained and the uncertainty estimations, which overlap for the finer grids, make sense. Figure 12 shows that the free surface also converges rapidly and that at $T_{rH} = 0.05$ it is very close to the finest threshold. Thus, varying $T_{rH}$ is an effective way to obtain grid convergence and the selected range for $T_{rH}$ is appropriate.

<figure><img src="image_12.png" /><figcaption>**Figure 11:** KCS resistance, grid convergence as a function of the number of cells, with estimated numerical uncertainty. The blue lines show the uncertainty interval of the experimental results [32].</figcaption></figure>

The results are compared with simulations on Hexpress meshes created with the C-Wizard automatic setup tool, which provides uniform cell sizes over the hull. Also, a combined Hexpress and systematic grid refinement (SGR) approach is shown, where the three coarsest grids are generated by Hexpress, using manual refinement at the bow and stern. The finest grids are created from these coarser grids by refining all cells once, using the adaptive refinement algorithm. Figure 11 compares the three results. The adapted series converges faster to its final value than the other two. A given numerical accuracy is therefore obtained on coarser grids: the same distance from the converged value is typically obtained with 1.5 to 2 times fewer cells than the C-Wizard series. Furthermore, the adapted and SGR series converge to roughly the same values for the resistance, which underlines the reliability of both series. The meshes of the C-Wizard series are not similar: in Hexpress, the boundary layer meshes become thinner as the outer mesh is refined. Therefore, the series converges to a lower value than the other two.
<figure><img src="image_13.png" /><figcaption>**Figure 12:** KCS, mesh convergence of the free-surface elevation. Top: result for $T_{rH} = 0.025$. Bottom: the difference of the results for $T_{rH} = 0.2L, 0.1L, 0.05L$ with respect to $0.025L$.</figcaption></figure>

## 6.2 DTMB 5415

The protocol (developed for the KCS) is applied unchanged to the DTMB 5415 model in towed condition at $Fr = 0.41$ and $Re = 1.74 \cdot 10^7$ following the measurements of [34]. This case represents a slender military vessel at full speed, where heavy wave breaking occurs both at the bow and the stern (figure 13).

The adapted-grid resistance is again compared with C-Wizard uniform grids and meshes with manual refinement and SGR (figure 14). The adapted series converges quickly and produces good uncertainty estimations; however, the converged resistance is smaller than for the experimental results, likely due to the difficulty of modelling the breaking wave behaviour. The non-adapted series seem to agree better with the experiments, but they converge to the same underestimated solution – albeit more slowly. Thus, the faster convergence and smaller uncertainty mean that the adapted-mesh results are the most reliable.

<figure><img src="image_14.png" /><figcaption>**Figure 13:** DTMB 5415, wave pattern for $T_{rH} = 0.025L$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14:** DTMB 5415, grid convergence and uncertainties for three series of meshes. The blue line represents the experimental result from [34].</figcaption></figure>

## 6.3 KCS in self-propulsion

To test the versatility of the ship refinement protocol, it is applied to a case for which it was never intended: the KCS in self-propulsion condition [35] with a body-force model to represent the propeller. This is a challenge since the mesh refinement has never been tested for a force field before. Also, a ship propeller is tiny with respect to the ship itself. Can something so small be captured by mesh adaptation that is configured for the flow around the ship?

<figure><img src="image_16.png" /><figcaption>**Figure 15:** KCS self-propulsion: axial velocity in the y-symmetry plane for $T_{rH} = 0.1L$ and location of the body force (pink box).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 16:** KCS self-propulsion: grid convergence of the resistance; the blue line represents the experimental result from [35].</figcaption></figure>

The initial Hexpress mesh is equal to the one from section 6.1; no particular refinement is applied around the actuator disk. Still, the adapted grid (figure 15) captures the propeller flow well, even at the coarse threshold $T_{rH} = 0.1L$. Thus, refinement based on the ship length scale is appropriate for a detail like the actuator disk. However, the mesh size throughout the disk region is equal to the minimum cell size. While this mesh is fine enough for the actuator disk, it is possible that the minimum cell size prevents the convergence to a mesh-independent solution. All resistance results (figure 16) agree well with the experiments; the difference between the series is similar to the DTMB 5415 for example. These force results confirm that the mesh adaptation protocol can be used for self-propulsion.

## 6.4 Flexible META hydrofoil

Finally, the hydrofoil protocol of [16] is tested. Like for the ships, this protocol has been established for rigid geometries. However, carbon-fiber foils undergo significant deformation under load and a consensus is developing among naval architects that this deformation should be taken into account for the evaluation of the forces. Therefore, the hydrofoil protocol is applied here with fluid-structure interaction, accomplished via a modal approach [30]: a structural solver is used to determine the deformation eigenmodes of the foil, then the flow solver computes the deformation as a linear combination of the first eigenmodes. Thus, no direct coupling of the structure and fluid solvers is needed.

The test case is a hydrofoil for an ocean-going catamaran, with a chord $c = 0.90m$ and span of about 2.2m, simulated at 10m/s. The foil has never been constructed so there are no experiments available. Since no structure was ever designed, the foil is considered to be made of a homogeneous isotropic material with a

Young's modulus of 60GPa. This simplification affects only the results, not the simulation procedure.

<figure><img src="image_18.png" /><figcaption>**Figure 17:** Deformed (red) and rigid (green) shape of the META hydrofoil.</figcaption></figure>

The fluid-structure interaction leads to significant deformation, as shown in figure 17 where the shape is compared with the undeformed foil shape. Both have been simulated with the same hydrofoil protocol. The resistance (figure 18) shows smooth convergence and similar behaviour for the rigid and flexible foils, with good uncertainty estimations. The uncertainty on the finest (23M cell) grid is about 1.6%. All deformed cases have about 200k more cells than the equivalent rigid case; this is likely because the deformed mesh is not perfectly aligned with the free surface, so it requires more isotropic cells to capture the surface.

Finally, figure 19 shows the difference in the drag between the flexible and rigid case. Although the convergence of this difference is not monotonic, the prediction is acceptable even on the coarsest grids. The more general conclusion is that the simulations allow to consistently determine the difference in performance between very similar geometries, even on coarse grids.
<figure><img src="image_19.png" /><figcaption>**Figure 18:** Convergence of the resistance for the rigid and flexible hydrofoils, with uncertainty.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19:** Convergence of the difference in resistance between the rigid and flexible hydrofoil.</figcaption></figure>

This is crucial for automatic shape optimisation and it demonstrates the reliability of the simulations.

# 7. CONCLUSIONS AND DISCUSSION

In preparing the industrial application of combined Hessian-based and free-surface mesh adaptation with ISIS-CFD, the first question is if such computations are possible today. For the adaptation algorithm itself, this is probably true since mesh adaptation is already in industrial use. The mesh deformation is a standard industrial technique in ISIS-CFD and has long since proven itself. This means that the generality and robustness of the user guidelines are the crucial issue. For calm-water resistance of displacement hulls, the tests in [12] show that the protocol can be applied unmodified to hull shapes ranging from slender to full, and to speeds from slow steaming to military top speeds. Furthermore, the protocol is applied successfully to body-force propeller models, a case for which it was never intended. Larger variations of the test case, such as fast planing hulls, were not tried. However, the hydrofoil test shows that a very similar protocol is adequate for these radically different geometries. Most calm-water cases are probably close enough to what was tested here, that even if the protocols in section 5 have to be adapted for them, it is

likely that suitable protocols exist. Thus, routine mesh adaptation for resistance simulations appears realistic today.

Does the proposed mesh adaptation offer advantages with respect to existing meshing approaches? The tests show that adapted meshes provide excellent grid convergence, often obtaining 1% difference with the converged resistance on grids of about 1M cells, which is much faster than on non-adapted meshes. Furthermore, mesh adaptation makes the simulation more reliable, since it automatically ensures that the mesh is suitable for the flow. While the correspondence of the resistance with experiments is often worse than for non-adapted meshes, the differences were traced to problems with the non-adapted meshes, such as insufficient resolution of breaking waves. Thus, this is not a numerical problem but a modelling issue. Finally, the grid convergence in all cases presented here is smooth, which leads to reliable convergence-based uncertainty estimations; this was also observed by [15] for mono-fluid flows and by [26] for wings and aircraft using compressible flow and tetrahedral mesh adaptation. Thus, the good uncertainty estimations are likely a property of the metric-based refinement. This is an important advantage for unstructured grids.

Open questions remain and the simulation protocol will probably be improved and extended. Still, mesh adaptation has proved its generality and the potential advantages are obvious, so it appears ready for industrial use. Furthermore, hydrofoil simulation in particular has the advantage of being a relatively new application. Many naval architects do not have established simulation practices and may therefore be more likely to adopt mesh adaptation as a standard practice. Thus, the hydrofoil may be the break-through application which will launch Hessian-based adaptation for hydrodynamic simulation. The future will tell.

# ACKNOWLEDGEMENTS

We thank the Institute Carnot MERS (ORUP project) and the research directorate of Centrale Nantes for funding the hydrofoil investigation.

# References

[1] Slotnick J., Khodadoust A., Alonso J., Darmofal D., Gropp W., Lurie E., Mavriplis D. "CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences." Tech. Rep. NASA/CR-2014-218178, NASA, 2014

[2] Balan A., Park M., Anderson K., Kamenetskiy D., Krakos J., Michal T., Alauzet F. "Verification of Anisotropic Mesh Adaptation for Turbu-
lent Simulations over ONERA M6 Wing.” *AIAA J*, vol. 58, 1–16, 2020

[3] Park M., Balan A., Clerici F., Alauzet F., Loseille A., Kamenetskiy D., Krakos J., Michal T., Galbraith M. “Verification of Viscous Goal-Based Anisotropic Mesh Adaptation.” *AIAA SciTech Forum*. Anaheim, CA, 11-21 January 2021

[4] Mercier R., Benard P., Lartigue G., Moureau V. “Dynamic adaptation of tetrahedral-based meshes for the simulation of turbulent premixed flames.” *17th International Conference on Numerical Combustion*. Aachen, Germany, 2019

[5] Nardoni C., Bordeu F., Cortia J. “Body-fitted discretization, accurate interface tracking and geometrical constraints in topology optimization.” *Adaptive Modelling and Simulation 2019*. El Campello, Spain, 2019

[6] Eskilsson C., Bensow R.E. “A mesh adaptive compressible Euler model for the simulation of cavitating flow.” *MARINE 2011*. Lisbon, Portugal, 2012

[7] Wang Z., Li L., Cheng H., Ji B. “Numerical investigation of unsteady cloud cavitating flow around the Clark-Y hydrofoil with adaptive mesh refinement using OpenFOAM.” *Ocean Eng*, vol. 206, 107349, 2020

[8] Windt J., Klaij C. “Adaptive mesh refinement in MARIN’s viscous flow solver ReFRESCO: implementation and application to steady flow.” *MARINE 2011*. Lisbon, Portugal, 2012

[9] Yilmaz N., Atlar M., Khorasanchi M. “An improved Mesh Adaption and Refinement approach to Cavitation Simulation (MARCS) of propellers.” *Ocean Eng*, vol. 171, 139–150, 2019

[10] Hildebrandt T., Reyer M. “Business and technical adaptivity in marine CFD simulations: Bridging the gap.” *COMPIT ’15*, pp. 394–405. Ulrichshusen, Germany, 2015

[11] Yvin C., Muller P. “Tip vortex cavitation inception without a cavitation model.” *19th Numerical Towing Tank Symposium*. St Pierre d’Oléron, France, 2016

[12] Wackers J., Deng G., Raymond C., Guilmineau E., Leroyer A., Queutey P., Visonneau M. “Adaptive grid refinement for ship resistance computations.” *Ocean Eng*, vol. 250, 110969, 2022

[13] Wackers J. “Industry application of adaptive grid refinement – the case of FINE/Marine.” *Adaptive Modelling and Simulation 2021*. Online, 2021

[14] Wackers J., Deng G.B., Guilmineau E., Leroyer A., Queutey P., Visonneau M. “Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation.” *Comput Fluids*, vol. 92, 209–222, 2014

[15] Wackers J., Deng G.B., Guilmineau E., Leroyer A., Queutey P., Visonneau M., Palmieri A., Livrerani A. “Can adaptive grid refinement produce grid-independent solutions for incompressible flows?” *J Comput Phys*, vol. 344, 364–380, 2017

[16] Richeux J. *Towards the automation of adaptive RANS simulations for hydrofoils*. M.Sc. thesis, Ecole Centrale de Nantes, France, 2022

[17] Queutey P., Visonneau M. “An Interface Capturing Method for Free-Surface Hydrodynamic Flows.” *Comput Fluids*, vol. 36, no. 9, 1481–1510, November 2007

[18] Rhie C.M., Chow W.L. “A numerical Study of the Turbulent Flow Past an Isolated Airfoil With Trailing Edge Separation.” *AIAA J*, vol. 17, 1525–1532, 1983

[19] Menter F. “Two-equation eddy-viscosity turbulence models for engineering applications.” *AIAA J*, vol. 32, no. 8, 1598–1605, 1994

[20] Duvigneau R., Visonneau M. “On the role played by turbulence closures in hull shape optimization at model and full scale.” *J Mar Sci Tech*, vol. 8, 11–25, 2003

[21] Loseille A., Dervieux A., Frey P., Alauzet F. “Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes.” *18th AIAA Conference*. Miami, FL, 2007

[22] Wackers J., Leroyer A., Deng G.B., Queutey P., Visonneau M. “Adaptive grid refinement for hydrodynamic flows.” *Comput Fluids*, vol. 55, 85–100, 2012

[23] Wackers J. *Adaptivity for complex flows*. HDR thesis, Université de Nantes, France, 2019

[24] George P.L., Hecht F., Vallet M.G. “Creation of internal points in Voronoi’s type method. Control adaptation.” *Adv Eng Software*, vol. 13, no. 5/6, 303–312, 1991

[25] George P.L., Borouchaki H. *Delaunay Triangulation and Meshing - Application to Finite Elements*. Hermes, 1998
[26] Alauzet F., Frazza L. “Feature-based and goal-oriented anisotropic mesh adaptation for RANS applications in aeronautics and aerospace.” *J Comput Phys*, vol. 439, 110340, 2021

[27] Leroyer A., Visonneau M. “Numerical methods for RANSE simulations of a self-propelled fish-like body.” *J Fluid Struct*, vol. 20, no. 3, 975–991, 2005

[28] Durand M. *Light and flexible Fluid/Structure Interaction, application to sailing boats*. Ph.D. thesis, Ecole Centrale de Nantes, France, 2012

[29] Durand M., Leroyer A., Lothodé C., Hauville F., Visonneau M., Floch R., Guillaume L. “FSI investigation on stability of downwind sails with an automatic dynamic trimming.” *Ocean Eng*, vol. 90, 129–139, 2014

[30] Mouton L., Leroyer A., Deng G.B., Queutey P., Soler T., Ward B. “Towards unsteady approach for future flutter calculations.” *J Sailing Tech*, vol. 7, 2018

[31] Robin P., Leroyer A., Richeux J., Prémorel D.d., Wackers J. “Starting off the right foot with foil sock approach and AGR criterion.” 24th *Numerical Towing Tank Symposium (NuTTS 2022)*. Zagreb, Croatia, 2022

[32] Kim W.J., Van D.H., Kim D.H. “Measurement of flows around modern commercial ship models.” *Exp Fluids*, vol. 31, 567–578, 2001

[33] Eça L., Hoekstra M. “A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies.” *J Comput Phys*, vol. 262, 104–130, 2014

[34] Olivieri A., Pistani F., Avanzini A., Stern F., Penna R. “Towing tank experiments of resistance, sinkage and trim, boundary layer, wake, and free surface flow around a naval combatant INSEAN 2340 model.” Tech. Rep. 421, IIHR, Iowa, 2001

[35] Hino T., editor. *CFD Workshop Tokyo 2005*. National Maritime Research Institute, Tokyo, Japan, 2005