High-resolution Numerical Methods Research Articles

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

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.

Intermittency of gravity wave turbulence on the surface of an infinitely deep fluid: Numerical experiment

The intermittency of gravity wave turbulence on the surface of an infinitely deep fluid has been studied numerically. We estimated the scaling exponent of surface elevation structure-functions directly from the equations of motion, e.g., by numerical solution of the underlying equations of motion, without any additional closure assumptions (e.g., phase randomisation) and reduction to the kinetic equations as in the conventional framework of wave turbulence theory. With our high-resolution numerical methods for accurately modelling the nonlinear water surface, we were able to evaluate the scaling exponents of the structure-function up to relatively high orders (p=15), which has not been previously reported. We investigated the effect of wave steepness and forcing on intermittency and compared our results with other analytical and numerical studies, including results on the intermittency of wave turbulence of a different nature (turbulence of bending waves on a plate and magnetohydrodynamic turbulence). Our results support the conjecture of universality for intermittency phenomena in wave turbulence.

Consistent high resolution interface-capturing finite volume method for compressible multi-material flows

Compressible multi-material flows are characterized by complex flow structures with a broad range of length scales and discontinuities associated with material interfaces and shock waves. High order and high resolution numerical methods are required to capture material interface as sharply as possible and to increase the resolution of complex structures along the material interface. This paper will present the application of the high resolution finite volume (FV) method based on the minimized dispersion and controllable dissipation (MDCD) reconstruction to compressible multi-material flows. The MDCD scheme with independent dispersion and dissipation provides a flexible mechanism to control the numerical dissipation. The adjustment of dissipation of the MDCD reconstruction will not affect the consistency, which is required for interface-diffusion capturing methods to prevent spurious oscillations. Several one- and two dimensional multi-material numerical simulations have indicated that the high resolution FV method based on MDCD reconstruction can capture material interfaces free of spurious oscillations and provide better resolved small-scale features than other numerical schemes on the same grid resolution.

High-resolution numerical investigation into the effects of winglet on the aerodynamic performance for a three-dimensional vertical axis wind turbine

Vertical axis wind turbines can harvest kinetic energy in the air stream from all directions and therefore have distinct advantages in the effective utilization of wind resource. Its aerodynamic performance highly depends on the blade geometrical configuration, especially on the tip shapes due to the negative effects of the shedding vortices. In this paper, spatially high-resolution numerical methods are adopted to study the evolution and mutual interaction process for the vortex system formed behind a vertical axis wind turbine with and without blade tip winglet, aiming at further understand the mechanism of the winglet effects on the turbine aerodynamic performance. For the cases without tip winglet, it has been shown that the high-resolution methods can properly capture the evolution and interaction of the omega-shaped vortices system and revealed the differences of the vortex mutual induction in the near wake region under three tip speed ratios. Then, the region in spanwise dominated by the effects of tip vortex is determined by proper orthogonal decomposition survey on the flow field of vertical profiles and based on which the winglet is installed hoping to maximize its three-dimensional effects on the spanwise flow field. The comparison of power coefficients for cases with and without winglet under the same swept area show that the winglet has significantly improved the effectiveness of energy harvesting in the studied tip speed ratios. Finally, the analysis of tip vortex circulation evolution demonstrates that the enhancement of power coefficient for the winglet case comes from the reduction of the spanwise region affected by tip vortex, instead of the strength decrease of the tip vortex itself. The conclusions can further reveal the flow mechanism of winglet effects on the vertical axis wind turbine aerodynamic performance, and therefore are helpful for guiding the optimized blade winglet design to achieve higher efficiency.

Investigating tool engagement in groundwood pulping: finite element modelling and in-situ observations at the microscale

Abstract With industrial groundwood pulping processes relying on carefully designed grit surfaces being developed for commercial use, it is increasingly important to understand the mechanisms occurring in the contact between wood and tool. We present a methodology to experimentally and numerically analyse the effect of different tool geometries on the groundwood pulping defibration process. Using a combination of high-resolution experimental and numerical methods, including finite element (FE) models, digital volume correlation (DVC) of synchrotron radiation-based X-ray computed tomography (CT) of initial grinding and lab-scale grinding experiments, this paper aims to study such mechanisms. Three different asperity geometries were studied in FE simulations and in grinding of wood from Norway spruce. We found a good correlation between strains obtained from FE models and strains calculated using DVC from stacks of CT images of initial grinding. We also correlate the strains obtained from numerical models to the integrity of the separated fibres in lab-scale grinding experiments. In conclusion, we found that, by modifying the asperity geometries, it is, to some extent, possible to control the underlying mechanisms, enabling development of better tools in terms of efficiency, quality of the fibres and stability of the groundwood pulping process.

Phenomenology of bubble-collapse-driven penetration of biomaterial-surrogate liquid-liquid interfaces

The paper presents phenomenology of interaction and penetration of liquid-liquid material interfaces initiated by shock-driven collapse of single and multiple microbubbles situated near the material interface. Previous experimental studies have established such a generic setting as relevant for the investigation of sonoporation, i.e., the perforation of live cells by microbubble collapses. We consider a planar or spherical, single- or dual-layer, material interface between a gelatin material and water. A single or several ideal-gas microbubbles are positioned near the interface. Bubble collapse is initiated by a shock wave with a pressure profile specific to laser generation and is flat when hitting the gas-water interface. The interfacial acoustic impedance match singles out the collapse-induced re-entrant jet as main event. High-resolution sharp-interface numerical methods are employed to ensure that wave dynamics, hydrodynamics, and interface transporting are accurately resolved. Bubble configurations are varied between single and double and between attached and with standoff distance. Parameters varied are shock-wave peak pressure and viscosity ratio between single and double layers of gelatin and the surrounding water. For inertia-dominated cases, two regimes are observed, the first characterized by linear growth of the penetration depth and the second by a t(2/3) scaling. The latter range is affected by viscosity which reduces penetration speed. The results show that process parameters, in particular shock overpressure, control not only penetration depth but also the size of the interface perforation, indicating means to steer processes in biomedical applications.

Simulation of turbulent flow in MATIS-H rod bundle with split-type mixing vanes

Turbulent flow simulations of 5 × 5 rod bundle with the split-type mixing grid are presented. The simulations were performed by the Jožef Stefan Institute (JSI) in the frame of its participation in the OECD/NEA CFD Benchmark exercise. The experimental data from the MATIS-H benchmark experiment located at the Korean Atomic Energy Research Institute (KAERI) were used for posteriori comparison with simulation results. To predict the flow between the rods and mixing vanes, a complex, fully hexagonal flow-oriented mesh was constructed, taking into account the geometry of mixing vanes. High resolution numerical methods and a scale-resolving turbulence model was used in ANSYS CFX simulation. Time-averaged velocity profiles, profiles of velocity fluctuations and integral velocity circulation at the specified downstream locations show good agreement with experimental data. In addition, the simulated flow patterns over the entire rod bundle cross-section are presented and discussed.

Simulation of antisolvent crystallization in impinging jets with coupled multiphase flow-micromixing-PBE

The KT (Kurganov and Tadmor, 2000) finite-volume central scheme is one of the most promising high-resolution numerical methods to solve the widely-used population balance equation (PBE) in crystallization and other areas. To meet the practical purpose of geometric-type particle size grid, the primary KT scheme was extended to a general form and validated for pure growth in homogeneous systems. Based on the extended KT scheme, a solver was developed that couples the general discretized PBE and amicromixing and a CFD mixture model in OpenFOAM (open-source field operation and manipulation). The simulation uses published parameters from and is compared to experimental antisolvent crystallization of lovastatin from a methanol-water mixture in an impinging jet. The effect of the existing solid crystals on some crystal properties is investigated in this work for the first time. The shapes of crystal size distribution (CSD) at various jet velocities are consistent with experimental observations. The geometric particle size partition is shown to be capable of improving the accuracy of simulation in divisions of the highest particle number densities or steep gradients in the number density. The existing solid crystals are also shown to have a non-negligible effect on the slurry flow crystallization systems once the mean crystal size reaches 20µm.

Immiscible thermo-viscous fingering in Hele-Shaw cells

We investigate immiscible radial displacement in a Hele-Shaw cell with a temperature dependent viscosity using two coupled high resolution numerical methods. Thermal gradients created in the domain through the injection of a low viscosity fluid at a different temperature to the resident high viscosity fluid can lead to the formation of unstable thermo-viscous fingers, which we explore in the context of immiscible flows. The transient, multi-zone heat transfer is evaluated using a newly developed auxiliary radial basis function-finite collocation (RBF-FC) method, which locally captures variation in flux and field variable over the moving interface, without the need for ghost node extrapolation. The viscosity couples the transient heat transfer to the Darcy pressure/velocity field, which is solved using a boundary element - RBF-FC method, providing an accurate and robust interface tracking scheme for the full thermo-viscous problem. We explore the thermo-viscous problem space using systematic numerical experiments, revealing that the early stage finger growth is controlled by the pressure gradient induced by the varying temperature and mobility field. In hot injection regimes, negative temperature gradients normal to the interface act to accelerate the interface, promoting finger bifurcation and enhancing the viscous fingering instability. Correspondingly, cold injection regimes stabilise the flow compared to isothermal cases, hindering finger formation. The interfacial mobility distribution controls the late stage bifurcation mode, with non-uniformities induced by the thermal diffusivity creating alternate bifurcation modes. Further numerical experiments reveal the neutral stability of the thermal effects on the fingering evolution, with classical viscous fingering dynamics eventually dominating the evolution. We conclude the paper with a mechanistic summary of the immiscible thermo-viscous fingering regime, providing the first detailed analysis of the thermal problem in immiscible flows.

MPI to Coarray Fortran: Experiences with a CFD Solver for Unstructured Meshes

High-resolution numerical methods and unstructured meshes are required in many applications of Computational Fluid Dynamics (CFD). These methods are quite computationally expensive and hence benefit from being parallelized. Message Passing Interface (MPI) has been utilized traditionally as a parallelization strategy. However, the inherent complexity of MPI contributes further to the existing complexity of the CFD scientific codes. The Partitioned Global Address Space (PGAS) parallelization paradigm was introduced in an attempt to improve the clarity of the parallel implementation. We present our experiences of converting an unstructured high-resolution compressible Navier-Stokes CFD solver from MPI to PGAS Coarray Fortran. We present the challenges, methodology, and performance measurements of our approach using Coarray Fortran. With the Cray compiler, we observe Coarray Fortran as a viable alternative to MPI. We are hopeful that Intel and open-source implementations could be utilized in the future.

Maintaining interface equilibrium conditions in compressible multiphase flows using interface capturing

An accurate treatment of material interfaces in compressible multiphase flows poses important challenges for high-resolution numerical methods. Although high-order interface-capturing schemes have been used to accurately simulate gas/liquid interfaces with the Euler equations, these methods can result in temperature spikes at material discontinuities. While this phenomenon is not problematic for Euler simulations, it gives rise to numerical errors when heat conduction is included. In this work, we identify the source of these errors and propose a methodology to prevent their occurrence for various models used to represent gas/liquid interfaces in compressible flows based on a “single-fluid” formulation, in which interfaces are represented by discontinuities in the material properties. Our focus lies in materials (gases and liquids primarily, but also solids) that can be described by a stiffened equation of state, though our approach is generalizable to other equations. We show that numerical approaches that prevent pressure oscillations at interfaces may generate temperature errors, which affect the energy (and pressure) through the heat conduction term. We demonstrate that the material properties entering the equation of state must be computed according to suitable transport equations in conservative or non-conservative forms; the pressure and temperature must be calculated based on the appropriate properties. To verify the analysis and compute problems with gas/liquid interfaces of relevance, we develop a three-dimensional, high-order accurate, solution-adaptive finite difference framework. In particular, we show that temperatures and pressures may be significantly overestimated in calculations of shock-induced bubble collapse in water if temperature errors are not prevented.

Analysis of finite-difference schemes based on exact and approximate solution of Riemann problem

The Riemann problem of one-dimensional arbitrary discontinuity breakdown for parameters of unsteady gas flow is considered as applied to the design of Godunov-type numerical methods. The problem is solved in exact and approximate statements (Osher-Solomon difference scheme used in shock capturing numerical methods): the intensities (the ratio of static pressures) and flow velocities on the sides of the resulting breakdowns and waves are determined, and then the other parameters are calculated in all regions of the flow. Comparison of calculation results for model flows by exact and approximate solutions is performed. The concept of velocity function is introduced. The dependence of the velocity function on the breakdown intensity is investigated. A special intensity at which isentropic wave creates the same flow rate as the shock wave is discovered. In the vicinity of this singular intensity approximate methods provide the highest accuracy. The domain of applicability for the approximate Osher-Solomon solution is defined by performing test calculations. The results are presented in a form suitable for usage in the numerical methods. The results obtained can be used in the high-resolution numerical methods.

Implementation and assessment of high-resolution numerical methods in TRACE

The 1st-order upwind differencing numerical scheme is widely employed to discretize the convective terms of the two-phase flow transport equations in reactor systems analysis codes such as TRACE and RELAP. While very robust and efficient, 1st-order upwinding leads to excessive numerical diffusion. Standard 2nd-order numerical methods (e.g., Lax–Wendroff and Beam–Warming) can effectively reduce numerical diffusion but often produce spurious oscillations for steep gradients. To overcome the difficulties with the standard higher-order schemes, high-resolution schemes such as nonlinear flux limiters have been developed and successfully applied in numerical simulation of fluid-flow problems in recent years. The present work contains a detailed study on the implementation and assessment of six nonlinear flux limiters in TRACE. These flux limiters selected are MUSCL, Van Leer (VL), OSPRE, Van Albada (VA), ENO, and Van Albada 2 (VA2). The assessment is focused on numerical stability, convergence, and accuracy of the flux limiters and their applicability for boiling water reactor (BWR) stability analysis. It is found that VA and MUSCL work best among of the six flux limiters. Both of them not only have better numerical accuracy than the 1st-order upwind scheme but also preserve great robustness and efficiency.

Simulations of the transient flow generated from a started flat plate

Transient operations are commonly founded in fluid machineries such as the starting, stopping, and variations of rotor speeds, etc. Flow generated from a started flat plate is of fundamental importance. Experiments have been done to observe the flow evolution in current researches. And in order to explore the flow in more detailed scale, some vortex methods with high resolution and other numerical methods were developed to solve various related problems by some researchers. But the promotion of vortex method to engineering application is rare due to its complexity and difficulty in specifying the boundary conditions. In order to build up a method of numerical study for such problems, a simplified model is built up with a flat plate. The development of two-dimensional viscous incompressible flow generated from an impulsively started and uniformly accelerated infinitesimally thin flat plate is simulated numerically. A dynamic mesh(DM) method based on the spring analogue and local remeshing is applied to realize the mesh motion caused by the started plate. Researches show that the mesh quality will decline under large grid shear force during the updating process. To conquer this problem, a region near the plate is separated to guarantee the mesh quality at location of interest which is the innovation of the present paper. All computations at least cover a period during which the plate translates 6 times its length. The simulated instantaneous velocity profiles, flow structures and drag coefficients under several Reynolds numbers (20⩽Re⩽126) and accelerations (20 m/s2⩽a⩽152 m/s2) are presented and compared with existing results in literatures. Comparisons are found to be satisfactory, confirming the validity of the current proposed method(region separated DM). The proposed DM method is firstly used to study the transient flow generated from a started flat plate and can be used in further study of transient characteristics during transient operations of turbo machineries.

Human impact and landscape utilization from the Mesolithic to medieval time traced by high spatial resolution pollen analysis and numerical methods

Pollen analysis in relation to archaeological excavations has been carried out by the coast of central Norway to investigate for the first time long-term human impact on the vegetation in this region. A total of 297 samples from eleven archaeological sites and connected bogs/peat profiles reveal the vegetation development from the Mesolithic until today. Principal Component Analysis (PCA) is used to identify and elucidate the development through time and space, and the results of pollen diagrams from bogs are compared to the results of on-site pollen data. Human impact has transformed the vegetation to a varying degree in different time periods. A dense population connected to the shoreline is slightly visible as opening-up of the forest in the Mesolithic, whereas re-use of Mesolithic settlement sites for cultivation and grazing in the Late Neolithic and Early Bronze Age (from c. 2300 cal BC) is clearly documented by palaeobotany but weakly recorded in the archaeological record. Areas for settlement and cultivation are found from the Late Bronze and Iron Ages and pollen analysis contributes to a detailed picture of utilization of the landscape. The investigation shows the potential of combining archaeological and botanical data and analysing several small pollen sequences within an excavation area to reveal spatial patterns in vegetation development and human impact through time.

Numerical simulation of cloud rise phenomena associated with nuclear bursts

We present numerical simulations of cloud evolution from nuclear explosions using high-resolution numerical methods. Our numerical approach includes a fluid mechanical model that is a combination of a compressible code (GEODYN) and a low Mach code (LMC). Early stages of nuclear explosions that are characterized by the blust wave propagation are simulated with an explicit code (GEODYN) that solves the compressible Navier–Stokes equations via a high-order Godunov scheme. As soon as the blust wave weakens (≈10 s) the subsequent cloud rise due to buoyancy forces can be effectively simulated by the LMC code. LMC is an implicit code based on a pressure projection technique, and derived from the compressible Navier–Stokes equations using an asymptotic analysis in Mach number. It analytically eliminates time step restrictions based on sound wave propagation and it is computationally efficient and accurate for simulations of cloud rise dynamics at later stages. We perform a series of cloud rise scenarios ranging from an idealized bubble rise problem to realistic air bursts. We analyze effects of compressible dynamics on cloud evolution at different stages. It is found that in a realistic configuration, interaction of the reflected shock wave from the ground with the fireball significantly affects cloud evolution, in contrast to the equivalent idealized bubble rise simulations. We validate the code predictions against available experimental data. It is demonstrated that, by providing the initial source from the compressible GEODYN code, late time flow evolution can be successfully simulated with the fast, efficient and accurate LMC code.

A Finite Element Solution of Structure-Borne Sound Attenuation for a Lightweight Timber Floor

Structure-borne sound attenuation of lightweight timber floor structures has been investigated by high-resolution numerical calculation methods. Two models have been studied both experimentally and numerically. The first considers a continuous chipboard plate and a supporting beam. The second considers two plates that are fastened on one beam. Damping effects, which arise from the two different discontinuities, are included. The real-valued modulus of elasticity is replaced with a complex value. In the calculations, the orthotropic material characteristics of the building material have also been taken into account. The results from the calculations and laboratory measurement confirm that the attenuation rate is highly dependent on the excitation location, if all the other parameters remain the same.

Simulation of transition and turbulence decay in the Taylor–Green vortex

Conventional large-eddy simulation (LES) and monotone integrated LES (MILES) are tested in emulating the dynamics of transition to turbulence in the Taylor–Green vortex (TGV). A variety of subgrid scale (SGS) models and high-resolution numerical methods are implemented in the framework of both incompressible and compressible fluid flow equations. Comparisons of the evolution of characteristic TGV integral measures are made with previously reported and new direct numerical simulation (DNS) data. The computations demonstrate that the convective numerical diffusion effects in the MILES methods can consistently capture the physics of flow transition and turbulence decay without resorting to an explicit SGS model, while providing accurate prediction of established theoretical findings for the kinetic energy dissipation, energy spectra, enstrophy and kinetic energy decay. All approaches tested provided fairly robust computational frameworks.

Inertial Motions during the Transient Adjustment of a Density Anomaly in the Equatorial Ocean with Application to the Western Pacific Warm Pool

Abstract This paper is focused on the spontaneous transient adjustment of a buoyant lens of water with uniform density, initially at rest in the vicinity of the equator. For parameters typical of the western Pacific warm pool, the adjustment is shown to generate finite-amplitude wave motions with period ∼8 days, which are not covered by the standard theory of linear equatorial waves. This mechanism may be at the origin of inertial motions at the early stages of ENSO events in the western Pacific Ocean. The lens adjustment is studied within the 11/2-layer reduced-gravity approximation on the equatorial β plane, using the high-resolution finite-volume numerical methods that are specially designed to handle outcropping isopycnals. Under the reduced-gravity approximation, a buoyant region of light water with outcropping boundaries in the vicinity of the equator is described by two parameters: the meridional-to-zonal scale aspect ratio δ and the ratio γ of the Coriolis force to the pressure force on its meridional boundary. For realistic parameters (δ ∼ 10−1; γ ∼ 1), the lens, initially at rest, spreads eastward in accord with nonrotating gravity current dynamics, whereas its westward extrusion is arrested so that the western edge splits into two anticyclonic vortices. Meanwhile finite-amplitude westward-propagating inertial wave motions develop at the interface between the spreading current and the ambient fluid. The inertial wave structure is shown to be consistent with the structure of stable wave modes predicted by linear analysis of small amplitude perturbations superimposed on a zonally symmetric equatorial current with outcropping isopycnals. A Wentzel–Kramers–Brillouin–Jeffreys ray-tracing analysis indicates that the inertial wave is emitted during the early stage of the gravity current evolution and then dispersed on the spreading current.

Numerical simulations of the interaction of accretion disks with young star magnetospheres

Aims. We present the results of the numerical simulations of the interaction between a magnetized star and an imperfectly conducting accretion disk. The star is rotating with constant angular velocity. The differentially rotating Keplerian disk is treated as a boundary condition. We are interested in the magnetic field topology dependence on the electrical conductivity of the disk. Methods. To analyze the “star–disk” interaction we numerically investigate the MHD equations using Godunov-type high resolution numerical methods. Results. It was found that in our model the “star–disk” interaction occurs with a quasi-periodic reconnection of the magnetic field coronal loops and plasmoid ejections. In the case of the perfect disk conductivity, the evolution of the coronal magnetic field leads to the periodic outflow of angular momentum from the disk. In the case of an imperfectly conducting disk, the configuration of the magnetic field is formed such that the disk angular momentum carried by the magnetic field is balanced by angular momentum carried by matter. It should be noted that we used the ideal MHD equation to obtain the solutions. The reconnection process in the disk corona depends on the numerical diffusivity that exists in our numerical code. Our simulations treat reconnection as occurring in current sheets. The thickness of the current sheet is broadened by numerical resistivity. Nevertheless, we suppose that the reconnection and plasmoid ejection takes place as well for real magnetic diffusivity. To verify the method and results we also used several more detailed grids to estimate the numerical diffusivity of the scheme. It is turned out that the setup model presented in the paper quite reasonable satisfies the goal of this paper, i.e., to investigate the regime of interaction between the magnetized star and the disk.