How to Calculate Mass Conservation and Energy Balance

Fabio Bocchi October 14, 2015
Share this on Facebook Share this on Twitter Share this on LinkedIn

As a technical support engineer, one of the most common technical questions I receive is: “How can I compute the mass conservation of a fluid flow simulation or the energy balance of a conjugate heat transfer simulation?” This is often requested to investigate and ensure a simulation’s accuracy. Here, I will demonstrate how to perform these calculations in COMSOL Multiphysics and introduce some of the predefined variables available for postprocessing the energy rate terms of the energy balance equation.

Let’s Start with Mass Conservation

To demonstrate the different topics covered in this blog post, I’ll use the example of an aluminum heat sink, which is often used to cool off electrical devices by dissipating heat. This tutorial model is available in its steady-state version in the COMSOL Multiphysics Application Library if you have the Heat Transfer Module or CFD Module.

The heat sink is made of aluminum, shaped with a cluster of pillars for cooling, and mounted on a chip that is made of a silica glass material. In the model setup, the heat sink sits inside of a rectangular channel with an inlet and outlet for air flow. The chip acts as a heat source and generates 1 watt.

Image showing the heat sink geometry.
The geometry of a basic heat sink.

In fluid dynamics, mass conservation leads to a well-known local continuity equation:

\frac{\partial\rho}{\partial t}+\nabla \cdot (\rho {\bf{u}}) = 0

By integrating this equation over the fluid domain and applying the divergence theorem, we obtain the global formulation of the mass conservation.

\int_\Omega \frac{\partial\rho}{\partial t}+\nabla \cdot (\rho {\bf{u}}) \ dV = \int_\Omega \frac{\partial\rho}{\partial t} \ dV +\int_{\partial\Omega} \rho{\bf{u}}\cdot {\bf{n}} \ dS


\int_\Omega \frac{\partial\rho}{\partial t} \ dV +\int_{\partial\Omega} \rho{\bf{u}}\cdot {\bf{n}} \ dS = 0

Let’s have a more specific look at the equation above. Whenever you are modeling a fluid flow simulation, you can evaluate this equation in order to check the mass conservation accuracy of your model. In any steady-state analysis, this equation reduces to \int_{\partial\Omega} \rho{\bf{u}} \cdot {\bf{n}} \ dS = 0 and states that the rate at which mass enters a system is equal to the rate at which mass leaves the system. In other words, the mass flow at inlets and outlets must balance.

A common mistake is to assume that mass conservation reduces to the volumetric flow rate conservation \int_{\partial\Omega} {\bf{u}}\cdot {\bf{n}} \ dS = 0. If the fluid density is constant, as for incompressible flow, the continuity equation simplifies to \nabla \cdot {\bf{u}} = 0, which states that the divergence of the flow velocity vanishes. In the case of incompressible flow, the assumption is then right. However, in most engineering problems, this assumption cannot be considered. Roughly speaking, in the case of motionless boundaries, we can consider it as long as the density does not depend on any variable (such as absolute pressure, temperature, concentration, etc.) that can lead to a change in density along the streamlines.

That said, let’s check the mass conservation accuracy of a fluid flow simulation within COMSOL Multiphysics. For this purpose, we can use the derived integration values — a useful postprocessing feature for calculating integrated quantities for each solution in a data set. To create a derived integration value, we go to Results > Derived Values > Integration > Surface Integration or Volume Integration. Then, we select the boundaries or volumes on which the integral should be carried out. Finally, we type the expression to be integrated.

Here, the expression for both inlet and outlet mass flow rates is written as: spf.rho*(u*nx+v*ny+w*nz), as shown in the figure below. spf.rho denotes the variable used to define the density of the fluid; (u,v,w) is the velocity field vector; and (nx,ny,nz) is the normal vector at the selected boundaries automatically computed by COMSOL Multiphysics. The volume mass flow rate can be expressed as d(spf.rho,t). The d operator is intended to perform a derivation of a variable with respect to another variable. Here, we take the time derivative of the fluid density.

Screenshot displaying the Derived Values functionality in COMSOL Multiphysics.
Step calculation for each mass flow rate term from the mass conservation equation.

In the static case of the heat sink simulation, the mass flow rates at the inlet and outlet boundaries are depicted below. The relative mass difference between the inlet and outlet is around 1e-5, which is less than the relative tolerance setting for the solver, which is set to 0.001. The mass conservation is therefore pretty accurate.

Screen capture showing the mass balance results of a stationary study.
Mass balance results for a stationary study.

We can also perform a time-dependent analysis of the same model. The total time of the simulation is set to be 20 minutes and the base of the heat sink experiences 1 watt of heat flux during the entire simulation time. The results of this simulation can be seen in the following animation. We set the simulation time to be long enough to reach both the thermal and fluid flow steady state.

The fluid flow streamlines are colored with the velocity norm and the temperature is plotted on the boundaries.

The mass balance results for the transient analysis are given in the following figure.

Image showing the mass balance results of a time-dependent study to calculate mass conservation.
Mass balance results for a time-dependent study.

The accuracy is still acceptable. Indeed, the mass balance values are well below 5e-7. The absolute solver tolerance specified in this simulation is 1e-6.

Calculating the Energy Balance

From the first law of thermodynamics and mechanical laws, the well-known global heat balance equation can be derived.


Here, Q_\mathrm{{exchange}} refers to the exchanged heat rate, taking into account conductive heat flux, -k\nabla T; radiative heat flux, {\bf{q}}_\mathrm{{rad}}; and additional heat sources, Q; such as electromagnetic heat sources (Joule heating), induction heating, or any user-defined heat source. P_\mathrm{{stress}} represents the stress power induced by mechanical stresses and E_\mathrm{{system}}=\int_\Omega E \ dm is the internal energy.

The stress power involved in this equation is converted into heat dissipation. The stress power expression is derived from the theory of continuum mechanics and can be written as

P_\mathrm{{stress}}=\int_\Omega (\sigma: {\bf{D}}) \ dV

where \sigma is the Cauchy stress tensor and {\bf{D}}= \frac{1}{2}(\nabla {\bf{u}} + \nabla {\bf{u}}^T) is the strain rate tensor.

In the case of a fluid, the stress tensor can be divided into a pressure part and a viscous part, \sigma=-p{\bf{I}}+\tau. The stress tensor then becomes the sum of work through pressure change and a viscous dissipation term as follows.

P_\mathrm{{stress}}=-\int_\Omega p(\nabla \cdot {\bf{u}}) \ dV+\int_\Omega (\tau : \nabla {\bf{u}}) \ dV

In COMSOL Multiphysics, we can either choose to add one, both, or none of these effects. Through the Non-Isothermal Flow multiphysics node, a check box for each effect is available and unchecked by default.

Screenshot showing how to access flow heating features in COMSOL Multiphysics.
The Include work done by pressure changes and Include viscous dissipation features.

In case of a conjugate heat transfer analysis, which is when the heat transfer equation is solved together with the Navier-Stokes and continuity equations, the following total energy flux becomes the conserved quantity.

{\bf{e}}_\mathrm{tot}=\rho{\bf{u}}E_0-k{\bf{\nabla}} T+{\bf{q}}_\mathrm{rad}-\sigma{\bf{u}}

where E_0=E+\frac{1}{2}{\bf{u}}\cdot {\bf{u}} is the total internal energy.

The total energy flux accounts for convective, conductive, and radiation heat fluxes. It contains additional terms that represent the convective kinetic energy, \rho\frac{{\bf{u}}}{2}({\bf{u}}\cdot {\bf{u}}), and the convective stress energy, \sigma{\bf{u}}. The energy balance equation then takes the following form:

\frac{d}{dt}\int_\Omega \rho E_0 \ dV +\int_{\partial\Omega} {\bf{e}}_\mathrm{tot} \cdot {\bf{n}} \ dS = \int_\Omega Q \ dV

In a stationary study, this expression reduces to

\int_{\partial\Omega} {\bf{e}}_\mathrm{tot} \cdot {\bf{n}} \ dS = \int_\Omega Q_ \ dV

For each quantity of this equation, a predefined global variable is available for postprocessing. The derived global evaluation values can be used to evaluate the variables. The table below summarizes the names of the different predefined variables of interest.

Total Accumulated Energy Rate \frac{d}{dt}\int_\Omega \rho E_0 \ dV \mathrm{ht.dEi0Int}
Total Net Energy Rate \int_{\partial\Omega} (\rho{\bf{u}}E_0-k{\bf{\nabla}} T+{\bf{q}}_\mathrm{rad}-\sigma{\bf{u}}) \cdot {\bf{n}} \ dS \mathrm{ht.ntefluxInt}
Total Heat Source \int_\Omega Q \ dV \mathrm{ht.QInt}

Therefore, the energy balance equation, written using the predefined variables in COMSOL Multiphysics, becomes

\mathrm{ht.dEi0Int} + \mathrm{ht.ntefluxInt} = \mathrm{ht.QInt}

Image showing the global evaluation settings.
Use of derived global evaluation values to evaluate the energy rate of predefined variables.

At steady state, the total accumulated energy rate vanishes. The total net energy rate and the total heat source must balance. The results for the stationary study are shown below. The relative error is again much less than the relative solver tolerance.

Screen capture displaying the energy balance.
Energy balance of a stationary analysis. The total net energy rate and total heat source must balance.

The energy rates are plotted below for the transient analysis. The total net energy rate increases progressively to finally reach its steady-state value, which balances the applied flux, 1 W, on the heat sink base. On the other hand, the total accumulated energy rate initially balances the total heat source and vanishes once the steady state is reached. Furthermore, the pink line denotes the energy balance absolute error, i.e., \mathrm{ht.dEi0Int} + \mathrm{ht.ntefluxInt}-\mathrm{ht.QInt}, which should be, in the best case, zero. The results show good agreement.

Plot comparing energy rate and time.
Energy rate versus time.

Concluding Remarks

We have discussed the theory of mass and energy conservation for static as well as transient conjugate heat transfer problems. We have also taken a look at how to calculate energy and mass balances with COMSOL Multiphysics in order to check the accuracy of simulation results. For this purpose, some useful derived values features have been introduced. Predefined energy rate variables are easy to use, and you can avoid handling the calculation of the energy rate expressions by yourself.

Although we have used a specific example to demonstrate the topics covered in this blog post, the demonstrated method can be extended to any conjugate heat transfer problem. For further reading about energy balances in COMSOL Multiphysics, please have a look at the Heat Transfer user guide.


  1. Vishwas Nesarikar October 14, 2015   11:16 pm

    Thank you for detailed explanation !
    Hope to see similar detailed blogs on other topics as well.

  2. Fabio Bocchi October 15, 2015   12:32 pm

    Hi Vishwas,

    Thanks for your feedback. We appreciate it.


  3. Nawaz Ahmad December 2, 2015   3:00 pm


    According to my observation based on working with COMSOL: COMSOL does not take into account the dispersive fluxes in mass balance, that is a major drawback. If you increase the dispersivity you have to face serious consequences in the form of mass balance error.

    I hope COMSOL experts will give attention to mass balance errors related to dispersion.


    Nawaz Ahmad

  4. Ali Mehrabifard February 3, 2017   2:45 am

    Thank you for your practical note. I am using surface integral technique to compute the mass flux in my doublet system in which there are two cylindrical shape wellbores, one of which is injection and the other one is production in a porous media. I sat a constant injection rate for the injection wellbore, and I interested in calculating the production rate on the production wellbore. However, in such a calculation, there is a point that when I am calculating the mass flux for injection well it is far different from the imposed value which I already sat. Accordingly, the other calculation including production mass flux (rate) is not reliable and indeed logically they are far from the expected values. I found some solutions, since it was the problem that many users are being faced. One solution was to remesh the model, I did that for a simple 2D it was working, but for the 3D model, the RAM memory is not supporting the calculations (I am using a 32 GB RAM memory). The other solution is to impose the injection flux (rate) by Weak Constraint, in this sense, based on my best knowledge, Lagrangian Multiplier method which is using here is a method by which maximum and minimum of an objective function is being calculated when there is a constraint function and boundary conditions. I got confused, how Weak Constraint can help me out?
    Is there any one who already modeled the wellbore in 3D with a line, if yes, how did u assign pressure value which can vary by depth?
    Thanks in advance,

  5. Katherine Sanchez May 15, 2017   3:57 pm

    Hi Fabio,

    I did not understand how they got the mass balance results for the transient analysis, You could share the expression of mass balance



  6. Caty Fairclough May 30, 2017   11:02 am

    Hi Katherine,

    Thanks for your comment.

    The mass balance is obtained through the “Derived values” in the image captioned “Step calculation for each mass flow rate term…”. Here, the expression for the normal mass flow (rho*u,n), which is the projection of the mass flow rate (rho*u) on the normal vector (n) to a boundary, is integrated both at the inlet and outlet boundaries. In addition, the accumulation of mass, d(rho)/dt, is integrated over the volume. The mass balance states that “in-out + accumulated = 0”.

    So the three derived values: “Inlet mass flow rate” (in, surface integral), “Outlet mass flow rate” (out, surface integral), and “Volume mass flow rate” (accumulation, volume integral) are the three derived values required to calculate the error in the total mass conservation (mass balance).

  7. Sunku Prasad J January 1, 2018   6:26 am

    thanks for the details of energy balance.
    I have one small doubt. How to perform energy balance COMSOL Multiphysics 4.3a ?.

  8. chenchen lee April 17, 2018   4:27 am

    Hi, Fabio Bocchi

    The value I calculated from inlet and outlet are available, but for the volume mass flow rate for every steps, it is “zero”, I do not know why. Thank you very much!

Loading Comments...