Third-order Runge-Kutta Method Research Articles

Numerical simulation of two-phase flows in 2-D petroleum reservoirs using a very high-order CPR method coupled to the MPFA-D finite volume scheme

The focus of the present work is to model the two-phase flow of oil and water in 2-D anisotropic petroleum reservoirs using a novel appropriate combination of a non-orthodox cell-centered Multipoint Flux Approximation with a Diamond Stencil (MPFA-D) finite volume method to discretize the pressure equation coupled to the very high-resolution Correction Procedure via Reconstruction (CPR) scheme for the discretization of the saturation equation. To suppress numerical oscillations (under/overshoots) near shocks, that are typical in higher-order schemes, and to deliver high accuracy in smooth regions of the solution, a hierarchical Multidimensional Limiting Process (MLP) is used in the reconstruction stage. The integration in time is carried out using a third-order Runge–Kutta method. The coupling of the pressure-saturation set of equations is carried out using a classical Implicit Pressure Explicit Saturation (IMPES) procedure. To properly couple the MPFA-D method with the CPR formulation, it is necessary to obtain an adequate velocity reconstruction throughout the control volumes of the mesh. Because the cell-centered finite volume method naturally delivers fluxes across control surfaces of the primal grid, then, we use a reconstruction operator based on the lowest order Raviart–Thomas interpolation functions and the Piola transformation, to get the complete knowledge of the velocity field throughout the domain. Finally, some representative problems were solved to evaluate the accuracy, efficiency, and shock-capturing capabilities of our new numerical methodology.

High-order maximum-principle-preserving and positivity-preserving weighted compact nonlinear schemes for hyperbolic conservation laws

In this paper, the maximum-principle-preserving (MPP) and positivity-preserving (PP) flux limiting technique will be generalized to a class of high-order weighted compact nonlinear schemes (WCNSs) for scalar conservation laws and the compressible Euler systems in both one and two dimensions. The main idea of the present method is to rewrite the scheme in a conservative form, and then define the local limiting parameters via case-by-case discussion. Smooth test problems are presented to demonstrate that the proposed MPP/PP WCNSs incorporating a third-order Runge-Kutta method can attain the desired order of accuracy. Other test problems with strong shocks and high pressure and density ratios are also conducted to testify the performance of the schemes.

A multi-material HLLC Riemann solver with both elastic and plastic waves for 1D elastic-plastic flows

A multi-material HLLC-type approximate Riemann solver with both elastic and plastic waves (MHLLCEP) is constructed for 1D elastic-plastic flows with the hypo-elastic model and the von Mises yielding condition. Although Cheng in 2016 [1] introduced a HLLC Riemann solver with elastic waves (HLLCE) for 1D elastic-plastic flows, Cheng assumed that pressure is continuous across a contact wave. This assumption may lead to numerical errors, especially for multi-material elastic-plastic flows. In our MHLLCEP, this assumption is not used again. Correspondingly, the errors introduced by the assumption are deleted, describing and evaluating the plastic waves are more accurate than that in the HLLCE. For the multi-material system, in this paper, a ghost cell method is used to eliminate the numerical oscillations and keep a high-order spatial reconstruction across the interface. Based on the MHLLCEP, combining with the third-order WENO reconstruction method and the third-order Runge–Kutta method in time, a high-order cell-centered Lagrangian scheme for 1D multi-material elastic-plastic flows is built in this paper. A number of numerical experiments are carried out. Numerical experiments show that the third-order scheme is robust, essentially non-oscillatory and, as suggested by numerical experiments, may also be convergent. Moreover, for multi-material elastic-plastic flows, the scheme with the MHLLCEP is more accurate and reasonable in resolving the multi-material interface than the scheme with the HLLCE.

Tsunami Intrusion and River Ice Movement

A two-dimensional wave model coupled with ice dynamics is developed to evaluate ice effects on shallow water wave propagation on a beach and in a channel. The nonlinear Boussinesq equations with ice effects are derived and solved by the hybrid technique of the Godunov-type finite volume method and finite difference method with the third-order Runge–Kutta method for time integration. The shock capturing method enables the model to simulate complex flows over irregular topography. The model is capable of simulating wave propagations accurately, including non-hydrostatic water pressure and wave dispersions. The ice dynamic module utilizes a Lagrangian discrete parcel method, based on smoothed particle hydrodynamics. The Boussinesq wave model is validated with an analytical solution of water surface oscillation in a parabolic container, an analytical solitary wave propagation in a flat channel, and experimental data on tsunami wave propagations. The validated model is then applied to investigate the interaction between ice and tsunami wave propagation, in terms of ice attenuation on tsunami wave propagations over a beach, ice deposition on the beach driven by the tsunami wave, and ice jam formation and release in a coastal channel with the intrusion of the tsunami wave. The simulated results demonstrated the interactions between tsunami waves and surface ice, including the maximum run up, ice movement along the beach, and ice jamming in a channel.

Symplectic simulation of dark solitons motion for nonlinear Schrödinger equation

In this paper, we study symplectic simulation of dark solitons motion of nonlinear Schrodinger equation (NLSE). The Ablowitz-Ladik model (A-L model) of NLSE can be expressed as a non-canonical Hamiltonian system. By using splitting technique, we construct explicit splitting K-symplectic methods for the A-L model. On the other hand, the A-L model can be transformed into a canonical system and standard symplectic methods can be employed to perform numerical simulation. A second order K-symplectic method and a second order symplectic method are employed to simulate one dark soliton and two dark solitons motion for the A-L model and its canonicalized system respectively. By comparing with a third-order non-symplectic Runge-Kutta method, we show the superiorities of the two symplectic methods in long-term tracking the motion of dark solitons and preserving the invariants. We also compare the CPU times of K-symplectic methods and standard symplectic methods and show that the former ones are more efficient. The energy-preserving scheme is also applied for non-canonical Hamiltonian systems. The numerical results demonstrate that the K-symplectic methods can nearly preserve the energy, the discrete invariants of A-L model and conserved quantities of NLSE, but the energy-preserving scheme can only exactly preserve the energy.

Constraint-consistent Runge–Kutta methods for one-dimensional incompressible multiphase flow

New time integration methods are proposed for simulating incompressible multiphase flow in pipelines described by the one-dimensional two-fluid model. The methodology is based on ‘half-explicit’ Runge–Kutta methods, being explicit for the mass and momentum equations and implicit for the volume constraint. These half-explicit methods are constraint-consistent, i.e., they satisfy the hidden constraints of the two-fluid model, namely the volumetric flow (incompressibility) constraint and the Poisson equation for the pressure. A novel analysis shows that these hidden constraints are present in the continuous, semi-discrete, and fully discrete equations.Next to constraint-consistency, the new methods are conservative: the original mass and momentum equations are solved, and the proper shock conditions are satisfied; efficient: the implicit constraint is rewritten into a pressure Poisson equation, and the time step for the explicit part is restricted by a CFL condition based on the convective wave speeds; and accurate: achieving high order temporal accuracy for all solution components (masses, velocities, and pressure). High-order accuracy is obtained by constructing a new third-order Runge–Kutta method that satisfies the additional order conditions arising from the presence of the constraint in combination with time-dependent boundary conditions.Two test cases (Kelvin–Helmholtz instabilities in a pipeline and liquid sloshing in a cylindrical tank) show that for time-independent boundary conditions the half-explicit formulation with a classic fourth-order Runge–Kutta method accurately integrates the two-fluid model equations in time while preserving all constraints. A third test case (ramp-up of gas production in a multiphase pipeline) shows that our new third-order method is preferred for cases featuring time-dependent boundary conditions.

Combined Compact Difference Scheme for Solving Unsteady Convection-Diffusion Equations

By using a fifth-order upwind discrete scheme for the convection term, a fourth-order central difference scheme for the diffusion term, we construct a boundary scheme that matches the interior point and obtain a stable combination semi-discrete compact scheme which is discretized by the third-order Runge-Kutta method in the time direction. Compared with the other fourth-order scheme, the error of the proposed method is small and the accuracy is high.

Trajectory errors of different numerical integration schemes diagnosed with the MPTRAC advection module driven by ECMWF operational analyses

Abstract. The accuracy of trajectory calculations performed by Lagrangian particle dispersion models (LPDMs) depends on various factors. The optimization of numerical integration schemes used to solve the trajectory equation helps to maximize the computational efficiency of large-scale LPDM simulations. We analyzed global truncation errors of six explicit integration schemes of the Runge–Kutta family, which we implemented in the Massive-Parallel Trajectory Calculations (MPTRAC) advection module. The simulations were driven by wind fields from operational analysis and forecasts of the European Centre for Medium-Range Weather Forecasts (ECMWF) at T1279L137 spatial resolution and 3 h temporal sampling. We defined separate test cases for 15 distinct regions of the atmosphere, covering the polar regions, the midlatitudes, and the tropics in the free troposphere, in the upper troposphere and lower stratosphere (UT/LS) region, and in the middle stratosphere. In total, more than 5000 different transport simulations were performed, covering the months of January, April, July, and October for the years 2014 and 2015. We quantified the accuracy of the trajectories by calculating transport deviations with respect to reference simulations using a fourth-order Runge–Kutta integration scheme with a sufficiently fine time step. Transport deviations were assessed with respect to error limits based on turbulent diffusion. Independent of the numerical scheme, the global truncation errors vary significantly between the different regions. Horizontal transport deviations in the stratosphere are typically an order of magnitude smaller compared with the free troposphere. We found that the truncation errors of the six numerical schemes fall into three distinct groups, which mostly depend on the numerical order of the scheme. Schemes of the same order differ little in accuracy, but some methods need less computational time, which gives them an advantage in efficiency. The selection of the integration scheme and the appropriate time step should possibly take into account the typical altitude ranges as well as the total length of the simulations to achieve the most efficient simulations. However, trying to summarize, we recommend the third-order Runge–Kutta method with a time step of 170 s or the midpoint scheme with a time step of 100 s for efficient simulations of up to 10 days of simulation time for the specific ECMWF high-resolution data set considered in this study. Purely stratospheric simulations can use significantly larger time steps of 800 and 1100 s for the midpoint scheme and the third-order Runge–Kutta method, respectively.

Nonstandard explicit third-order Runge-Kutta method with positivity property

When one solves differential equations, modeling physical phenomena, it is of great importance to take physical constraints into account. More precisely, numerical schemes have to be designed such that discrete solutions satisfy the same constraints as exact solutions. Based on general theory for positivity, with an explicit third-order Runge-Kutta method (we will refer to it as RK3 method) positivity is not ensured when applied to the inhomogeneous linear systems and the same result is regained on nonlinear positivity for this method. Here we mean by positivity that the nonnegativity of the components of the initial vector is preserved. Nonstandard finite differences (NSFDs) schemes can improve the accuracy and reduce computational costs of traditional finite difference schemes. In addition to NSFDs produce numerical solutions which also exhibit essential properties of solution. In this paper, we investigate the positivity property for nonstandard RK3 method when applied to the numerical solution of special nonlinear initial value problems (IVPs) for ordinary differential equations (ODEs). We obtain new results for positivity which are important in practical applications. We provide some numerical examples to illustrate our results.

Numerical Analysis on the Two-Dimensional Unsteady Magnetohydrodynamic Compressible Flow through a Porous Medium

In the present study, the unsteady magnetrohydrodynamic (MHD) flow of compressible fluid with variable thermal properties has been numerically investigated. The electrically conducting fluid flows through a porous media channel. The uniform magnetic field is applied perpendicular to the direction of the flow. The wall is assumed to be non-conducting and maintained at two different temperatures. The thermal conductivity and viscosity of the fluid change with temperature. Sixth - Order Accurate Compact Finite Difference scheme together with the Third-order Runge-Kutta method is used to solve a set of non-linear equations. The results of the calculation are expressed in the form of the velocity and temperature at different values of the magnetic field and porosity. The proposed mathematical model and numerical methods have been validated by comparing with the results of previously published studies that the compared results reveal the same trends. The difference is due to the compressibility and property variation effects. The results showed that the magnetic field and variable properties considerably influences the flows that is compressible thereby affecting the heat transfer as well as the wall shear stress.

Numerical simulation of 1-D oil and water displacements in petroleum reservoirs using the spectral finite volume method

The study and development of high-order discretization techniques for the modeling and simulation of multiphase flows in petroleum reservoirs are still a challenge from the computational viewpoint due to the difficulty posed by some physical features such as heterogeneity and anisotropy of the medium, that are of paramount importance in this class of applications. In this paper, we present a preliminary investigation of the very high order Spectral Finite Volume Method (SFVM) for the discretization of the saturation equation that describes the 1-D Oil-Water displacement through porous media. In this equation, the hyperbolic term is discretized using a high order SFVM, the capillary term is approximated using a second order CVFD (Control Volume Finite Difference) method and the integration in time is carried out using a third-order Runge–Kutta method. On the other hand, to prevent numerical oscillations near shocks, that are typical from higher order schemes, a slope limiter procedure is used in the reconstruction stage. To evaluate the performance of the SFVM some simple model problems, including capillarity and gravity, are solved, yielding accurate and efficient results. The use of very high order formulations can be of utmost importance, particularly for multidimensional and real large-scale applications with complex wave structures.

A finite-volume HLLC-based scheme for compressible interfacial flows with surface tension

Shock waves are often used in experiments to create a shear flow across liquid droplets to study secondary atomization. Similar behavior occurs inside of supersonic combustors (scramjets) under startup conditions, but it is challenging to study these conditions experimentally. In order to investigate this phenomenon further, a numerical approach is developed to simulate compressible multiphase flows under the effects of surface tension forces. The flow field is solved via the compressible multicomponent Euler equations (i.e., the five equation model) discretized with the finite volume method on a uniform Cartesian grid. The solver utilizes a total variation diminishing (TVD) third-order Runge–Kutta method for time-marching and second order TVD spatial reconstruction. Surface tension is incorporated using the Continuum Surface Force (CSF) model. Fluxes are upwinded with a modified Harten–Lax–van Leer Contact (HLLC) approximate Riemann solver. An interface compression scheme is employed to counter numerical diffusion of the interface. The present work includes modifications to both the HLLC solver and the interface compression scheme to account for capillary force terms and the associated pressure jump across the gas–liquid interface. A simple method for numerically computing the interface curvature is developed and an acoustic scaling of the surface tension coefficient is proposed for the non-dimensionalization of the model. The model captures the surface tension induced pressure jump exactly if the exact curvature is known and is further verified with an oscillating elliptical droplet and Mach 1.47 and 3 shock-droplet interaction problems. The general characteristics of secondary atomization at a range of Weber numbers are also captured in a series of simulations.

A high-order multi-zone cut-stencil method for numerical simulations of high-speed flows over complex geometries

In this paper, we present a method for performing uniformly high-order direct numerical simulations of high-speed flows over arbitrary geometries. The method was developed with the goal of simulating and studying the effects of complex isolated roughness elements on the stability of hypersonic boundary layers. The simulations are carried out on Cartesian grids with the geometries imposed by a third-order cut-stencil method. A fifth-order hybrid weighted essentially non-oscillatory scheme was implemented to capture any steep gradients in the flow created by the geometries and a third-order Runge–Kutta method is used for time advancement. A multi-zone refinement method was also utilized to provide extra resolution at locations with expected complex physics. The combination results in a globally fourth-order scheme in space and third order in time. Results confirming the method's high order of convergence are shown. Two-dimensional and three-dimensional test cases are presented and show good agreement with previous results. A simulation of Mach 3 flow over the logo of the Ubuntu Linux distribution is shown to demonstrate the method's capabilities for handling complex geometries. Results for Mach 6 wall-bounded flow over a three-dimensional cylindrical roughness element are also presented. The results demonstrate that the method is a promising tool for the study of hypersonic roughness-induced transition.

Development of DETAC and its application to the hydrogen detonation analysis

In this study, a numerical analysis code (DETAC, Detonation Analysis Code) for hydrogen detonation during the reactor severe accident was developed using Fortran 90 language, and the simulation was performed for the hydrogen detonation. A global-chemistry model was adopted to simulate the chemical reaction. The Euler equations were solved using third-order Runge-Kutta method with fifth-order weighted essentially non-oscillatory scheme handling the convection flux. Afterward, the hydrodynamics solver was verified by comparison of predicted results and exact solutions of four cases of shock tube problems. A hydrogen detonation in a pipe was simulated to verify this code by comparing the results with the classical C-J theory. Furthermore, this code was applied to the hydrogen detonation analysis in the compartment of BWR building. Two cases with different ignition locations were analyzed in this paper and the maximum pressure of these cases were 7.5 MPa and 8.0 MPa, respectively. The pressure and the temperature during detonation were affected by the ignition location. The results indicated that the possibility of reactor building destruction exists if the hydrogen detonation occurs.

Fourier Finite-Element Method with Linear Basis Functions on a Sphere: Application to Elliptic and Transport Equations

Abstract The Fourier finite-element method (FFEM) on the sphere, which performs with an operation count of O(N2 log2N) for 2N × N grids in spherical coordinates, was developed using linear basis functions. Dependent field variables are expanded with the Fourier series in the longitude, and the Fourier coefficients are represented with a series of first-order finite elements. Different types of pole conditions were incorporated into the Fourier coefficients of the scalar and vector variables in order to avoid discontinuity at the poles. For the Laplacian operator, the linear element was defined as a function of the sine of latitude instead of the latitude. The FFEM was applied to the derivatives of the first- and second-order elliptic equations and the transport equations. The scale-selective high-order Laplacian-type filter was implemented as a hyperviscosity. For the first-order derivative the fourth-order convergence rate of the accuracy, as is expected from the theoretical analysis, was achieved. Elliptic equations were found to be solved accurately without pole discontinuity, and the convergence rate turned out to be second order. The cosine bell advection, time-differenced with a third-order Runge–Kutta method, showed that the squared-norm error convergence rate was slightly above second order. Both the Gaussian bell advection and the deformational flow produced the theoretical convergence rate of fourth order. The high-order filter was found to be effective in maintaining a quasi-uniform resolution over the sphere, and thus allowed a large time step size. Sensitivity experiments of cosine bell advection over the poles revealed that the CFL number, as defined with the maximum grid size on the global domain, can be taken to be as large as unity.

Modelo matemático del transporte de una toxina en una red trófica marina

Algunos casos de intoxicación por consumo de peces contaminados, como la intoxicación por ciguatera, ocurren inesperadamente y no son fáciles de detectar previamente, dado que los peces que portan la toxina no presentan aspecto y/o síntomas de enfermos. En este trabajo proponemos un modelo matemático para el transporte y acumulación de una toxina a través de una cadena alimentaria. El modelo se plantea mediante un sistema rígido de ecuaciones diferenciales que describen la dinámica. Se analiza la estabilidad local de la solución de equilibrio. Se discuten diferentes escenarios de aparición de brotes de una toxina a partir de simulaciones numéricas obtenidas mediante un esquema de discretización que combina un método de Runge-Kutta de tercer orden y la regla del trapecio, evitando la rigidez del sistema. Los resultados muestran que el tiempo que tarda en desaparecer la toxina en la red trófica depende del estado en que se encuentra la dinámica poblacional al momento del brote. Esta información puede emplearse para establecer un tiempo de veda en la pesca de tal manera que la toxina se reduzca a niveles inocuos para la salud humana.

Linear Unsteady CFD of Vibrating Blades of Turbomachinery

This paper presents the numerical calculation of unsteady aerodynamic forces on the vibrating cascade blade row using the implicit time integration method. It is confirmed that first and second-order Euler methods, which are the implicit methods, require less computational time than the third-order Runge-Kutta method, which is the explicit method, by two orders of magnitude.

A coarse-constrained multiscale method for accelerating incompressible flow computations

We present a coarse-constrained multiscale (CCM) method for accelerating incompressible flow computations. Reducing the number of degrees of freedom of the Poisson solver by powers of two in the primitive variable fractional-step method, or the vorticity-stream function formulation of the problem accelerates these computations while, for the first level of coarsening, retaining the same level of accuracy in the fine-resolution velocity field variables. CCM is a modular approach that facilitates data transfer with simple interpolations and uses black-box solvers for the Poisson and advection-diffusion equations in the flow solver. Here, we implement a particular CCM method, CCMRK3, that uses the full weighting operation for mapping from fine to coarse grids, the third-order Runge-Kutta method for time stepping, and finite differences for the spatial discretization. After solving the Poisson equation on a coarsened grid, bilinear interpolation is used to obtain the fine data for consequent time stepping on the full grid. We compute several benchmark flows: the Taylor-Green vortex, a vortex pair merging, a double shear layer, decaying turbulence, the Taylor-Green vortex on a distorted grid, and laminar flow over a circular cylinder. In all cases we use either FFT-based or V-cycle multigrid linear-cost Poisson solvers. We favor the FFT-based Poisson solver for Cartesian grid applications due to its robustness. A linear acceleration rate is obtained for all the cases we consider, with reduction factors in computational time between 2 and 42. We also find that the computational savings increases with increasing distortion ratio on non-Cartesian grids, which makes the CCM method a useful tool for generalized curvilinear incompressible flow solvers.

Relationship between vortex and sound in under-expanded supersonic impinging jets

SUMMARY The study of an under-expanded supersonic jet impinging on a flat plate by using large-eddy simulation is reported. A third-order upwind compact difference and a fourth-order symmetric compact scheme are employed to discretize the nondimensional axisymmetric compressible Favre-filtered Navier–Stokes equations in space, whereas the third-order Runge–Kutta method with the total variation diminishing property is adopted to deal with the temporal discretization. The numerical simulation successfully captures the shock wave and vortex structures with different scales in the flow field. Waves with high and low frequencies traveling forward and reflecting back, and sound sources in different locations can be observed. By comparison with the frequency of the impinging tone from the experiment, it can be deduced that the change of pressure and swirling strength in the shear layer, pressure change on the impinging plate, and vortex merging in the jet shear layer are interdependent with the impinging tone. The effects of nozzle lip thickness on the impinging jet flow field have been investigated. The results show that the values of pressure fluctuation and vortex swirling strength in the shear layer near the nozzle have an extremum with the variation of the nozzle lip thickness. The results provide a theoretical foundation for the design of supersonic nozzles. Copyright © 2011 John Wiley & Sons, Ltd.

Strong-stability-preserving 3-stage Hermite–Birkhoff time-discretization methods

Strong-stability-preserving (SSP) time-discretization methods have a nonlinear stability property that makes them particularly suitable for the integration of hyperbolic conservation laws. A collection of SSP explicit 3-stage Hermite–Birkhoff methods of orders 3 to 7 with nonnegative coefficients are constructed as k-step analogues of third-order Runge–Kutta methods, incorporating a function evaluation at two off-step points. Generally, these new methods have larger effective CFL coefficients than the hybrid methods of Huang with the same step number k. They have larger maximum scaled step sizes than hybrid methods on Burgers' equations.