Volume Penalization Research Articles

Direct aeroacoutic simulation of a ducted axial fan with acoustic liners using a volume penalization lattice Boltzmann solver on heterogeneous massively parallel platforms

Aeroacoustic noise of an axial fan is directly solved to understand its mechanism, using our in-house Lattice Boltzmann Method (LBM) solver. To make the LBM solver industrially successful, we must reduce calculation time while maintaining its accuracy. We adopt Volume Penalization method for moving boundaries to achieve independent parallelism, and to utilize GPU (Graphical Processing Units) for parallel computation. In this study, we present a comparison between numerical and experimental results of an axial fan in a lined duct. The acoustic liner is designed to reduce 2nd BPF noise, and the noise reduction of 9.8dB is observed in the experiment. With our LBM code, the noise reduction was reproduced within an accuracy of 2.5dB, and profound insights are provided through the detailed investigation of the numerical results.

Adjoint-based shape optimization for compressible flow based on volume penalization method

AbstractReducing the resistance of compressible flow around a blunt body is of great interest in engineering applications, while an efficient shape optimization method for compressible flows remains far from well established, especially for high Mach numbers. To this end, a volume penalization method for simulating compressible flows past a no-slip and isothermal solid is established by introducing an artificial body force and a heat sink into the governing equations. The level-set functions are introduced as design variables, and the cost functional is defined as the total drag acting on the solid. Then, a continuous adjoint-based shape optimization algorithm for drag reduction is developed by deriving the adjoint equations, the adjoint boundary conditions, and the shape update formula. Both the forward and adjoint simulations are verified by existing studies. The results show that the relative deviations of the drag coefficients obtained in the present study from those reported in the reference studies are around 5% at most, and also a comparable drag reduction rate and also optimal shapes can be reproduced by the present optimization scheme for benchmark problems at relatively low Mach numbers considered in previous studies. Finally, the present method is applied to shape optimization of an initially two-dimensional cylinder and also a three-dimensional sphere in the transonic regime of Ma∞ = 1.2. The drag reduction of over 20% is achieved for both two-dimensional and three-dimensional cases.

Theoretical considerations of the volume penalization immersed boundary method for turbulent flows

This Letter explores the volume penalization immersed boundary method for turbulent flows from a more physical perspective. The volume penalization approach consists of introducing a penalty source into the governing equations, resulting in a flow akin to a porous medium with low permeability. Although penalizing the turbulent equations conventionally involves adding a similar penalty source as in the original equations, this work reveals an alternative formulation that includes an additional term with physical meaning. The novelty of this Letter is to consider the penalized flow with an additional property, the fluid resistance, establishing a cross correlation with fluctuating velocity for further modeling.

The effect of shear-induced segregation on bidisperse suspension settling

Shear is a common phenomenon in practical mineral processes, particularly in pipeline transportation. When applied to realistic polydisperse systems consisting of particles of different sizes, shear can lead to distinct interactions of each size component with the carrier fluid, resulting in segregation between sizes. This paper focuses on the effect of size segregation induced by shear on the settling process of suspension flows containing different-sized spherical particles. In this regard, the bidisperse system is considered as the simplest example of polydisperse flows, ensuring the investigation retains the essential features of the system of interest without added complexity. The Volume Penalization Immersed Boundary Method, as a resolved Eulerian–Lagrangian CFD-DEM approach, is employed to gain detailed information on individual particles. A simple cross shear is generated by the moving side walls, and its effect is initially evaluated on two simulation cases containing single spherical particles of different sizes. The results from the single particle simulation cases reveal a higher acceleration of the larger particle towards the centreline due to the higher lateral force proportional to its size. The resolved approach is then used for the bidisperse multiparticle systems with diameter ratios of l=2 and l=5, both subjected to a similar sheared field. In the suspension with l=2, the scalar dispersion function, which measures the degree of migration of each size component between the side walls, indicates segregation by size as the shear rate increases. This is evident in the motion of large particles towards the centre and the accumulation of smaller particles near the walls, which arises from higher proportional force on the larger particles observed in the single particle simulation cases. Although the increase in the diameter ratio to l=5 does not significantly change the migration behaviour of large particles, it provides a more uniform distribution of smaller ones. In this case, the small particles find it more easily to fit into the gaps between the large particles for l=5. Overall, this accumulation of each size component at specific regions in the sheared domain affects the local interaction of settling particles with the fluid flow, resulting in an increase in the settling rate of large particles up to 40% in the considered bidisperse cases. This highlights the impact of shear on the interaction of spherical particles in a multiple sized system, specifically in the practical sheared applications including mineral slurries transportation.

Flat flow solution to the mean curvature flow with volume constraint

Abstract In this paper I will revisit the construction of a global weak solution to the volume preserving mean curvature flow via discrete minimizing movement scheme by Mugnai, Seis and Spadaro [L. Mugnai, C. Seis and E. Spadaro, Global solutions to the volume-preserving mean-curvature flow, Calc. Var. Partial Differential Equations 55 2016, 1, Article ID 18]. This method is based on the gradient flow approach due to Almgren, Taylor and Wang [F. Almgren, J. E. Taylor and L. Wang, Curvature-driven flows: a variational approach, SIAM J. Control Optim. 31 1993, 2, 387–438] and Luckhaus and Sturzenhecker [S. Luckhaus and T. Sturzenhecker, Implicit time discretization for the mean curvature flow equation, Calc. Var. Partial Differential Equations 3 1995, 2, 253–271] and my aim is to replace the volume penalization with the volume constraint directly in the discrete scheme, which from practical point of view is perhaps more natural. A technical novelty is the proof of the density estimate which is based on second variation argument.

Adjoint-based shape optimization for radiative transfer in porous structure for volumetric solar receiver

Due to their large surface areas, porous structures are often used for volumetric solar receivers for efficient utilization of the solar energy. However, the precise simulation and optimization of radiative heat transfer inside porous structures are still challenging due to their complex geometries. In the present study, an adjoint-based shape optimization framework for radiative transfer in porous structures is first proposed. Specifically, an arbitrary complex geometry of a porous structure is represented by a level-set function embedded in a Cartesian grid system. The volume penalization method combined with the discrete ordinate method is employed to simulate a forward radiative transfer problem in a porous structure. Then, a cost functional for an adjoint analysis is defined by the volume average of the divergence of the radiative heat flux, which corresponds to the total heat absorbed by the porous structure with a negative sign. The adjoint equations for minimizing the cost functional, and the corresponding shape update formula are derived and implemented in an open-source software, OpenFOAM. The developed adjoint-based shape optimization framework is applied to three types of initial porous structures, i.e., SC Kelvin, BCC Kelvin, and square honeycomb models. For all the initial porous models, the front surface commonly extends towards the direction of the incident radiation, and consequently notable improvements of the absorption efficiency can be achieved through the optimization. In addition, a systematic parametric survey for the initial shape of the SC Kelvin model is conducted to reveal the impacts of the number of unit cells, the solid emissivity, and the solid temperature on the optimal shapes and the resulting absorption efficiencies.

Numerical Simulation of Cavitation Control around a Circular Cylinder Using Porous Surface by Volume Penalized Method

In this work, we conducted a numerical study on the cavitation flow around a circular cylinder with Re=200 and σ=1, through the implementation of a porous coating. The primary objective addressed the effectiveness of utilizing a porous surface to control cavitation. We analyzed the cavitation dynamics around the cylinder and the hydrodynamic performance at different permeability levels of the porous surfaces (K=10−12−10−10). The flow was governed by the density-based homogeneous mixture model, and the volume penalization method was used to deal with the porous layer. A high-order compact numerical method was adopted for the simulation of the cavitating flow through solving the preconditioned multiphase equations. The hydrodynamic findings demonstrated that the fluctuations in the lift coefficient decreased when the porous layer was applied. However, it is not possible to precisely express an opinion about drag because the drag coefficient may vary, either increasing or decreasing, depending on the permeability within a constant thickness of the porous layer. The results revealed that the application of a porous layer led to the effective suppression of cavitation vortex shedding. In addition, a reduction of the shedding frequency was obtained, which was accompanied by thinner and elongated vortices in the wake region of the cylinder. With the proper porous layer, the inception of the cavitation on the cylinder was suppressed, and the amplitude of pressure pulsations due to the cavitation shedding mechanism was mitigated.

A comparison of force computations for modeling fluid - structure interaction problems with the lattice Boltzmann volume penalisation method

A comparison of force computations for modeling fluid - structure interaction problems with the lattice Boltzmann volume penalisation method

Implementation of a level-set-based volume penalization method for solving fluid flows around bluff bodies

A volume penalization-based immersed boundary technique is developed and thoroughly validated for fluid flow problems, specifically flow over bluff bodies. The proposed algorithm has been implemented in an open source field operation and manipulation (OpenFOAM), a computational fluid dynamics solver. The immersed boundary method offers the advantage of inserting a complex solid object inside a Cartesian grid system, and therefore, the governing equations can be applied to such a simpler grid arrangement. For capturing the fluid–solid interface more accurately, the grid is refined near the solid surface using topoSetDict and refineMeshDict utilities in OpenFOAM. In order to avoid any numerical oscillation and to compute the gradients accurately near the interface, the present volume penalization method (VPM) is integrated with a signed distance function, which is also referred to as a level-set function. Benchmark problems, such as flows around a cylinder and a sphere, are considered and thoroughly validated with the results available in the literature. For the flow over a stationary cylinder, the Reynolds number is varied so that it covers from a steady two-dimensional flow to an unsteady three-dimensional flow. The capability of the present solver has been further verified by considering the flow past a vibrating cylinder in the cross-stream direction. In addition, a flow over a sphere, which is inherently three-dimensional due to its geometrical shape, is validated in both steady and unsteady regimes. The results obtained by the present VPM show good agreement with those obtained by a body-fitted grid using the same numerical scheme as that of the VPM, and also with those reported in the literature. The present results indicate that the VPM-based immersed boundary technique can be widely applicable to scientific and engineering problems involving flow past stationary and moving bluff bodies of arbitrary geometry.

Multi-objective topology optimization of heat transfer surface using level-set method and adaptive mesh refinement in OpenFOAM

The present study proposes a new efficient and robust algorithm for multi-objectives topology optimization of heat transfer surfaces to achieve heat transfer enhancement with a less pressure drop penalty based on a continuous adjoint approach. It is achieved with a customized OpenFOAM solver, which is based on a volume penalization method for solving a steady and laminar flow around iso-thermal solid objects with arbitrary geometries. The fluid-solid interface is captured by a level-set function combined with a newly proposed robust reinitialization scheme ensuring that the interface diffusion is always kept within a single local grid spacing. Adaptive mesh refinement is applied in near-wall regions automatically detected by the level-set function to keep high resolution locally, thereby reduces the overall computational cost for the forward and adjoint analyses. The developed solver is first validated in a drag reduction problem of a flow around a two-dimensional cylinder at the Reynolds numbers of 10 and 40 by comparing reference data. Then, the proposed scheme is extended to heat transfer problems in a two-dimensional flow at the Prandtl number of 0.7 and 6.9. Finally, three-dimensional topology optimization for multi-objective problems is considered for cost functionals with different weights for the total drag and heat transfer. Among various solutions obtained on the Pareto front, 4.0% of heat transfer enhancement with 12.6% drag reduction is achieved at the Reynolds number of 10 and the Prandtl number of 6.9. Moreover, the optimization of a staggered pin-fin array demonstrates that the optimal shapes and arrangement of the fins strongly depend on the number of rows from the inlet. Specifically, the pin-fins in the first and third rows extended in the upstream direction further enhance heat transfer, while the fins in the second row vanish to reduce pressure loss.

A combined volume penalization/selective frequency damping approach for immersed boundary methods: Application to moving geometries

High-order numerical techniques and immersed boundary methods (IBMs) are gaining popularity to avoid constructing body-fitted meshes while accurately resolving complex flows on Cartesian grids. Recently, we have presented a new treatment for the immersed boundary method based on the combination of volume penalization and selective frequency damping (SFD) [J. Kou and E. Ferrer, “A combined volume penalization/selective frequency damping approach for immersed boundary methods applied to high-order schemes,” J. Comput. Phys. 472, 111678 (2023).], that offers improved accuracy for nonmoving geometries. The objective of SFD is to remove nonphysical, high-frequency oscillations inside the solid body and by doing so enhance the accuracy in the fluid region. The present paper extends the new immersed boundary treatment to moving geometries. The convergence of this approach is first validated by the method of manufactured solutions, where we design a one-dimensional advection–diffusion case, with a moving interface, to validate the numerical accuracy. Second, we simulate an unsteady flow past a plunging circular cylinder (Navier–Stokes solver). In this case, we show that the combination of volume penalization and SFD provides improved accuracy for moving geometries.

Generalisation of the Penalised Wall Function Method for the Simulation of Turbulent flows With Unfavourable Pressure Gradients

The penalized wall function method for simulation of compressible near-wall turbulent flow regions in the numerical modeling of viscous compressible flows is developed. The method is formulated as a differential condition to match the outer and the wall function solutions and is based on a generalized characteristic-based volume penalization method to transfer shear stress from the outer region of the boundary layer to the wall. The method is modified to extend its applicability to turbulent flows with adverse pressure gradient, when separation and reattachment zones are formed, as well as to use computational meshes with coarser near-wall resolution. These advantages are demonstrated for two test problems, namely, the flow over a flat plate with zero and adverse pressure gradients.

Impact of parametric variation to achieve extreme mechanical metamaterials through topology optimization

Metamaterials are manufactured structures having extreme properties which do not exist in nature. Topology optimization (TO) provides an excellent alternative to the intuition-based homogenization approach used in commercial packages for automating the design of metamaterials. However, multiple topological parameters must be decided to get the tailored effective properties. This research aims to study the optimum parameters to be used in TO for developing incompressible material. To achieve this, we design microstructures to maximize the bulk modulus. A number of parametric studies are performed by varying topological control parameters such as volume fraction, penalization power, filter radius, mesh size and initial design. In addition to this, the combinations of these parameters are also investigated. These studies will be helpful for the design of incompressible metamaterials. The bulk modulus for topologically optimized designs approaches the Hashin–Shtrikman upper bound (HSUB). The results show that appropriate parameter selection can significantly improve the microstructure’s incompressibility. Such materials can have potential applications in designing aerospace and energy-absorbing components.

Drag reduction effect of distributed roughness on the transitional flow state using direct numerical simulation

We analyzed the details of Tollmien–Schlichting wave transitions on a sandy rough surface through direct numerical simulations by resolving each roughness level using the volume penalization and zonal methods. Although wall roughness has generally been discussed in terms of increased drag in fully developed turbulent flows, we found that the sandy roughness on a flat plate suppresses the growth of both the turbulent kinetic energy and friction drag coefficient during the laminar-turbulent transitional flow state. The roughness size is on the microscale order, assuming similarity based on the Reynolds number and considering the case of the surface of a transonic airplane. Therefore, we call such three-dimensional and complex-shaped roughness with this effect distributed microroughness (DMR). On the sandy rough surface, the Tollmien–Schlichting wave collapses into a three-dimensional flow structure in the initial stage; however, the turbulent kinetic energy remains small. Generally, the turbulent kinetic energy of the subsequent three-dimensional turbulent component increases in correlation with the size of the primary T-S vortex. On the rough surface, the T-S roll vortex collapses before it maximally grows due to the spanwise disturbance caused by the three-dimensionality of the sandy roughness shape, resulting in the three-dimensional turbulence associated with the roll remaining small. Near the wall surface, the velocity fluctuates due to the distributed roughness, and the statistical peaks of the kinetic energy budget tend to be away from a wall. Although the details of the mechanism by which changes in the energy balance cause a decrease in the productive turbulent energy have not been elucidated, we will clarify the changes due to roughness parameters in a concurrent study.

Assessment of volume penalization and immersed boundary methods for compressible flows with various thermal boundary conditions

A continuous forcing immersed boundary method (IBM) and a volume penalization method are investigated to simulate compressible viscous flows past an isothermal or adiabatic solid obstacle using a body non-conformal Cartesian grid. The present methods are validated on a large range of flow configurations such as incompressible and compressible flows (steady or unsteady) past isothermal or adiabatic solid obstacles. An adiabatic condition with implicit treatment of the source term is derived for the penalization method, taking advantage of the framework of the ghost-cell methods. Accuracy and computational cost of the methods are discussed. Three-dimensional supersonic flow configurations are carried out with the penalization method. The coupling between the porosity parameter and the different boundary conditions of the present penalization method is investigated in order to properly reproduce the wave reflection on a solid wall. Finally, a new penalization model is designed to improve the shock wave reflection on a solid surface while significantly reducing the computational cost compared to the previous porosity model.

Application of immersed boundary methods to non-Newtonian yield-pseudoplastic flows

The interaction force between fluid and particles in particulate flows is highly dependent on the rheology of fluid. In some engineering applications, like mineral processing, the rheological behaviour of the carrier fluid is characterized by a yield-pseudoplastic non-Newtonian model. The immersed boundary method (IBM) is a numerical strategy to simulate the interaction between phases in the particulate systems using a Cartesian fixed grid. It imposes the solid boundary condition in the computational domain through a modification (e.g. a forcing term) in the governing equations. However, there is still little study in the literature discussing the applicability of different IBMs to yield-pseudoplastic flows. In this paper, the numerical methodology in several versions of immersed boundary methods, including the Volume Penalization IBM (VP-IBM), Indirect Imposition of Discrete Forcing IBM (ID-IBM), and Direct Imposition of Discrete Forcing (DD-IBM), is described. The numerical implementation of the methods is first validated using Newtonian benchmark cases. Then the solvers are utilized for the simulation of yield-pseudoplastic cases with stationary or moving particles. The investigation shows that in the explicit forcing methods, like VP-IBM and ID-IBM, the forcing procedure is conducted using an intermediate solution over the entire computational domain. The dependence of viscosity to strain rate in non-Newtonian cases leads to an inevitable sharp change in viscosity near the solid surface, giving rise to stiffness of the intermediate solution in the interface region. This makes the explicit forcing approaches based on intermediate solution incompatible with the flows with non-Newtonian yield-pseudoplastic rheology. On the other hand, the implicit forcing used in DD-IBM directly implements the boundary conditions on the boundary cells, making this method free from the perturbations related to the highly viscous solid region. In this regard, the simulation with the Direct Imposition IBMs like ghost-cell shows a good agreement between the predicted results and the experimental measurements for the stationary/moving cases. Therefore, it is shown that the DD-IBM is a reliable approach to modelling non-Newtonian flows.

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 κ) 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 κ are considered. Second-order pointwise accuracy is achieved for both velocity and pressure solutions for reasonably small values of κ. Error saturation occurs when κ is extremely small, but the convergence rate of the solver does not degrade. The solver converges faster as κ 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 κ 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.

Characteristic-Based Volume Penalization-Imposed Wall Function Method for Turbulent Boundary Layer Modeling

A method to approximate near-wall boundary conditions for the compressible Reynolds-Averaged Navier–Stokes equations is proposed. The differential formulation to match the external and the wall function solutions is reformulated in a form of the generalized characteristic-based volume penalization method to model the transfer of the shear stress from the outer region of the boundary layer to the wall. The exchange location is specified implicitly in terms of a localized source term in the boundary layer equation written as a function of the distance from the wall normalized by the viscous length scale. The shear stress on the wall is determined by solving an auxiliary equation for the wall-stress imposing the analytical wall function solution through the characteristic-based volume penalization method. The proposed method noticeably reduces the near-wall mesh resolution requirements without a significant modification of the numerical algorithm and completely eliminates the ill-defined explicit solution matching procedure. The developed approach is numerically implemented using the vertex-centered control volume method on structured meshes. Its effectiveness is demonstrated by solving two test problems: the two-dimensional channel flow and turbulent flow over an infinitely thin plate.

Adjoint-based shape optimization for radiative transfer using level-set function and volume penalization method

A inverse design of a radiative heat transfer surface is established based on an adjoint analysis combined with a level-set function and a volume penalization method. The level-set function is employed to represent an arbitrary complex geometry embedded in a Cartesian grid system, whereas the volume penalization method is applied to simulate a forward radiative transfer problem. Then, the sensitivity of a current shape with respect to a prescribed cost functional is obtained from adjoint analysis, and it is used to update the shape. The entire inverse design procedures including the forward and adjoint analyses as well as the shape update based on the adjoint analysis are implemented to an open-source CFD solver, OpenFOAM. Two pure radiation problems from the existing literature are considered for verifying the proposed method. The resultant optimal geometries in both problems show good agreement with the reference data. Finally, the present method is applied to a more complex three-dimensional problem to validate the strength of its applicability to complex geometries.

A modified equation analysis for immersed boundary methods based on volume penalization: Applications to linear advection–diffusion equations and high-order discontinuous Galerkin schemes

The Immersed Boundary Method (IBM) is a popular numerical approach to impose boundary conditions without relying on body-fitted grids, thus reducing the costly effort of mesh generation. To obtain enhanced accuracy, IBM can be combined with high-order methods (e.g., discontinuous Galerkin). For this combination to be effective, an analysis of the numerical errors is essential. In this work, we apply, for the first time, a modified equation analysis to the combination of IBM (based on volume penalization) and high-order methods (based on nodal discontinuous Galerkin methods) to analyze a priori numerical errors and obtain practical guidelines on the selection of IBM parameters. The analysis is performed on a linear advection–diffusion equation with Dirichlet boundary conditions. Three ways to penalize the immersed boundary are considered, the first penalizes the solution inside the IBM region (classic approach), whilst the second and third penalize the first and second derivatives of the solution. We find optimal combinations of the penalization parameters, including the first and second penalizing derivatives, resulting in minimum errors. We validate the theoretical analysis with numerical experiments for one- and two-dimensional advection–diffusion equations.