# How to write this non-linear PDE with the finite volume method?

I wish to solve a coupled system of non-linear equation of this form,

$$u_t = -(\mathcal{F})_x + f(x,u,w) \\ w_t = \mathcal{F} + g(x,u,w)$$

by stepping the equations forward in time. The first equation is the advection-diffusion equation with non-linear reaction term, and where the flux has the definition,

$$\mathcal{F}=au-du_x$$

First equation

$$\int_{\Omega}u_t~dx = -\int_{\Omega}(\mathcal{F})_x~dx + \int_{\Omega}f(x,u,w)~dx \\ \tilde{u}_t = \frac{1}{h_j}\left( -\mathcal{F}_{j+\frac{1}{2}} + \mathcal{F}_{j-\frac{1}{2}} \right) + \tilde{f}(x_j,u_j,w_j)$$

where $h_j$ is the width (not volume because I am considering 1D only) of the $j$-th cell.

Second equation

However, how should the second equation be written, can I do the following?

$$\int_{\Omega}w_t~dx = \int_{\Omega}\mathcal{F}~dx + \int_{\Omega}g(x,u,w)~dx \\ \tilde{w}_t = \frac{1}{h_j}\mathcal{F}_jh_j + \tilde{g}(x_j,u_j,w_j)$$

Moreover, if we assume that the flux in the $j$-th cell is constant across the cell then the integral becomes,

$$\int_{\Omega}\mathcal{F}~dx = \mathcal{F}_jh_j$$

This is simply base $\times$ height of the cell, but because I want to know the cell average values I divide again by the cell width. Finally I get the following for the second equation,

$$\tilde{w}_t = \mathcal{F}_j + \tilde{g}(x_j,u_j,w_j)$$

Q1. Is this correct?

In terms of implementation I know 3 of the terms in $\mathcal{F}_j$, they are:

• $a_j$, the cell advection velocity,
• $d_j$, the cell diffusion coefficient,
• $u_j$, the cell mass concentration,

However, $u_x$ is not known.

Q2. Can I simply write $u_x$ using the values of concentration interpolated to the cell faces? For example,

$$u_x = \frac{1}{h_j}\left( u_{j+\frac{1}{2}} - u_{j-\frac{1}{2}} \right)$$

This doesn't seem quite right because I was taught that in the finite-volume method, concentrations are defined in the cells and fluxes are defined on the faces.

Update

Regarding definition of fluxes.

For the first equation I define the fluxes as,

$$\mathcal{F}_{j+\frac{1}{2}} = a_{j+\frac{1}{2}}\left( \frac{h_{j+1}}{2h_{+}} u_j + \frac{h_j}{2h_{+}} u_{j+1} \right) - d_{j+\frac{1}{2}} \frac{u_{j+1}-u_j}{h_{+}}$$

$$\mathcal{F}_{j-\frac{1}{2}} = a_{j-\frac{1}{2}}\left( \frac{h_{j}}{2h_{-}} u_{j-1} + \frac{h_{j-1}}{2h_{-}} u_{j} \right) - d_{j-\frac{1}{2}} \frac{u_{j}-u_{j-1}}{h_{-}}$$

where $h_j$ is the cell width, and $h_{\pm}$ are the centroid distances in the forward and backward direction.

Here linear interpolation is used to define the variable at the cell faces $u(x_{j-1/2})$ and $u(x_{j+1/2})$, and forward difference and backward difference is used to define $u_x(x_{j+1/2})$ and $u_x(x_{j-1/2})$. This gives a central discretisation for both the advection and diffusion terms (however, upwinding can still be introduced easily).

However for the second equation the flux is not defined on the cell's left and right face, it is defined as the cell average, so I need to use a different definition for,

$$\mathcal{F}_j = a(x_j)u(x_j) - d(x_j)u_x(x_j)$$

Clearly, $a(x_j)$, $u(x_j)$, and $d(x_j)$ are already known. The only term than needs defining for the second equation is $u_x(x_j)$. So maybe I could use,

$$u_x(x_j) = \frac{1}{h_j}\left( u_{j+1/2} - u_{j-1/2} \right)$$

with linear interpolation, which as you pointed out would give a second order accurate central difference.

Could I also use?

$$u_x(x_j) = \frac{1}{h_{+}}\left( u_{j+1} - u_{j} \right)$$

In my comment below, this is what I was getting at. How do define the $u_x(x_j)$ for the second equation, what are the pros. and cons. of using a central or upwinding scheme for this term.

Second equation

However, how should the second equation be written, can I do the following?

$$\int_{\Omega}w_t~dx = \int_{\Omega}\mathcal{F}~dx + \int_{\Omega}g(x,u,w)~dx \\ \tilde{w}_t = \frac{1}{h_j}\mathcal{F}_jh_j + \tilde{g}(x_j,u_j,w_j)$$

Moreover, if we assume that the flux in the $j$-th cell is constant across the cell then the integral becomes,

$$\int_{\Omega}\mathcal{F}~dx = \mathcal{F}_jh_j$$

This is also true if the flux is linear in the cell, not just constant.

Q1. Is this correct?

Yes.

However, $u_x$ is not known.

Q2. Can I simply write $u_x$ using the values of concentration interpolated to the cell faces? For example,

$$u_x = \frac{1}{h_j}\left( u_{j+\frac{1}{2}} - u_{j-\frac{1}{2}} \right)$$

Yes, in finite difference/volume you generally have to compute the derivatives from the conserved variables. Note that to find the derivative $u_x$ at $j$ you can use a central difference (which may be unstable depending on the form of the diffusion and you might need to apply upwinding), which would be

$$u_{x,j} = \frac{1}{2h} (u_{j+1} - u_{j-1}) + O(h^2)$$

assuming constant $h$ cell size. If you were using a linear interpolation to find $u_{j-1/2} = 1/2(u_{j-1}+u_{j})$ you can see that the above expression for central differencing using $j+1$ and $j-1$ is exactly the same as your $u_x = \frac{1}{h_j}\left( u_{j+\frac{1}{2}} - u_{j-\frac{1}{2}} \right)$

• Thank you! The diffusion component to the flux is $du_x$ where $d>0$ so is it correct that upwinding can be introduced by choosing, $$du_x = \frac{d}{h}(u_{j+1} - u_{j})$$ (assuming constant step size)? The advantages of upwind would be that is prevents oscillation when $u_x$ is very large (sharp gradients). What would be the other benefits and how do I know which scheme (central or upwind) is more stable? – boyfarrell Jul 10 '13 at 23:32
• Not quite. You don't upwind the diffusive operator, but rather the advective one. Physically you can understand this as diffusion happening in all directions (elliptical) and advection having a direction (hyperbolic) - you can't smell things downwind from you. You'll find this exact question in any intro CFD text - look for von Neumann stability analysis. And first look at strictly the advection equation and be able to solve that before adding diffusion. Character-limited answer: Central is more accurate (2nd order), but unstable. Simple upwinding is stable, but only first order accurate. – Aurelius Jul 11 '13 at 2:42
• You picked up all the important points in less that 500 characters, impressive! What I was asking (it's slightly unclear) is how to define the flux term $F_j$. I have made an update which hopefully clarifies what I mean. – boyfarrell Jul 11 '13 at 8:28
• There's a few methods, but one of the huge achievements in cfd history is this en.m.wikipedia.org/wiki/MUSCL_scheme follow that scheme exactly for determining face fluxes and you'll be good to go. – Aurelius Jul 11 '13 at 10:32
• The MUSCL seems to be a good solution, in that quantities will all remain positive in the equation (I need the solution to remain positive for a physical valid simulation). However, I would first like to implement the simplest approach. Am I correct in saying that if I use an upwind discretisation for $u_x$ (see the updated above), then this second equation will remain positive? A central discretisation would cause oscillations. – boyfarrell Jul 11 '13 at 15:18