First-order Upwind Research Articles

Comparisons of experimentally and numerically determined statistics for predicting low-occurrence wind speeds around a 1:1:2 block model

Accurate prediction of low-occurrence wind speeds around urban structures is crucial for safe building design. Although Large-eddy simulation (LES) is commonly used as a high-fidelity model as compared with the Reynolds-Averaged Navier–Stokes (RANS) simulations, the present validation process relies on the comparison of fundamental statistics of the mean and standard deviations. The discrepancies in LESs and wind-tunnel experiments (WTEs) are unclear in terms of physical quantities characterizing the unsteadiness of the simulated turbulent flow such as probability density and power spectral densities, and low-occurrence winds speeds. Therefore, this study aims to evaluate the effectiveness of LES in predicting unsteady wind behavior around a 1:1:2 block model. The study identifies prominent differences to improve the accuracy of unsteady numerical simulations especially for the purpose of predicting low-occurrence wind speeds. Various advection schemes in LESs were investigated, including first-order upwind, second-order linear, and dynamic interpolation schemes. The results show significant discrepancies, particularly in higher-order statistics and low-occurrence wind speeds, with WTE consistently exhibiting higher energy levels across all frequencies. These findings highlight the need to refine advection schemes to enhance their predictive accuracy. LESs with minimal numerical errors from discretization schemes can substantially improve urban wind assessments and contribute to the design of safer structures.

A novel multi-fidelity optimization framework for high-altitude propellers

A novel optimization methodology is presented for optimizing high-altitude propellers operating in the lower stratosphere. The proposed methodology combines a Bayesian optimization approach for most design variables and Vortex Theory for the optimization of twist distribution and rotational velocity. For the Bayesian optimization part of the problem, a multi-fidelity approach is employed with three levels of modeling fidelity. The considered fidelity levels, from lower to higher, are: Vortex Theory, and 3D Reynolds-Averaged Navier-Stokes (RANS) with the use of γ−Reθ transition model, converged with first-order upwind, and second-order upwind for the momentum equations. Additionally, the optimization problem involves two objectives: maximizing the cruise efficiency and minimizing the total volume of the propeller. To enhance the performance of multi-fidelity, multi-objective (MFMO) acquisition functions, a novel algorithm is proposed. The proposed optimization methodology demonstrates its effectiveness in generating several high-performing propeller designs and proves to be more efficient than an equivalent problem formulation relying solely on Bayesian optimization. Post-processing of optimization data reveals new design variable bounds that correspond to more efficient blade shapes. These findings provide insights into some important geometric features of high-altitude propellers and establish guidelines for future optimization efforts within the specified thrust and power consumption ranges. Finally, simulations of the best-performing designs across a range of advance ratios and altitudes confirm their high performance in various operating conditions.

Simulation of two-phase flow in 3D fractured reservoirs using a projection-based Embedded Discrete Fracture Model on Unstructured tetrahedral grids (pEDFM-U)

The modeling of fluid flow in heterogeneous, anisotropic and fractured porous media is relevant in many applications, including hydrocarbon and groundwater extraction, dispersion of contaminants, hydrogen or carbon dioxide (CO2) storage. Thus, accurate and scalable simulation of fluid flow through these formations continues to be a great challenge. The presence of fractures that are explicitly modeled, ranging from flow barriers to highly conductive channels, significantly increases the complexity of the numerical simulation. The Embedded Discrete Fracture Model (EDFM) produces results that can be as accurate as those obtained using equidimensional Discrete Fracture Models (DFM) at a much lower computational cost for highly conductive fractures, but it is not adequate for impermeable barriers. In contrast, the pEDFM (projection-based Embedded Discrete Fracture Model) produces accurate results for both, channels and barriers by using additional matrix-fracture connectivities. In its current versions, it is restricted to cartesian and corner-point grid geometries and to the classical Two Point Flux Approximation (TPFA) scheme. In this work, for the first time, the pEDFM is extended to handle unstructured tetrahedral meshes (pEDFM-U). In addition, interface fluxes are approximated by using the Multipoint Flux Approximation method with a Diamond stencil (MPFA-D). This allows for simulation of highly heterogeneous and anisotropic geo-models. The developed method is robust and can deal with full permeability tensors on arbitrary tetrahedral meshes. The advective terms are discretized considering an implicit First Order Upwind (FOU) method. Our method is implemented within the DARSim (Delft Advanced Reservoir Simulation) open-source simulator framework. Through several test cases, the proposed method showed to be robust and capable to accurately capture the effects of both high and low permeability fractures for general tetrahedral meshes, under arbitrary heterogeneous and anisotropic permeability tensors for the rock matrix.

A very high order Flux Reconstruction (FR) method for the numerical simulation of 1-D compositional fluid flow model in petroleum reservoirs

Compositional reservoir simulation is an extremely important tool for modeling the fluid flow in complex petroleum reservoirs, as in cases where the reservoir fluid is volatile and composed of several pseudo-components with distinct characteristics, or reservoirs that require the use of enhanced oil recovery process. For these cases, simpler black-oil models are not suitable. The compositional model involves the solution of a large system of partial differential equations, that comprises the mass conservation, the Darcy’s law and the fugacity constraints. However, the high number of equations and constraints render the compositional problems extremely complex and computationally demanding solutions. To improve accuracy and reduce the high computational costs, one can use higher than first or even second-order methods to approximate the advective flux terms in the hyperbolic conservation laws that describe the multicomponent transport in the reservoir. Those very high-order schemes can replace low-order approaches, such as the First Order Upwind (FOU) method, which is traditionally used in commercial petroleum reservoir simulators (CMG, 2019; ECLIPSE, 2022), obtaining more accurate solutions with reduced computational cost. In this work, we consider a three-phase and isothermal fluid flow of water, oil and gas and assume that there is no mass transfer between the water and the hydrocarbon phases. Physical dispersion and capillary effects are neglected. To solve the system of Partial Differential Equations (PDEs), we used the classical IMPEC (Implicit Pressure Explicit Composition) approach and, to model the complex phase behavior, we used the Peng–Robinson equation of state (PR-EOS). The classical Two Point Flux Approximation (TPFA) method is used to spatially discretize the diffusion terms in the pressure equation. For the first time in literature the very high order Flux Reconstruction (FR) method is adapted for the numerical simulation of the compositional model in petroleum reservoirs. The FR method is a numerical scheme employed to discretize the hyperbolic conservation laws using general grids. To avoid spurious oscillations in the vicinity of solution discontinuities, the h-MLP (hierarchical Multi-dimensional Limiting Process) is applied in the reconstruction stage. For the time integration, we have used the third-order Runge–Kutta approach. To appraise the robustness of our formulation, we have solved some one-dimensional benchmark problems found in literature and compared the FR solution with the FOU method and with the second-order Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) method. From our results, the FR method has shown improved accuracy, efficiency and capability of detecting sharper shocks when compared to the other methods that we have considered for the numerical modeling of the compositional flow in the present paper.

Spatiotemporal velocities, aggregation, and false dispersion: numerical considerations for the modeling of colonial and motile harmful algae

PurposeThe purpose of this study is to investigate the errors arising from the numerical treatment of model processes, paying particular attention to the impact of key system features including widely variable dispersion coefficients, spatiotemporal velocities of algal cells, and the aggregation of algae from single cells to large colonies. An advection–dispersion model has been presented to describe the vertical transport of colonial and motile harmful algae in a lake environment.Design/methodology/approachModel performance is examined for two different numerical treatments of the advective term: first-order upwind and quadratic upwind with a stability-preserving flux limiter (SMART). To determine how these schemes impact predictions, comparisons are made across a sequence of models with increasing complexity.FindingsUsing first-order upwinding for advection–dispersion calculations with a time oscillating velocity field leads to oscillatory numerical dispersion. Subjecting an initially uniform distribution of large-sized algal colonies to a spatiotemporal velocity creates a concentration pulse, which reaches a steady-state width at high-grid Peclet numbers when using the SMART scheme; the pulse exhibits contraction–expansion behavior throughout a velocity cycle at all Peclet numbers when using first-order upwinding. When aggregation dynamics are included with advection-dominated spatiotemporal transport, results indicate the SMART scheme predicts larger peak concentration values than those predicted by first-order upwind, but peak location and the time to large colony appearance remain largely unchanged between the two advective schemes.Originality/valueTo the best of the authors’ knowledge, this study is the first numerical investigation of a novel advection–dispersion model of vertical algal transport. In addition, a generalized expression for the effective dispersion coefficient of temporally variable flow fields is presented.

Improved fast iterative method for higher calculation accuracy of traveltimes

Fast Iterative Method (FIM), a grid-based ray-tracing algorithm, has a computational efficiency far exceeding that of conventional fast marching method (FMM) in serial computing, because it adopts an iterative algorithm and can be calculated in parallel. However, the finite difference scheme of the FIM only employs the first-order Godunov upwind finite difference scheme, resulting in low calculation accuracy of arrival traveltimes. This paper proposes a new algorithm of FIM, which utilizes second-order and mixed-order difference schemes under the parallel computing structure, effectively improves the calculation accuracy of traveltimes, and basically maintains the original computational efficiency. Especially, improved finite difference calculation in the diagonal direction of square grids and the double-grid technique near the source are applied to further reduce overall calculation error. Numerical experiments show that for the given model, the calculation error can be reduced from 2 ms of the original FIM calculation to 0.05 ms, and the calculation accuracy is increased by 40 times, while the time consumption increases only 27%.

An Innovative Approach for Numerical Solution of the Unsteady Convection-Dominated Flow Problems

In this study, convection-diffusion equation is solved numerically using four different space discretization methods namely first-order upwinding, second-order central difference, cubic (partially upwinded) and cubic-TVD (Total Variation Diminishing) techniques. All methods are compared with the analytical solution. The first-order method is not close to the analytical solution due to the numerical dispersion. The higher-order techniques reduce numerical dispersion. However, they cause another numerical error, unphysical oscillation. This study proposes an innovative approach on cubic-TVD method to eliminate undesired oscillations. Proposed model decreases numerical errors significantly compared to previously developed techniques. Moreover, numerical results of presented model quite close to the analytical solution. Finally, all Matlab codes of numerical and analytical solutions for convection-diffusion equation are added to Appendix in order to facilitate other researchers’ work.

Restoring the conservativity of characteristic-based segregated models: Application to the hybrid lattice Boltzmann method

A general methodology is introduced to build conservative numerical models for fluid simulations based on segregated schemes, where mass, momentum, and energy equations are solved by different methods. It is especially designed here for developing new numerical discretizations of the total energy equation and adapted to a thermal coupling with the lattice Boltzmann method (LBM). The proposed methodology is based on a linear equivalence with standard discretizations of the entropy equation, which, as a characteristic variable of the Euler system, allows efficiently decoupling the energy equation with the LBM. To this extent, any LBM scheme is equivalently written under a finite-volume formulation involving fluxes, which are further included in the total energy equation as numerical corrections. The viscous heat production is implicitly considered thanks to the knowledge of the LBM momentum flux. Three models are subsequently derived: a first-order upwind, a Lax–Wendroff, and a third-order Godunov-type schemes. They are assessed on standard academic test cases: a Couette flow, entropy spot and vortex convections, a Sod shock tube, several two-dimensional Riemann problems, and a shock–vortex interaction. Three key features are then exhibited: (1) the models are conservative by construction, recovering correct jump relations across shock waves; (2) the stability and accuracy of entropy modes can be explicitly controlled; and (3) the low dissipation of the LBM for isentropic phenomena is preserved.

A very high-order flux reconstruction approach coupled to the MPFA-QL finite volume method for the numerical simulation of oil-water flows in 2D petroleum reservoirs

In recent years, high-order methods for the discretization of the advection dominated saturation equation, and that results from the modelling of the 2D Oil-Water displacement through porous formations, have received renewed attention from the CFD community, as a suitable alternative to the traditional first-order upwind (FOU) finite volume (FV) scheme, which is the standard in industry. A key aspect of high-order methods applied to numerical simulation of petroleum reservoirs is the lower computational cost for a comparable degree of accuracy against their lower-order counterparts. In this context, in the present paper, for the first time in literature, we propose the use of the very high order Flux Reconstruction (FR) formulation for the discretization of the advective terms in the saturation equation using unstructured quadrilateral meshes. To solve the elliptic pressure equation, we use the recently developed Multipoint Flux Approximation method with a Quasi-Local Stencil (MPFA-QL). This method is very robust and able to handle highly heterogeneous and anisotropic permeability tensors. The system of equations is solved in a segregated manner by using the classical Implicit Pressure Explicit Saturation (IMPES) procedure. To properly couple the cell-centered MPFA-QL method with the FR formulation, it is necessary to obtain an adequate velocity reconstruction throughout the control volumes of the mesh. As 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 first-order Brezzi–Douglas–Marini (BDM) space, and the Piola transformation, to get the complete knowledge of the velocity field throughout the domain. To eliminate spurious oscillations around the discontinuities, which are typical in high-order methods, we have used a hierarchical multidimensional limiting process (hMLP). To compute fluxes over the control surfaces we have used the Roe-E approximate Riemann solver. Finally, through some representative examples, we show the accuracy, efficiency and robustness of our new high order numerical methodology to model oil-water displacements in highly heterogeneous and anisotropic petroleum reservoirs.

A high resolution numerical method for solving atmospheric models

We develop a numerical method to solve atmospheric model equations, which is a hybrid of the exponential time differencing (ETD) and the equation solving solution gradient (ESSG) method. The Coriolis force term due to geostrophic effects is handled by the ETD method and the convective flux balance terms are implemented by the ESSG method. In the ESSG method, gradients of solution variables are solved by their corresponding difference equations instead of interpolating from nearby solution values as in the conventional finite volume method. In this way, significant features of the model equations can be preserved and correctly considered to provide a reliable simulation. Feasibility of the proposed method is verified in a number of applications and shows it is more accurate than other common schemes, including the first-order upwind (UD) scheme and the second-order minimum modulus (minmod) scheme. We apply the method to a one-dimensional shock tube, a two-dimensional Riemann problem, and an atmospheric vortex flow. Our studies indicate that the exponential time differencing together with the equation solving solution gradient method provides very accurate simulations for the atmospheric flow model equations.

Proposals for Sewage Management at Luruaco Lake, Colombia

This study presents numerical simulations of faecal coliforms dynamics in Luruaco lake, located in Atl\'antico Department, Colombia. The velocity field is obtained through a two-dimensional horizontal (2DH) model of Navier-Stokes equations system. The transport equation of faecal coliforms concentration is provided from a convective-diffusive-reactive equation. The lake's geometry is built through cubic spline and multi-block methods. The discretization method by Finite Differences and the First Order Upwind (FOU) are applied to the 2DH model. The Mark and Cell (MAC) method is used to determine numerically the velocity field of water flow. Numerical simulations are carried out for a 72-hour period in order to understand the influence of faecal coliforms injections from each tributary. From the qualitative and quantitative analysis of the factors that influence faecal coliforms dynamics, proposals are presented, which aim to reduce contamination in some regions of Lake Luruaco. The numerical simulations show that the best option to improve water quality in the lake is the implementation of two actions, the diversion of the Lim\'on stream to the Negro stream and the installation of a sewage treatment plant at the mouth of the Negro stream. Other less expensive proposals are also presented.

Large-eddy simulation of flow around an isolated building: A step-by-step analysis of influencing factors on turbulent statistics

We designed and exclusively performed a wind tunnel experiment to obtain high-quality and high-resolution measurement results for three components of mean and fluctuating velocities around a 1: 1: 2 isolated building to establish a benchmark case for large-eddy simulation (LES) as a part of the activities of the Architectural Institute of Japan. Subsequently, we conducted LES to investigate the effect of computational grid arrangement on turbulent statistics by changing three grid arrangements. Finally, we performed a comparison by changing spatial discretization schemes, sub-grid scale (SGS) turbulence models, and convergence criteria. A grid in which the width of a building was discretized into 20 grids yielded sufficient resolution to reproduce the mean flow and second-order turbulent statistics. Under this appropriate grid arrangement, in the condition that the blending ratio of the first-order upwind interpolation scheme with the second-order linear interpolation scheme from 5 to 10%, the appropriate second-order statistics of wind velocity around the isolated building was obtained and numerical oscillation was effectively removed. Then, with this grid, any SGS model yielded appropriate results. In addition, the convergence criteria had little effects on the LES results if stable solution was obtained. Because the convergence criteria greatly affect the total calculation cost, weaker appropriate criteria should be set.

Two-Dimensional Numerical Study of the Transient Flow Conditions in Complete Shock Tunnel

In the current research, an axisymmetric model is developed to study high-speed unsteady flow in the test section of a 7 meter-long shock tunnel. The computational calculations of the shock tunnel are conducted using the Fluent CFD solver. The Finite Volume Method (FVM) is used to discretize the governing equations of mass, momentum, and energy. The accuracy of the numerical model is investigated with first-order upwind, second-order upwind, and third-order MUSCL schemes. Adaptive mesh refinement is implemented to resolve the shock wave and contact surface regions accurately. The numerical results are compared with theoretical calculations and experimental data from experimental tests and the comparison shows good agreement. Different test gases of Helium, Air and CO2, are utilized in the current study. The results show that steady test conditions are maintained for a longer test time by adjusting the pressure ratio and gas combination across the diaphragm. The highest shock wave speed and strength are achieved for a gas combination of Helium-CO2, but a longer test duration is observed when using Air as the test gas.

A well-balanced and positivity-preserving SPH method for shallow water flows in open channels

A well-balanced and positivity-preserving meshless method based on smoothed particle hydrodynamics (SPH) is developed to simulate one-dimensional (1D) and two-dimensional (2D) shallow water (SW) flows in open channels with irregular geometries. A new form of the characteristic equations that govern the water-surface level and water velocity is introduced to specify the numerical inflow/outflow boundary conditions. An additional condition, derived from temporal discretization to determine the time-step size, forces the water depth to be positive. A 1D finite volume shallow water (FVSW) model based on the first-order Godunov upwind method is built to conduct a comparison of the 1D meshless-based and mesh-based SW models. Six benchmark cases – still water, single trapezoidal and rectangular and prismatic and non-prismatic channels, and a dendritic channel network – are employed to validate the proposed models and compared with the exact and mesh-based numerical solutions. A real-world case of the Chicago Area Waterways System (CAWS) is investigated to highlight the performance of the proposed 1D model for a practical hydraulic system.

Physics-Capturing Discretizations for Spectral Wind-Wave Models

This paper discusses the discretization methods that have been commonly employed to solve the wave action balance equation, and that have gained a renewed interest with the widespread use of unstructured grids for third-generation spectral wind-wave models. These methods are the first-order upwind finite difference and first-order vertex-centered upwind finite volume schemes for the transport of wave action in geographical space. The discussion addresses the derivation of these schemes from a different perspective. A mathematical framework for mimetic discretizations based on discrete calculus is utilized herein. A key feature of this algebraic approach is that the process of exact discretization is segregated from the process of interpolation, the latter typically involved in constitutive relations. This can help gain insight into the performance characteristics of the discretization method. On this basis, we conclude that the upwind finite difference scheme captures the wave action flux conservation exactly, which is a plus for wave shoaling. In addition, we provide a justification for the intrinsic low accuracy of the vertex-centred upwind finite volume scheme, due to the physically inaccurate but common flux constitutive relation, and we propose an improvement to overcome this drawback. Finally, by way of a comparative demonstration, a few test cases is introduced to establish the ability of the considered methods to capture the relevant physics on unstructured triangular meshes.

Assessing Numerical Aspects of Transitional Flow Simulations Using the RANS Equations

ABSTRACT This paper addresses the influence of the scheme used in the convective terms of turbulence and transition transport equations on the numerical accuracy of transitional flow simulations. Three transition models are combined with the Shear Stress Transport turbulence model: the model, the γ model and the Amplification Factor Transport model. Two airfoils and a 6:1 prolate spheroid are used as test cases. Using first-order upwind in the convective terms of the k and ω transport equations has a weak effect on the numerical accuracy of the simulations when there is no transition model. However, when combined with the transition models, there is a significant improvement if second-order schemes are used. The iterative convergence obtained for the first-order approach was similar to that observed with the second-order schemes. Therefore, common ‘rules of thumb’ from RANS simulations at high Reynolds numbers may fail when transition models are included.

Implicit fast sweeping method for hyperbolic systems of conservation laws

Implicit time-accurate methods are often used to integrate stiff problems where explicit schemes impose severe time step restrictions. This paper presents an efficient numerical framework based on the fast sweeping method (FSM) for solving linear and nonlinear hyperbolic systems of conservation laws. The solution at each discrete location is computed by sweeping the numerical domain in several predetermined directions that follow the causality of the characteristic families. The use of a fractional step strategy eliminates the need for a solution selection criterion while one-sided stencils limit the number of sweeps to at most 2d for d space dimensions. This work focuses on the first-order implicit upwind method since it constitutes the building block for high-order conservative schemes. For problems where the degree of stiffness evolves over time, implicit-explicit hybridization can be accomplished with the same algorithm by simply switching the stencil at each time level. As opposed to traditional implicit solvers, the sweeping method does not require a local time linearization of the fluxes thereby preserving the nonlinear stability properties of the original implicit scheme. It also avoids the large computational and memory requirements associated with solving large block-diagonal systems of equations. A series of one- and two-dimensional test cases are presented for the inviscid Burgers' equation and the reactive Euler equations. The results indicate that the implicit FSM can allow for a major reduction in the number of time steps even in the presence of discontinuous solution profiles.

The importance of compatible sub-grid scale and spatial discretization models on the simulation of large-scale pool fire

Using discretization technique is one of the challenges associated with numerical modeling. In the Large Eddy Simulation (LES) method, this issue becomes more critical because it affects the sub-grid scale (SGS) model and creates errors. Using an LES method, we investigated the numerical errors and compatibility of different discretization methods with SGS in pool fire modeling. First-order, central, and second-order upwind and linear upwind stabilized transport (LUST) schemes were examined in Smagorinsky (S_SGS) and k-equation (K_SGS) sub-grid scales. First-order upwind and second-order upwind methods estimated the velocity field, temperature, and perturbation with a significant error, so that in the plume area, there was 90% error with the first-order upwind and 50% error with the second-order upwind concerning the results of the mean square horizontal velocity perturbation. With an error of 20%, the LUST method had a better agreement with the experimental results. Besides, SGS did not have a significant effect on the results of this method. However, the second-order upwind method was more consistent with the S_SGS model. To improve the results of upwind methods, Quadratic Upstream Interpolation for Convective Kinematics and Monotonic Upwind Scheme for Conservation Laws discretization methods are recommended.

Performance analysis of three iteration-free numerical methods for fast and accurate simulation of thermal dynamics in district heating pipeline

Numerical methods can provide sufficient information for thermal dynamics of district heating (DH) network. However, high computation burden restricts its application, especially for large-scale DH networks. This dilemma can be overcome by establishing fast numerical calculation procedure and properly selecting the calculation steps without reduction of numerical accuracy. The first-order implicit upwind method is an iteration-free numerical approach for fast simulation of DH pipeline thermal dynamics, while the first-order precision restricts its numerical performance. To further improve the simulation performance, two iteration-free numerical methods with high-order precision, the second-order implicit upwind method and third-order semi-implicit QUICK method are developed. Validations of the three numerical methods are conducted with the measured data of a real DH pipeline. Preferred calculation steps of the three methods are studied via comprehensive numerical experiments. Based on the preferred step size, simulation comparisons of the numerical methods are performed. Results show that the outlet temperature fluctuations of DH pipeline can be feasibly predicted by these methods with satisfying accuracy. The second-order implicit upwind method perform best considering least computation burden (0.003 s) within acceptable error. With the application of the second-order implicit upwind method, simulation analyses on DH pipeline under varying flow velocities are performed.

Locally bound-preserving enriched Galerkin methods for the linear advection equation

In this work, we introduce algebraic flux correction schemes for enriched (P1⊕P0 and Q1⊕P0) Galerkin discretizations of the linear advection equation. The piecewise-constant component stabilizes the continuous Galerkin approximation without introducing free parameters. However, violations of discrete maximum principles are possible in the neighborhood of discontinuities and steep fronts. To keep the cell averages and the degrees of freedom of the continuous P1/Q1 component in the admissible range, we limit the fluxes and element contributions, the complete removal of which would correspond to first-order upwinding. The first limiting procedure that we consider in this paper is based on the flux-corrected transport (FCT) paradigm. It belongs to the family of predictor-corrector algorithms and requires the use of small time steps. The second limiting strategy is monolithic and produces nonlinear problems with well-defined residuals. This kind of limiting is well suited for stationary and time-dependent problems alike. The need for inverting consistent mass matrices in explicit strong stability preserving Runge–Kutta time integrators is avoided by reconstructing nodal time derivatives from cell averages. Numerical studies are performed for standard 2D test problems.