Analysis of the performance of two moving mesh methods applied to semi-submerged structures with a rotational DOF

2.1. Mathematical model

The governing equations of the flow are composed of the continuity and the momentum equations (Schilichting, 1979), in which the water and air are fluids considered to be incompressible: (1)

where i, j, k =1, 2, ui are the velocity components, t is the time, p is the pressure, ρ is the specific mass, and gi the gravity acceleration components. τij is the viscous stress tensor and μ is the fluid viscosity. The RANS equations (ReynoldsAveraged Navier-Stokes) are obtained by the decomposition of the instantaneous velocity and the pressure of the NavierStokes equations into average and fluctuating components and integrating these equations over time. The turbulence of the flow is considered by means of the Reynolds stress tensor according to the Boussinesq hypothesis. In this study, the turbulence model, which allows closing the equation system and relating the Reynolds stress tensor with the variables of the flow, is the k–ω SST (Shear Stress Transport) turbulence model (Menter, 1965). The free surface movement is defined by using the VoF (Volume of Fluid) technique proposed by Hirt and Nichols (1981), which consists of identifying the free surface position through the volume fraction, f, which is null for the air, 1 for the water and 0.5 for the free surface position. The volume fraction transport equation is given by 

The numerical model of ANSYS Fluent® (2016) is used to develop the 2D wave flume and the semi-submerged moving body, whose sketch is shown in Figure 1. The volume finite technique (Períc and Ferzinger, 1997), in which the variables are defined at the center of the mesh cells, is applied to solve the RANS and VoF equations. The time integration is carried out by using an implicit formulation and a first-order time discretization.

Figure 1. Sketch of cases with semi-submerged bodies with a rotational DOF. DOF: Degree Of Freedom.

The SIMPLEC (Semi-Implicit Method for Pressure Linked Equations-Consistent) algorithm is used for the coupling of pressure and velocity. The convective terms of the momentum equations are determined by the third-order MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws) scheme and the secondorder UPWIND is used for the equations of the k–ω SST turbulence model. The pressure is determined by PRESTO! (PREssure STaggering Option) scheme, which is classically used to model the wave propagation in ANSYS Fluent® (2016). The relaxation factors are 0.3 for the pressure, 0.7 for the momentum, and 0.8 for the turbulence kinetic energy (k) e the specific dissipation rate (ω). In the VoF method, the explicit formulation is used for the volume fraction and the Geo-Reconstruct scheme is used for the air-water interface.

In the SM method, the interfaces between the adjacent zones of the mesh slip each other along the time. Therefore, the intersections between the interfaces are updated each time step. In this study, the moving zone is circular and follows the rotation of the circumference around its center, which slips over the circumference that belongs to the fixed external zone of the domain. Therefore, in this interface, the mesh is non-conform.

The following boundary conditions are imposed: atmospheric pressure on the top of the flume; walls and bottom with no-slip conditions; and k = 10-6 m2/s2 and ω=1 s-1 on the wave maker, the top of the flume, and the end of the flume (Lin e Liu, 1998; Elhanafi et al., 2016). The free surface level at rest, null velocity components, the hydrostatic pressure on the water, and the atmospheric pressure on the air are initial conditions.

In the DM method, the smoothing method is used to deform the mesh along the time in a zone around the moving body. The mesh nodes are displaced to follow the body movement and, at the same time, provide a better geometric quality of the cells. The diffusion-based smoothing technique is employed to move the nodes, which is governed by the diffusion equation given by:

2.2. The rotational DOF equation

A structure with a rotational DOF in a plane subjected to a moment is governed by the equation (Fig. 1):

where u is the mesh node velocity and γ is the diffusion coefficient, which is calculated as a function of the distance between boundaries, as follows

where Ӫ is the angular acceleration of the structure, I is the mass moment of inertia in relation to the hinge, Mf is the moment caused by the reaction forces of the system coupled to the structure due to the friction and/or forces imposed by external devices, Mp is the moment due to the structure weight, and Mh is the hydrodynamic moment imposed by the flow, caused by the pressure, the viscous stress and the buoyance forces applied to the surface of the structure.

2.3. Moving mesh methods

In this study, two different methods of the moving mesh are used: the Sliding Mesh (SM) and the Dynamic Mesh (DM), which are briefly described below.

The movement of the moving mesh zone is governed by the rigid body rotation of the moving boundary around the circumference center. A six DOF system is solved at each instant according to the moving body Law. The moving body equation is imposed by a specific function named DEFINE_ SDOF_PROPERTIES in ANSYS Fluent®, which is wrote in an UDF (User Defined Functions). The moving body equation of a system is calculated by using an implicit updated.

where d is the normalized distance between boundaries and α is the diffusion parameter, whose value used in this study is 1.5, which is recommended by the ANSYS Fluent® manual for body rotations. Equation (6) is discretized by the finite element method and allows obtainig the mesh node velocity. Therefore, the node positions are updated at each instant as follows: 

The mesh cell quality can be evaluated by an orthogonal quality index which varies from 0 to 1, in which the higher the value, the better the cell quality (ANSYS, 2015).

Figure 2. Sketch of the computational domain to the buoy case.

3.1. Buoy

This case is based on the Jung (2004) experiments, which investigated the behavior of a semi-submerged hexahedral buoy in a flume 0.9 m deep and 0.9 m wide. The buoy has only the rotational DOF around the gravity center (mass moment of inertia I=0.236 kg.m2), which is at the free surface at the rest. The buoy has a height and length of 0.1 and 0.3 m, respectively, which fills all the width of the flume and, consequently, allows considering a 2D flow behavior. In this study, two cases are analyzed: the free rotation decay from an inclination of 15°, which identifies its natural oscillation period; and the interaction between the buoy and an incident regular wave with a period T=0.93 s and a height H=0.027 m. Figure 2 shows a sketch of the computational domain, in which the wave flume is 6L + 0.3 m long (L=1.35 m), with the buoy located at the center of the flume. The friction of the moving parts due to the buoy rotation is considered as a viscous friction, in which the damping coefficient of c = 0.3056 Nms/rad is adopted (Calderer et al., 2014). In the case of the free decay test, the left and right boundaries are considered vertical walls, whereas, in the case of the incident wave over the buoy, there is a static wave maker

on the left boundary and a radiation boundary on the right. At the wave maker, velocity components and the volume fraction are imposed, corresponding to the characteristics of the incident wave, by using the active absorption technique (Didier et al., 2017; Teixeira et al., 2017; Didier and Teixeira, 2022). At the radiation boundary, the active absorption technique is also applied to avoid the wave reflection into the flume. Figure 3 shows the mesh, which is composed of a wave propagation zone downstream and upstream of the buoy and a central zone around the buoy. In the wave propagation zone, the mesh is structured and the cell sizes are L/70 in the horizontal direction. There are two transition zones to adapt the diminishing of cells in L/50 near the wave maker/flume and the end boundary, and e (cell size at the buoy) at the adjacent zone around the buoy. In the vertical direction of the wave propagation zone, there is a region 2H high around the free surface at the rest with cell sizes of H/20. From this region to the bottom and the top, there is a stretching of the cell size. This configuration follows previously validated and applied ones described in several works, such as Teixeira et al. (2013), Mendonça et al. (2018), Teixeira and Didier (2023), Güths et al., (2022), and Wiener et al. (2022). The central zone around the buoy has a nonstructured mesh, composed of quadrilaterals, which includes a

Figure 4. Time series of the angle of rotation of the buoy in the case of buoy decay test. circular domain defined by a ray of 0.4 m. When the DM method is used, the mesh is deformable in this zone. In the case of the SM method, this circumference is adopted as a slip boundary. The cell size inside the circle, in which the buoy is inserted, is e=1.23 m, defined after a mesh convergence analysis. At the external part of the circumference, cells around the free surface at the rest has sizes equal to e. The cell sizes increase in direction to the bottom (around 2.5 times e) and the top, with higher increasing, since the fluid is the air in this region. This mesh is used in both studied cases. In the case of the buoy decay, the inclination of 15° is the initial condition about the horizontal direction and the buoy is released to oscillate due to the action of the water buoyancy. Figure 4 shows the time series of the angle of rotation of the buoy by using both SM and DM methods and a comparison with the experimental results of Jung (2004). The amplitude of the angle of rotation diminishes along the time due to the damping caused by both the viscous water forces and the friction on the moving buoy. It can be noticed that the results obtained by both methods are very similar and very close to the experimental one. The oscillation period is the natural period of the system, which is TN=0.93 s, obtained by the numerical simulations and the experiments. The second case deals to the incident regular wave, with T=0.93 s and H=0.027 m, over the buoy. Figure 5 shows the time series of the angle of rotation after the stabilization of the flow obtained by both moving methods. It may be emphasized that the buoy oscillates with the same period of the fundamental period of the incident wave. Besides, the clockwise and counterclockwise oscillations have similar behaviors, which is expected since the wave propagates in deep water. Results are compared with experimental ones developed by Jung (2004). It can be observed that both methods present similar results and agree well with the results obtained by experiments. The maximum and minimum amplitudes of the angle are 11.6° and -11.7°, respectively, for the SM method and 11.4° and -11.4° for DM method. These values are very similar to the experimental ones, which are 11.5° e -11.7°, respectively. Therefore, in terms of accuracy, both methods show a very good performance.

Figure 5. Time series of the angle of rotation of the buoy subjected to an incident regular wave with T=0.93 s and H=0.027 m. In the case of the DM method, it is very important to analyze the mesh quality. Figure 6 shows the orthogonal quality index to the cells that belongs the zone in which the mesh deforms. It can be noticed that the orthogonal quality index before the mesh deformation (Fig. 6a) is very good (0.70–0.95) to excellent (0.95–1.00), according to ANSYS (2015). It indicates that distortions in the cells are not significant. Besides, even in the maximum amplitude of the angle of rotation of the buoy (Fig. 6b), the distortions of cells are not significant and orthogonal quality index is also very good (0.70–0.95) to excellent (0.95–1.00). Another important criterion to evaluate the performance of each type of moving method is related to the computation cost. The simulation by using the DM method demands 1.44 and 1.32 times the computational time obtained by the SM method for the free decay test and the incident regular wave over the buoy, respectively. Therefore, considering that both methods present similar accuracy, the SM method is the most suitable for cases with similar characteristics to the buoy one, due to its lower computational cost. It is important to emphasize that, in these cases, the use of the remeshing technique was not necessary, since the distortions of the mesh due to the movement of the buoy were not significant. In

Figure 6. Orthogonal quality index of the cells at (a) the initial condition and (b) the maximum angle of rotation for the incident wave over the buoy using the DM method. DM: Dynamic mesh.

Figure 7. Sketch of the computational domain for the incident wave over the Pendulor device. cases in which the remeshing technique must be used, the computational time increases and, consequently, this criterion must be re-evaluated.

3.2. Pendulor device

Figure 7 shows the sketch of the computational domain for the case of the incident regular wave over the Pendulor device. The flap of the device with thickness of 0.5 m, and length lp=7 m, in which 2 m are immersed (ap), is located at the center of the flume 6L long and h=10 m deep. The hinge is located at the emerged tip of the flap, which has a specific mass of 1450 kg/m3 and a mass moment of inertia of I=81667 kg.m2. In the wave maker, a second order Stokes regular wave is generated with T=12 s and H=1.5 m. Generally, in this type of device, a hydraulic system is responsible to absorb the mechanical energy (Power takeoff, PTO), in which reaction forces of cylinders cause a damping moment (MPTO). In this work, the viscous damping model is used to represent this moment, as follows: 

where cPTO is the viscous damping coefficient of the hydraulic system, which is 500 kNms/rad in this study. The mean time mechanical power absorbed by the device is given by 

Figure 8 shows the mesh whose configuration is similar to the buoy case, with two zones of wave propagation and a central zone that includes the flap. The mesh of the central zone around the flap is a non-structured one, composed of quadrilaterals with a minimal edge size of 0.062 m (chosen after a mesh convergence analysis) located at the plate boundary and in the region around the free surface. The cell size increases in direction to the bottom and the top flume, with a similar configuration to the previous case. In the case of the SM method, there is a moving circumference that follows the movement of the device and slips on a circumference coupled to the static mesh. When the DM method is used, the mesh is deformable in the central zone. Figure 9 shows the velocity magnitude and the stream lines in a region around the flap for a stabilized periodic flow at four

Figure 9. Velocity magnitude and stream lines around the flap using SM method. (a) In the vertical position and moving counterclockwise, (b) in the maximum angle of rotation, (c) in the vertical position and moving clockwise, and (d) in the minimum angle of rotation. SM: Sliding mesh.

instants by using the SM method: a) in the vertical position and moving counterclockwise, b) in the maximum angle of rotation, c) in the vertical position and moving clockwise, and d) in the minimum angle of rotation. It can be observed that the fluid dynamics is significant around the submerged tip of the flap at every instant. At the instant b), a high velocity intensity can be noted near the tip of the flap and also over a wide area near the bottom of the flume at the rear of the flap. Besides, the emerged tip of the flap is high enough to avoid overtopping due to this incident wave. These results are very similar to the ones obtained using DM method. Figure 10 shows comparisons of the time series of the angle of rotation, angular velocity, and angular acceleration obtained by both methods. It can be noticed that the time

series of both moving methods are very similar. The higher difference is observed in the angular acceleration, which presents a significant non-linear behavior. The angle of rotation oscillates with the predominant frequency equal to the fundamental frequency of the incident wave. The amplitudes counterclockwise, corresponding to the horizontal displacement in the wave propagation direction, are 16.85° and 17.02° for SM and DM, respectively. These values are higher than the amplitudes clockwise, which are 5.26° and 4.61°, respectively. It is expected since the wave propagates on the flume in intermediate water. The angular velocity also has higher amplitudes counterclockwise and shows a higher non-linearity than the angle of rotation. The angular velocity amplitudes counterclockwise and clockwise are 0.130 and

Figure 10. Time series of the (a) angle of rotation, (b) angular velocity and (c) angular acceleration for the Pendulor device.

Figure 11. Time series of (a) damping moment and (b) absorbed power for the Pendulor device. 0.091 rad/s, respectively for the DM method and 0.127 and 0.082 rad/s for the SM method. The angular acceleration has also its predominant frequency equal to the one of the incident wave, although other frequencies are observed significantly. Although small differences between both methods are observed in these time series, the general behaviors are similar. Figure 11 shows the time series of the damping moment (MPTO) and the absorbed power (PPTO) of the device. The results obtained by using SM and DM are similar, with little differences, mainly around the maximum amplitudes clockwise. The damping moment, which directly depends on the angular velocity of the flap, shows a behavior with oscillation predominantly in the fundamental frequency of the incident wave. The amplitudes counterclockwise are 64.91 and 63.37 kNm/m for SM and DM, respectively,

which are higher than those clockwise, 45.81 and 40.88 kNm/m, respectively. The absorbed power is higher during the counterclockwise than the clockwise, following the behavior of the angular velocity and the damping moment. The power amplitudes counterclockwise are 8.43 and 8.04 kW/m for SM and DM, respectively, and 4.20 and 3.35 kW/m clockwise. The temporal mean absorbed powers are

2.69. and 2.55 kW/m for SM and DM, respectively, whose

differences are about 5%. Figure 12 shows the orthogonal quality index around the flap by using the DM method before and after the mesh deformation at the maximum angle of rotation counterclockwise. It can be observed that the initial non-structured mesh has a very good index. When the flap is positioned at the maximum angle of rotation, the orthogonal quality remains good. However,

Figure 12. Orthogonal quality index in (a) the initial condition and (b) in the maximum angle of rotation for the case of the Pendulor device using the DM method. DM: Dynamic mesh.

some distortions are noticed between the tip of the flap and the bottom since the bottom imposes a fixed movement condition on the mesh. It can be emphasized that excessive cell distortions can cause a loss of accuracy. The computational time demanded by the DM method is 1.29 times the one of the SM method, which confirms the same tendency of the buoy case. This point is a very important criterion for deciding the best moving method to be used. In this case, as in the buoy one, the cell distortions were not large enough to use the remeshing method. If remeshing method was used, probably, the computational time of the DM method would be higher.