1
$\begingroup$

I am trying to model a 1-d advection-convection numerically, using an upwind scheme. I'm using the following equation to calculate the value of internal cells:

$$C_x^{t+1} = C_x^{t} + D\frac{\Delta t}{\Delta x^2} ( C_{x-1}^t-2 C_{x}^t+C_{x+1}^t) + u\frac{\Delta t}{\Delta x}(C_{x+1}^t-C_{x}^t)$$

and I'm using the Robin boundary condition (a.k.a. insulating boundary). In other words, I'm trying to guarantee no mass will leave the model environment. The boundary value is calculated by setting the leftmost and right most cells as:

$ C_1^{t+1} = D \frac{C_2^{t}}{u+D} $ and $ C_n^{t+1} = C_{n-1}^t\frac{(u\Delta x+D)}{D} $

Where $C$ is the concentration, $D$ the diffusion coefficient, $u$ is the velocity and $x$ is the current cell x-coordinate which goes from $1$ to $n$.

I understand this scheme should be numerically stable and accurate when both the following conditions are met:

$ u\frac{\Delta t}{\Delta x} + 2D\frac{\Delta t}{\Delta x^2} \leq1 $ and $\frac{u\Delta x}{D} < \frac{2}{1- u\frac{\Delta t}{\Delta x}} $

However, if I run this model with parameters $D=0.05$,$u=0.1$,$\Delta t=0.001$,$\Delta x=1$, there is some mass loss after 10,000 time steps, even tough the conditions above are met. It may seem a lot, but since my timesteps are very small it means less $t<10$ when that starts happening.

Any idea, of what could cause this loss of mass?

Below there's code to conduct this simulation in R.

transfun = function(time,state,parms) {
    D = parms$D
u = parms$u
    dt = parms$dt
dx = parms$dx
    n  = length(state)


    new = numeric(n)

    #Boundary cells
    state[1]= state[2]* D/(u+D)
    state[n]= state[n-1]* (u*dx+D)/D


    # Internal cells
    for(a  in 2:(n-1)) {
        state[a] = state[a] +D*dt/dx^2 * ( state[a-1] + -2*state[a] + state[a+1]) + 
                 u * dt/dx * (-state[a] + state[a+1])
    }

    #Boundary cells
    state[1]= state[2]* D/(u+D)
    state[n]= state[n-1]* (u*dx+D)/D


    return(list(state))
}

## Hand made loop 
# Ensure parameters follow the stability rule: 0<= C^2 <= 2s <= 1 
# where C = u*dt/dx and s = D*dt/dx^2
D=0.05
u=0.1
dt=0.001
dx=1
C = u*dt/dx
s = D*dt/(dx^2)
n = 15
T = 10000
print(C+2*s < 1)
print(u*dx/D<2/(1-C))



resul = matrix(NA,T,n)
resul[1,] = rep(0,n)
resul[1,ceiling(15/2)] = 1
for(a in seq(2,T,1)) {
    resul[a,] = transfun(1,state=resul[a-1,],parms=list(D=D,u=u,dt=dt,dx=dx))[[1]]
}

,

$\endgroup$
4
  • $\begingroup$ Finite difference methods are not conservative. If you want to conserve the physical properties of the model you should use an FVM or a DG method. Even then there will be some violation due to finite precision arithmetic. However, it is a scale issue. With FDM, the violation is dependent on the mesh size, while with FVM/DG it is usually at the order of $10^{-13}$ or less and it depends on the machine epsilon and the smallest normalized positive number. $\endgroup$ – Abdullah Ali Sivas Apr 21 at 17:58
  • $\begingroup$ This may be slight off-topic. If finite difference methods are not conservative, why do they exist? Do they have any practical use or they have been completed replaced by FVM/DG? Is it possible to say, for example, that altough mass is being loss the relative density in each cell remain the same? That's the only use I can imagine for a non-conservative method. $\endgroup$ – JMenezes Apr 22 at 10:28
  • 1
    $\begingroup$ They are mostly obsolete, but they still have some uses both pedagogical and otherwise. Also, most of the methods we have available now are not conservative. Even DG and FVM in space are not conservative if they are coupled with a non-conservative time stepping methods like some explicit RK methods. It is more of a matter of awareness and the level of "error"/"violation of conservation" that you are willing to accept. If it is above the tolerance, you can always refine the mesh (of course, this increases the computational cost) uniformly or adaptively. $\endgroup$ – Abdullah Ali Sivas Apr 22 at 14:19
  • $\begingroup$ @JMenezes Well, one advantage of finite difference method in comparison to FVM or generally FEM is that: it is much easier to understand the math behind the method and also implementation is much easier than FEM or FVM. Finite difference still is a good choice for solving PDEs in structured grids even beyond pedagogy and for serious research. $\endgroup$ – Alone Programmer Apr 22 at 15:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.