Additional Source Term Research Articles

A dissipative extension to ideal hydrodynamics

Abstract We present a formulation of special relativistic, dissipative hydrodynamics (SRDHD) derived from the well-established Müller-Israel-Stewart (MIS) formalism using an expansion in deviations from ideal behaviour. By re-summing the non-ideal terms, our approach extends the Euler equations of motion for an ideal fluid through a series of additional source terms that capture the effects of bulk viscosity, shear viscosity and heat flux. For efficiency these additional terms are built from purely spatial derivatives of the primitive fluid variables. The series expansion is parametrized by the dissipation strength and timescale coefficients, and is therefore rapidly convergent near the ideal limit. We show, using numerical simulations, that our model reproduces the dissipative fluid behaviour of other formulations. As our formulation is designed to avoid the numerical stiffness issues that arise in the traditional MIS formalism for fast relaxation timescales, it is roughly an order of magnitude faster than standard methods near the ideal limit.

Induced gravitational waves: the effect of first order tensor perturbations

Scalar induced gravitational waves contribute to the cosmological gravitational wave background. They can be related to the primordial density power spectrum produced towards the end of inflation and therefore are a convenient new tool to constrain models of inflation. These waves are sourced by terms quadratic in perturbations and hence appear at second order in cosmological perturbation theory. While the focus of research so far was on purely scalar source terms we also study the effect of including first order tensor perturbations as an additional source. This gives rise to two additional source terms: a term quadratic in the tensor perturbations and a cross term involving mixed scalar and tensor perturbations. We present full analytical expressions for the spectral density of these new source terms and discuss their general behaviour. To illustrate the generation mechanism we study two toy models containing a peak on small scales. For these models we show that the scalar-tensor contribution becomes non-negligible compared to the scalar-scalar contribution on smaller scales. We also consider implications for future gravitational wave surveys.

Uniformly high-order bound-preserving OEDG schemes for two-phase flows

This paper proposes and rigorously analyzes novel high-order bound-preserving oscillation-eliminating discontinuous Galerkin (BP-OEDG) schemes with the Harten–Lax–van Leer (HLL) numerical flux for the Kapila five-equation two-phase flow model. The evolution equation for the volume fraction is formulated as a conservative advection equation with an additional non-conservative source term. Utilizing a technical splitting approach, we design a uniformly high-order discretization for this source term, incorporating an “upwind” discretization of the non-conservative product at cell interfaces based on the HLL wave speeds. This method inherently ensures the symmetry of the volume fractions, which is crucial for establishing the bound-preserving property. The positivity of partial densities and internal energy is rigorously proven using the geometric quasilinearization (GQL) approach, which transforms the nonlinear pressure positivity constraint into equivalent linear constraints. To suppress potential spurious oscillations, we incorporate a scale-invariant and linearity-invariant oscillation elimination (OE) procedure that damps the DG modal coefficients after each Runge–Kutta stage, as proposed in [M. Peng, Z. Sun and K. Wu, OEDG: Oscillation-eliminating discontinuous Galerkin method for hyperbolic conservation laws, Math. Comput. (2024), https://doi.org/10.1090/mcom/3998]. This OE procedure, acting as a post-processing filter based on the jumps of the DG solution at cell interfaces, is easy to implement, maintains the Abgrall equilibrium condition around an isolated material interface, and preserves the high-order accuracy of the DG schemes. The effectiveness and robustness of the proposed high-order BP-OEDG schemes are demonstrated through several benchmark numerical experiments.

Computational fluid dynamics simulation of cryogenic vertical upflow boiling under Earth gravity

This study presents numerical simulations and their validation for flow boiling of liquid nitrogen (LN2) in a vertical upflow orientation, with a primary aim to understand the complex two-phase flow and heat transfer phenomena important to space applications. The computational fluid dynamics (CFD) model utilized the coupled level set volume of fluid (CLSVOF) method, incorporating additional source terms for bubble collision dispersion force and shear lift force in the momentum conservation equation to enhance simulation accuracy. The simulations were conducted for two mass velocities (G = 526 and 804 kg/m2s) and three different heat flux levels (approximately 10%, 30%, and 70% of critical heat flux (CHF) under Earth gravity. The model was validated against measured wall temperature data acquired from the authors’ previous experimental studies, demonstrating average deviations of less than 2.8 K across all operating conditions. The simulated two-phase flow contours illustrated various flow patterns, including bubbly, slug, churn, and annular. Both mass velocity and heat flux were observed to impact the onset of nucleate boiling (ONB), bubble nucleation, growth, and coalescence, and overall vapor structure. The simulations also offered insight into axial and radial void fraction and velocity profiles, revealing local flow acceleration trends synchronized with void fraction development. A comparison between predicted and measured bulk fluid temperature profiles showed excellent agreement, further validating the CFD model’s accuracy and practical usefulness for two-phase cryogenic flow boiling simulations in space applications.

Suppression of blowup by slightly superlinear degradation in a parabolic–elliptic Keller–Segel system with signal-dependent motility

In this paper, we consider an initial–Neumann boundary value problem for a parabolic–elliptic Keller–Segel system with signal-dependent motility and a source term. Previous research has rigorously shown that the source-free version of this system exhibits an infinite-time blowup phenomenon when dimension N≥2. In the current work, when N≤3, we establish uniform boundedness of global classical solutions with an additional source term that involves slightly super-linear degradation effect on the density, of a maximum growth order slogs, unveiling a sufficient blowup suppression mechanism. The motility function considered in our work takes a rather general form compared with recent works (Fujie and Jiang, 2020; Lyu and Wang, 2023) which were restricted to the monotone non-increasing case. The cornerstone of our proof lies in deriving an upper bound for the second component of the system and an entropy-like estimate, which are achieved through tricky comparison skills and energy methods, respectively.

Tilt effect on squeeze-film damping in electrostatic micro-actuators

In this study, we theoretically explore the impact of tilt on the switching dynamics of a MEMS device operated by a pair of electrostatic actuators and subject to pronounced squeeze-film damping in the transition and free molecular flow regimes. Our investigation reveals that inducing tilt, achieved through the application of uneven voltages on the electrodes, introduces additional source terms in the modified Reynolds equation governing the pressure distribution in the squeeze-film gap. These added terms result in an accelerated decay of the squeeze film force, facilitating faster device switching between the initial and target configurations.To analyze the device dynamics, we develop a coupled numerical model and employ evolutionary multi-objective optimization to determine an optimal actuation scheme. This optimization is carried out for scenarios where initially different voltages are applied to the electrodes during a certain pre-defined time interval before transitioning to the steady voltage corresponding to the target gap size at equilibrium. Our findings highlight that the suggested actuation strategy leads to substantial improvements in the switching time of the device. Furthermore, we provide insights into the conditions favoring tilt actuation over parallel-motion actuation where the initially applied voltages are equal.

Spherically symmetric electrovac spacetimes in conformal Killing gravity

Recently Harada proposed a third-order gravitational theory for the derivatives of the metric. Mantica and Molinari showed that Harada’s theory may be recast into the form of Einstein’s field equations (EFEs) with an additional source term which is a second-order conformal Killing tensor. Accordingly they named the theory conformal Killing gravity. However, they overlooked the fact that all solutions of the new theory (except those satisfying EFEs) admit a non-trivial second-order Killing tensor. Harada derived an analog of the Schwarzschild solution. Recently Tarciso et al obtained a generalization of Harada’s vacuum solution analogous to the Reissner–Nordström solution. However, like Harada they assumed a restricted form for a static spherically symmetric metric. In this study the most general spherically symmetric static vacuum solution of Harada’s theory and its generalization with a Maxwell electromagnetic field as source were obtained. The validity of Birkhoff’s theorem for static spherically symmetric electrovac fields in conformal Killing gravity is investigated.

GQL-based bound-preserving and locally divergence-free central discontinuous Galerkin schemes for relativistic magnetohydrodynamics

This paper develops novel and robust central discontinuous Galerkin (CDG) schemes of arbitrarily high-order accuracy for special relativistic magnetohydrodynamics (RMHD) with a general equation of state (EOS). These schemes are provably bound-preserving (BP), meaning they consistently preserve the upper bound for subluminal fluid velocity and the positivity of density and pressure, while also (locally) maintaining the divergence-free (DF) constraint for the magnetic field. For 1D RMHD, the standard CDG method is exactly DF, and its BP property is proven under a condition achievable by the BP limiter. For 2D RMHD, we design provably BP and locally DF CDG schemes based on the suitable discretization of a modified RMHD system, which is the relativistic analogue of Godunov's symmetrizable form of the non-relativistic MHD system (Godunov, 1972 [19]). A key novelty in our schemes is the meticulous discretization of additional source terms in the modified RMHD equations, so as to precisely counteract the influence of divergence errors on the BP property across overlapping meshes. Notably, we provide rigorous proofs of the BP property for our CDG schemes and first establish the theoretical connection between BP and discrete DF properties on overlapping meshes for RMHD. Owing to the absence of explicit expressions for primitive variables in terms of conserved variables, the constraints of physical bounds are strongly nonlinear, making the BP proofs highly nontrivial. We overcome these challenges through technical estimates within the geometric quasilinearization (GQL) framework (Wu and Shu, 2023 [49]), which equivalently converts the nonlinear constraints into linear ones. Furthermore, we introduce a new 2D cell average decomposition on overlapping meshes, which relaxes the theoretical BP CFL constraint and reduces the number of internal nodes, thereby enhancing the efficiency of the 2D BP CDG method. Finally, we implement the proposed CDG schemes for extensive RMHD problems with various EOSs, demonstrating their robustness and effectiveness in challenging scenarios like ultra-relativistic blasts and jets in strongly magnetized environments.

A consistent methodology to transport a passive scalar with the geometric Volume-of-Fluid method isoAdvector

This paper presents a new methodology allowing the discretisation of phase specific transport equations within the Volume-of-Fluid interface capturing method framework. The method uses a sharp interface algorithm to compute the transport of species concentration. The interface reconstruction and advection is provided by the geometric advection scheme isoAdvector[1] implemented in the OpenFOAM® library. When discretising the transport equation for species, lack of consistency with the free surface advection scheme can lead to numerical errors, causing conservation or boundedness issues. This work addresses the issue of consistency in convective transport of species and is divided in two parts. First, a new interpolation procedure is used to compute face values from cell-centered values. Then, the diffusive operator of the transport equation is corrected. Finally, a set of test cases are presented to validate the transport equation's consistency with the free surface advection. Species transport across the interface is not part of the scope of this article, however, this methodology can further be used to study mass transfer at the gas-liquid interface using additional mass source terms that are not discussed here.

Fundamental solution of the time-space bi-fractional diffusion equation with a kinetic source term for anomalous transport

The purpose of this paper is to study the fundamental solution of the time-space bi-fractional diffusion equation incorporating an additional kinetic source term in semi-infinite space. The equation is a generalization of the integer-order model ∂tρ(x,t)=∂x2ρ(x,t)-ρ(x,t)\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${\\partial _{{t}} {\\rho (x,t)}} = {\\partial ^2_{{x}} {\\rho (x,t)}} - { \\rho (x,t)}$$\\end{document} (also known as the Debye–Falkenhagen equation) by replacing the first-order time derivative with the Caputo fractional derivative of order 0<α<1\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$0<\\alpha < 1$$\\end{document}, and the second-order space derivative with the Riesz-Feller fractional derivative of order 0<β<2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$0< \\beta <2$$\\end{document}. Using the Laplace-Fourier transforms method, it is shown that the parametric solutions are expressed in terms of the Fox’s H-function that we evaluate for different values of α\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\alpha$$\\end{document} and β\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\beta$$\\end{document}.

CFD investigation on start-up transient of HPR1000 secondary side passive residual heat removal heat exchanger

In advanced pressurized water reactor HPR1000, the secondary side passive residual heat removal heat exchanger (SPRS HX) is a crucial component. It removes the decay heat using boiling and convection heat transfer in the event of accidents. The robust functioning of the SPRS HX significantly influences overall reactor safety. To investigate the thermal hydraulic attributes of SPRS HX, computational fluid dynamics (CFD) is used in this study. The tube region is simulated using the porous media method (PMM). To model the two-phase flow in the tank, the drift flux model (DFM) is employed. The assessment of heat transfer rate from the primary to the secondary sides is conducted by empirical correlations. The governing equations are addressed through the utilization of FLUENT. Specifically, the momentum and energy equations are augmented with additional source terms that account for flow resistance and heat transfer in the tube region. These modifications are enacted through User Defined Functions (UDF). Three-dimensional profiles of the void fraction, fluid temperature and velocity are acquired, allowing for a comprehensive analysis of the heat transfer attributes of tube region.

Models of quasi-discontinuous solar-wind streams

Context. The modeling of the solar-wind outflow patterns is addressed in terms of local transient distortions of the flow, temperature, and density profiles due to the presence of local energy sources. A recently introduced related new class of analytically derived quasi-discontinuous solar-wind solutions is numerically approached. Aims. The analytical discontinuous solutions can asymptotically obtained from steady-state and time-dependent models in the limit of very localized external heating. The aim of the current study is to develop a numerical confirmation for the presence of quasi-discontinuous distortions of the wind profiles by mimicking the local energy sources with additional source terms in the governing equations of the numerical models. Methods. Corresponding systems of ordinary and partial differential equations, respectively, are formulated employing prescribed heating functions. After a comparison of sequences of numerically obtained steady-state solutions with the analytical one, the stability of the former is tested with a time-dependent simulation. Results. The analytical discontinuous solutions are asymptotically reproduced with the quasi-discontinuous steady-state and time-dependent numerical solutions in the limit of vanishingly small width (compared to the other characteristic length scales of the system) of the heating function. Conclusions. The interpretation that such solutions result from strongly localized heating has been confirmed both qualitatively and quantitatively. The applied numerical approach enables the building of more complex, multidimensional counterpart models and local profiles of typical local energy sources that are presumably responsible for the dynamical properties of the solar-wind patterns found.

Building a reduced order model from CFD data on leaves to evaluate optimal climate conditions

Climate control in vertical farms or greenhouses is crucial in order to have the optimal climate for the plants. This climate can be simulated using Computational Fluid Dynamics (CFD). This way, different ventilation methods can for example be compared. Previously, plants were modeled as a porous zone, in which additional source terms are added to the energy or transpiration transport equation. In this study, ten different leaves under different orientations are studied. The heat and mass balance for each leaf is solved and the boundary layer is adequately resolved. Different leaf orientations and positions to the flow are used, so that each leaf has its own distinct boundary layer. Leaf temperature and humidity are a function of the incoming flow, but also of the incoming temperature, humidity and stomatal resistance. In order to be able to quickly assess the effects of changing boundary conditions, on the temperature and humidity of the leaves, a Reduced Order Model (ROM) is made from the CFD-data. This method allows to quickly assess the effects of a stomatal resistance, the incoming radiation or a changing inlet temperature or humidity. A ROM was built using the data from one simulation. In order to see whether the ROM is able to predict other configurations, two additional cases were solved in CFD and were compared to the ROM. Leaf temperatures and humidities only differed maximally 3.14% and were quite good at predicting this. Predicted leaf transpiration rates differed slightly more, as these differed maximally 10.5%. Overall, the ROM has proven to be a valuable and numerically cost-effective tool to predict overall leaf climate and can be a good tool when later the climate of a full plant is simulated.

GO-MELT: GPU-optimized multilevel execution of LPBF thermal simulations

Computationally modeling the laser powder bed fusion process can be expensive despite the fact that the laser covers a small fraction of the entire domain. Hence, we developed GO-MELT, a multilevel approach inspired by the variational multiscale method. GO-MELT calculates the full temperature solution using three coupled thermal solvers, each solving one of three overlapping domains, subsequently referred to as levels. Level 1’s solver obtains a global, coarse-scale solution. Level 2’s solver computes a meso-scale thermal field spanning the melt pool region. Level 3’s solver obtains a fine-scale temperature field spanning the laser’s immediate region. Coarser solutions provide boundary conditions for finer solvers. Finer solutions compute subgrid scales that influence coarser solvers through additional source terms. Being independently meshed, Levels 2 and 3 can track a laser along its tool path without remeshing. Moreover, fixed-sized, structured meshes allow for GPU acceleration using Google’s JAX library with just-in-time compilation. Five case studies were conducted to demonstrate proper convergence of GO-MELT, to show that GO-MELT can reach the same accuracy as a uniform mesh with fewer degrees of freedom, and to quantify the computational load from each level’s solver. A simulated production run averaged 1.64ms per time step after taking 9.36h to complete 20.5 million time steps, which is approximately 678× faster than a uniform mesh solver with identical resolution and GPU acceleration. Future work will implement different solvers into GO-MELT to improve fidelity and speed.

Revisiting RANS turbulence modelling for bubble-induced turbulence: Effects of surfactants

Recent experiments of dispersed bubbly flow (Ma et al., 2023) show that contamination of the carrier phase plays a considerable role in bubble-induced turbulence (BIT): although the bubbles rise slower with adding surfactants, the BIT increases. This trend, however, cannot be reproduced by any current BIT models in the framework of Euler–Euler Reynolds-averaged Navier–Stokes two-equation turbulence model, including the model proposed by Ma et al. (2017). The present work revisits this widely used model (Ma et al., 2017) as a representative, which considers the BIT with an additional source term in the k- and ɛ (or ω) - equation, respectively. We have systematically investigated the reason why such a BIT model is unable to capture the BIT modulation caused by different degrees of contamination and discuss the new possibility for the prefactor in the source term of the ɛ-equation that may improve the results.

Adaptive Turbulence Model for Leading Edge Vortex Flows Preconditioned by a Hybrid Neural Network

Eddy-viscosity-based turbulence models provide the most commonly used modeling approach for computational fluid dynamics simulations in the aerospace industry. These models are very accurate at a relatively low cost for many cases but lack accuracy in the case of highly rotational leading edge vortex flows for mid to low aspect-ratio wings. An enhanced adaptive turbulence model based on the one-equation Spalart–Allmaras turbulence model is fundamental to this work. This model employs several additional coefficients and source terms, specifically targeting vortex-dominated flow regions, where these coefficients can be calibrated by an optimization procedure based on experimental or high-fidelity numerical data. To extend the usability of the model from single or cluster-wise calibrated cases, this work presents a preconditioning approach of the turbulence model via a neural network. The neural network provides a case-unspecific calibration approach, enabling the use of the model for many known or unknown cases. This extension enables aircraft design teams to perform low-cost Reynolds-averaged Navier–Stokes simulations with increased accuracy instead of complex and costly high-fidelity simulations.

Application of reactor dosimetry techniques for source term validation in F-18 production with a medical cyclotron

Imaging techniques such as positron emission tomography need radionuclides that can be produced with cyclotrons. The involved nuclear reactions mostly produce neutron and gamma radiation, which must be shielded. In order to be able to determine the required thickness of the shielding for a newly bought medical cyclotron, corresponding radiation protection calculations were carried out. The necessary neutron source term was supplied by the manufacturer. To verify this source term, additional source terms were calculated with the MCNP6 and the FLUKA code packages. The results with the source, on base of code internal nuclear models, produced comparable results, but the neutron yield with the source supplied by the manufacturer turned out to be lower by a factor of about 5. To validate the results, experiments were carried out on an already existing cyclotron. Neutron fluence was measured with standard monitors which are used in reactor dosimetry. The experiments were performed during regular 18F production. Activities of the nuclides were measured by gamma spectroscopy and compared with the calculated activities.

Study on the effects of manifold structure on the gas flow distribution uniformity of anode of PEMFC stack with 140-cell

The gas maldistribution of the large proton exchange membrane fuel cell (PEMFC) stack leads to the deterioration of power performance. This numerical investigation studies the gas flow rate distribution in anode of two U-type stacks (a 15-cell stack and a 140-cell stack) which are composed of the same type of unit fuel cells. The anode gas channels of the two stacks are regarded as porous medium whose characteristic parameters are determined from pressure drop test data. To consider the nitrogen permeation from cathode to anode gas channel during the purging process, an expression of an additional momentum source term is proposed which helps to better fit the simulation results with pressure drop test data. Results show that the anode of the 140-cell stack has a quite good flow distribution uniformity. However, its asymmetric feed header of the manifold structure still has some negative effect. The same 140-cell stack installed with a rotational symmetric feed header is investigated and the simulation result shows that it can reduce the original maximum non-uniformity by about 50%.

Dispersive hydrodynamics in non-Hermitian nonlinear Schrödinger equation with complex external potential

In this paper dispersive hydrodynamics associated with the non-Hermitian nonlinear Schrödinger (NLS) equation with generic complex external potential is studied. In particular, a set of dispersive hydrodynamic equations are obtained. They differ from their classical counterparts (without an external potential), by the presence of additional source terms that alter the density and momentum equations. When restricted to a class of Wadati-type complex potentials, the resulting hydrodynamic system conserves a modified momentum and admits constant intensity/density solutions. This motivates the construction and study of an initial value problem (IVP) comprised of a centred (or non-centred) step-like initial condition that connects two constant intensity/density states. Interestingly, this IVP is shown to be related to a Riemann problem posed for the hydrodynamic system in an appropriate traveling reference frame. The study of such IVPs allows one to interpret the underlying non-Hermitian Riemann problem in terms of an ‘optical flow’ over an obstacle. A broad class of non-Hermitian potentials that lead to modulationally stable constant intensity states are identified. They are subsequently used to numerically solve the associated Riemann problem for various initial conditions. Due to the lack of translation symmetry, the resulting long-time dynamics show a dependence on the location of the step relative to the potential. This is in sharp contrast to the NLS case without potential, where the dynamics are independent of the step location. This fact leads to the formation of diverse nonlinear wave patterns that are otherwise absent. In particular, various gain-loss generated near-field features are present, which in turn drive the optical flow in the far-field which could be comprised of various rich nonlinear wave structures, including DSW-DSW, DSW-rarefaction, and soliton-DSW interactions.

Modeling of Wind Turbine Interactions and Wind Farm Losses Using the Velocity-Dependent Actuator Disc Model

This paper analyzes the interaction of wind turbines and losses in wind farms using computational fluid dynamics (CFD). The mathematical model used consisted of three-dimensional Reynolds-averaged Navier–Stokes (RANS) equations, while the presence of wind turbines in the flow was simulated as additional source terms. The novelty of the research is the definition of the source term as a velocity-dependent actuator disc model (ADM). This allowed for modeling the operation of a wind farm consisting of real wind turbines, characterized by power coefficients Cp and thrust force coefficients CT, which are a function of atmospheric wind speed. The calculations presented used a real 5 MW Gamesa turbine. Two different turbine spacings, 5D and 10D, where D is the diameter of the turbine, and two different locations corresponding to the offshore and onshore conditions were examined. The proposed model can be used to analyze wind farm losses not only in terms of the geometric distribution of individual turbines but also in terms of a specific type of wind turbine and in the entire wind speed spectrum.