Showing papers in "Computers & Mathematics With Applications in 2020"

Lump, mixed lump-stripe and rogue wave-stripe solutions of a (3+1)-dimensional nonlinear wave equation for a liquid with gas bubbles

Abstract: Liquids with gas bubbles are very common in science, engineering, physics, nature and life. Under investigation in this paper is a (3+1)-dimensional nonlinear wave equation for a liquid with gas bubbles. With respect to the velocity of the liquid-gas bubble mixture, we obtain the lump, rogue wave, mixed lump-stripe soliton and mixed rogue wave-stripe soliton solutions via the symbolic computation. Based on the mixed lump-stripe soliton solutions, we construct the rogue wave and investigate the fusion and fission phenomena between a lump and the one-stripe soliton. We graphically study the mixed lump-stripe soliton under the influence of the parameters β , γ , δ and ξ , which represent the dispersion, perturbed effect, disturbed wave velocities along the y and z (i.e., the two transverse) directions, respectively. With the decreasing value of β to -1, the graph from a lump and one-stripe soliton shows a soliton; with the increasing value of γ to -3, location of the lump moves along the negative x axis; with the value of δ increasing to 0.5, location of the lump moves along the positive x axis; with the increasing value of ξ to -3, location and range of the lump soliton keep unchanged. With respect to the velocity of the mixture, we obtain the interaction between a rogue wave and a pair of stripe solitons according to mixed rogue wave-stripe soliton solutions. A lump is provided.

A two-grid MMOC finite element method for nonlinear variable-order time-fractional mobile/immobile advection–diffusion equations

Abstract: A fully discrete two-grid modified method of characteristics (MMOC) scheme is proposed for nonlinear variable-order time-fractional advection–diffusion equations in two space dimensions. The MMOC is used to handle the advection-dominated transport and the two-grid method is designed for efficiently solving the resulting nonlinear system. Optimal L 2 error estimates are derived for both the MMOC scheme and the corresponding two-grid MMOC scheme. Numerical experiments are presented to demonstrate the accuracy and the efficiency of the proposed method.

A unified view of some numerical methods for fractional diffusion

Abstract: In recent years, a number of numerical methods for the solution of fractional Laplace and, more generally, fractional diffusion problems have been proposed. The approaches are quite diverse and include, among others, the use of best uniform rational approximations, quadrature for Dunford–Taylor-like integrals, finite element approaches for a localized elliptic extension into a space of increased dimensions, and time stepping methods for a parabolic reformulation of the fractional differential equation. A systematic comparison, both theoretical and experimental, of these approaches has thus far been lacking. A main contribution of the present work is the observation that all approaches mentioned above can, in fact, be interpreted as realizing different rational approximations of a univariate function over the spectrum of the original (non-fractional) diffusion operator. While this is obvious for some of the methods, it is a new result in particular for extension-based and time stepping approaches. This observation allows us to cast all described methods into a unified theoretical and computational framework, which has a number of benefits. Theoretically, it enables us to develop new convergence proofs for several of the studied methods, clarifies similarities and differences between the approaches, suggests how to design new and improved methods, and allows a direct comparison of the relative performance of the various methods. Practically, it provides a single, simple to implement, efficient and fully parallel algorithm for the realization of all studied methods; for instance, this immediately yields a fast and memory-efficient way of realizing all tensor product extension methods and lets us parallelize the otherwise inherently sequential time stepping approach. Finally, we present a detailed numerical study comparing all investigated methods for various fractional exponents and draw conclusions from the results. The comparison is made fair by the central insight that the computational effort of all these methods depends only on a single parameter, the degree of the underlying rational approximation. As a point of comparison, we also test a simple rational approximation method based on a black-box direct rational approximation algorithm which performs very well in practice.

Analysis of an augmented moving least squares approximation and the associated localized method of fundamental solutions

Abstract: The localized method of fundamental solutions (LMFS) is an efficient meshless collocation method that combines the concept of localization and the method of fundamental solutions (MFS). The resultant system of linear algebraic equations in the LMFS is sparse and banded and thus, drastically reduces the storage and computational burden of the MFS. In the LMFS, the moving least square (MLS) approximation, based on fundamental solutions, is used to construct approximate solution at each node. In this paper, this fundamental solutions-based MLS approximation, named as an augmented MLS (AMLS) approximation, is generalized to any point in the computational domain. Computational formulas, theoretical properties and error estimates of the AMLS approximation are derived. Then, taking Laplace equation as an example, this paper sets up a framework for the theoretical error analysis of the LMFS. Finally, numerical results are presented to verify the efficiency and theoretical results of the AMLS approximation and the LMFS. Convergence and comparison researches are conducted to validate the accuracy, convergence and efficiency.

Free vibration analysis of FG-CNTRC shell structures using the meshfree radial point interpolation method

Abstract: This study conducts the first-known free vibration analysis of functionally graded carbon nanotubes-reinforced (FG-CNTRC) shell structures using the meshfree radial point interpolation method (RPIM). The modified first-order shear deformation theory (modified FSDT) is implemented to get the realistic effect of the transverse shear deformation with its parabolic distribution. Numerical examples are carried out to examine the convergence and accuracy of the element-free RPIM method in its application to the free vibration FG-CNTRC analysis of shell structures. Results obtained using the proposed meshfree method, demonstrate that the improved FSDT is very successful compared to closed-form solutions and finite element results using different shell theories.

The p- and hp-versions of the virtual element method for elliptic eigenvalue problems

Abstract: We discuss the p - and h p -versions of the virtual element method for the approximation of eigenpairs of elliptic operators with a potential term on polygonal meshes. An application of this model is provided by the Schrodinger equation with a pseudo-potential term. As an interesting byproduct, we present for the first time in literature an explicit construction of the stabilization of the mass matrix. We present in detail the analysis of the p -version of the method, proving exponential convergence in the case of analytic eigenfunctions. The theoretical results are supplied with a wide set of experiments. We also show numerically that, in the case of eigenfunctions with finite Sobolev regularity, an exponential approximation of the eigenvalues in terms of the cubic root of the number of degrees of freedom can be obtained by employing h p -refinements. Importantly, the geometric flexibility of polygonal meshes is exploited in the construction of the h p -spaces.

The conforming virtual element method for polyharmonic problems

Abstract: In this work, we exploit the capability of virtual element methods in accommodating approximation spaces featuring high-order continuity to numerically approximate differential problems of the form ( − Δ ) p u = f , p ≥ 1 . More specifically, we develop and analyze the conforming virtual element method for the numerical approximation of polyharmonic boundary value problems, and prove an abstract result that states the convergence of the method in suitable norms.

Multi-dimensional spectral tau methods for distributed-order fractional diffusion equations

Abstract: The distributed-order fractional diffusion equation is a generalization of the standard fractional diffusion equation that can model processes lacking power-law scaling over the whole time-domain. An important application of distributed-order diffusions is to model ultraslow diffusion where a plume of particles spreads at a logarithmic rate. To broaden the range of applicability of distributed-order fractional diffusion models, efficient numerical methods are needed to solve the model equation. In this work, we develop spectral tau schemes to discretize the fractional diffusion equation with distributed-order fractional derivative in time and Dirichlet boundary conditions. The model solution is expanded in multi-dimensions in terms of Legendre polynomials and the discrete equations are obtained with the tau method. Numerical examples are provided to highlight the convergence rate and the flexibility of this approach. The proposed spectral tau methods yield an exponential rate of convergence when the solution is smooth. Our results confirm that nonlocal numerical methods are best suited to discretize distributed-order fractional differential equations as they naturally take the global behavior of the solution into account.

Variational multiscale interpolating element-free Galerkin method for the nonlinear Darcy–Forchheimer model

Abstract: In this paper, the variational multiscale interpolating element-free Galerkin (VMIEFG) method is developed to obtain the numerical solution of the nonlinear Darcy–Forchheimer model. We use the interpolating moving least squares method instead of the moving least squares approximation to construct meshless shape functions with delta function properties. Then the flux boundary condition of the Darcy–Forchheimer model can be handled easily. Hughes’ variational multiscale (HVM) method is applied to overcome the numerical oscillation caused by equal-order basis for the velocity and pressure. Moreover, the HVM ensures that the resultant formulation in the VMIEFG method is consistent and the stabilization parameter (or tensor) appears naturally. Consequently, the stabilization parameter is free of user-defined. The fixed point iteration method is used to deal with the nonlinear term. Some numerical examples are provided to illustrate the stability and performance of the proposed method for solving the nonlinear Darcy–Forchheimer model.

Direct meshless local Petrov–Galerkin (DMLPG) method for time-fractional fourth-order reaction–diffusion problem on complex domains

Abstract: A new numerical scheme has been developed based on the fast and efficient meshless local weak form i.e direct meshless local Petrov–Galerkin (DMLPG) method for solving the fractional fourth-order partial differential equation on computational domains with complex shape. The fractional derivative is the Riemann–Liouville fractional derivative. At first, a finite difference scheme with the second-order accuracy has been employed to discrete the time variable. Then, the DMLPG technique is employed to achieve a full-discrete scheme. The time-discrete scheme has been studied in terms of unconditional stability and convergence order by the energy method in the L 2 space. Also, some numerical results are presented to show the efficiency and accuracy of the proposed technique on the simple and complex domains with the irregular and non-regular grid points.

A (2+1)-dimensional time-dependent Date–Jimbo–Kashiwara–Miwa equation: Painlevé integrability and multiple soliton solutions

Abstract: In this work, we extend the (2+1)-dimensional Date- Date–Jimbo–Kashiwara–Miwa equation to a new version with time-dependent coefficients. The integrability property of this developed equation is demonstrated by exhibiting the Painleve test. We show that this time-dependent equation admits multiple real and multiple complex soliton solutions. Moreover, we derive other solitonic and periodic solutions.

Projection-based reduced order models for a cut finite element method in parametrized domains

Abstract: This work presents a reduced order modeling technique built on a high fidelity embedded mesh finite element method. Such methods, and in particular the CutFEM method, are attractive in the generation of projection-based reduced order models thanks to their capabilities to seamlessly handle large deformations of parametrized domains and in general to handle topological changes. The combination of embedded methods and reduced order models allows us to obtain fast evaluation of parametrized problems, avoiding remeshing as well as the reference domain formulation, often used in the reduced order modeling for boundary fitted finite element formulations. The resulting novel methodology is presented on linear elliptic and Stokes problems, together with several test cases to assess its capability. The role of a proper extension and transport of embedded solutions to a common background is analyzed in detail.

A fast method for variable-order Caputo fractional derivative with applications to time-fractional diffusion equations

Abstract: In this paper, we propose a fast algorithm for the variable-order (VO) Caputo fractional derivative based on a shifted binary block partition and uniform polynomial approximations of degree r . Compared with the general direct method, the proposed algorithm can reduce the memory requirement from O ( n ) to O ( r log n ) storage and the complexity from O ( n 2 ) to O ( r n log n ) operations, where n is the number of time steps. As an application, we develop a fast finite difference method for solving a class of VO time-fractional diffusion equations. The computational workload is of O ( r m n log n ) and the active memory requirement is of O ( r m log n ) , where m denotes the size of spatial grids. Theoretically, the unconditional stability and error analysis for the proposed fast finite difference method are given. Numerical results of one and two dimensional problems are presented to demonstrate the well performance of the proposed method.

Heat transfer enhancement of Water-Al2O3 nanofluid in an oval channel equipped with two rows of twisted conical strip inserts in various directions: A two-phase approach

Abstract: This is a numerical study of convective heat transfer of the Water-Al2O3 nanofluid in an oval channel using two-phase mixture model. The channel is fitted with two rows of twisted conical strip inserts with various directions relative to each other leading to three different combinations of the mentioned inserts, namely inward Co-Conical inserts (CCI-inward), Counter-Conical inserts (CoCI), and outward Co-Conical inserts (CCI-outward) in which its lower wall is exposed to a constant heat flux. The effect of Reynolds number ranging from 250 to 1000, nanofluid volume fraction ranging from 1 to 3 % and conical strip insert combinations are examined on the fluid flow patterns and heat transfer characteristics. The results showed that among the three combinations of the twisted conical strip insert, CCI-inward locally presents the highest values of heat transfer coefficient, as about 17% higher than plain tube, considering the nature of the secondary flow created in this case. It is also found that the effect of increasing nanofluid concentration on the channel thermal performance is more significant at higher values of R e number; however, the pressure drop difference between the three models is subtle.

Temporal second order difference schemes for the multi-dimensional variable-order time fractional sub-diffusion equations

Abstract: A special point on each time interval is found for the approximation of the variable-order time Caputo derivative, which makes at least second order approximation accuracy be obtained. On this basis, two difference schemes are proposed for the multi-dimensional variable-order time fractional sub-diffusion equations, which have second order accuracy in time, second order and fourth order accuracy in space, respectively. The obtained difference schemes are proved to be uniquely solvable. The convergence and stability of the schemes in the discrete H 1 -norm are analyzed by utilizing the energy method. Some numerical examples are presented to verify the theoretical results.

A two-stage adaptive scheme based on RBF collocation for solving elliptic PDEs

Abstract: In this paper we present a new adaptive two-stage algorithm for solving elliptic partial differential equations via a radial basis function collocation method. Our adaptive meshless scheme is based at first on the use of a leave-one-out cross validation technique, and then on a residual subsampling method. Each of phases is characterized by different error indicators and refinement strategies. The combination of these computational approaches allows us to detect the areas that need to be refined, also including the chance to further add or remove adaptively any points. The resulting algorithm turns out to be flexible and effective through a good interaction between error indicators and refinement procedures. Several numerical experiments support our study by illustrating the performance of our two-stage scheme. Finally, the latter is also compared with an efficient adaptive finite element method.

Two meshless methods based on local radial basis function and barycentric rational interpolation for solving 2D viscoelastic wave equation

Abstract: In this paper, 2D viscoelastic wave equation is solved numerically both on regular and irregular domains. For spatial approximation of viscoelastic wave equation two meshless methods based on local radial basis function and barycentric rational interpolation are proposed. Both of the spatial approximation methods do not need mesh, node connectivity or integration process on local subdomains so they are truly meshless. For local radial basis function method we used an existing algorithm in literature to choose an acceptable shape parameter. Time marching is performed with fourth order Runge Kutta method. L 2 and L ∞ error norms for some test problems are reckoned to indicate efficiency and performance of the proposed two methods. Also, stability of the methods is discussed. Acquired results confirm the applicability of the proposed methods for 2D viscoelastic wave equation. We have performed some comparisons between the proposed two methods in the sense of accuracy and computational cost. From the comparisons, we have observed that performance of the barycentric rational interpolation in the sense of accuracy is slightly better than the performance of local radial basis function however computational cost of the local radial basis function is less than the computational cost of barycentric rational interpolation.

Linearized ADI schemes for two-dimensional space-fractional nonlinear Ginzburg–Landau equation

Abstract: Space and time approximations for two-dimensional space fractional complex Ginzburg–Landau equation are examined The schemes under consideration are discreted by the second-order backward differential formula (BDF2) in time and two classes of the fractional centered finite difference methods in space A linearized technique is employed by the extrapolation We prove the unique solvability and stability for both numerical methods The convergence of both numerical methods is analyzed at length utilizing the energy argument, and the convergence orders under the optimal step size ratio are O ( τ 2 + h 2 ) and O ( τ 2 + h 4 ) in the sense of the discrete L 2 -norm, where τ is the time step size, h = max { h x , h y } , and h x , h y are spatial grid sizes in the x -direction and y -direction, respectively In addition, we construct a multistep alternating direction implicit (ADI) scheme and a multistep compact ADI scheme based on BDF2 for the efficiently numerical implementation Finally, numerical examples are carried out to verify our theoretical results

Generalized finite difference method for solving stationary 2D and 3D Stokes equations with a mixed boundary condition

Abstract: In the present work, a generalized finite difference method (GFDM), a meshless method based on Taylor-series approximations, is proposed to solve stationary 2D and 3D Stokes equations. To overcome the troublesome pressure oscillation in the Stokes problem, a new simple formulation of boundary condition for the Stokes problem is proposed. This numerical approach only adds a mixed boundary condition, the projections of the momentum equation on the boundary outward normal vector, to the Stokes equations, without any other change to the governing equations. The proposed formulation can be easily discretized by the GFDM. The GFDM is evolved from the Taylor series expansions and moving-least squares approximation, and the derivative expressed of unknown variables as linear combinations of function values of neighboring nodes. Numerical examples are utilized to verify the feasibility of the proposed GFDM scheme not only for the Stokes problem, but also for more involved and general problems, such as the Poiseuille flow, the Couette flow and the Navier–Stokes equations in low-Reynolds-number regime. Moreover, numerical results and comparisons show that using the GFDM to solve the proposed formulation of the Stokes equations is more accurate than the classical formulation of the pressure Poisson equation.

An implicit finite difference analysis of magnetic swimmers propelling through non-Newtonian liquid in a complex wavy channel

Abstract: There are many unicellular tiny organisms which can self-propel collectively through non-Newtonian fluids by means of producing undulating deformation. Example includes nematodes, rod shaped bacteria and spermatozoa. Here we use Taylor’s swimming sheet model, with non-Newtonian fluid bounded with in a complex wavy walls of a two-dimensional channel. Oldroyd-4 constant fluid is approximated as cervical mucus and MHD effects are also considered. After utilizing lubrication and creeping flow assumption the reduced non-linear differential equation is solved (by implicit finite difference technique) so that it will satisfy the dynamic equilibrium condition for steady propulsion. For a special (Newtonian) case the expressions of swimming speed and flow rate are also presented. We also demonstrate that the rheological properties of non-Newtonian fluid can assist or resist the pack of micro-organisms (swimming sheet), while the larger undulation amplitude in swimmer’s body and magnetic field in downward direction can enhance the propulsion speed. The solution obtained via implicit finite difference method is also validated by a built in MATLAB routine bvp-4c. This built in function is based on collocation technique. Moreover an excellent correlation is achieved for both numerical methods.

Primal–dual weak Galerkin finite element methods for elliptic Cauchy problems

Abstract: The authors propose and analyze a well-posed numerical scheme for a type of ill-posed elliptic Cauchy problem by using a constrained minimization approach combined with the weak Galerkin finite element method. The resulting Euler–Lagrange formulation yields a system of equations involving the original equation for the primal variable and its adjoint for the dual variable, and is thus an example of the primal–dual weak Galerkin finite element method. This new primal–dual weak Galerkin algorithm is consistent in the sense that the system is symmetric, well-posed, and is satisfied by the exact solution. A certain stability and error estimates were derived in discrete Sobolev norms, including one in a weak L 2 topology. Some numerical results are reported to illustrate and validate the theory developed in the paper.

A mixed MoL–TMoL for the numerical solution of the 2D Richards’ equation in layered soils

Abstract: Water infiltration into layered soils is studied, considering a two dimensional spatial domain. The focus is on the treatment of discontinuity at the intersection of non-overlapping soils. The novelty of this paper is based on a mixed MoL–TMoL, which merges desirable features of both Method of Lines (MoL) and Transversal Method of Lines (TMoL); such a numerical approach allows us the numerical treatment of the solution at the discontinuous interfaces by means of Filippov theory for dynamical systems. Numerical simulations, based on implicit and semi-implicit schemes of low accuracy, are provided for validating this approach.

Fully-Discrete Explicit Locally Entropy-Stable Schemes for the Compressible Euler and Navier-Stokes Equations

Abstract: Recently, relaxation methods have been developed to guarantee the preservation of a single global functional of the solution of an ordinary differential equation. Here, we generalize this approach to guarantee local entropy inequalities for finitely many convex functionals (entropies) and apply the resulting methods to the compressible Euler and Navier–Stokes equations. Based on the unstructured h p -adaptive SSDC framework of entropy conservative or dissipative semidiscretizations using summation-by-parts and simultaneous-approximation-term operators, we develop the first discretizations for compressible computational fluid dynamics that are primary conservative, locally entropy stable in the fully discrete sense under a usual CFL condition, explicit except for the parallelizable solution of a single scalar equation per element, and arbitrarily high-order accurate in space and time. We demonstrate the accuracy and the robustness of the fully-discrete explicit locally entropy-stable solver for a set of test cases of increasing complexity.

An alternating direction implicit Galerkin finite element method for the distributed-order time-fractional mobile–immobile equation in two dimensions

Abstract: In this paper, we shall present the alternating direction implicit (ADI) Galerkin finite element method (FEM) for solving the distributed-order time-fractional mobile–immobile equation in two dimensions. In the time direction, the backward Euler method is used to deal with the temporal first-order derivative, and the weighted and shifted Grunwald formula is employed to discretize the distributed-order time-fractional derivative. Galerkin FEM is used for discretization of the spatial direction, and then an ADI Galerkin finite element scheme is constructed. The stability and L 2 -norm convergence are proved. Several numerical experiments are provided to verify our theoretical analysis.

A C1 Virtual Element Method on polyhedral meshes

Abstract: The purpose of the present paper is to develop C 1 Virtual Elements in three dimensions for linear elliptic fourth order problems, motivated by the difficulties that standard conforming Finite Elements encounter in this framework. We focus the presentation on the lowest order case, the generalization to higher orders being briefly provided in the Appendix. The degrees of freedom of the proposed scheme are only 4 per mesh vertex, representing function values and gradient values. Interpolation error estimates for the proposed space are provided, together with a set of numerical tests to validate the method at the practical level.

Stabilized reduced basis methods for parametrized steady Stokes and Navier–Stokes equations

Abstract: It is well known in the Reduced Basis approximation of saddle point problems that the Galerkin projection on the reduced space does not guarantee the inf–sup approximation stability even if a stable high fidelity method was used to generate snapshots. For problems in computational fluid dynamics, the lack of inf–sup stability is reflected by the inability to accurately approximate the pressure field. In this context, inf–sup stability is usually recovered through the enrichment of the velocity space with suitable supremizer functions. The main goal of this work is to propose an alternative approach, which relies on the residual based stabilization techniques customarily employed in the Finite Element literature, such as Brezzi–Pitkaranta, Franca–Hughes, streamline upwind Petrov–Galerkin, Galerkin Least Square. In the spirit of offline–online reduced basis computational splitting, two such options are proposed, namely offline-only stabilization and offline–online stabilization. These approaches are then compared to (and combined with) the state of the art supremizer enrichment approach. Numerical results are discussed, highlighting that the proposed methodology allows to obtain smaller reduced basis spaces (i.e., neglecting supremizer enrichment) for which a modified inf–sup stability is still preserved at the reduced order level.

Asymptotic stability of a stochastic May mutualism system

Abstract: This paper deals with a mutualism system in random environments, in which the cooperation of two species is guaranteed depending on increasing the carrying capacity of each other. Mathematically, we obtain the existence and uniqueness of a stable stationary distribution by means of Markov semigroup theory and Fokker–Planck equation. In addition, we prove that densities of the distributions for the solutions converge in L 1 to an invariant density under appropriate conditions. Numerical simulations are also presented to illustrate the theoretical results.

Unconditional superconvergence analysis of a linearized Crank–Nicolson Galerkin FEM for generalized Ginzburg–Landau equation

Abstract: In this paper, we study an effective numerical method for the generalized Ginzburg–Landau equation (GLE). Based on the linearized Crank–Nicolson difference method in time and the standard Galerkin finite element method with bilinear element in space, the time-discrete and space–time discrete systems are both constructed. We focus on a rigorous analysis and consideration of unconditional superconvergence error estimates of the discrete schemes. Firstly, by virtue of the temporal error results, the regularity for the time-discrete system is presented. Secondly, the classical Ritz projection is used to obtain the spatial error with order O ( h 2 ) in the sense of L 2 -norm. Thanks to the relationship between the Ritz projection and the interpolated projection, the superclose estimate with order O ( τ 2 + h 2 ) in the sense of H 1 -norm is derived. Thirdly, it follows from the interpolated postprocessing technique that the global superconvergence result is deduced. Finally, some numerical results are provided to confirm the theoretical analysis.

On the application of isogeometric finite volume method in numerical analysis of wet-steam flow through turbine cascades

Abstract: The isogeometric finite volume analysis is utilized in this research to numerically simulate the two-dimensional viscous wet-steam flow between stationary cascades of a steam turbine for the first time. In this approach, the analysis-suitable computational mesh with “curved” boundaries is generated for the fluid flow by employing a non-uniform rational B-spline (NURBS) surface that describes the cascade geometry, and the governing equations are then discretized by the NURBS representation. Thanks to smooth and accurate geometry representation of the NURBS formulation, the employed isogeometric framework not only resolves issues concerning the conventional mesh generation techniques of the finite volume method in steam turbine problems, but also, as validated against well-established experiments, significantly improves the accuracy of the numerical solution. In addition, the shock location in the cascade is predicted and tracked with a sufficient accuracy.

A Hybrid High-Order method for the incompressible Navier–Stokes problem robust for large irrotational body forces

Abstract: We develop a novel Hybrid High-Order method for the incompressible Navier–Stokes problem robust for large irrotational body forces. The key ingredients of the method are discrete versions of the body force and convective contributions in the momentum equation formulated in terms of a globally divergence-free velocity reconstruction. Two key properties are mimicked at the discrete level, namely the invariance of the velocity with respect to irrotational body forces and the non-dissipativity of the convective term. A full convergence analysis is carried out, showing optimal orders of convergence under a smallness condition involving only the solenoidal part of the body force. The performance of the method is illustrated by a complete panel of numerical tests, including comparisons that highlight the benefits with respect to more standard formulations.