publications

2024

1. Coupling finite elements of class C1 on composite curved meshes for second order elliptic problems

   Ashish Bhole, Hervé Guillard, Boniface Nkonga, and 1 more author

   International Journal for Numerical Methods in Fluids, 2024

   Abs

   Summary Finite elements of class 𝒞1 are suitable for the computation of magnetohydrodynamics instabilities in tokamak plasmas. In addition, isoparametric approximations allow for a precise alignment of the mesh with the magnetic field line. Mesh alignment is crucial to achieve axisymmetric equilibria accurately. It is also helpful to deal with the anisotropy nature of magnetized plasma flows. In this numerical framework, several practical simulations are now available. They help to understand better the operation of existing devices and predict the optimal strategies for using the international ITER tokamak under construction. However, a mesh-aligned isoparametric representation suffers from the presence of critical points of the magnetic field (magnetic axis, X-point). We here explore a strategy that combines aligned mesh out of the critical points with non-aligned unstructured mesh in a region containing these points. By this strategy, we can avoid highly stretched elements and the numerical difficulties that come with them. The mesh-aligned interpolation uses bi-cubic Hemite-Bézier polynomials on a structured mesh of curved quadrangular elements. On the other hand, we assume reduced cubic Hsieh-Clough-Tocher finite elements on an unstructured triangular mesh. Both meshes overlap, and the resulting formulation is a coupled discrete problem solved iteratively by a suitable one-level Schwarz algorithm. In this paper, we will focus on the Poisson problem on a two-dimensional bounded regular domain. This elliptic equation is a simplified version of the axisymmetric tokamak equilibrium one at the asymptotic limit of infinite major radius (large aspect ratio).

2023

1. Stabilized bi-cubic Hermite Bézier finite element method with application to gas-plasma interactions occurring during massive material injection in tokamaks

   Ashish Bhole, Boniface Nkonga, José Costa, and 3 more authors

   Computers & Mathematics with Applications, 2023

   Abs HTML

   Development of a numerical tool based upon the high-order, high-resolution Galerkin finite element method (FEM) often encounters two challenges: First, the Galerkin FEMs give central approximations to the differential operators and their use in the simulation of the convection-dominated flows may lead to the dispersion errors yielding entirely wrong numerical solutions. Secondly, high-order, high-resolution numerical methods are known to produce high wave-number oscillations in the vicinity of shocks/discontinuities in the numerical solution adversely affecting the stability of the method. We present the stabilized finite element method for plasma fluid models to address the two challenges. The numerical stabilization is based on two strategies: Variational Multiscale (VMS) and the shock-capturing approach. The former strategy takes into account (the approximation of) the effect of the unresolved scales onto resolved scales to introduce upwinding in the Galerkin FEM. The latter adaptively adds the artificial viscosity only in the vicinity of shocks. These numerical stabilization strategies are applied to stabilize the bi-cubic Hermite Bézier FEM in the computational framework of the nonlinear magnetohydrodynamics (MHD) code JOREK. The application of the stabilized FEM to the challenging simulation of Shattered Pellet Injection (SPI) in JET-like plasma is presented. It is shown that the developed numerical stabilization model improves the stability of the underlying numerical algorithm and the computational cost required to reveal the complex physics is reduced. The physical and numerical models presented can be used to perform expensive simulations of the plasma applications in large computational domains such as JET, and ITER.

2022

1. Treatment of polar grid singularities in the bi-cubic Hermite-Bézier approximations: Isoparametric finite element framework

   Ashish Bhole, Boniface Nkonga, Stanislas Pamela, and 2 more authors

   Journal of Computational Physics, 2022

   Abs HTML

   We propose a numerical treatment for the geometric singularity at the polar grid center encountered in the application of the isoparametric bi-cubic Hermite Bézier finite element method. The treatment applies a set of new basis functions at the polar grid center in the numerical algorithm where the new basis functions are simply the linear transformations of the original basis functions. The linear transformation comes out naturally by analyzing the interpolation formula at the polar grid center. The proposed polar treatment enforces the C1 regularity in the physical space and preserves the order of the accuracy of the interpolation. The treatment is applied in the nonlinear MHD code JOREK. With the help of a range of numerical tests, it is demonstrated that the polar treatment improves the stability and accuracy of the numerical approximation near the polar grid center. The polar treatment presented can be applied to the grid center of circular or non-circular polar grids and is also applicable for the bi-cubic Hermite finite element method.

2020

1. Extended full-MHD simulation of non-linear instabilities in tokamak plasmas

   S. J. P. Pamela, A. Bhole, G. T. A. Huijsmans, and 5 more authors

   Physics of Plasmas, Oct 2020

   Abs HTML

   Non-linear magnetohydrodynamic (MHD) simulations play an essential role in active research and understanding of tokamak plasmas for the realization of a fusion power plant. The development of MHD codes such as JOREK is a key aspect of this research effort. In this paper, we present an operational version of the full-MHD model implemented in JOREK, a significant advancement from the reduced-MHD model used for previous studies, where assumptions were made on the perpendicular dynamics and the toroidal magnetic field. The final model is presented in detail, and benchmarks are performed using both linear and non-linear simulations, including comparisons between the new full-MHD model of JOREK and the previously extensively studied reduced-MHD model, as well as results from the linear full-MHD code CASTOR3D. For the cases presented, this new JOREK full-MHD model is numerically and physically reliable, even without the use of numerical stabilization methods. Non-linear modeling results of typical tokamak instabilities are presented, including disruption and edge-localized-mode physics, most relevant to current open issues concerning future tokamaks such as ITER and DEMO.

2. A path conservative finite volume method for a shear shallow water model

   Praveen Chandrashekar, Boniface Nkonga, Asha Kumari Meena, and 1 more author

   Journal of Computational Physics, Oct 2020

   Abs HTML

   The shear shallow water model provides an approximation for shallow water flows by including the effect of vertical shear in the model. This model can be derived from the depth averaging process by including the second order velocity fluctuations which are neglected in the classical shallow water approximation. The resulting model has a non-conservative structure which resembles the 10-moment equations from gas dynamics. This structure facilitates the development of path conservative schemes and we construct HLL, 3-wave and 5-wave HLLC-type solvers. An explicit and semi-implicit MUSCL-Hancock type second order scheme is proposed for the time integration. Several test cases including roll waves show the performance of the proposed modeling and numerical strategy.

2019

1. A discontinuous Galerkin method for a two dimensional reduced resistive MHD model

   Praveen Chandrashekar, Boniface Nkonga, and Ashish Bhole

   Computers & Fluids, Oct 2019

   Abs HTML

   We are concerned with the numerical approximation of an incompressible ionized gas (plasma) flowing in a toroidal geometry. We also assume that the flow is independent of the toroidal coordinate and the resulting model is thus 2-D. We consider a symmetric formulation, the so-called Reduced Resistive MHD model, where the governing equation gives the evolution of the axial current and vorticity (scalar variables). These equations are written in a quasi conservative form and, using the Discontinuous Galerkin (DG) framework, an approximation strategy is proposed and analyzed for triangular meshes. This approach combines a Galerkin projection of the velocity stream function and magnetic flux to obtain divergence-free approximations, together with a DG approximation of the evolutionary equations for current and vorticity, while the integration is performed using Crank–Nicholson scheme. The designed methodology is validated on the kink-mode and the tilting MHD instabilities.

2. Fluctuation splitting Riemann solver for a non-conservative modeling of shear shallow water flow

   Ashish Bhole, Boniface Nkonga, Sergey Gavrilyuk, and 1 more author

   Journal of Computational Physics, Oct 2019

   Abs HTML

   In this paper we propose a fluctuation splitting finite volume scheme for a non-conservative modeling of shear shallow water flow (SSWF). This model was originally proposed by Teshukov (2007) in [14] and was extended to include modeling of friction by Gavrilyuk et al. (2018) in [7]. The directional splitting scheme proposed by Gavrilyuk et al. (2018) in [7] is tricky to apply on unstructured grids. Our scheme is based on the physical splitting in which we separate the characteristic waves of the model to form two different hyperbolic sub-systems. The fluctuations associated with each sub-systems are computed by developing Riemann solvers for these sub-systems in a local coordinate system. These fluctuations enable us to develop a Godunov-type scheme that can be easily applied on mixed/unstructured grids. While the equation of energy conservation is solved along with the SSWF model in [7], in this paper we solve only SSWF model equations. We develop a cell-centered finite volume code to validate the proposed scheme with the help of some numerical tests. As expected, the scheme shows first order convergence. The numerical simulation of 1D roll waves shows a good agreement with the experimental results. The numerical simulations of 2D roll waves show similar transverse wave structures as observed in [7].

2016

1. Roles of bulk viscosity on Rayleigh-Taylor instability: Non-equilibrium thermodynamics due to spatio-temporal pressure fronts

   Tapan K. Sengupta, Aditi Sengupta, Nidhi Sharma, and 3 more authors

   Physics of Fluids, Sep 2016

   Abs HTML

   Direct numerical simulations of Rayleigh-Taylor instability (RTI) between two air masses with a temperature difference of 70 K is presented using compressible Navier-Stokes formulation in a non-equilibrium thermodynamic framework. The two-dimensional flow is studied in an isolated box with non-periodic walls in both vertical and horizontal directions. The non-conducting interface separating the two air masses is impulsively removed at t = 0 (depicting a heaviside function). No external perturbation has been used at the interface to instigate the instability at the onset. Computations have been carried out for rectangular and square cross sections. The formulation is free of Boussinesq approximation commonly used in many Navier-Stokes formulations for RTI. Effect of Stokes’ hypothesis is quantified, by using models from acoustic attenuation measurement for the second coefficient of viscosity from two experiments. Effects of Stokes’ hypothesis on growth of mixing layer and evolution of total entropy for the Rayleigh-Taylor system are reported. The initial rate of growth is observed to be independent of Stokes’ hypothesis and the geometry of the box. Following this stage, growth rate is dependent on the geometry of the box and is sensitive to the model used. As a consequence of compressible formulation, we capture pressure wave-packets with associated reflection and rarefaction from the non-periodic walls. The pattern and frequency of reflections of pressure waves noted specifically at the initial stages are reflected in entropy variation of the system.

2. Non-equilibrium Thermodynamics of Rayleigh-Taylor instability

   Tapan K Sengupta, Aditi Sengupta, K S Shruti, and 2 more authors

   Journal of Physics: Conference Series, Oct 2016

   Abs

   Rayleigh-Taylor instability (RTI) has been studied here as a non-equilibrium thermodynamics problem. Air masses with temperature difference of 70K, initially with heavier air resting on lighter air isolated by a partition, are allowed to mix by impulsively removing the partition. This results in interface instabilities, which are traced here by solving two dimensional (2D) compressible Navier-Stokes equation (NSE), without using Boussinesq approximation (BA henceforth). The non-periodic isolated system is studied by solving NSE by high accuracy, dispersion relation preserving (DRP) numerical methods described in Sengupta T.K.: High Accuracy Computing Method (Camb. Univ. Press, USA, 2013). The instability onset is due to misaligned pressure and density gradients and is evident via creation and evolution of spikes and bubbles (when lighter fluid penetrates heavier fluid and vice versa, associated with pressure waves). Assumptions inherent in compressible formulation are: (i) Stokes’ hypothesis that uses zero bulk viscosity assumption and (ii) the equation of state for perfect gas which is a consequence of equilibrium thermodynamics. Present computations for a non-equilibrium thermodynamic process do not show monotonic rise of entropy with time, as one expects from equilibrium thermodynamics. This is investigated with respect to the thought-experiment. First, we replace Stokes’ hypothesis, with another approach where non-zero bulk viscosity of air is taken from an experiment. Entropy of the isolated system is traced, with and without the use of Stokes’ hypothesis. Without Stokes’ hypothesis, one notes the rate of increase in entropy to be higher as compared to results with Stokes’ hypothesis. We show this using the total entropy production for the thermodynamically isolated system. The entropy increase from the zero datum is due to mixing in general; punctuated by fluctuating entropy due to creation of compression and rarefaction fronts originating at the interface and reflecting from the walls.

2014

1. Error dynamics of diffusion equation: Effects of numerical diffusion and dispersive diffusion

   Tapan K. Sengupta, and Ashish Bhole

   Journal of Computational Physics, Oct 2014

   Abs HTML

   An error propagation equation has been obtained for the numerical solution of an initial-boundary value problem governed by linear diffusion equation. The error propagation equation is analyzed to identify different types of errors as those due to sub- or super-diffusion and dispersive diffusion errors. Quantification of errors has also been used for numerical analysis of some central difference methods to solve this equation.

2013

1. Direct numerical simulation of 2D transonic flows around airfoils

   Tapan K. Sengupta, Ashish Bhole, and N.A. Sreejith

   Computers & Fluids, Oct 2013

   Abs HTML

   High-accuracy, time-accurate compressible Navier–Stokes solvers have been developed for transonic flows. These solvers use optimized upwind compact schemes (OUCS) and four-stage, fourth order explicit Runge–Kutta (RK4) time integration scheme, details of which can be obtained in Sengupta [Sengupta TK. High Accuracy Computing Methods, fluid flows and wave phenomena. UK: Cambridge Univ. Press; 2013]. Although these compact schemes have been developed originally for direct simulation of incompressible flows, it is shown here that the same can be used for compressible flows, with shock-boundary layer interactions clearly captured for flow past NACA 0012 and NLF airfoils. Numerical higher order diffusion terms which are used for incompressible flows, have been replaced here by the pressure-based artificial diffusions proposed by Jameson et al. [Jameson A, Schmidt W, Turkel E. Numerical solution of the Euler equations by finite volume methods using Ruge–Kutta time stepping schemes. AIAA Paper 1981-1259. AIAA 14th fluid and plasma dynamics conference. Palo Alto, CA; 1981]. Such second and fourth order diffusion terms are used adaptively at selective points, located by the pressure switch. Developed computational methods used here are validated for cases with and without shocks, for which experimental results are available. Apart from surface pressure coefficient, contours of physical quantities are presented to explain the time-accurate results. Presented methods are robust and the results can be gainfully used to study shock formation, drag divergence and buffet onset of flow over airfoils.