Temporal Discretization Schemes Research Articles

An efficient numerical method for pricing a Russian option with a finite time horizon

In this paper, we present a finite difference scheme for a linear complementarity problem with a mixed boundary condition arising from pricing a Russian option with a finite time horizon. An implicit Euler method for the temporal discretization and second-order difference schemes on a piecewise uniform mesh for the spatial discretization are used to solve the linear complementarity problem with a mixed boundary condition. It is shown that the transformed discrete operator satisfies a maximum principle, which is used to derive the error estimate. It is proved that the scheme is first- and second-order convergent with respect to the temporal and spatial variables, respectively. Numerical experiments verify the validity of the theoretical results.

Implicit large eddy simulation of passive scalar transfer in compressible planar jet

AbstractImplicit large eddy simulation (ILES) of passive scalar transfer in compressible turbulence is evaluated for subsonic and supersonic turbulent planar jets. The ILES used in this study relies on fully explicit numerical schemes for spatial and temporal discretization and low‐pass and shock‐capturing filters used as an implicit subgrid‐scale (SGS) model. The ILES results are compared with the direct numerical simulation (DNS) database of the same flows. The ILES results exhibit good agreements with the DNS for first‐ and second‐order statistics of velocity and passive scalar. The scalar transport by turbulent velocity fluctuations is well captured by the ILES. The temporal evolution of the jet strongly depends on the jet Mach number, where a higher Mach number results in the delay of jet development. The Mach number dependence of velocity and passive scalar fields is consistent between the ILES and DNS. The low‐pass filters used as the implicit SGS model contribute to the dissipation of turbulent kinetic energy and scalar variance. Under the present numerical conditions, the filters account for about 50% of the dissipation in a fully developed turbulent jet. The dissipation rate in the ILES, which is the sum of the grid‐scale and SGS dissipation rates, is very close to the dissipation rate in the DNS, and the amount of the SGS dissipation is well controlled by the low‐pass filters. The filters also dump numerical oscillations in the velocity field caused by strong pressure waves outside the supersonic jet at the high Mach number.

Multiphysics for nuclear energy applications using a cohesive computational framework

With the recent development of advanced numerical algorithms, software design, and low-cost high-performance computer hardware, reliance on coupled multiphysics to predict the behavior of complex physical systems is beginning to become standard practice. This is especially true in nuclear energy applications where strong nonlinear interdependencies exist between reactor physics, radiation transport, multi-scale nuclear fuels performance, thermal fluids, etc. Resolving these nonlinear dependencies requires choices in multiphysics software approaches. Two main multiphysics modeling and simulation approaches have emerged. The first is based upon “code coupling” where disparate physics codes of different software design, code languages, and spatial and temporal integration schemes are coupled together with relatively complex data passing interfaces. The second multiphysics software approach is to employ a “cohesive” framework where all physics applications are developed with a common software design, i.e., data structures, syntax, input format, integrated spatial and temporal discretization schemes, etc. Here we present the Multiphysics Object-Oriented Simulation Environment (MOOSE) development and runtime framework and describe the framework’s cohesive modeling and simulation multiphysics approach. Then, a “cohesive-like” extension of the MOOSE framework is presented where MOOSE-based physics software applications are efficiently coupled to non-MOOSE (external) physics codes to form multiphysics applications using MOOSE’s unique interface capabilities. Finally, several examples of MOOSE’s cohesive and cohesive-like multiphysics applications will be demonstrated. These multiphysics demonstrations will incorporate both MOOSE-based applications and external codes, including Nek5000, RELAP-7, TRACE, BISON, and Pronghorn.

Mixed variational one-dimensional dynamic thermo-viscoplasticity for wave propagation

In view of recent developments in advanced manufacturing and nanomaterial technologies, which have enabled realization of novel material architectures at small-scale, a mixed variational temperature-dependent rate-dependent thermo-viscoplastic continuum element is proposed and developed for elastic, plastic, and thermal wave propagation problems. The variational scheme is based on a Hamiltonian approach, while for the temporal and spatial discretization scheme a discrete calculus of variations approach is adopted. Within this framework, an exponential temperature dependence for material properties and second sound effects is included. The developed framework sets the foundations that can lead to new class of temperature dependent models for capturing thermoplastic behavior and wave propagation of different materials. The robustness of the proposed variational approach, and aspects of coupled dynamic temperature-dependent thermoplastic response are demonstrated through several case studies. Applications include elastic, plastic, and thermal wave propagation and reflection in semi-infinite and finite length space at small temporal scales.

Non-linearly stable reduced-order models for incompressible flow with energy-conserving finite volume methods

A novel reduced-order model (ROM) formulation for incompressible flows is presented with the key property that it exhibits non-linearly stability, independent of the mesh (of the full order model), the time step, the viscosity, and the number of modes. The two essential elements to non-linear stability are: (1) first discretise the full order model, and then project the discretised equations, and (2) use spatial and temporal discretisation schemes for the full order model that are globally energy-conserving (in the limit of vanishing viscosity). For this purpose, as full order model a staggered-grid finite volume method in conjunction with an implicit Runge-Kutta method is employed. In addition, a constrained singular value decomposition is employed which enforces global momentum conservation. The resulting ‘velocity-only’ ROM is thus globally conserving mass, momentum and kinetic energy. For non-homogeneous boundary conditions, a (one-time) Poisson equation is solved that accounts for the boundary contribution. The stability of the proposed ROM is demonstrated in several test cases. Furthermore, it is shown that explicit Runge-Kutta methods can be used as a practical alternative to implicit time integration at a slight loss in energy conservation.

Fully resolved simulation of dense suspensions of freely evolving buoyant particles using an improved immersed boundary method

Fully resolved simulation of flows with buoyant particles is a challenging problem since buoyant particles are lighter than the surrounding fluid, and as a result, the two phases are strongly coupled together. In this work, the virtual force stabilization technique introduced by Schwarz et al. [Schwarz, S., Kempe, T., & Fröhlich, J. (2015). A temporal discretization scheme to compute the motion of light particles in viscous flows by an immersed boundary method. J. Comput. Phys., 281, 591–613] is extended to simulate buoyant particle suspensions with high volume fractions (up to 40%). It is concluded that the dimensionless numerical model constant Cv in the virtual force technique should increase with volume fraction. The behavior of a single rising particle, two in-line rising particles, and buoyant particle suspensions are studied. In each case, results are compared with experimental works on bubbly flows to highlight the differences and similarities between buoyant particles and bubbles. Finally, the drag coefficient is extracted from simulations of buoyant particle suspensions at different volume fractions, and based on that a drag correlation is presented.

Determination of Numerical Errors in Under-Resolved DNS of Turbulent Non-isothermal Flows

The methodology to quantify the numerical dissipation in Under-resolved DNS (UDNS) based on the balance of the kinetic energy equation by Schranner et al. (Comput Fluids 114:84–97, 2015), has been examined in this study. Furthermore, this methodology has been extended for considering active scalars, based on both balance of the kinetic energy and thermal variance equations for assessing the quality of UDNS. As a first step towards the analysis of reactive flows, where temperature is actively coupled to the momentum equation, the turbulent Rayleigh–Bénard convection in a cubical enclosure is considered here. The simulations have been carried out for different Prandtl numbers (Pr=0.1, 1.0, 10) to analyse different levels of momentum-temperature coupling for a representative value of nominal Rayleigh number (Ra=10^7) in order to investigate the effects of grid resolution and discretisation schemes on the numerical dissipation. It has been found that the numerical dissipation for both kinetic energy (R_{k}) and thermal variance (R_{t}) can range between positive or negative values depending on the combined effects of temporal and spatial discretisation schemes and grid spacing. The positive values of R_{k} and R_{t} indicate a loss of the kinetic energy and thermal variance due to numerical errors, whereas the negative values of R_{k} and R_{t} result in a moderate attenuation of these quantities. Accordingly, the numerical dissipation results obtained here (i.e. R_{k} and R_{t}) have been utilised for evaluating the effective values of the Rayleigh and Prandtl number (i.e. Ra_{eff}, Pr_{eff}) for UDNS. Using the scaling of the mean Nusselt number (Nu sim {Ra^a} {Pr^b}), it has been shown that the error relative to exact mean Nusselt number can be estimated based on R_{k}, R_{t} and the effective mean Nusselt number of the UDNS (Nu_{UDNS}). Overall, this methodology has been found a promising tool for the quality assessment of UDNS for the applications of non-isothermal flows (e.g. Rayleigh–Bénard convection).

MPI-AMRVAC: A parallel, grid-adaptive PDE toolkit

We report on the latest additions to our open-source, block-grid adaptive framework MPI-AMRVAC, which is a general toolkit for especially hyperbolic/parabolic partial differential equations (PDEs). Applications traditionally focused on shock-dominated, magnetized plasma dynamics described by either Newtonian or special relativistic (magneto)hydrodynamics, but its versatile design easily extends to different PDE systems. Here, we demonstrate applications covering any-dimensional scalar to system PDEs, with e.g. Korteweg–de Vries solutions generalizing early findings on soliton behavior, shallow water applications in round or square pools, hydrodynamic convergence tests as well as challenging computational fluid and plasma dynamics applications. The recent addition of a parallel multigrid solver opens up new avenues where also elliptic constraints or stiff source terms play a central role. This is illustrated here by solving several multi-dimensional reaction–diffusion-type equations. We document the minimal requirements for adding a new physics module governed by any nonlinear PDE system, such that it can directly benefit from the code flexibility in combining various temporal and spatial discretization schemes. Distributed through GitHub, MPI-AMRVAC can be used to perform 1D, 1.5D, 2D, 2.5D or 3D simulations in Cartesian, cylindrical or spherical coordinate systems, using parallel domain-decomposition, or exploiting fully dynamic block quadtree-octree grids.

Scalability of OpenFOAM Density-Based Solver with Runge–Kutta Temporal Discretization Scheme

Compressible density-based solvers are widely used in OpenFOAM, and the parallel scalability of these solvers is crucial for large-scale simulations. In this paper, we report our experiences with the scalability of OpenFOAM’s native rhoCentralFoam solver, and by making a small number of modifications to it, we show the degree to which the scalability of the solver can be improved. The main modification made is to replace the first-order accurate Euler scheme in rhoCentralFoam with a third-order accurate, four-stage Runge-Kutta or RK4 scheme for the time integration. The scaling test we used is the transonic flow over the ONERA M6 wing. This is a common validation test for compressible flows solvers in aerospace and other engineering applications. Numerical experiments show that our modified solver, referred to as rhoCentralRK4Foam, for the same spatial discretization, achieves as much as a 123.2% improvement in scalability over the rhoCentralFoam solver. As expected, the better time resolution of the Runge–Kutta scheme makes it more suitable for unsteady problems such as the Taylor–Green vortex decay where the new solver showed a 50% decrease in the overall time-to-solution compared to rhoCentralFoam to get to the final solution with the same numerical accuracy. Finally, the improved scalability can be traced to the improvement of the computation to communication ratio obtained by substituting the RK4 scheme in place of the Euler scheme. All numerical tests were conducted on a Cray XC40 parallel system, Theta, at Argonne National Laboratory.

Superconvergence analysis of an energy stable scheme for nonlinear reaction-diffusion equation with BDF mixed FEM

A step-2 backward differential formula (BDF) temporal discretization scheme is constructed for nonlinear reaction-diffusion equation and superconvergence results are studied by mixed finite element method (FEM) with the elements Q11 and Q01×Q10 unconditionally. In particular, we apply an artificial regularization term to guarantee the energy stability of the step-2 BDF scheme. Splitting technique is utilized to get rid of the ratio between the time step size τ and the subdivision parameter h. Temporal error estimates in H2-norm are derived by use of the function's monotonicity, which leads to the regularities of the solutions for the time-discrete equations. Spatial error estimates in L2-norm are deduced to bound the numerical solution in L∞-norm. Unconditional superconvergence estimates of un in H1-norm and q→n=∇un in (L2)2-norm with order O(h2+τ2) are obtained. The global superconvergent properties are deduced through above results. Two numerical examples testify the theoretical analysis.

Strong convergence of a stochastic Rosenbrock-type scheme for the finite element discretization of semilinear SPDEs driven by multiplicative and additive noise

This paper aims to investigate the numerical approximation of a general second order parabolic stochastic partial differential equation(SPDE) driven by multiplicative and additive noise. Our main interest is on such SPDEs where the nonlinear part is stronger than the linear part, usually called stochastic dominated transport equations. Most standard numerical schemes lose their good stability properties on such equations, including the current linear implicit Euler method. We discretize the SPDE in space by the finite element method and propose a novel scheme called stochastic Rosenbrock-type scheme for temporal discretization. Our scheme is based on the local linearization of the semi-discrete problem obtained after space discretization and is more appropriate for such equations. We provide a strong convergence of the new fully discrete scheme toward the exact solution for multiplicative and additive noise and obtain optimal rates of convergence. Numerical experiments to sustain our theoretical results are provided.

Simulating Turbulent Thermal Plasma Flows for Nanopowder Fabrication

This article presents descriptions of theoretical models and numerical methods for simulating turbulent thermal plasma flow with nanopowder growth. Turbulence models must express turbulent and laminar states because both states co-exist with thermal plasmas showing large density variation and transport properties. Time-dependent 3D simulations are conducted based on Large Eddy Simulation using a dynamic Smagorinsky model. Results show significant difference depending on temporal and spatial discretization schemes and velocity–pressure coupling algorithms. Simulation results demonstrate that advanced numerical methods with high-order accuracy should be used for long and robust computations capturing steep gradients of nanopowder concentration and plasma temperature and 3D dynamic motions of multiscale vortices, which are turbulent features of thermal plasma flows with low Mach numbers. A thermal plasma jet generates a double-layer structure of inner high-temperature thick vortex rings and outer low-temperature thin vortex rings near the nozzle exit. Flowing downstream, these vortices interact, deform, and break up. Consequently, plasma transits to a complex thermal flow. The widely spreading distribution of multiscale vortices agrees with experimental observations, which are not simulated using conventional methods. Nanopowder is generated from material vapour by nucleation and condensation at interfacial regions between plasma and cold gas. Those regions include numerous vortices. Therefore, the vortices convey the nanopowder, producing a complex distribution of nanopowder. Simultaneously, the nanopowder diffuses and increases in size, decreasing in number by interparticle coagulation. Cross-correlation analysis suggests that a nanopowder distribution distant from a plasma jet can be controlled through temperature fluctuation control at the upstream plasma fringe.

A simplified finite volume lattice Boltzmann method for simulations of fluid flows from laminar to turbulent regime, Part I: Numerical framework and its application to laminar flow simulation

In this paper, an unstructured grid based finite volume lattice Boltzmann method (FVLBM) that can be used for the simulation of incompressible laminar flows is presented and studied in detail. This method is derived from a simple modification of the cell-vertex unstructured grid based FVLBM proposed by Stiebler et al. (2006). Compared with other complex flux reconstruction methods, the present scheme has a low computational cost and can also achieve second-order spatial accuracy. Furthermore, depending on the use of the different temporal discretization schemes, the temporal accuracy can be adjusted for both steady and unsteady flows. Besides, some comparisons of the computational cost and accuracy with another FVLBM scheme are also presented. Meanwhile, different boundary conditions are illustrated that are easy to implement on complex geometries. To validate the present method, four cases are carried out, including a Couette flow driven by one plate for an accuracy test, flow in a square cavity, flow around a single circular cylinder and a more complex flow around double circular cylinders. Numerical experiments show that the present scheme can simulate steady and unsteady flows at relatively high Reynolds number with relatively few grid cells, thus demonstrating the good capability of the present method.

Numerical Solution of Viscoelastic Fluid‐Structure‐Diffusion Systems with Applications in Ophthalmology

AbstractThe research of fluid‐structure interaction problems is a continuously growing field, especially regarding applications in medicine and biology. We present the coupling of a potentially viscoelastic fluid with multiple hyperelastic structures incorporating chemical processes in the arbitrary Lagrangian Eulerian framework. This monolithic formulation allows a robust numerical solution with Newton's method. The discretization is based on the backward Euler scheme for temporal discretization and the Galerkin finite element method for spatial discretization. This fluid‐structure interaction problem is applied to ophthalmology in order to improve the medical treatment of retinal diseases. The physiological processes include the elastic response of various structures like the sclera, lens and iris coupled to the fluid‐like vitreous which is modeled by a viscoelastic Burgers type model for the healthy case and by the Newtonian Navier‐Stokes equations for the pathological case. Since most medical treatments are based on the injection of medicine we furthermore study the drug distribution, which is modeled by convection‐diffusion‐reaction equations, in the whole eye for healthy and non‐healthy pathologies.

CFD Julia: A Learning Module Structuring an Introductory Course on Computational Fluid Dynamics

CFD Julia is a programming module developed for senior undergraduate or graduate-level coursework which teaches the foundations of computational fluid dynamics (CFD). The module comprises several programs written in general-purpose programming language Julia designed for high-performance numerical analysis and computational science. The paper explains various concepts related to spatial and temporal discretization, explicit and implicit numerical schemes, multi-step numerical schemes, higher-order shock-capturing numerical methods, and iterative solvers in CFD. These concepts are illustrated using the linear convection equation, the inviscid Burgers equation, and the two-dimensional Poisson equation. The paper covers finite difference implementation for equations in both conservative and non-conservative form. The paper also includes the development of one-dimensional solver for Euler equations and demonstrate it for the Sod shock tube problem. We show the application of finite difference schemes for developing two-dimensional incompressible Navier-Stokes solvers with different boundary conditions applied to the lid-driven cavity and vortex-merger problems. At the end of this paper, we develop hybrid Arakawa-spectral solver and pseudo-spectral solver for two-dimensional incompressible Navier-Stokes equations. Additionally, we compare the computational performance of these minimalist fashion Navier-Stokes solvers written in Julia and Python.

Convergence and error analysis of an automatically differentiated finite volume based heat conduction code

Purpose This paper aims to investigate the convergence and error properties of a finite volume-based heat conduction code that uses automatic differentiation to evaluate derivatives of solutions outputs with respect to arbitrary solution input(s). A problem involving conduction in a plane wall with convection at its surfaces is used as a test problem, as it has an analytical solution, and the error can be evaluated directly. Design/methodology/approach The finite volume method is used to discretize the transient heat diffusion equation with constant thermophysical properties. The discretized problem is then linearized, which results in two linear systems; one for the primary solution field and one for the secondary field, representing the derivative of the primary field with respect to the selected input(s). Derivatives required in the formation of the secondary linear system are obtained by automatic differentiation using an operator overloading and templating approach in C++. Findings The temporal and spatial discretization error for the derivative solution follows the same order of accuracy as the primary solution. Second-order accuracy of the spatial and temporal discretization schemes is confirmed for both primary and secondary problems using both orthogonal and non-orthogonal grids. However, it has been found that for non-orthogonal cases, there is a limit to the error reduction, which is concluded to be a result of errors in the Gauss-based gradient reconstruction method. Originality/value The convergence and error properties of derivative solutions obtained by forward mode automatic differentiation of finite volume-based codes have not been previously investigated.

Consistent upwinding for sequential fully implicit multiscale compositional simulation

There is strong interest to design sequential fully implicit (SFI) methods for compositional flow simulations with convergence properties that are comparable to fully implicit (FI) methods. SFI methods decompose the fully coupled system into a pressure equation and a transport system of the components. During the pressure update, the compositions are frozen, and during the transport calculations, both the pressure and total velocity are kept constant. The two systems are solved sequentially, and the solution, which is a fully implicit one, is obtained by controlling the splitting errors due to the decoupling. Having an SFI scheme that enjoys a convergence rate similar to FI makes it possible to design specialized numerical methods optimized for the different parabolic and the hyperbolic operators, as well as the use of high-order spatial and temporal discretization schemes. Here, we use the multiscale restriction-smoothed basis (MsRSB) method for the parabolic operator. We also show that phase-potential upwinding is incompatible with the total velocity formulation of the fluxes, which is common in SFI schemes. We observe that in cases with strong gravity or capillary pressure, it is possible to have flow reversals. These reversals can strongly affect the convergence rate of SFI methods. In this work, we employ phase upwinding (PU) as well as implicit hybrid upwinding (IHU) with a SFI method. IHU determines the upwinding direction differently for the viscous, buoyancy, and capillary pressure terms in the phase velocity expressions. The use of IHU leads to a consistent SFI scheme in terms of both pressure and compositions, and it improves the SFI convergence significantly in settings with strong buoyancy or capillarity. We demonstrate the robustness of the IHU-based SFI algorithm across a wide parameter range. Realistic compositional models with gas and water injection are presented and discussed.

LES and RANS modelling of under-expanded jets with application to gaseous fuel direct injection for advanced propulsion systems

A density-based solver with the classical fourth-order accurate Runge-Kutta temporal discretization scheme was developed and applied to study under-expanded jets issued through millimetre-size nozzles for applications in high-pressure direct-injection (DI) gaseous-fuelled propulsion systems. Both large eddy simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) turbulence modelling techniques were used to evaluate the performance of the new code. The computational results were compared both quantitatively and qualitatively against available data from the literature. After initial evaluation of the code, the computational framework was used in conjunction with RANS modelling (k-ω SST) to investigate the effect of nozzle exit geometry on the characteristics of gaseous jets issued from millimetre-size nozzles. Cylindrical nozzles with various length to diameter ratios, namely 5, 10 and 20, in addition to a diverging conical nozzle, were studied. This study is believed to be the first to provide a direct comparison between RANS and LES within the context of nozzle exit profiling for advanced high-pressure injection systems with the formation of under-expanded jets. It was found that reducing the length of the straight section of the nozzle by 50% resulted in a slightly higher level of under-expansion (∼2.6% higher pressure at the nozzle exit) and ∼1% higher mass flow rate. It was also found that a nozzle with 50% shorter length resulted in ∼6% longer jet penetration length. At a constant nozzle pressure ratio (NPR), a lower nozzle length to diameter ratio resulted in a noticeably higher jet penetration. It was found that with a diverging conical nozzle, a fairly higher penetration length could be achieved if an under-expanded jet formed downstream of the nozzle exit compared to a jet issued from a straight nozzle with the same NPR. This was attributed to the radial restriction of the flow and consequently formation of a relatively smaller reflected shock angle. With the conical nozzle used in this study and a 30 bar injection pressure, an under-expanded hydrogen jet exhibited ∼60% higher penetration length compared to an under-expanded nitrogen jet at 100 μs after start of injection. Moreover, the former jet exhibited ∼22% higher penetration compared to a nitrogen jet issued through the conical profile with 150 bar injection pressure.

Pizza: an open-source pseudo-spectral code for spherical quasi-geostrophic convection

We present a new pseudo-spectral open-source code nicknamed pizza. It is dedicated to the study of rapidly-rotating Boussinesq convection under the 2-D spherical quasi-geostrophic approximation, a physical hypothesis that is appropriate to model the turbulent convection that develops in planetary interiors. The code uses a Fourier decomposition in the azimuthal direction and supports both a Chebyshev collocation method and a sparse Chebyshev integration formulation in the cylindrically-radial direction. It supports several temporal discretisation schemes encompassing multi-step time steppers as well as diagonally-implicit Runge-Kutta schemes. The code has been tested and validated by comparing weakly-nonlinear convection with the eigenmodes from a linear solver. The comparison of the two radial discretisation schemes has revealed the superiority of the Chebyshev integration method over the classical collocation approach both in terms of memory requirements and operation counts. The good parallelisation efficiency enables the computation of large problem sizes with $\mathcal{O}(10^4\times 10^4)$ grid points using several thousands of ranks. This allows the computation of numerical models in the turbulent regime of quasi-geostrophic convection characterised by large Reynolds $Re$ and yet small Rossby numbers $Ro$. A preliminary result obtained for a strongly supercritical numerical model with a small Ekman number of $10^{-9}$ and a Prandtl number of unity yields $Re\simeq 10^5$ and $Ro \simeq 10^{-4}$. pizza is hence an efficient tool to study spherical quasi-geostrophic convection in a parameter regime inaccessible to current global 3-D spherical shell models.

Analysis of Robustness of Variational Multiscale Error Estimators for the Forward Propagation Study

Error estimation and control play a crucial role in the accuracy and reliability of numerical simulations. Variational multiscale (VMS) methods have emerged as effective techniques for error estimation in forward propagation studies. The robustness and accuracy of VMS error estimators under different conditions and parameters need to be thoroughly investigated. This research paper focuses on the analysis of the robustness of variational multiscale error estimators for forward propagation studies. The study explores the performance and limitations of VMS error estimators under various scenarios, including different mesh resolutions, material properties, and boundary conditions. The research aims to provide insights into the effectiveness and reliability of VMS error estimators, contributing to the improvement of forward propagation simulations. The analysis of robustness of variational multiscale error estimators for the forward propagation study addresses the important challenge of accurately assessing and quantifying the errors in computational models used for forward propagation. Variational multiscale error estimators have gained significant attention as they provide a reliable framework for error estimation in complex physical simulations. However, the robustness of these estimators, particularly in the context of forward propagation, remains an open question. To analyze the robustness of variational multiscale error estimators for the forward propagation study. We begin by formulating the problem and defining the mathematical framework for error estimation. Next, we investigate the performance and reliability of variational multiscale error estimators under different scenarios, including variations in model parameters, mesh resolution, and boundary conditions. We evaluate the accuracy and consistency of the estimators by comparing the estimated errors with reference solutions or known analytical solutions.
 We examine the influence of different factors, such as spatial and temporal discretization schemes, on the robustness of variational multiscale error estimators. We assess the sensitivity of the estimators to modeling assumptions and potential sources of error, such as numerical approximation and solution regularization. The effect of different physical phenomena, such as material nonlinearity or complex boundary conditions, on the performance of the estimators. The findings of this provide valuable insights into the robustness and limitations of variational multiscale error estimators for the forward propagation study. Understanding the performance characteristics and potential weaknesses of these estimators is crucial for accurately assessing the reliability and uncertainty of computational models used in forward propagation. The analysis of robustness contributes to the development of improved error estimation techniques and enhances the confidence in the results obtained from computational simulations. This focuses on the analysis of robustness of variational multiscale error estimators for the forward propagation study. By investigating their performance under different scenarios and evaluating their sensitivity to various factors, we provide a comprehensive understanding of the reliability and limitations of these estimators. The insights gained from this analysis have implications for improving the accuracy and confidence in computational models used for forward propagation in diverse fields, including engineering, physics, and computational sciences.