Volume Penalization Research Articles

Approximate acoustic boundary conditions in the time-domain using volume penalization.

This paper presents an immersed boundary method for modeling complex impedance boundary conditions in wave-based finite-difference time-domain simulations. The fully parallelizable and physically motivated Brinkman method allows for the representation of complex geometries on simple Cartesian grids as porous material by introducing a friction term and an effective volume. The parameters are specified using blending functions, enabling impedance boundary conditions without the need for grid fitting or special boundary treatment. Representative acoustic configurations are analyzed to assess the method. In detail, acoustic materials on and in front of a rigid wall, a reacting surface as well as fully reflecting walls are examined. Comparison with analytical solutions shows satisfactory agreement of the resulting impedances in the range from 20 Hz up to 4 kHz. The method is derived for the (non-)linear Euler equations and the acoustic wave equation. An extensive stability analysis is carried out.

Computational aerodynamics of insect flight using volume penalization

The state-of-the-art of insect flight research using advanced computational fluid dynamics techniques on supercomputers is reviewed, focusing mostly on the work of the present authors. We present a brief historical overview, discuss numerical challenges and introduce the governing model equations. Two open source codes, one based on Fourier, the other based on wavelet representation, are succinctly presented and a mass-spring flexible wing model is described. Various illustrations of numerical simulations of flapping insects at low, intermediate and high Reynolds numbers are presented. The role of flexible wings, data-driven modeling and fluid–structure interaction issues are likewise discussed.

Numerical analysis of heat and mass transfer through beds of spherical and non-spherical elements

Many issues related to mass and heat transfer through beds of granular materials are still not fully understood. In this work, non-isothermal turbulent flow is analysed within granular layers of spherical and non-spherical elements. We apply a volume penalization (VP) approach formulated in the framework of an immersed boundary technique (IB) on Cartesian computational meshes. It allows modelling flows around solid objects with almost arbitrarily complex shapes and in any form of contact. The validation of the solution accuracy is performed against ANSYS Fluent simulations using body-fitted meshes and experimental literature data. It shows the capability of the IB-VP approach for the simulations of flows in complex geometries. The main research focuses on the comparison of the influence of various types of particles and their temperature on vorticity, turbulence level and pressure drop inside and behind the granular bed. In particular, we analyse how the shape of the solid particles affects the efficiency of heat transfer in different flow conditions. The obtained results reveal the occurrence of very complex flow structures (recirculation and stagnation regions) inside beds. Comparison of results also point out preferred configurations of the beds.

Comparison of one-dimensional and three-dimensional glottal flow models in left-right asymmetric vocal fold conditions.

While the glottal flow is often simplified as one-dimensional (1D) in computational models of phonation to reduce computational costs, the 1D flow model has not been validated in left-right asymmetric vocal fold conditions, as often occur in both normal and pathological voice production. In this study, we performed three-dimensional (3D) and 1D flow simulations coupled to a two-mass model of adult male vocal folds and compared voice production at different degrees of left-right stiffness asymmetry. The flow and acoustic fields in 3D were obtained by solving the compressible Navier-Stokes equations using the volume penalization method with the moving vocal fold wall as an immersed boundary. Despite differences in the predicted flow pressure on vocal fold surface between the 1D and 3D flow models, the results showed reasonable agreement in vocal fold vibration patterns and selected voice outcome measures between the 1D and 3D models for the range of left-right asymmetric conditions investigated. This indicates that vocal fold properties play a larger role than the glottal flow in determining the overall pattern of vocal fold vibration and the produced voice, and the 1D flow simplification is sufficient in modeling phonation, at least for the simplified glottal geometry of this study.

Penalized Wall Function Method for Turbulent Flow Modeling

A novel penalized wall function method for simulations of wall-bounded compressible turbulent flows is proposed. The new approach is based on the Reynolds-averaged Navier–Stokes (RANS) equations to model the outer region of the turbulent boundary layer, while the inner part is approximated by the equilibrium wall function model. 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 and to impose the wall-stress boundary conditions on the RANS solution. The exchange location is specified implicitly through a localized source term in the boundary layer equation, written as a function of the normalized distance from the wall. The wall-stress condition is determined by solving an auxiliary equation for the wall-stress, ensuring the correct matching of the RANS and the wall function solutions at the exchange layer. The proposed method noticeably reduces the near-wall mesh resolution requirements without significant modification of the RANS solver and removes the ill-defined explicit matching procedure, commonly used by traditional wall function-based methods. The penalized wall function approach is implemented using the vertex-centered control volume method on unstructured computational grids. The effectiveness of the developed penalized wall function method is demonstrated for twodimensional bump-in-channel flow for the Spalart–Allmaras turbulence model.

Numerical modeling of a self-propelled dolphin jump out of water

A computational model is developed to investigate the jump of a self-propelled dolphin out of water. This model relies on the Navier–Stokes equations, where a fictitious domain approach with the volume penalization method is used for fluid-structure coupling, and the continuous surface force approach is used to model the water–air interface, the latter being tracked in a level-set framework. The dolphin's geometry is based on freely available data from the literature. While body deformation is imposed, the leading linear and angular displacements are computed from Newton’s laws. Numerical simulations show that it is necessary to generate large propulsives forces to allow the jump out of water. When the dolphin is out of water, its trajectory follows a purely ballistic one.

A combined volume penalization / selective frequency damping approach for immersed boundary methods applied to high-order schemes

There has been an increasing interest in developing efficient immersed boundary method (IBM) based on Cartesian grids, recently in the context of high-order methods. IBM based on volume penalization is a robust and easy to implement method to avoid body-fitted meshes and has been recently adapted to high order discretizations [1]. This work proposes an improvement over the classic penalty formulation for high-order solvers based on flux reconstruction. We include a selective frequency damping (SFD) approach [2] acting only inside solid body defined through the immersed boundary masking, to damp spurious oscillations. An encapsulated formulation for the SFD method is implemented, which can be used as a wrapper around an existing time-stepping code. The numerical properties have been studied through eigensolution analysis based on the advection equation. These studies not only show the advantages of using the SFD method as an alternative of the traditional volume penalization, but also show the favorable properties of combining both approaches. This new approach is then applied to the Navier-Stokes equation to simulate steady flow past an airfoil and unsteady flow past a circular cylinder. The advantages of the SFD method in providing improved accuracy are reported.

Volume penalization method for solving coupled radiative-conductive heat transfer problems in complex geometries

A novel immersed boundary method for simulating coupled radiative-conductive heat transfer based on a volume penalization method is proposed and verified. A smooth source/sink term is introduced to the governing equations for radiative and conductive heat transfer in order to take into account the absorption, reflection and emission at a solid surface. The proposed algorithm is successfully implemented into an open-source CFD solver, OpenFOAM with adaptive mesh refinement to guarantee sufficiently fine grid resolution near solid surfaces. Benchmark problems of pure radiation and combined conduction-radiation are considered for verifying the proposed method. The present simulation results under various absorption and emission conditions show good agreement with analytical solutions as well as existing literature data within around 5% errors. In addition to its superior performance in terms of accuracy, the present approach does not require special treatments around the solid surface with arbitrary complex geometry to impose boundary conditions for radiative heat transfer, and thereby wide applicability in scientific and engineering problems.

Improved Gulf Stream separation through Brinkman penalization

The advantage of a smooth representation of bathymetry in terrain-following σ-coordinate ocean models is compromised by the need to avoid numerical errors on steep slopes associated with pressure gradient discretization or spurious diapycnal diffusion. Geopotential z-coordinate models avoid these errors, but greatly underrepresent the interaction of flow with a topographic slope, especially when the bathymetry is underresolved. Hybrid coordinate models are also deficient because it is difficult to find a satisfactory compromise between z and σ coordinates. More general vertical coordinates (not just combinations of z and σ) can also be used, in particular for ocean interior, but without solving the problems associated with the representation of bathymetry. With volume penalization, we do not seek a compromise, but rather a correction to the usual coordinate systems that realistically recovers continuous and steep bathymetry. The Brinkman volume penalization method studied here is a modified version of the one introduced in Debreu et al.(2020) that simplifies the numerical implementation of the penalization, increases robustness and improves its computational performance for realistic long-term simulations, while preserving accuracy. We apply this penalization method to the Gulf Stream separation problem that has puzzled modelers for decades. The method improves the representation of the flow-topography interaction and achieves realistic separation of the Gulf Stream at resolutions as coarse as 1/8°. In addition, it provides a tool to separate the effect of eddy activity and topographic slope when changing grid resolution. This has never before been possible because at coarse resolution none of the usual coordinate systems can properly represent a steep continental slope. Our results show that realistic bathymetry is more important than eddy activity in ensuring realistic Gulf Stream separation, even though many recent studies tend to focus on the eddy activity. A steep slope can exert a stabilizing influence that promotes a strong mean slope current with strong inertia that helps it separate from the coast at the topographic curvature of Cape Hatteras. We anticipate that a successful topographic slope correction will be very valuable for climate models, as their current resolution is far from sufficient to represent western boundary currents (WBCs) using traditional coordinate systems. Our results suggest that a climate model with a 1/4° resolution using volume penalization — and perhaps also some parameterization of the eddy-mean flow interaction to energize the WBCs — could represent ocean circulation much more realistically than a model at the same resolution, but without volume penalization.

A particle resolved simulation approach for studying shock interactions with moving, colliding solid particles

This work applies a new combination of techniques for the fully resolved simulation of compressible, gas–particle multiphase flows. The adaptive wavelet collocation method is used to dynamically, and efficiently, adapt the computational grid to localized flow features and the particles. A characteristic-based volume penalization method that imposes arbitrary Dirichlet, Neumann, or Robin-type immersed boundary conditions, is used to enforce the no-slip condition at particle surfaces. A hard-sphere collision model is applied to capture the particle–particle collisions. Proof of concept test cases are presented, showcasing the dynamic grid adaptation and fully resolved two-way coupling between the phases that is possible with this approach. Results for a shock-driven single cylinder under viscous and inviscid conditions are presented along with a demonstration of a shock interacting with a cloud of randomly distributed cylinders and spheres.

A nonperiodic Chebyshev spectral method avoiding penalization techniques for a class of nonlinear peridynamic models

AbstractIn the framework of elastodynamics, peridynamics is a nonlocal theory able to capture singularities without using partial derivatives. The governing equation is a second order in time partial integro‐differential equation. In this article, we focus on a one‐dimensional nonlinear model of peridynamics and propose a spectral method based on the Chebyshev polynomials to discretize in space. The main capability of the method is that it avoids the assumption of periodic boundary condition in the solution and can benefit of the use of the fast Fourier transform. We show its convergence and find that the method results to be very efficient in terms of accuracy and execution time with respect to spectral methods based on the Fourier trigonometric polynomials associated to a volume penalization technique.

Comparison of one-dimensional and three-dimensional glottal flow models in left-right asymmetric vocal fold conditions

While the glottal flow is often simplified as one-dimensional (1D) in computational models of phonation to reduce computational costs, the 1D flow model has not been validated in left-right asymmetric vocal fold conditions, as often occur in both normal and pathological voice production. It is unclear to what extent the 1D model approximates the effect of three-dimensional (3D) flow phenomena and fluid-structure interaction. In this study, we performed 1D and 3D flow simulations coupled with the two-mass vocal fold model and compared vocal fold vibration patterns at different degrees of left-right stiffness asymmetry. The flow and acoustic fields in 3D were predicted by solving the compressible Navier-Stokes equations using the volume penalization method and considering the moving vocal fold wall as an immersed boundary. The results showed that vocal fold vibration amplitudes and left-right phase differences in the 3D flow were predicted by the 1D flow model under conditions of small left-right asymmetry, while vocal fold vibration amplitudes were underestimated at conditions of large left-right asymmetry. This indicates that 1D flow models may be sufficient in modeling phonation under left-right asymmetric conditions, although the performance can be further improved by more accurately predicting air pressure on vocal fold surface.

A pseudo-spectral based efficient volume penalization scheme for Cahn–Hilliard equation in complex geometries

In this paper, we have developed an efficient volume penalization based diffuse-filtering scheme to solve variable mobility based Cahn–Hilliard equation in complex geometries. The main novelty of the work is that a dealiased pseudo-spectral scheme with immersed interface method (IIM) is proposed for solving any generalized concentration-dependent mobility function-based Cahn–Hilliard (CH) equation in complicated computational domains. An indicator function based mobility parameter is introduced to perform simulation of binary spinodal decomposition problem at a lower computational expense in complex geometries by solving a single phase field equation. Due to the smooth removal of high frequency Fourier components, the solution of the present RK4 based diffuse-filtering scheme does not display spurious currents when suitable low-pass filtering strategy and adequately resolved mobility indicator are incorporated. The traditional and memory optimized zero padding schemes are also implemented to show the comparative performance of different dealiasing schemes for the variable mobility based Cahn–Hilliard equation. It is found that the diffuse-filtering scheme displays reasonable accuracy similar to the zero padding based schemes but its average CPU time is significantly lower, which indicates better computational performance of the scheme for the variable mobility Cahn–Hilliard equation. Time variation of the characteristic length scale during spinodal decomposition of a binary mixture agrees well with the analytical prediction. The optimal three stage SSPRK3 temporal scheme is employed and it is found that time step size can be increased approximately 1.4 times than the classical RK4 scheme reducing total CPU time. Oscillation free numerical solution and conservation of order parameter are obtained for the complex geometry based spinodal decomposition problem. A radially-averaged structure factor is introduced to quantify resolution issues of the dealiasing schemes for the spinodal decomposition problem in different complex geometries.

Computational aeroacoustics of aerofoil leading edge noise using the volume penalization-based immersed boundary methods

Broadband noise due to the turbulence-aerofoil interaction, which is also called the leading edge noise, is one of the major noise sources of aircraft (including the engine). To study the noise properties numerically is a popular approach with the increasing power of computers. Conventional approaches of using body-fitted grids at the boundaries would be convoluted due to the complex geometries, which can constrain the efficiency of parametric studies. A promising approach to tackle this issue is to use the immersed boundary method (IBM). Among various IBM variants, the volume penalization (VP) approach employs a masking function to identify the immersed solid boundary, and continuous forcing terms are added to the original flow governing equations to account for the boundary conditions. It is, therefore, efficient and easy to implement into the existing computational aeroacoustics solvers. In this work, the VP-based IBM is used to simulate the leading edge noise by combining with the advanced synthetic turbulence method. The simulations are conducted for both the isolated aerofoils and cascade, and the results are compared with the well-validated body-fitted grid solutions. The viscosity effect is also highlighted by comparing the results obtained by solving both Euler and Navier–Stokes equations.

Pressure-Tight and Non-stiff Volume Penalization for Compressible Flows

Embedding geometries in structured grids allows a simple treatment of complex objects in fluid simulations. Various methods for embedding geometries are available. The commonly used Brinkman-volume-penalization models geometries as porous media, and approximates a solid object in the limit of vanishing porosity. In its simplest form, the momentum equations are augmented by a term penalizing the fluid velocity, yielding good results in many applications. However, it induces numerical stiffness, especially if high-pressure gradients need to be balanced. Here, we focus on the effect of the reduced effective volume (commonly called porosity) of the porous medium. An approach is derived, which allows reducing the flux through objects to practically zero with little increase of numerical stiffness. Also, non-slip boundary conditions and adiabatic boundary conditions are easily constructed. The porosity terms allow keeping the skew symmetry of the underlying numerical scheme, by which the numerical stability is improved. Furthermore, very good conservation of mass and energy in the non-penalized domain can be achieved, for which the boundary smoothing introduces a small ambiguity in its definition. The scheme is tested for acoustic scenarios, for near incompressible and strongly compressible flows.

A coupled Two-relaxation-time Lattice Boltzmann-Volume penalization method for flows past obstacles

In this article, a coupled Two-relaxation-time Lattice Boltzmann-Volume penalization (TRT-LBM-VP) method is presented to simulate flows past obstacles. Two relaxation times are used in the collision operator, of which one is related to the fluid viscosity and the other one is related to the numerical stability and accuracy. The volume penalization method is introduced into the TRT-LBM by an external forcing term. In the procedure of the TRT-LBM-VP, the processes of interpolating velocities on the boundaries points and distributing the force density to the Eulerian points are unneeded. Performing the TRT-LBM-VP on a certain point, only the variables of this point are needed. As a consequence, the TRT-LBM-VP can be conducted parallelly. From the comparison between the result of the cylindrical Couette flow solved by the TRT-LBM-VP and that solved by the Single-relaxation-time LBM-VP (SRT-LBM-VP), the accuracy of the TRT-LBM-VP is higher than that of the SRT-LBM-VP. Flows past a single circular cylinder, a pair of cylinders in tandem and side-by-side arrangements, two counter-rotating cylinders and a NACA-0012 airfoil are chosen as numerical experiments to verify the present method further. Good agreements between the present results and those in the previous literatures are achieved.

A volume penalization immersed boundary method for flow interactions with aquatic vegetation

A volume-penalization immersed boundary (VPIB) method was developed to study flow interactions with aquatic vegetation. The model has been validated with data from laboratory experiments and previous high-fidelity models with satisfactory results. Sensitivity analyzes on both penalty parameter and thickness parameter were conducted, and optimal values for these parameters are recommended. The validated model has been applied to study the effects of swaying motion of vegetation stems on the flow dynamics at both vegetate-stem scale and patch scale. The swaying motion of the vegetation stem is prescribed following a cubic law that peaks at the top and decreases to zero at the bottom. At stem-scale, the hydrodynamics depend on the Keulegan Carpenter number (KC), which is defined as the maximum excursion of the vegetation stem to the diameter of the stem. Simulations with three KC values were carried out. For KC≥1, the flow turbulence is significantly enhanced by the swaying motion of the stem, and turbulence becomes more isotropic in the wake. The swaying motion of vegetation stems caused a 5% increase of the bottom shear stress at the shoulders of the stem, and the effect is negligible in the wake. At patch-scale, the hydrodynamics depend on the effective Keulegan Carpenter number based on the patch size of the vegetation patch, and the solid volume fraction for dense vegetation canopy. Solid volume fraction was varied while maintaining the same effective Keulegan Carpenter in the simulations. When the effective Keulgen Carpenter number is small (KC<1), effects of the swaying motion of vegetation stems on the large patch-scale dynamics are not significant, including both the turbulence statistics and the bottom stress.

波状壁の上を通過する粘弾性流体のせん断流れの数値解析

A two-dimensional shear flow of a viscoelastic fluid over a wavy wall was numerically analyzed by the finite difference method. The upper convected Maxwell (UCM) model was used as the constitutive equation, and the unsteady flow problem was solved numerically using the stream function-vorticity method. Variables are discretized using a staggered grid with uniform mesh width. The volume penalization method was used to represent the moving upper wall and the fixed lower wavy wall that induces a vorticity perturbation. The flow in an infinitely long channel is computed using periodic boundary conditions in the flow direction. The results show that the vorticity intensifies over the wavy surface with positive and negative vorticity stripes, at which extreme wall-normal velocity fluctuations are also observed. At relatively high Deborah numbers, the lower ends of these extreme fluctuation regions are close to the critical line associated with a change of type of the governing partial differential equations in a Couette flow. In contrast, the upper ends of the extreme fluctuation region are close to the critical layer of the shear flow, where the flow velocity coincides with the viscoelastic Alfvén wave speed.

Volume Penalization Method for Solving Coupled Radiative-Conductive Heat Transfer Problems in Complex Geometries

Volume Penalization Method for Solving Coupled Radiative-Conductive Heat Transfer Problems in Complex Geometries

Eigensolution analysis of immersed boundary method based on volume penalization: Applications to high-order schemes

This paper presents eigensolution and non-modal analyses for immersed boundary methods (IBMs) based on volume penalization for the linear advection equation. This approach is used to analyze the behavior of flux reconstruction (FR) discretization, including the influence of polynomial order and penalization parameter on numerical errors and stability. Through a semi-discrete analysis, we find that the inclusion of IBM adds additional dissipation without changing significantly the dispersion of the original numerical discretization. This agrees with the physical intuition that in this type of approach, the solid wall is modelled as a porous medium with vanishing viscosity. From a stability point of view, the variation of penalty parameter can be analyzed based on a fully-discrete analysis, which leads to practical guidelines on the selection of penalization parameter. Numerical experiments indicate that the penalization term needs to be increased to damp oscillations inside the solid (i.e. porous region), which leads to undesirable time step restrictions. As an alternative, we propose to include a second-order term in the solid for the no-slip wall boundary condition. Results show that by adding a second-order term we improve the overall accuracy with relaxed time step restriction. This indicates that the optimal value of the penalization parameter and the second-order damping can be carefully chosen to obtain a more accurate scheme. Finally, the approximated relationship between these two parameters is obtained and used as a guideline to select the optimum penalty terms in a Navier-Stokes solver, to simulate flow past a cylinder.