An effective preconditioning strategy for volume penalized incompressible/low Mach multiphase flow solvers

The volume penalization (VP) or the Brinkman penalization (BP) method is a diffuse interface method for simulating multiphase fluid-structure interaction (FSI) problems in ocean engineering and/or phase change problems in thermal sciences and engineering. The method relies on a penalty factor (which is inversely related to body's permeability K) that must be large to enforce rigid body velocity in the solid domain. When the penalty factor is large, the discrete system of equations becomes stiff and difficult to solve numerically. In this paper, we propose a projection method-based preconditioning strategy for solving volume penalized (VP) incompressible and low-Mach Navier-Stokes equations. The projection preconditioner enables the monolithic solution of the coupled velocity -pressure system in both single phase (uniform density and viscosity) and multiphase (variable density and viscosity) flow settings. In this approach, the penalty force is treated implicitly, which is allowed to take arbitrary large values without affecting the solver's convergence rate or causing numerical stiffness. It is made possible by including the penalty term in the pressure Poisson equation (PPE), which was not included in previous works that solved VP incompressible Navier-Stokes equations using the projection method. We show how and where the Brinkman penalty term enters the PPE by re-deriving the projection algorithm for the VP method. Solver scalability under grid refinement is demonstrated, i.e., convergence is achieved with the same number of iterations regardless of the problem size. A manufactured solution in a single phase setting is used to determine the spatial accuracy of the penalized solution. Various values of body's permeability K are considered. Second-order pointwise accuracy is achieved for both velocity and pressure solutions for reasonably small values of K. Error saturation occurs when K is extremely small, but the convergence rate of the solver does not degrade. The solver converges faster as K decreases, contrary to prior experience. Two multiphase fluid-structure interaction (FSI) problems from the ocean engineering literature are also simulated to evaluate the solver's robustness and performance (in terms of its number of iterations). The proposed solver also allows us to investigate the effect of K on the motion of the contact line over the surface of the immersed body. It also allows us to investigate the dynamics of the free surface of a solidifying metal. & COPY; 2023 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons .org /licenses /by /4 .0/).