Pub Date : 2025-10-25DOI: 10.1016/j.compfluid.2025.106885
Janina Bender , Thomas Izgin , Philipp Öffner , Davide Torlo
For hyperbolic conservation laws, the famous Lax–Wendroff theorem delivers sufficient conditions for the limit of a convergent numerical method to be a weak (entropy) solution. This theorem is a fundamental result, and many investigations have been done to verify its validity for finite difference, finite volume, and finite element schemes, using either explicit or implicit linear time-integration methods.
Recently, the use of modified Patankar (MP) schemes as time-integration methods for the discretization of hyperbolic conservation laws has gained increasing interest. These schemes are unconditionally conservative and positivity-preserving and only require the solution of a linear system. However, MP schemes are by construction nonlinear, which is why the theoretical investigation of these schemes is more involved. We prove an extension of the Lax–Wendroff theorem for the class of MP methods. This is the first extension of the Lax–Wendroff theorem to nonlinear time integration methods with just an additional hypothesis on the total time variation boundedness of the numerical solutions. We provide some numerical simulations that validate the theoretical observations.
{"title":"The Lax–Wendroff theorem for Patankar-type methods applied to hyperbolic conservation laws","authors":"Janina Bender , Thomas Izgin , Philipp Öffner , Davide Torlo","doi":"10.1016/j.compfluid.2025.106885","DOIUrl":"10.1016/j.compfluid.2025.106885","url":null,"abstract":"<div><div>For hyperbolic conservation laws, the famous Lax–Wendroff theorem delivers sufficient conditions for the limit of a convergent numerical method to be a weak (entropy) solution. This theorem is a fundamental result, and many investigations have been done to verify its validity for finite difference, finite volume, and finite element schemes, using either explicit or implicit linear time-integration methods.</div><div>Recently, the use of modified Patankar (MP) schemes as time-integration methods for the discretization of hyperbolic conservation laws has gained increasing interest. These schemes are unconditionally conservative and positivity-preserving and only require the solution of a linear system. However, MP schemes are by construction nonlinear, which is why the theoretical investigation of these schemes is more involved. We prove an extension of the Lax–Wendroff theorem for the class of MP methods. This is the first extension of the Lax–Wendroff theorem to nonlinear time integration methods with just an additional hypothesis on the total time variation boundedness of the numerical solutions. We provide some numerical simulations that validate the theoretical observations.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"304 ","pages":"Article 106885"},"PeriodicalIF":3.0,"publicationDate":"2025-10-25","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145414644","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-24DOI: 10.1016/j.compfluid.2025.106877
Spencer H. Bryngelson
Ensemble-averaged polydisperse bubbly flow models require statistical moments of the evolving bubble size distribution. Under step forcing, these moments reach statistical equilibrium in finite time. However, the transitional phase before equilibrium and cases with time-dependent forcing are required to predict flow in engineering applications. Computing these moments is expensive because the integrands are highly oscillatory, even when the bubble dynamics are linear. Ensemble-averaged models compute these moments at each grid point and time step, making cost reduction important for large-scale bubbly flow simulations. Traditional methods evaluate the integrals via traditional quadrature rules. This approach requires a large number of quadrature nodes in the equilibrium bubble size, each equipped with its own advection partial differential equation (PDE), resulting in significant computational expense. We formulate a Levin collocation method to reduce this cost. Given the differential equation associated with the integrand, or moment, the method approximates it by evaluating its derivative via polynomial collocation. The differential matrix and amplitude function are well-suited to numerical differentiation via collocation, and so the computation is comparatively cheap. For an example excited polydisperse bubble population, the first moment is computed with the presented method at relative error with 100 times fewer quadrature nodes than the trapezoidal rule. The gap increases for smaller target relative errors: the Levin method requires times fewer points for a relative error of . The formulated method maintains constant cost as the integrands become more oscillatory with time, making it particularly attractive for long-time simulations. Mechanistically, the transient behavior of the moments is set by two effects: resonance detuning across bubble sizes, which causes phase mixing of oscillations, and viscous damping, which removes radial kinetic energy. The proposed formulation isolates the oscillations while keeping the remaining terms smooth, so accuracy does not deteriorate at late times.
{"title":"Fast integration method for averaging polydisperse bubble population dynamics","authors":"Spencer H. Bryngelson","doi":"10.1016/j.compfluid.2025.106877","DOIUrl":"10.1016/j.compfluid.2025.106877","url":null,"abstract":"<div><div>Ensemble-averaged polydisperse bubbly flow models require statistical moments of the evolving bubble size distribution. Under step forcing, these moments reach statistical equilibrium in finite time. However, the transitional phase before equilibrium and cases with time-dependent forcing are required to predict flow in engineering applications. Computing these moments is expensive because the integrands are highly oscillatory, even when the bubble dynamics are linear. Ensemble-averaged models compute these moments at each grid point and time step, making cost reduction important for large-scale bubbly flow simulations. Traditional methods evaluate the integrals via traditional quadrature rules. This approach requires a large number of quadrature nodes in the equilibrium bubble size, each equipped with its own advection partial differential equation (PDE), resulting in significant computational expense. We formulate a Levin collocation method to reduce this cost. Given the differential equation associated with the integrand, or moment, the method approximates it by evaluating its derivative via polynomial collocation. The differential matrix and amplitude function are well-suited to numerical differentiation via collocation, and so the computation is comparatively cheap. For an example excited polydisperse bubble population, the first moment is computed with the presented method at <span><math><msup><mn>10</mn><mrow><mo>−</mo><mn>3</mn></mrow></msup></math></span> relative error with 100 times fewer quadrature nodes than the trapezoidal rule. The gap increases for smaller target relative errors: the Levin method requires <span><math><msup><mn>10</mn><mn>4</mn></msup></math></span> times fewer points for a relative error of <span><math><msup><mn>10</mn><mrow><mo>−</mo><mn>8</mn></mrow></msup></math></span>. The formulated method maintains constant cost as the integrands become more oscillatory with time, making it particularly attractive for long-time simulations. Mechanistically, the transient behavior of the moments is set by two effects: resonance detuning across bubble sizes, which causes phase mixing of oscillations, and viscous damping, which removes radial kinetic energy. The proposed formulation isolates the oscillations while keeping the remaining terms smooth, so accuracy does not deteriorate at late times.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"304 ","pages":"Article 106877"},"PeriodicalIF":3.0,"publicationDate":"2025-10-24","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145414646","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-24DOI: 10.1016/j.compfluid.2025.106883
S. Bennie , M. Fossati
Presented in the following work is a comprehensive analysis of wake vortex encounters with residential structures. From the results of high fidelity LES simulations, the dynamics and underlying flow structures which govern these potentially damaging encounters have been identified. Through evaluation of the pressure loads transmitted to the roof surface, the potential for damage to occur to a residential structure as a result of wake vortex exposure has been evaluated for a variety of cases. Regarding the building’s design, structures possessing larger pitch angles and thus steeper roofs have been found to sustain the largest peak loads for their encounter with an identical wake vortex system as compared to their flatter roofed counterparts. Similarly, upon assessing the effect of the environmental conditions it was observed that for increasingly turbulent atmospheres, the wake vortex encounter would occur sooner and with a reduced intensity compared to more neutral conditions. These behaviours have been attributed to the effects of secondary flow structures formed from the shedding of vorticity from the building surface or from wake vortex interactions with the eddies that comprise the atmospheric environment. These secondary flow structures energise wake vortex instability mechanisms thus leading to the variations in pressure loads sustained by the roof. With respect to the impact orientation, we note that there exists a minimal difference on the pressure loads generated during a wake vortex encounter for small angular offsets up to .
{"title":"Aircraft wake vortex encounters with residential structures","authors":"S. Bennie , M. Fossati","doi":"10.1016/j.compfluid.2025.106883","DOIUrl":"10.1016/j.compfluid.2025.106883","url":null,"abstract":"<div><div>Presented in the following work is a comprehensive analysis of wake vortex encounters with residential structures. From the results of high fidelity LES simulations, the dynamics and underlying flow structures which govern these potentially damaging encounters have been identified. Through evaluation of the pressure loads transmitted to the roof surface, the potential for damage to occur to a residential structure as a result of wake vortex exposure has been evaluated for a variety of cases. Regarding the building’s design, structures possessing larger pitch angles and thus steeper roofs have been found to sustain the largest peak loads for their encounter with an identical wake vortex system as compared to their flatter roofed counterparts. Similarly, upon assessing the effect of the environmental conditions it was observed that for increasingly turbulent atmospheres, the wake vortex encounter would occur sooner and with a reduced intensity compared to more neutral conditions. These behaviours have been attributed to the effects of secondary flow structures formed from the shedding of vorticity from the building surface or from wake vortex interactions with the eddies that comprise the atmospheric environment. These secondary flow structures energise wake vortex instability mechanisms thus leading to the variations in pressure loads sustained by the roof. With respect to the impact orientation, we note that there exists a minimal difference on the pressure loads generated during a wake vortex encounter for small angular offsets up to <span><math><msup><mn>20</mn><mo>∘</mo></msup></math></span>.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"304 ","pages":"Article 106883"},"PeriodicalIF":3.0,"publicationDate":"2025-10-24","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145464296","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-22DOI: 10.1016/j.compfluid.2025.106882
Alex Kleb, Krzysztof J. Fidkowski, Joaquim R.R.A. Martins
Computational fluid dynamics is essential for designing aircraft, turbines, and other engineering systems. However, generating suitable computational meshes for complex geometries remains the primary bottleneck in analysis workflows, often requiring days of manual effort. Traditional boundary-conforming meshes excel at capturing near-wall physics in viscous flows but demand specialized expertise and extensive preprocessing time. Cartesian cut-cell methods provide automatic mesh generation for complex geometries in minutes, yet they struggle with high Reynolds number viscous flows where boundary layers exhibit rapid velocity changes that require prohibitively fine resolution for isotropic elements. The fundamental challenge is accurately modeling boundary layer physics on automatically generated meshes without sacrificing the computational efficiency that makes such methods attractive. In this work, we show that an ordinary differential equation (ODE) wall function incorporating pressure-momentum balance enables accurate high Reynolds number viscous flow predictions on coarse Cartesian cut-cell meshes. Our approach solves a one-dimensional boundary value problem at each wall boundary face that accounts for the transition from the viscous dominated near-wall region to the inviscid wake region, allowing forcing points to operate effectively at . Unlike traditional wall functions, the ODE is not limited to the logarithmic layer and maintains accuracy in strong pressure gradient environments typical of aerodynamic applications. The ODE can achieve correct skin friction predictions on meshes more than four times coarser than analytical wall functions require. The ODE wall functions are solved with a robust Newton–Krylov implementation that utilizes adaptive mesh refinement. It converges reliably across diverse flow conditions while solving hundreds of degrees of freedom in fewer than ten linear iterations. These results demonstrate that automatic high-fidelity viscous flow analysis is achievable without manual mesh generation expertise.
{"title":"Solving high Reynolds number flows on Cartesian cut-cell meshes using an ODE wall function with momentum balance","authors":"Alex Kleb, Krzysztof J. Fidkowski, Joaquim R.R.A. Martins","doi":"10.1016/j.compfluid.2025.106882","DOIUrl":"10.1016/j.compfluid.2025.106882","url":null,"abstract":"<div><div>Computational fluid dynamics is essential for designing aircraft, turbines, and other engineering systems. However, generating suitable computational meshes for complex geometries remains the primary bottleneck in analysis workflows, often requiring days of manual effort. Traditional boundary-conforming meshes excel at capturing near-wall physics in viscous flows but demand specialized expertise and extensive preprocessing time. Cartesian cut-cell methods provide automatic mesh generation for complex geometries in minutes, yet they struggle with high Reynolds number viscous flows where boundary layers exhibit rapid velocity changes that require prohibitively fine resolution for isotropic elements. The fundamental challenge is accurately modeling boundary layer physics on automatically generated meshes without sacrificing the computational efficiency that makes such methods attractive. In this work, we show that an ordinary differential equation (ODE) wall function incorporating pressure-momentum balance enables accurate high Reynolds number viscous flow predictions on coarse Cartesian cut-cell meshes. Our approach solves a one-dimensional boundary value problem at each wall boundary face that accounts for the transition from the viscous dominated near-wall region to the inviscid wake region, allowing forcing points to operate effectively at <span><math><mrow><msup><mi>y</mi><mo>+</mo></msup><mo>></mo><mn>600</mn></mrow></math></span>. Unlike traditional wall functions, the ODE is not limited to the logarithmic layer and maintains accuracy in strong pressure gradient environments typical of aerodynamic applications. The ODE can achieve correct skin friction predictions on meshes more than four times coarser than analytical wall functions require. The ODE wall functions are solved with a robust Newton–Krylov implementation that utilizes adaptive mesh refinement. It converges reliably across diverse flow conditions while solving hundreds of degrees of freedom in fewer than ten linear iterations. These results demonstrate that automatic high-fidelity viscous flow analysis is achievable without manual mesh generation expertise.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"304 ","pages":"Article 106882"},"PeriodicalIF":3.0,"publicationDate":"2025-10-22","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145414647","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-20DOI: 10.1016/j.compfluid.2025.106880
Niccolò Tonicello , Guido Lodato , Matthias Ihme
The present work focuses on the extension of the Spectral Difference (SD) scheme to the five-equation Baer-Nunziato model for the simulation of immiscible compressible fluids. This five-equation model is augmented with the Allen-Cahn regularisation to avoid both over-diffusion and over-thinning of the phase field representing the interface. In order to preserve contact discontinuities, in the reconstruction step of the SD scheme, a change of variables from conservative to primitive is used. This approach is shown to be beneficial in avoiding pressure oscillations at material interfaces. An extensive series of numerical tests, considering both two- and three-dimensional problems, are performed to assess accuracy and robustness of the present method. Specifically, both laminar and turbulent flows, as well as low-Mach and highly compressible flows, are considered, including cases with and without surface tension. The proposed change of variables is shown to improve the stability of the scheme, significantly reducing pressure oscillations at the material interfaces. This improved robustness enables the method to achieve accurate and stable solutions across a broad range of flow conditions.
{"title":"Extension of a spectral difference method for the diffused-interface five-equation model","authors":"Niccolò Tonicello , Guido Lodato , Matthias Ihme","doi":"10.1016/j.compfluid.2025.106880","DOIUrl":"10.1016/j.compfluid.2025.106880","url":null,"abstract":"<div><div>The present work focuses on the extension of the Spectral Difference (SD) scheme to the five-equation Baer-Nunziato model for the simulation of immiscible compressible fluids. This five-equation model is augmented with the Allen-Cahn regularisation to avoid both over-diffusion and over-thinning of the phase field representing the interface. In order to preserve contact discontinuities, in the reconstruction step of the SD scheme, a change of variables from conservative to primitive is used. This approach is shown to be beneficial in avoiding pressure oscillations at material interfaces. An extensive series of numerical tests, considering both two- and three-dimensional problems, are performed to assess accuracy and robustness of the present method. Specifically, both laminar and turbulent flows, as well as low-Mach and highly compressible flows, are considered, including cases with and without surface tension. The proposed change of variables is shown to improve the stability of the scheme, significantly reducing pressure oscillations at the material interfaces. This improved robustness enables the method to achieve accurate and stable solutions across a broad range of flow conditions.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"304 ","pages":"Article 106880"},"PeriodicalIF":3.0,"publicationDate":"2025-10-20","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145414645","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-18DOI: 10.1016/j.compfluid.2025.106879
Jiali Tu , Haijian Yang , Shanlin Qin , Weifeng Guo , Rongliang Chen
Aortic dissection is a serious clinical condition characterized by a tear in the intima of the aortic wall. To better understand and treat this complex condition, researchers increasingly use hemodynamic simulations based on fluid-structure interaction (FSI). However, time-dependent three-dimensional FSI simulations in aortic dissection are complex and often inefficient in terms of computational time and memory. In this paper, we present a highly scalable parallel numerical method to simulate the blood flow in a full-sized aorta with dissection. The blood flow is modeled using the unsteady incompressible Navier-Stokes equations in an arbitrary Lagrangian-Eulerian framework, and the aortic wall is modeled as linear elastic material in a Lagrangian description. The entire system is discretized monolithically by a stabilized finite element method in space and a fully implicit scheme in time, and the resulting algebraic system is then solved using the Newton-Krylov-Schwarz method. Hemodynamic parameters within the aortic dissection are examined, revealing differences from simulations that rely solely on computational fluid dynamics. Scalability tests on a supercomputer demonstrate a parallel efficiency of 44.1 % with up to 2304 processor cores, reducing the simulation time for an entire cardiac cycle to 0.36 h.
{"title":"A highly scalable parallel simulation of blood flow with fluid-structure interaction in patient-specific aortic dissection","authors":"Jiali Tu , Haijian Yang , Shanlin Qin , Weifeng Guo , Rongliang Chen","doi":"10.1016/j.compfluid.2025.106879","DOIUrl":"10.1016/j.compfluid.2025.106879","url":null,"abstract":"<div><div>Aortic dissection is a serious clinical condition characterized by a tear in the intima of the aortic wall. To better understand and treat this complex condition, researchers increasingly use hemodynamic simulations based on fluid-structure interaction (FSI). However, time-dependent three-dimensional FSI simulations in aortic dissection are complex and often inefficient in terms of computational time and memory. In this paper, we present a highly scalable parallel numerical method to simulate the blood flow in a full-sized aorta with dissection. The blood flow is modeled using the unsteady incompressible Navier-Stokes equations in an arbitrary Lagrangian-Eulerian framework, and the aortic wall is modeled as linear elastic material in a Lagrangian description. The entire system is discretized monolithically by a stabilized finite element method in space and a fully implicit scheme in time, and the resulting algebraic system is then solved using the Newton-Krylov-Schwarz method. Hemodynamic parameters within the aortic dissection are examined, revealing differences from simulations that rely solely on computational fluid dynamics. Scalability tests on a supercomputer demonstrate a parallel efficiency of 44.1 % with up to 2304 processor cores, reducing the simulation time for an entire cardiac cycle to 0.36 h.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"303 ","pages":"Article 106879"},"PeriodicalIF":3.0,"publicationDate":"2025-10-18","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145359266","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-17DOI: 10.1016/j.compfluid.2025.106869
Gino I. Montecinos , Eleuterio F. Toro
The WENO-DK reconstruction [M. Dumbser and M. Käser, JCOMP.221:693-723, (2007)] is a type of WENO procedure in which, for the one-dimensional case only the leftmost, centered and rightmost stencils are involved. For even orders the central stencil contains more elements than degrees of freedom and an overdetermined system is solved by means of a least-squares approach. Here, it is numerically investigated the impact of choosing the smallest and largest central stencil around the cell of interests and proposed two variants to obtain the central polynomial where the solution of overdetermined systems is not needed. Implementations of the proposed approaches in the framework of fully discrete ADER schemes for the linear advection equation and the Euler equations of gas dynamics are reported. Comparisons with conventional WENO and conventional WENO-DK confirm that the proposed variants of WENO-DK are a suitable compromise between simplicity and accuracy in the context of ADER schemes, implemented up to the tenth order of accuracy in space and time.
{"title":"Variants for the WENO-DK reconstruction of even orders in the framework of ADER methods for very high orders of accuracy","authors":"Gino I. Montecinos , Eleuterio F. Toro","doi":"10.1016/j.compfluid.2025.106869","DOIUrl":"10.1016/j.compfluid.2025.106869","url":null,"abstract":"<div><div>The WENO-DK reconstruction [M. Dumbser and M. Käser, <em>JCOMP.</em> <strong>221</strong>:693-723, (2007)] is a type of WENO procedure in which, for the one-dimensional case only the leftmost, centered and rightmost stencils are involved. For even orders the central stencil contains more elements than degrees of freedom and an overdetermined system is solved by means of a least-squares approach. Here, it is numerically investigated the impact of choosing the smallest and largest central stencil around the cell of interests and proposed two variants to obtain the central polynomial where the solution of overdetermined systems is not needed. Implementations of the proposed approaches in the framework of fully discrete ADER schemes for the linear advection equation and the Euler equations of gas dynamics are reported. Comparisons with conventional WENO and conventional WENO-DK confirm that the proposed variants of WENO-DK are a suitable compromise between simplicity and accuracy in the context of ADER schemes, implemented up to the tenth order of accuracy in space and time.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"303 ","pages":"Article 106869"},"PeriodicalIF":3.0,"publicationDate":"2025-10-17","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145325680","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-17DOI: 10.1016/j.compfluid.2025.106881
A.R. Kocharina , D.V. Chirkov
The incompressible Navier-Stokes equations are solved using the finite-volume artificial compressibility method. A Godunov-type scheme with an exact Riemann solver is developed for the evaluation of inviscid fluxes across cell faces. To this end, the exact solution of the one-dimensional Riemann problem for the artificial compressibility equations is obtained using the method of -diagrams. The uniqueness of the solution is rigorously proven. The method is then extended to the multidimensional case. Two approaches for evaluation of the tangential velocity component are examined and discussed. A high-order variant of the Godunov scheme based on third-order MUSCL interpolation is proposed. At that, non-uniformity of the grid is taken into account. An implicit formulation of the scheme is developed, and the linearization process is described in detail. The proposed scheme is compared with the well-established Roe scheme through a series of steady-state two-dimensional benchmark problems, including inviscid and viscous flows around a circular cylinder and the 2D lid-driven cavity flow. The performance of the schemes on non-orthogonal grids is also investigated. Finally, both Roe and Godunov schemes are applied to the simulation of a three-dimensional turbulent flow in a hydraulic turbine flow passage. The results show that while both schemes exhibit comparable accuracy and convergence, the Godunov scheme offers advantages for inviscid simulations on highly non-orthogonal grids.
{"title":"Godunov scheme for numerical solution of incompressible Navier-Stokes equations","authors":"A.R. Kocharina , D.V. Chirkov","doi":"10.1016/j.compfluid.2025.106881","DOIUrl":"10.1016/j.compfluid.2025.106881","url":null,"abstract":"<div><div>The incompressible Navier-Stokes equations are solved using the finite-volume artificial compressibility method. A Godunov-type scheme with an exact Riemann solver is developed for the evaluation of inviscid fluxes across cell faces. To this end, the exact solution of the one-dimensional Riemann problem for the artificial compressibility equations is obtained using the method of <span><math><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>p</mi><mo>)</mo></mrow></math></span>-diagrams. The uniqueness of the solution is rigorously proven. The method is then extended to the multidimensional case. Two approaches for evaluation of the tangential velocity component are examined and discussed. A high-order variant of the Godunov scheme based on third-order MUSCL interpolation is proposed. At that, non-uniformity of the grid is taken into account. An implicit formulation of the scheme is developed, and the linearization process is described in detail. The proposed scheme is compared with the well-established Roe scheme through a series of steady-state two-dimensional benchmark problems, including inviscid and viscous flows around a circular cylinder and the 2D lid-driven cavity flow. The performance of the schemes on non-orthogonal grids is also investigated. Finally, both Roe and Godunov schemes are applied to the simulation of a three-dimensional turbulent flow in a hydraulic turbine flow passage. The results show that while both schemes exhibit comparable accuracy and convergence, the Godunov scheme offers advantages for inviscid simulations on highly non-orthogonal grids.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"304 ","pages":"Article 106881"},"PeriodicalIF":3.0,"publicationDate":"2025-10-17","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145360932","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-16DOI: 10.1016/j.compfluid.2025.106878
A.A. Gavrilov , A.V. Shebelev , A.V. Minakov
The paper presents the results of testing of a Eulerian model of two-phase turbulent non-Newtonian flow with coarse particles, proposed by the authors. The model includes equations for two-phase flow with rheological correlations and an equation for particle concentration transfer taking into account interfacial slip. The turbulence model takes into account the modulation of turbulence by particles. The proposed model has been validated on the problems of steady turbulent flow of shear thinning viscoplastic fluid with heavy particles in a horizontal pipe. The impact of Reynolds number and rheological parameters on the reliability of numerical simulations was examined. A comparison of experimental data with DNS-DEM simulation data has demonstrated that the proposed model is capable of accurately predicting the distribution of particle concentration, particle velocity, as well as carrier flow and pressure drop in the channel.
{"title":"Statistical model of turbulent flow of shear-thinning viscoplastic fluid with solid particles","authors":"A.A. Gavrilov , A.V. Shebelev , A.V. Minakov","doi":"10.1016/j.compfluid.2025.106878","DOIUrl":"10.1016/j.compfluid.2025.106878","url":null,"abstract":"<div><div>The paper presents the results of testing of a Eulerian model of two-phase turbulent non-Newtonian flow with coarse particles, proposed by the authors. The model includes equations for two-phase flow with rheological correlations and an equation for particle concentration transfer taking into account interfacial slip. The turbulence model takes into account the modulation of turbulence by particles. The proposed model has been validated on the problems of steady turbulent flow of shear thinning viscoplastic fluid with heavy particles in a horizontal pipe. The impact of Reynolds number and rheological parameters on the reliability of numerical simulations was examined. A comparison of experimental data with DNS-DEM simulation data has demonstrated that the proposed model is capable of accurately predicting the distribution of particle concentration, particle velocity, as well as carrier flow and pressure drop in the channel.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"304 ","pages":"Article 106878"},"PeriodicalIF":3.0,"publicationDate":"2025-10-16","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145360933","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Pub Date : 2025-10-16DOI: 10.1016/j.compfluid.2025.106874
Anna Schwarz , Daniel Kempf , Jens Keim , Patrick Kopper , Christian Rohde , Andrea Beck
High-order methods are well-suited for the numerical simulation of complex compressible turbulent flows, but require additional stabilization techniques to capture instabilities arising from the underlying non-linear hyperbolic equations. This paper provides a detailed comparison of the effectiveness of entropy stable discontinuous Galerkin methods for the stabilization of compressible (wall-bounded) turbulent flows. For this investigation, an entropy stable discontinuous Galerkin spectral element method is applied on Gauss–Legendre and Gauss–Lobatto nodes. In the compressible regime, an additional stabilization technique for shock capturing based on a convex blending of a low-order finite volume with the high-order discontinuous Galerkin operator is utilized. The present investigation provides a systematic study from convergence tests, to the Taylor–Green vortex and finally to a more intricate turbulent wall-bounded 3D diffuser flow, encompassing both weakly compressible and compressible regimes. The comparison demonstrates that the DGSEM on Gauss–Lobatto nodes is generally less accurate for an equal amount of degrees of freedom. Conversely, it is faster than the DGSEM on Gauss–Legendre nodes due to a less severe time step restriction and simpler numerical operator. A performance comparison reveals that the DGSEM on Gauss–Lobatto nodes generally outperforms the DGSEM on Gauss nodes for under-resolved turbulence in the subsonic regime on a periodic domain. Conversely, the opposite effect can be observed for wall-bounded flows as well as the supersonic regime, the latter depending of course on the chosen shock-capturing scheme. To the author’s knowledge, this is the first time for which a comparison of entropy stable DGSEM on Gauss–Lobatto and Gauss–Legendre has been performed for compressible, wall-bounded turbulent flows with separation.
{"title":"Comparison of entropy stable collocation high-order DG methods for compressible turbulent flows","authors":"Anna Schwarz , Daniel Kempf , Jens Keim , Patrick Kopper , Christian Rohde , Andrea Beck","doi":"10.1016/j.compfluid.2025.106874","DOIUrl":"10.1016/j.compfluid.2025.106874","url":null,"abstract":"<div><div>High-order methods are well-suited for the numerical simulation of complex compressible turbulent flows, but require additional stabilization techniques to capture instabilities arising from the underlying non-linear hyperbolic equations. This paper provides a detailed comparison of the effectiveness of entropy stable discontinuous Galerkin methods for the stabilization of compressible (wall-bounded) turbulent flows. For this investigation, an entropy stable discontinuous Galerkin spectral element method is applied on Gauss–Legendre and Gauss–Lobatto nodes. In the compressible regime, an additional stabilization technique for shock capturing based on a convex blending of a low-order finite volume with the high-order discontinuous Galerkin operator is utilized. The present investigation provides a systematic study from convergence tests, to the Taylor–Green vortex and finally to a more intricate turbulent wall-bounded 3D diffuser flow, encompassing both weakly compressible and compressible regimes. The comparison demonstrates that the DGSEM on Gauss–Lobatto nodes is generally less accurate for an equal amount of degrees of freedom. Conversely, it is faster than the DGSEM on Gauss–Legendre nodes due to a less severe time step restriction and simpler numerical operator. A performance comparison reveals that the DGSEM on Gauss–Lobatto nodes generally outperforms the DGSEM on Gauss nodes for under-resolved turbulence in the subsonic regime on a periodic domain. Conversely, the opposite effect can be observed for wall-bounded flows as well as the supersonic regime, the latter depending of course on the chosen shock-capturing scheme. To the author’s knowledge, this is the first time for which a comparison of entropy stable DGSEM on Gauss–Lobatto and Gauss–Legendre has been performed for compressible, wall-bounded turbulent flows with separation.</div></div>","PeriodicalId":287,"journal":{"name":"Computers & Fluids","volume":"303 ","pages":"Article 106874"},"PeriodicalIF":3.0,"publicationDate":"2025-10-16","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"145359265","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":3,"RegionCategory":"工程技术","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}