Simulation Flow Conditions of P5168 Propeller
It is well known that the compressible flow equations face difficulties at low Mach numbers due to the large ratio of the acoustic and convective time scales, which leads to an ill-conditioned system when solving low-speed or incompressible flows. The time-dependent system of the Euler and Navier-Stokes equations exhibits stiffness that is strongly dependent on the Mach number and the Reynolds number. In this regard, Briley et al. (Briley et al., 1983) first introduced the preconditioning method using a simple constant preconditioning matrix added to a non-dimensional form of the isoenergetic equations. This generally improved convergence for a test case with reference Mach number Mr = 0.05 using the ADI factorization scheme in primitive variables. However, when applying this preconditioning to rotating flows in either fixed or rotating frame formulation, Sheng, et al. (Sheng & Wang, 2006; Wang & Sheng, 2005) observed the instability of the scheme due to the large variation of rotating speeds across the computational domain. Furthermore, it was found that the convergence of the preconditioned equations is sensitive to the selection of the reference Mach number, especially in rotating flows with a wide range of radial speeds and physical time scales. It was later proved using the Fourier footprint analysis (Wang & Sheng, 2005) that the eigensystem of the compressible governing equations can be significantly affected by both free stream and rotating speeds in rotating flows. A modified preconditioning scheme was thus proposed (Sheng & Wang, 2006; Wang & Sheng, 2005), in which both the global reference Mach number and the rotating Mach number are considered in the formulation of the preconditioning matrix. In general, this modified preconditioning scheme has improved the convergence and accuracy of compressible flows in subsonic, transonic and supersonic Mach number regimes. In this study, the modified preconditioning is further investigated and validated for predicting incompressible viscous flows in rotating machinery, such as a marine propeller P5168.
One of the most important characters for marine propellers is the cavitation observed in high speed flows, which is of vital importance because of the damage of metal surfaces and degradation of performance of lifting surfaces. It is also a source of high-frequency noise in connection with acoustic detection of ships and submarines. Cavitation would take place when the local pressure drops to the vapor pressure. Therefore, accurate prediction of the velocity and pressure field is essential for understanding the process of cavitation inception and improving the hydrodynamic performance of the marine propeller. Since the condition for cavitation inception is related to the tip-vortex location, strength, convection, cavitation position on the propeller surface, and the minimum pressure location. The present study will investigate these important flow characteristics using modified preconditioning scheme, and compare the computed hydrodynamics characters of P5168 with the LDV measurement (Chesnakas & Jessup, 1998).
2. Preconditioned method
2.1. Governing equations
The governing equations for solving the hydrodynamic flow in the present study are based on an unsteady three-dimensional Reynolds-averaged Navier-Stokes equations, coupled with the one-equation turbulent model of Spalart and Allmaras (Spalart & Allmaras, 1991). The unsteady compressible equations are cast in a rotating Cartesian coordinates with a rotational speed of Ω, and are written in terms of absolute velocities. It represents a system of conservation laws for a control volume that relates the rate of change of a vector of node-based average state variables to the flux through the volume surface. The following reference quantities are chosen to normalize the flow variables and equations: density, ρ r ; velocity, U rr ; temperature, T r ; pressure, ρ r U r 2 ; length, L r ; time, L r /U r ; energy and enthalpy, C p T r . The non-dimensional conservative form of the governing equations can be written in integral form as
The above equation can also be written in terms of primitive variables by introducing a transformation matrix
The projected normal flux vector to the face of a control volume is given as
where E c = (γ−1)M r 2 is an Eckert number. Θ is the relative velocity in the normal direction of the control surface in a rotating frame.Let Ω×r be the speed of the control volume at position r, Θ is given as
The viscous stresses in the flux vector are given as
and the heat flux term is
where μ is the viscosity, Pr is the Prandt number, and Re is the Reynolds number based on the reference values.
Since the governing equations are cast in the rotating coordinate system, body forces (Coriolis and centrifugal forces) are added in the momentum equations as a source term due to the non-inertial reference frame, which can be expressed as
The formulation of the governing equations in a rotating frame provides advantages for faster convergence in unsteady flows and a more stable solution without the unsteady effects of the grid motion. Furthermore, the above formulation can be easily restored into the ones in the fixed frame by setting the rotational speed Ω to zero in Eq. (6), making the body force source term disappear from the governing equations. All other relations, including the equation of state, remain the same under the current formulation. Using the absolute velocity components to represent the formulation in rotating frame offers great flexibility for solving the governing equations in both unsteady mode (fixed frame) and steady mode (rotating frame).
2.2. Preconditioned equations
In order to solve flows from essentially incompressible regime to supersonic Mach number, a preconditioning matrix Γq -1 is introduced into the time derivative terms of the compressible governing equations, which can be expressed as
The selection of preconditioning parameter β(M r ) will be discussed in the following section.
The finite–volume scheme needs the approximation of both inviscid and viscous fluxes on the face of a control volume. In order to develop characteristic–based flux formulation, the eigensystem of the preconditioned equations needs to be derived.
The preconditioned Eq. (7) can also be expressed as
The eigenvalues of
2.3. Preconditioning parameter
From the above equations, it is seen that the eigenvalues and eginevectors of the system matrix have been modified due to the introduction of a preconditioned matrix. If the preconditioning parameter β is selected as one, no preconditioning is applied to the governing equations, and thus the original set of eigenvalues for compressible flowsis restored.
Briley, et al. (Briley et al., 2003) suggested that a global constant preconditioning parameter is evaluated based on the reference Mach number as:
where the reference Mach number should represent the global flow property for the problem being solved. For a steady and non–rotating flow, the above choice provided well–behaved eigenvalues in all characteristic waves, and thus improved the convergence and accuracy for solutions up to nearly incompressible Mach number (Briley et al., 2003). However, it was found that such a preconditioning may cause instability for low Mach number rotating flows, such as the case of a marine propeller P5168 presented here. In rotating machinery, the flow filed is complicated with secondary swirling flow structures due to the high-speed rotating elements. It is thus rational to consider the speed of rotating element in determining the preconditioning parameter. In this work the modified preconditioning parameter is employed for flows aimed at unsteady flows in rotating machinery:
where ΜΩ is evaluated based on the rotational reference speed, and α is a scalar factor to tune the rotational Mach number in the preconditioned scheme. Numerical tests conducted in the previous work (Sheng & Wang, 2006; Wang & Sheng, 2005) indicated that the value of α should be in the range of 0.5~1.0 in order to achieve the optimal convergence and accuracy in preconditioned solutions. In the following, an analysis of the Fourier footprint for the eigensystem of the two-dimensional Euler equations provides evidences for the modification of the preconditioning parameter in low Mach number rotating flows.
2.4. Fourier footprint analysis
Fourier analysis is a useful tool for inspecting analytical and numerical solutions of the partial differential equations of mathematical physics. It is used to investigate the accuracy related to the approximation of hyperbolic equations and characteristics of all frequency modes of the discretized equation system (Vichnevetsky & Bowles, 2005). Euler equation is a pure hyperbolic system in which Fourier analysis finds its full power. The Euler Fourier footprint reveals the characters of error propagation and damping of all modes produced by the numerical discretization. The Fourier footprint of a given spatial discretization scheme is obtained by inserting the discrete Fourier transforms in the semi-discrete equation of the finite-volume governing equations. Detailed analytical formulation of the eigenvalues of Fourier footprint of the preconditioned Euler equation can be found in the reference (Wang, 2005).
In terms of Von Neumann stability analysis (O’Brien et al., 1951), the convective part of numerical solutions of hyperbolic equations corresponds to the extent of the Fourier footprints along the imaginary axis, and the dissipation component corresponds to the extent of the negative real axis. In Fig. 1 (a) for a low Mach number non-rotating flow at Mr = 0.01, there exists large disparity among phase velocities of convective waves (in red and black colors) and acoustic waves (in blue and green colors), which is visually indicated by the extension of the imaginary part of Fourier footprints. The group velocity, which describes the propagation of transient errors, is approximately equal to the amplitude of the sinusoid that defines the imaginary component in the plots. Thus, the damping and propagation of transient errors are much less effective for the convective waves, especially for low frequency modes. For this non-rotating flow, the preconditioned system with β = Μ r 2 equalized the real parts of the complex eigenvalues among all wave families, and rescaled the real and imaginary components to the same order of magnitude as indicated in Fig. 1 (b).
Figure 2 presents the Euler Fourier footprint for the same low Mach number flow, but with a rotating Mach number Μ Ω = 0.1. Without preconditioning, large disparity still remains among the convective and acoustic wave speeds in the same order of non-rotating flow, see Fig.2 (a). When the original preconditioning parameter is used ( β = Μ r 2 ) for Mach number relative to a rotating frame of Μ Ω = 0.1, as shown in Fig. 2(b), the real part extensions of Fourier footprints in all wave families are properly rescaled, and the damping effects of the high frequency error modes are well balanced between acoustic and convective waves. However, there are more footprints collapsed to the origin so that neither damping nor propagation of those low frequency modes is possible, and the system will not converge.
This may explain the reason why the original preconditioning showed instability for low Mach number rotating flows. Using a modified a modified preconditioning parameter as β Μ r 2 +Μ Ω 2 proposed here, the extensions along the negative real axis for all wave modes are clustered to the same order of magnitude, as indicated in Fig. 2(c). In addition, the imaginary parts of the complex eigenvalues are mostly away from the negative real axis. The propagation properties are thus maintained at the ordinary levels so that the low frequency error modes can rapidly propagate out of the computational domain. The convergence of the Euler system is thus improved.
3. Discretization scheme
The preconditioned governing equations are discretized using a finite volume formulation. This is done by dividing the computational domain into a finite number of elements, from which control volumes are formed that surround each vertex in the mesh. The flow variables are stored at the vertices of the element. The preconditioned governing Eq. (7) or (9) is then numerically integrated over the closed boundaries of control volumes surrounding each node, and the preconditioning matrix
3.1. Convective flux evaluation
The convective numerical flux of preconditioned governing equations is calculated using in a formula similar to Roe’s flux, as proposed by Briley, et al. (Briley et al., 2003). Consider the surface j of a control volume whose left and right states are denoted by L and R, the numerical flux projected onto the surface j can be written as:
Since the eigensystem and the flux are evaluated using the above averaged variables instead of Roe variables, the current approach is considered to be an extension of the Roe flux approximation.
In the above Eq. (16) quantities
where r is the vector that extends from the central node to the midpoint of each edge, and
3.2. Diffusive flux evaluation
Since the local state can be found from the reconstructed solutions at the vertex of each edge, the essence of computing the viscous fluxes is finding the gradients of the velocity vector and the temperature at the midpoints of the interface. A directional derivative approach for viscous gradients has been developed by Hyams (Hyams, 2002). In the present work, a normal derivatibbe approach is presented.
Introducing two unit tangential vectors
It is easy to verify that unit vectors are orthogonal to each other. The components of any gradient can be expressed in terms of derivatives in
The derivatives in the normal direction
3.3. Temporal discretization
The temporal discretization is considered with an implicit Euler method for the system of governing equations. A general expression is available for the temporal discretization, which is given as:
3.4. Newton’s method
The above nonlinear system of Eq. (26) is solved by Newton’s method, which yields a linear system of equations at each time step. The use of Newton’s method allows relatively large time step to be used in numerical computations, a considerable savings compared to the explicit method, in particular fro the unsteady time-accurate problems. Denote
The Newton’s method solves the following equation:
where m = 1,2,3... are the Newton steps, with
In the above Eq. (29), the first term on the right hand side is the contribution from the unsteady time derivative of q, the second term is the contribution from the steady state residual (both inviscid and viscous) of the governing equations, and the third comes from the source term in the rotating reference frame. The flux Jacobian of the steady state residual can be evaluated by an approximate method. Consider the matrix
A similar method is used to evaluate the viscous flux Jacobian by taking the derivatives of discretized viscous flux formulation.
3.5. Gauss-Seidel relaxation
The nonlinear system of equations is linearized by Newton’s method, which results in a sparse system of equations at each time. The solution of the sparse system of equations is obtained by a relaxation scheme in which
In the forward pass,
4. Turbulence model
Several turbulence closure models have been developed for high-Reynolds number flows. One of them is the Spalart-Allmaras one-equation turbulence model (Spalart & Allmaras, 1991), which formulates a transport equation for the turbulence Reynolds number, and is the only one used in the present study. From the original Spalart-Allmaras formulation, a transport equation for the turbulence Reynolds number
where the first and second terms on the right hand side of Eq. (34) is the source production, and the third term is the turbulence dissipation. Constants in Eq. (34) can be found in the original paper (Spalart & Allmaras, 1991).
5. Characteristic-based boundary conditions
In order to compute internal and external flows at arbitrary Mach numbers, four commonly used characteristic–based boundary conditions have been developed (Sheng & Wang, 2003), including far field boundary, inflow boundary based on total conditions or mass flow, outflow boundary based on back pressure, and impermeable wall. The characteristic variables are derived based on preconditioned system matrix in order to be consistent with interior point treatment. To solve the governing equations in a rotating frame, a source term is included which should be implicitly integrated into the characteristic boundary conditions. In the following, characteristic-based boundary conditions are provided for completeness.
5.1. Far field
The preconditioned system matrix can be diagonalized as
where subscript b denotes a point on the boundary, and subscript r denotes an appropriate reference point depending on the sign of the eigenvalues on the boundary. Five characteristic variables for the three-dimensional preconditioned system are
The direction of propagation of five characteristic waves is determined by the sign of each associated eigenvalue. Since all boundary faces have an positive outward pointing normal with the current unstructured grid topology, a negative sign of the eigenvalue means that the related characteristic wave propagates from the free stream to the boundary nodes, and the reference value is the free stream value; while a positive sign of the eigenvalue means that the characteristic wave propagates from the interior of the computational domain to the boundary, so the reference value should be taken from the interior node. The primitive variables at the boundary are then computed by solving the following equations:
For internal flows in rotating machineries, total conditions (total pressure Pt, total temperature Tt) are usually given at the inlet, along with the flow angle
Since there are six unknowns (ρ b , p b , u b , v b , w b , U b ) out of the above five equations, one more relation is required to solve the equations. For the subsonic inflow, all eigenvalues except
An alternate way to specify the inflow boundary conditions is to use the total temperature Tt, mass flow rate
All flow variables at the boundary are obtained by solving the above equations.
Characteristic variable boundary conditions are applied to outflows as well. For the subsonic outflow, the first four eigenvalues λ(1,2,3,4) are positive and the fifth eigenvalue λ(5) is negative. Flow variables at the boundary are connected with the interior values through the first four characteristic relations. The fifth relation is replaced by the static pressure at the exit
The boundary values (ρ b , p b , u b , v b , w b ) are solved by the above equations.
For a supersonic outflow case, the fifth eigenvalue λ(5) is positive, and no characteristic waves can propagate from downstream to the upstream. The last equation (57) is replaced by the fifth characteristic relation (58) where characteristic variables come from the interior field:
5.4 Impermeable wall
For impermeable surfaces, the specified condition is that of no flow past perpendicular to the surface, i.e.,
It is noticed that the first four characteristic relations (59) through (62) can be obtained by setting
To validate the new preconditioned algorithm for low Mach number flows in rotating machinery, application of a hydrodynamic flow about the David Taylor marine propeller P5168 is performed in this study. This is a five-bladed, controllable–pitch propeller (Fig. 3) with a design advance coefficient of J = 1.27. The diameter of the propeller is 0.4027m. The 3-D LDV tip-vortex flow field measurement was conducted by Chesnakas and Jessup (Chesnakas & Jessup, 1998) in the CDNSWC 36-inch water tunnel. Four different advance coefficients were measured with detailed velocity data behind the propeller. The numerical simulations were conducted under the four advance ratios corresponding to the experiment. The body-fitted mixed-element anisotropic unstructured mesh was generated using the Advancing-Front Local Reconnection (AFLR) grid generation technique (Marcum, 2001). The far field boundary of the computational domain extends to 3.5 blade diameters in both upstream and downstream, and 1.5 blade diameters in the radial direction. The unstructured mesh contains 3.825 million nodes and 12.182 million elements. For the viscous gird used in this work, the y+ distribution is less than 1.5 over the solid bodies to provide a good viscous sublayer resolution. In addition to the fine grid region over all the viscous surfaces for achieving good boundary layer resolution, a number of interior faces were created downstream in the domain along blade wake paths to refine the volume grid in order to capture the fine details of the flow field in the downstream of the propeller. Those interior faces are transparent to the flow solver with no physic boundary conditions required. As moving further beyond the shaft rear face, the computational mesh becomes coarse to reduce the computational costs. The spatial variation of grid resolution is expected to affect the accuracy of the numerical simulations and its influence will be assessed in the investigation of the simulation results. Computation domain was partitioned into 64 sub-domains for parallel computations.
The velocity measurements were performed by Chesnakas and Jessup (Chesnakas & Jessup, 1998) in the water tunnel, and is used in this work for code validation. Detailed flow conditions for the simulation cases are summarized in Table 1. The characteristic based far-field boundary condition is applied to the far-field boundary with implicit treatment of the source terms. The no-slip boundary conditions are applied to the propeller blade and shaft
|Advance Coefficient J=V/nD||Free Stream V(m/s)||Rotational Speed RPM||Tip Velocity V tip (m/s)||Tip Mach Number M tip||Reynolds Number Re (106)|
surfaces in the rotating frame. All solid surfaces are assumed to be adiabatic, i.e., no heat fluxes across the solid surfaces. The Reynolds number is based on the free stream speed and the diameter of the propeller. The reference Mach number Μ r = 0.1 is selected for all the cases based on the free stream condition. The preconditioning parameter is estimated by the modified scheme as β = Μ r 2 +Μ Ω 2 , in which Μ Ω is set equal to the propeller tip Mach number listed in Table 1 for each case. This seems to provide the best accuracy and convergence. Due to the adoption of the relative rotating frame, flows of the rotating propeller can be considered as steady problems. The local time stepping is used in all the simulations with CFL number of 10. One Newton’s iteration and eight Gauss-Seidel sub-iterations are used at each time step. Simulations are initiated from a uniform flow, and converged steady solutions are achieved in less than 1000 iterations. Fig. 4 shows the simulation history of the computed propeller thrust coefficient to demonstrate the convergence at each advance ratio of the propeller.
6.1. Propeller performance
The propeller has been open-water tested in the CDNSWC towing basin, where the cavitation inception conditions are related to either thrust loading (KT) or tunnel velocity (J). The numerical simulations are performed under the same four advance ratios J=0.98, 1.1, 1.27, 1.51 that have been tested in the water tunnel. The computed thrust and torque on the propeller are calculated by integrating the pressure over the propeller surfaces. The thrust and torque coefficients are defined as
where ρ is the flow density, n is the rate of revolutions in the unit of revolution per second, and D is the propeller diameter. Fig. 5 shows that computed thrust and torque coefficients for the P5168 propeller agree well with experiment in all advance ratios. Results obtained by the current preconditioning method seem to provide excellent predictions of the propeller thrust and torque comparing with the experimental data.
6.2. Circumferentially averaged velocity
The measured and computed circumferentially averaged velocities at two downstream locations (x/R=0.2386 and 0.8378) are shown in Fig. 6. The filled symbols are measured values and unfilled symbols are computed values. Predicted axial and tangential velocities match the experiment quietly well, although a slight underestimated radial velocity is seen in Fig. 6 (c). This is related to the prediction of tip-vortex flow in the downstream of the
propeller, which will be discussed in the following section. However, the overall trend in velocity distributions are well captured in the computation, including in the blade tip-vortex region at r/R=1.0. The computed results showed favorable agreement with the experimental data.
6.3. Mean velocity contours
To further investigate the tip-vortex flow generated by the propeller, the computed contours of the axial, radial, and tangential velocity components in the rotating frame at a stream-wise location of x/R = 0.2386 are presented in Fig. 7. The tip vortices generated by the blade tip are clearly visible in both the computation and the experiment, where the vortex core is represented as blue “spokes” in the plot of low velocity. There is a flow circulation around the tip vortex, as well as the radial velocity generated by the blade wake. The strong asymmetry of the vortex can be seen here. It also indicates that the current preconditioned scheme does capture the general feature of the tip-vortex flow. The axial and radial velocity components (Fig. 7(a) and (b)) showed a good agreement regarding the locations of tip vortices between the computation and measurement. However, the numerical diffusion was observed in the current second-order discretization scheme in the vortex core region, as predicted vortex strength appeared not as strong as what was measured in the experiment. To improve the prediction accuracy in the vortex core region, an adaptively refined mesh or a higher-order spatial discretization scheme (Sheng & Wang, 2009) is necessary. Nevertheless, the modified preconditioning method was proved to provide stable and reasonably accurate predictions to the hydrodynamic viscous flows in rotating machinery, which would otherwise impossible to achieve with the original preconditioning method, not to mention the compressible flow solver.
6.4. Tip-vortex trajectory
The tip vortex trajectory was predicted and compared with the experiment measurements. An approximate approach to locate the tip vortex core location is to limit the search for a minimum pressure within the region of the known vortex indicated by a swirl parameter. The physical interpretation for this criterion is given by Jeong and Hussain (Jeong & Hussain, 1995). Within the vortex core, pressure tends to have a local minimum on the axis of a circulatory or swirling motion when the centrifugal force is balanced by the pressure force, which is strictly true only in a steady inviscid planar flow. Using this method, the tip-vortex trajectories were identified from computed flows at the advance ratio of J=1.1. Fig. 8 compared the predicated radial locations of the tip vortex with the experimental observation. In the experiment, the tip vortex was tracked from the blade tip to an axial location of approximately one diameter downstream where x/R=2.0. The current simulation conducted in this study can predict the tip vortex radial trajectory accurately up to a limited distance downstream, due to the grid coarsening in the downstream of the propeller. As shown in Fig. 8, the predicted tip-vortex trajectory is tracked to the location about only one radius downstream at x/R=1.0 before dissipated into the flow field. The largest disagreement between the prediction and the experiment measurement occurs at the last station of x/R=1.0. Beyond that point it appeared to be difficult to trace the tip vortex as it has completely smeared out in the flow field. Examination of the grid resolution of the propeller indicated a finer grid resolution only coving a limited distance in the downstream of the propeller. The grid cells became very coarser as moving beyond the shaft rear surface at station x/R=0.69. Fig. 9 clearly demonstrated the tip vortex convection which is quickly dissipated due to the limited grid resolution. As discussed in the previous section, a refined computational mesh and a higher-order discretization scheme may reduce the numerical dissipation of this vortical flow, although the computational costs would be significantly higher.
6.5. Tip vortex convection and decay
Another way to trace the tip vortex of the propeller is to use a non-dimensional parameter termed intrinsic swirl parameter
6.6. Tip-vortex attachment
The water tunnel experiment also observed suction-side tip vortex attachment occurred aft of the blade tip over certain range of tested advance ratio. The locations of the tip vortex attachment on the five blades showed minor variation from blade to blade, and their average radial location has been estimated. In the open-water experiment at an advance ratio of J=0.98, this average location was estimated at the blade trailing edge of 0.99R radius, while at J=1.1, the attachment point was moved aft to 0.998R radius of the trailing edge. From the current computed results, the iso-surface of swirl parameter τ=1.0 is generated to visualize the tip vortex structure as shown in Fig. 12. The iso-surface is shaded by static
pressure. The propeller geometry is not plotted in the figure, but the tip region of the propeller blade is seen clearly, which is actually the swirl iso-surface in the blade boundary layer. At the blade tip, the attachment occurs aft a local minimum pressure center at the trailing edge. The iso-surface of the swirl parameter indicated that the tip vortex attachment occurs at the trail edge on the blade suction side between 0.994R-0.999R. The variation of the attachment points at J=0.98 and J=1.1 are not distinguishable except that the local minimum pressure center upstream the attachment location has lower pressure at J=0.98 and originates stronger tip vortex. The current simulations correctly captured the tip-vortex attachment observed in the experiment, which is important for the improved understanding of the cavitation inception of the marine propeller.
A modified preconditioning method was investigated and validated in the prediction of hydrodynamic viscous flows for marine propeller P5168. The preconditioning parameter is based on the reference Mach number and rotating Mach number to provide stable and accurate solutions for low Mach number flows in rotating machinery, while the original preconditioning method failed to provide converged solutions in the case presented in this study. The predicted overall propeller performance, circumferentially averaged velocities, mean velocity contours, tip-vortex trajectory and decay at blade downstream stations are validated against the experimental data. The comparison between the computation and the experiment indicates that the current preconditioned solver was able to capture the general features of the tip vortex generated by the propeller. In particular, the predicted thrust and torque coefficients at several advance ratios matched well with the open-water experimental data. However, the current numerical simulation showed a quick decay of the tip vortex in the propeller downstream, due to a lack of enough grid resolution. Future studies include using adaptive grid technique to locally refine the computational mesh in the vortex core region, and a higher-order spatial discretization schemes to further reduce the numerical diffusion in the solver.
The author wishes to thank Qiuying Zhao of The University of Toledo and Xiao Wang of High Performance Computing Collaboratory at Mississippi State University for helping processing the simulation data and formatting the manuscript.