I trying to do a simulation where there are two quantities, $\Delta\theta(x,t) = \theta(x,t) - \theta_{o}$ and $\Delta\nu(x,t) = \nu(x,t) - \nu_{o}$. These quantities are chemical concentrations resultant from reactions. There is no significant convective flow only diffusive flow.

The quantities $\Delta\theta $ $\Delta\nu$ are related to each other by a reaction reaction, this reaction is a source term for both dependent on the other.

Furthermore, both variables have source terms that depend on themselves, i.e. if $\Delta\nu$ or $\Delta\theta\neq0$ there are reactions that try to drive them to equilibrium.

1) Because of this dependence, when I'm setting up my equations, should I put my source terms within the matrix (and not on the b vector)?

2) Is FVM still conservative if I have variable dependent sources?


Functional form of the equation:

So we don't get confused, let: $\Delta \theta = \theta - \theta_{o} \rightarrow \theta$, because $\theta_o$ is the equilibrium value and we're only interested in quantities away from equilibrium. Then:

$\frac{\delta \theta(t,x)}{\delta t} = \Delta D \Delta \theta(t,x) + Source(\theta,\nu)$

Notice, that I explicitly show the Source depenence on $\theta$ and $\nu$. So when I write the $a_{p}$ term it contains, not just the diffusive terms from $a_E$ and $a_W$, but also from the source.

Also I solve it in phase space by replacing the time-domain $\theta$ with a frequency domain one.

  • 1
    $\begingroup$ I think it would be useful to state the equations you are considering. $\endgroup$ Commented Feb 10, 2016 at 19:33
  • $\begingroup$ @WolfgangBangerth Edited to show the form of the equations used. The form is the same for both theta and nu. Fundamentally, however, the question is, is FVM conservative when S depends on the variable, and if not what are the implications? Every book I look assumed S depends on the position or something else except the variable we're trying to solve for. $\endgroup$
    – user19247
    Commented Feb 10, 2016 at 23:57
  • $\begingroup$ Your first question is unclear because you introduce a variable $b$ that is never mentioned (much less defined) elsewhere. Your second question is unclear since we don't even know enough about your equation to determine whether the exact solution is conservative. In fact, you haven't even stated what quantity it is that you would like to conserve. $\endgroup$ Commented Feb 17, 2016 at 17:50

1 Answer 1


The finite volume method is said to be conservative because you use a relation between surface integrals on the faces of the elements and volume integrals on the elements to construct the matrix to solve. From wikipedia:

These terms are then evaluated as fluxes at the surfaces of each finite volume. Because the flux entering a given volume is identical to that leaving the adjacent volume, these methods are conservative

Did you write your own fvm solver? or are you using a package? If you wrote your own you should be able to check that indeed "the flux entering a given volume is identical to that leaving the adjacent volume" because you are familiar with how you build your matrix.

I did partially write a fvm solver for turbulent flow like 10 years ago. From what I remember, if everything is done properly, I don't see why in this case your problem should not be conservative.

The FVM method is said to be conservative because the way the system is build. It involves using the divergence theorem on a small volume. And when all the math is laid out, you realize that indeed the system is conservative. I don't see why whether the source term depends on the variable or not is important, because you'll integrate as before. But again, to be sure, the best is write down the development of fvm on your equation.

So my bet is that yes, it is still convervative. But to be sure we would have to see how you discretize the integrals and interpolate the average values on the faces of the volumes. If done properly, it should be conservative.

EDIT: Let's expand this with some math. I'm going to follow (even with some copy/paste) the wikipedia article on the FVM because is available to all. Specifically, the general conservation law using this PDE:

$ \dfrac{\partial \mathbf{u}}{\partial t} + \nabla\cdot \mathbf{f}(\mathbf{u}) = 0$

Here, $\bf u$ represents a vector of states and $\bf f $ represents the corresponding flux tensor. Again we can sub-divide the spatial domain into finite volumes or cells. For a particular cell, $i$ , we take the volume integral over the total volume of the cell, $v _{i} $ , which gives,

$ \int_{v_i} \dfrac{\partial \mathbf{u}}{\partial t} dv + \int_{v_i} \nabla\cdot \mathbf{f}(\mathbf{u}) dv = 0$

now let's write the first integral in terms of the average value $\bar{u}$ of $u$ in that cell. That is, $ \int_{v_i} \dfrac{\partial \mathbf{u}}{\partial t} = v_i \dfrac{\partial \mathbf{\bar{u}}}{\partial t}$

to solve the second integral, we use the divergence theorem as follows:

$ \int_{v_i} \nabla\cdot \mathbf{f}(\mathbf{u}) dv = \oint_{s_i} \mathbf{f}(\mathbf{u}) \cdot \mathbf{\bar{n}} ds $

where $s_i$ is the surface that covers the element $i$. That way we get to:

$ v_i \dfrac{\partial \mathbf{\bar{u}}}{\partial t} + \oint_{s_i} \mathbf{f}(\mathbf{u}) \cdot \mathbf{\bar{n}} ds = 0$

Here is when it gets tricky, and how you do it in your equation it may affect your result.

now, imagine that you discretize a volume in triangles. Assign to the center the average value. You can use finite differences to get rid of the time derivative, and get the next value of $u$ using an explicit scheme. To calculate the surface integral, you can interpolate the value between the different points (where you have the average of $u$) at the surface. Wikipedia puts it like this:

Again, values for the edge fluxes can be reconstructed by interpolation or extrapolation of the cell averages. The actual numerical scheme will depend upon problem geometry and mesh construction. ... Finite volume schemes are conservative as cell averages change through the edge fluxes. In other words, one cell's loss is another cell's gain!

from your equation, you have to figure out how to integrate the source term. From your comment, note that the second integral depends too on $u$. The issue is that you have a time derivative, meaning that in your linear system Ax=b, A will have the coefficients of the next value of $u$, while b will be computed using the current values of $u$. This is of course implementation dependent but I don't see other way to do it. Going back to your equation, (from this point on all is speculation because I don't know what your function is), let's assume that

$\int_{v_i} \mathbf{S}(\mathbf{u},v) = v_i \mathbf{S}(\mathbf{\bar{u}},v)$

if this is correct, then you have

$ \dfrac{\partial \mathbf{\bar{u}}}{\partial t} = \dfrac{1}{v_i}\oint_{s_i} \mathbf{f}(\mathbf{u}) \cdot \mathbf{\bar{n}} ds + \mathbf{S}(\mathbf{\bar{u}},v)$

which can be solve numerically. How? it depends on how you solve the time derivative, how you make the interpolation between the nodes to calculate the surface integrals. Note that you practically have the same problem as before.

Is this conservative? well, how is your function source? is the above assumption valid? how do you construct your system? are you writing your code or using a library? But like I said, if things are done properly, and source is "well behaved", my bet is yes.

Hope this helps. Like I said, you have to go through your code. Take a good FVM book and go through the method. I hope this info has helped and pointed you in the right direction.

  • 1
    $\begingroup$ Div = 0, and therefore Gauss' theorem is applicable, when the there are no sources. From what I understand, when there are sources FVM puts them in the b vector, to recover the relation ap = aw + ae + an+ as. But I can only put the source in the b vector if it is a constant wrt to the variable I'm trying to solve. $\endgroup$
    – user19247
    Commented Feb 12, 2016 at 14:12
  • $\begingroup$ I'll try to elaborate and put some equations to explain myself better later on the afternoon. Maybe that will help you. At the end you have to go through your code. $\endgroup$
    – cauchi
    Commented Feb 12, 2016 at 14:18
  • $\begingroup$ Done. Hope it helps $\endgroup$
    – cauchi
    Commented Feb 12, 2016 at 17:32
  • $\begingroup$ Thank you, I have to read it in detail though. Since I'm not concerned with transients and I have a sinusoidal excitation, I transform the system into phase space and write $\theta \rightarrow i \omega \Theta$. That is I use phasors to solve it. $\endgroup$
    – user19247
    Commented Feb 14, 2016 at 13:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.