I like your summary, Youssef Hafez.

To add: one could argue how relatively important local conservation is compared to higher-order accuracy, for various applications, such that there are some flow problems for which FV methods would be preferred. The latter can indeed occur in some (extremely) high-Re flow over terrain in atmospheric \ wind applications.

Thank you, Mark — I agree with your point.

Indeed, the trade-off between local conservation and higher-order accuracy is very application-dependent. In cases such as high-Reynolds number flow over complex terrain, particularly in atmospheric and wind engineering applications, strict local conservation of mass, momentum, and energy — as naturally enforced by the Finite Volume Method — can be crucial to maintaining numerical stability and physical realism over long simulation times.

That said, FEM’s strength lies in its ability to minimize the residual in a variational sense and flexibly handle irregular geometries with higher-order basis functions, which can be very advantageous in environmental, hydraulic, and structural-fluid interaction problems.

In practice, I see FEM and FVM as complementary tools rather than competing ones — with the “better” choice driven by the physics of the problem, the computational resources, and the desired balance between accuracy, conservation, and robustness.

You are wrong… FVM is a FEM for a specific choice of the weight function.

I am not stating a new fact, this is the theory, see for example the textbook of Leveque.

Accuracy order is something you can get both in FEM and FVM as you want, with the due computational effort.

However, FEM is generally expensive, at present DG formulation has received more attention.

Mark Kelly

conservation in discrete sense is a fundamental property of numerical schemes since ensure proper velocity propagation of travelling waves.

Filippo Maria Denaro

Dear Prof. Denaro,

Thank you for sharing your thoughts on the similarities between FEM and FVM. While I appreciate the theoretical connections between the two, saying they are the same is—if I may—like saying water is the same as ice: chemically identical (H₂O), but physically with very different properties and behavior.

I am aware that FVM can be derived as a special case of FEM under specific choices of weighting functions. However, in actual numerical implementation and performance, FEM and FVM are distinct. These differences become especially important in applied computational work, where accuracy, convergence behavior, and computational cost can vary substantially.

A fair comparison should ideally involve the same problem with the same number of cells/elements. For example, in Kadia et al. (2024) (DOI: 10.1061/JHEND8.HYENG-13837), OpenFOAM (FVM) was used to study secondary currents, turbulence, and shear stresses in supercritical narrow open channel flows. For an aspect ratio of 1.0, their mesh had 62 × 62 cells, required 9,659 iterations to converge, and took 18.6 minutes of execution time. In contrast, using the Finite Element Method (Hafez, 1995) for a channel of the same aspect ratio, I achieved convergence with a coarser 20 × 20 element mesh in just 570 iterations, with normalized residuals well below the normalized velocity tolerance of 0.0001. Unfortunately, the OpenFOAM study did not report residuals, so a direct error comparison is not possible, but the differences in iteration count, grid size, and computational efficiency are clear.

For these reasons, while the mathematical relationship between FEM and FVM is well established, their practical numerical characteristics—and thus their performance in real-world hydrodynamic and turbulence simulations—can be markedly different.

Reference: Hafez, Y.I. (1995). A k–ε Turbulence Model for Predicting the Three-Dimensional Velocity Field and Boundary Shear in Open and Closed Channels. Civil Engineering Department, Colorado State University, Fort Collins, Colorado, USA.

Dear Youssef Hafez

a well assessed comparison requires to adopt more classical test cases in CFD. No complex physics and turbulence modelling.

Then, FEM and FVM must be compared using the same level of theoretical accuracy. During the last years DG method received more attention, however they are computationally more expensive than FVM.

I suggest to produce a comparison between FEM and FVM for the case of the Taylor solution, using a scheme with the same theorectical accuracy and comparing the computational cost.

Again, discrete conservation is mandatory.

Filippo Maria Denaro

Dear Prof. Denaro,

Thank you for your reply and for highlighting the importance of classical test cases in CFD comparisons.

My earlier comment—that FEM yields the minimum errors or residuals, particularly when the problem can be cast in a variational formulation—was referring to a theoretical property of the method. This result is well established in numerical analysis and does not require numerical computation to verify; it follows directly from the Galerkin orthogonality principle and is widely recognized in the literature.

Of course, I have no objection to comparing FEM and FVM for standard benchmark problems—such as duct flow, cavity flow, backward-facing step, flow around bluff bodies, or heat transport—provided they are implemented with the same order of accuracy and on the same mesh resolution. Such tests can be useful to quantify practical differences in accuracy, convergence rate, and computational cost.

Regarding discrete conservation, I fully agree it is necessary, but it is not by itself sufficient to guarantee accuracy. The fact that the flux leaving one cell face exactly equals the flux entering the adjacent cell is correct, but if that flux is a poor approximation of the true physical flux, then conserving it exactly is of little value. In practice, the accuracy of the approximated flux is far more important than the fact that it is conserved at the discrete level—particularly in simulations where small flux errors can accumulate into large physical errors.

For this reason, while conservation is an important requirement, it must go hand-in-hand with the minimization of residuals and the accurate representation of the underlying physics—areas where FEM has well-documented theoretical advantages in problems with a variational structure.

Best regards,

It is true that classical FEM can be more computationally expensive than FVM, and this has sometimes led practitioners to prefer faster, less costly methods such as FVM, often justifying the choice on the basis of its discrete conservation property. However, discrete conservation alone is not sufficient to guarantee accuracy. The fact that the flux leaving one cell face equals the flux entering the adjacent cell is fine, but if this flux is a poor approximation to the true physical flux, then conserving it perfectly is of little practical benefit. Accuracy in the flux computation is far more important, and this is where FEM has a theoretical advantage—particularly if the number of elements is increased to match the resolution used in typical FVM simulations.

If FEM is allowed the same number of cells/elements as FVM, the results are often superior in accuracy, stability, and convergence. The argument that FEM is “too slow” is becoming less relevant as computational power and speed continue to grow. Furthermore, FEM has developed conservative variants—such as the Mixed Finite Element Method (MFEM), Control Volume Finite Element Method (CVFEM), and Discontinuous Galerkin (DG) method—that maintain local discrete conservation while preserving FEM’s high-order accuracy and geometric flexibility.

For these reasons, I believe that while computational cost is a practical consideration, it should not, in itself, deter the use of FEM—particularly when accuracy and robustness are critical, and when modern hardware can mitigate its higher cost.

yes, I agree except the concept of accuracy... In FEM a lot of theoretical results are obtained with the constraint of the linear problem.

Furthermore, accuracy order I am discussing for FVM is related to the power p in O(hp) in the local truncation error of the numerical flux, that is not the number of cells. On a same grid size h, I can get different order of accuracy.

I do not consider OpenFOAM a proper reference for the performance of FVM.

Filippo Maria Denaro Dear Prof. Denaro,

Thank you for your thoughtful points. In addition to the important issues of flux conservation and formal accuracy order, I believe there is an even more fundamental matter to consider.

The governing differential equations for fluid flow—such as the Navier–Stokes equations—are derived for infinitesimal control volumes, where physical diffusion can always balance advective transport, no matter how large the latter is. However, when we discretize these equations numerically, our cells or elements are finite, not infinitesimal. This means that the discrete numerical equations are already an approximation of the true physics from the outset, even before any turbulence modeling or numerical scheme choices are made.

Because of this, evaluating numerical methods solely by comparing their final results against analytical solutions or measured data, while necessary, is far from sufficient. This is akin to dividing 64 by 16, incorrectly “canceling” the 6 in numerator and denominator, and still arriving at the correct answer of 4—right result, wrong reasoning. In numerical modeling, it is the validity of the intermediate steps and the internal consistency of the method that truly matter, not just whether the final number matches an experiment.

Judging a numerical method or turbulence model solely by how well it matches end results is, therefore, much like judging a book by its cover—it can be very misleading. Robust evaluation should include scrutiny of the underlying discretization and how faithfully it preserves the intended physics at the discrete level, not just whether the output “looks right.”

Best regards,

To tell the truth, I am not sure to understand what you mean in your comment.

Any finite volume, no matter of the size, fulfills the physics properties of the integral equation.

Finite element and finite volume methods give identical results. Finite element method needs book keeping of grids (coordinate as well as connectivity).

Filippo Maria Denaro Dear Prof. Denaro

Let me restate my point using the one-dimensional viscous diffusion term in the Navier–Stokes equations as an example, without relying on symbolic equations.

In the one-dimensional incompressible case, the viscous term is proportional to the viscosity multiplied by the second spatial derivative of velocity. This term comes from assuming a Newtonian fluid with constant viscosity (the so-called Boussinesq approximation) and is derived from the net viscous flux into an infinitesimal control volume.

1. How the “second derivative” actually appears

Imagine an element or control volume that starts at a position we call x and ends at x plus delta-x. The viscous flux into the element at the left face is proportional to the velocity gradient at x, while the flux out at the right face is proportional to the velocity gradient at x plus delta-x.

The net viscous effect inside the element is the difference between these two fluxes, divided by the element length delta-x. When we let delta-x shrink to an infinitesimal value, this difference-over-length becomes the second derivative of velocity. Importantly, this second derivative is really located at the center of the element, at x plus half delta-x. When the element is infinitesimally small, the center and the node location are essentially the same, so the standard Navier–Stokes form is exact.

2. Why element size matters

If the element length delta-x is not sufficiently small, the velocity curvature (second derivative) at the element center is not the same as at the node. In that case, the standard second-derivative term no longer exactly represents the true net viscous flux for that element.

This is a physical-model limitation, not just a numerical truncation error. The issue is not that the integral form of the equations fails—the integral form can hold for any volume—but that the differential form is only strictly valid in the limit of infinitesimally small volumes. In numerical simulations, our finite cells do not satisfy this limit, so the accuracy of the viscous term depends directly on element size.

3. Extension to other equations

The same reasoning applies to heat conduction in the heat equation, where the thermal diffusion term is proportional to the thermal conductivity multiplied by the second spatial derivative of temperature. The second derivative form is exact only in the infinitesimal-volume limit.

4. Adjusting for finite-size effects

Some numerical schemes attempt to compensate for this mismatch between physical diffusion scales and numerical cell size by introducing explicit time or length-scale corrections—examples include artificial diffusion, delayed-correction schemes, and certain forms of upwinding with a physical interpretation. These are essentially patches to recover accuracy lost when elements are not small enough.

5. On the shear stress model in Navier–Stokes

It is also important to remember that the Navier–Stokes equations assume a linear relationship between shear stress and velocity gradient. This is valid for Newtonian fluids at moderate strain rates.

However, non-linear stress–strain models can be essential in more complex situations. In particular, they are needed to correctly predict the formation and strength of secondary current vortex structures that arise in the cross-sections of non-circular ducts—such as square or rectangular ducts—and in open channels. These secondary currents are driven by anisotropy in the turbulence and cannot be fully captured by the linear Newtonian formulation.

In short: The integral form of the equations is valid for any finite volume, but the differential form—where viscosity multiplies the second derivative of velocity—is only exact in the infinitesimal limit. When applied to finite-sized numerical elements, the accuracy of this term depends on element size, regardless of whether the method is finite volume, finite element, or finite difference. Furthermore, when dealing with complex secondary flows such as vortex structures in non-circular ducts and open channels, adopting non-linear stress–strain models becomes critical for capturing the correct physics.

Dear Youssef Hafez

here you are my replies (R) about your topic.

_______

Let me restate my point using the one-dimensional viscous diffusion term in the Navier–Stokes equations as an example, without relying on symbolic equations.

In the one-dimensional incompressible case, the viscous term is proportional to the viscosity multiplied by the second spatial derivative of velocity. This term comes from assuming a Newtonian fluid with constant viscosity (the so-called Boussinesq approximation) and is derived from the net viscous flux into an infinitesimalcontrol volume.

1. How the “second derivative” actually appears

Imagine an element or control volume that starts at a position we call x and ends at x plus delta-x. The viscous flux into the element at the left face is proportional to the velocity gradient at x, while the flux out at the right face is proportional to the velocity gradient at x plus delta-x.

The net viscous effect inside the element is the difference between these two fluxes, divided by the element length delta-x. When we let delta-x shrink to an infinitesimal value, this difference-over-length becomes the second derivative of velocity. Importantly, this second derivative is really located at the centerof the element, at x plus half delta-x. When the element is infinitesimally small, the center and the node location are essentially the same, so the standard Navier–Stokes form is exact.

R. What you are addressing is just the mathematical approach to derive the pointwise differential form of the NSE. The physical conservation laws start from the Reynolds theorem for a finite control volume. The differential form requires strong assumption on regularity, the physical integral laws not. The diffusive flux 2*mu*S is the Newtonian model, there are modern proposals where that is not introduced and the equations remain hyperbolic, not parabolic. Infinitesimal concept is mathematics, in physics we have the continuum model that allows us to express the equations. But the variable we introduce are averaged on a small but finite volume.

2. Why element size matters

If the element length delta-x is not sufficiently small, the velocity curvature (second derivative) at the element center is not the same as at the node. In that case, the standard second-derivative term no longer exactly represents the true net viscous flux for that element.

This is a physical-model limitation, not just a numerical truncation error. The issue is not that the integral form of the equations fails—the integral form can hold for any volume—but that the differential form is only strictly valid in the limit of infinitesimally small volumes. In numerical simulations, our finite cells do not satisfy this limit, so the accuracy of the viscous term depends directly on element size.

R. Not at all, the diffusive flux is integrated over a face of the volume of finite size. This is exact! Forget the second derivative. Approximation appears only when we discretize the flux and the surface integral. That is only in the sense of the truncation error.

What is relevant, and you seem not to consider, is the action of the diffusion at physical scales. Diffusive flux starts to act at the Taylor microscale and ends up to the Kolmogorov length scale, a finite length. There is no infinitesimal physical diffusive scale.

If you are concerned that your pointwise PDE are approximated by a discrete set of equation on a finite grid size h, then just use h enough fine to be compared to the Kolmogorov length (DNS). Any other greater grid size will produce an implicit filtering (Nyquist grid cut-off plus effects of discretization) and you are in the field of the filtered (regularized) equations (LES). It not a coincidence that the integral form of the NSE are also theoretically linked to the volume-filtered LES equations.

3. Extension to other equations

The same reasoning applies to heat conduction in the heat equation, where the thermal diffusion term is proportional to the thermal conductivity multiplied by the second spatial derivative of temperature. The second derivative form is exact only in the infinitesimal-volume limit.

4. Adjusting for finite-size effects

Some numerical schemes attempt to compensate for this mismatch between physical diffusion scales and numerical cell size by introducing explicit time or length-scale corrections—examples include artificial diffusion, delayed-correction schemes, and certain forms of upwinding with a physical interpretation. These are essentially patches to recover accuracy lost when elements are not small enough.

R. there are theoretical and hystorical papers from Los Alamos researchers (see a review of L. Margolin) that show how artificial diffusion is a natural consequence of the finite length size, in the physical sense, not as an approach to compensate the numerical approximation. In other words, the regularized variable on a finite volume size is naturally smoothed at high wavenumbers.

All you are addressing for PDE is just the topic about the local truncation error analysis. If you convert this in the PDE on the finite length scale then you are discussing the meaning of the integral conservation laws.

In conclusion, you observations should be extended to a more general framework.

5. On the shear stress model in Navier–Stokes

It is also important to remember that the Navier–Stokes equations assume a linearrelationship between shear stress and velocity gradient. This is valid for Newtonian fluids at moderate strain rates.

However, non-linear stress–strain models can be essential in more complex situations. In particular, they are needed to correctly predict the formation and strength of secondary current vortex structures that arise in the cross-sections of non-circular ducts—such as square or rectangular ducts—and in open channels. These secondary currents are driven by anisotropy in the turbulence and cannot be fully captured by the linear Newtonian formulation.

R. Here you are discussing about a totally different topic, that is the turbulence modelling in RANS formulation. That adds nothing to the previous discussion about the physics since here the variables are statistically averaged.

In short: The integral form of the equations is valid for any finite volume, but the differential form—where viscosity multiplies the second derivative of velocity—is only exact in the infinitesimal limit. When applied to finite-sized numerical elements, the accuracy of this term depends on element size, regardless of whether the method is finite volume, finite element, or finite difference. Furthermore, when dealing with complex secondary flows such as vortex structures in non-circular ducts and open channels, adopting non-linear stress–strain models becomes critical for capturing the correct physics.

Filippo Maria Denaro

Dear Prof. Denaro,

Thank you for your detailed comments. I agree that starting from the Reynolds transport theorem and working in the integral form is the most rigorous physical basis for the Navier–Stokes equations, and that finite volume methods respect these conservation laws for any control volume.

However, my main point is not about whether the integral or differential form is the “right” starting point, but about how accuracy is achieved in the discrete representation of the equations.

In finite difference and finite volume methods, the approximation of derivatives or fluxes relies on local algebraic differences. The accuracy of these methods depends directly on grid size and on how well the discrete stencil mimics the continuous operator. Errors are controlled by truncation analysis, but there is no inherent mechanism ensuring that the discrete residual is minimized globally.

By contrast, the finite element method introduces the orthogonality condition through the weighted residual (Galerkin) formulation. This ensures that the error in the approximate solution is orthogonal to the chosen function space. In other words, FEM does not only conserve fluxes—it minimizes the residual in a mathematically optimal sense across the domain. This additional property is what often makes FEM more accurate than FVM or FDM, particularly when dealing with diffusion-dominated problems or complex geometries.

So while all three methods conserve fluxes when properly implemented, FEM achieves something extra: the orthogonality condition provides a stronger consistency framework, ensuring that the discrete equations approximate the continuous physics in a more systematic and accurate way.

This is the main argument I wanted to highlight, and why I believe FEM has an intrinsic accuracy advantage in certain cases.

Dear Youssef Hafez

what you are not understanding is that in FVM the residual is projected onto the space of piecewise-constant test functions.

In this sense, Finite Volume can be interpreted as a Petrov–Galerkin method, where the test functions are cell indicators.

Check my old paper on Numerical Heat Transfer.

There is no superiority or inferiorit in that.

Accuracy order depends on the discretiziation.

Filippo Maria Denaro

Dear Prof. Denaro

In the Finite Element Method, the residual is not left uncontrolled. Instead, it is forced to be orthogonal to the chosen set of weighting or test functions. In the classical Galerkin formulation, where the test functions are taken from the same space as the trial (solution) functions, this orthogonality condition ensures that the residual is minimized in a well-defined mathematical sense. For many problems, particularly elliptic ones, this leads to the so-called minimum residual property in an energy norm. Even in the Petrov–Galerkin case, where the test space differs from the trial space, the residual still satisfies an orthogonality condition with respect to that chosen test space, so its behavior is systematically controlled.

In contrast, the Finite Volume Method does not multiply the residual by a family of test functions. Instead, it integrates the governing equations over each control volume and enforces the balance of fluxes. This means that the residual is forced to average to zero over the control volume, which guarantees local conservation. However, this averaging is not equivalent to an orthogonality condition. The residual inside a control volume can have oscillations or variations that are not constrained, as long as their average is zero. Therefore, the residual is not minimized in any particular norm, and there is no guarantee of orthogonality beyond the trivial case of constant test functions.

This highlights a key conceptual difference. In the Finite Element Method, residual orthogonality leads to systematic control and, in many cases, minimum residual solutions. In the Finite Volume Method, conservation is guaranteed, but residual minimization and orthogonality are not. Algebraically, one can show that some Finite Volume schemes are equivalent to a Petrov–Galerkin Finite Element formulation with piecewise constant test functions. Yet, with such simple test functions, the minimum residual and orthogonality properties degenerate to a trivial case.

In summary, the Finite Element Method and the Finite Volume Method both stem from the same governing differential equations but enforce their discrete forms differently. FEM emphasizes residual orthogonality and minimization, while FVM emphasizes conservation through residual averaging. This distinction explains why FEM solutions can sometimes be more systematically accurate, whereas FVM solutions are strictly conservative but may not enjoy the same minimum residual property.

Filippo Maria Denaro

Dear Prof. Denaro,

I appreciate your clarification and I do acknowledge that FVM can be cast in a Petrov–Galerkin framework with piecewise-constant test functions, ensuring conservation by making the average residual vanish in each control volume. This is indeed a valid interpretation.

However, it is equally important to highlight that this is not equivalent to the Galerkin orthogonality principle as realized in the finite element framework. In FEM, the residual is enforced to be orthogonal to a chosen functional space of test functions, which implies a minimization property of the error in an appropriate norm. This orthogonality condition is what provides FEM with a rigorous error control mechanism and convergence guarantees.

By contrast, FVM’s use of constant test functions does not enforce such orthogonality. It guarantees conservation, but it does not minimize the residual in a functional sense. This is a fundamental difference that cannot be overlooked.

Moreover, modern FEM variants (such as mixed FEM, stabilized FEM, and discontinuous Galerkin methods) are capable of enforcing local conservation, so the classical argument that “FEM is not conservative” is not a fundamental limitation of the method. When conservation and residual minimization are combined, FEM clearly provides a mathematically superior framework compared to standard FVM.

With these distinctions clarified, I think we have reached the point where all the essential issues are on the table. Further debate would likely only circle back over the same ground. If you remain unconvinced, I propose that we allow the readers of this forum and the ResearchGate discussion to draw their own conclusions.

Respectfully, Youssef Hafez

In FVM, residual is orthogonal to the piecewice test functions.

There are people expert in FEM that can express the theory better than me. For example Antonella Abbà Antonella Abbà

Filippo Maria Denaro

Dear Prof. Denaro,

Thank you again for your comments. I believe there may have been some oversight, since in my earlier reply I already explained the fundamental difference: the Finite Volume Method, in its classical form, only applies a trivial test function (the constant function), whereas the Finite Element Method enforces orthogonality of the residual with respect to a whole functional space of test functions.

This distinction is not something that requires running an application case—it is already fully evident in the first principle equations of both methods:

• In the Finite Volume Method, as you also illustrated in the figure you attached (on CFD-Online under the same discussion topic), the residual is integrated over a control volume and forced to vanish only in the sense of its average. This is equivalent to testing the residual only against the constant function. No higher-order orthogonality is imposed, and therefore there is no minimization of the residual in any energy norm.
• In the Finite Element Method (Galerkin), the residual is multiplied by every function in the chosen test space (which, in the Galerkin case, is the same as the trial space). The resulting condition means that the residual is orthogonal to an entire finite-dimensional space of functions, not just to constants. This orthogonality is exactly what leads to systematic error minimization and the well-known convergence guarantees of FEM.

Thus, while one can algebraically reinterpret FVM as a Petrov–Galerkin method with piecewise-constant test functions, this immediately shows the limitation: the orthogonality condition degenerates to the trivial case. FEM, by contrast, enforces orthogonality in a much richer space, which provides true error control.

That is why saying that FVM has the same orthogonality property as FEM is misleading. FVM guarantees local conservation, but FEM guarantees both conservation (when designed accordingly) and residual minimization through genuine orthogonality. This is a fundamental difference that arises from the formulation itself, not from numerical experiments.

Respectfully, Youssef Hafez