Finite difference method help (bvp)

[beedle]

Hi, I'm currently writing a code to solve a steady-state boundary problem across multiple layers of a system. The system involves diffusion-reaction of various species in a porous medium. I am simply using central finite differences to model this setup, which says that essentially -div*J+Q = 0. (I take J to be mole flux of a species and Q to be the rate of creation of that species per volume). It is one dimensional, so basically, my equations typically look something like this:

$$J = vC - D_i\frac{dc}{dx}$$

$$-\frac{dJ}{dx} + Q_i = 0$$

$$-v\frac{dC}{dx}+D_i\frac{d^2C}{dx^2}+Q_i=0$$

or, in terms of central differences

$$-v\frac{C_{j+1}-C_{j-1}}{h_i}+D_i\frac{C_{j+1}-2C_{j}+C_{j-1}}{h_i ^2} +Q_i = 0$$

where the i's designate which phase/region I am in. So I have written this set of equations in various regions of my problem. At the interfaces of these regions, I know I am supposed to say that the fluxes and concentrations of the species are continuous. Take the interface in the following example to be located at some point j, and let phase 1 lie to the left and phase 2 lie to the right. I want to write something like this at point j:

$$A(J_{j-1/2} - J_{j+1/2}) + Q_{overall}dV = 0$$

So here is what I am thinking for the centered finite differences:

$$(v\frac{C_j+C_{j-1}}{2}-D_{i=1}\frac{C_j-C_{j-1}}{h_{i=1}}) -(v\frac{C_{j+1}+C_{j}}{2}-D_{i=2}\frac{C_{j+1}-C_{j}}{h_{i=2}})+Q_{i=1}\frac{h_{i=1}}{2}+Q_{i=2}\frac{h_{i=2}}{2}=0$$

Again, Q is rate of rxn per unit volume, so I am weighing according to the step size h between phases. I'd just like to check if this is the correct way of imposing the continuity of fluxes between phases. I bring this up because a code I am looking at from a former grad student in my research group has somewhat different looking fd's. He neglected the terms involving the v's, but his other terms still look different:

$$(D_{i=1}\frac{C_j-C_{j-1}}{h_{i=1} ^2}) -(D_{i=2}\frac{C_{j+1}-C_{j}}{h_{i=2} ^2})+\frac{Q_{i=1}+Q_{i=2}}{2}=0$$

the two h's are not necessarily equal, so I am not clear about this way of doing it and how it makes sense. If the two h's -are- equal, then my way reduces to his way. If my way is wrong, could someone please explain why? Thanks. I will be on for a while so if anything in my post needs clarification just ask.

 

[jvicens]

Hello,

Thank you for sharing your code and thoughts on how to impose the continuity of fluxes between phases. You are correct in your approach to include the v and Q terms in your central finite differences equations. These terms are important as they represent the transport and reaction of the species in the system.

In the code you are looking at from the former grad student, it seems that they have neglected the transport terms (v) and only included the diffusion and reaction terms (D and Q). This approach may work in certain cases where the transport is negligible, but it is not a general solution and may lead to incorrect results.

In order to understand the difference between your approach and the former grad student's, it is important to look at the physical meaning of the terms in the equations. The v term represents the transport of the species due to advection (flow) in the system. This is an important term to include, especially in porous media systems where advection can have a significant impact on the transport of species. The D term represents the diffusion of the species in the system, and the Q term represents the reaction rate of the species per unit volume.

In the former grad student's approach, the transport term (v) is neglected and the diffusion and reaction terms are averaged at the interface. This may work in certain cases, but it is not a general solution. In your approach, the transport term is included and the diffusion and reaction terms are weighted according to the step size between phases. This is a more accurate representation of the physical processes happening in the system.

Therefore, I would recommend using your approach in your code as it includes all the important terms and is a more accurate representation of the system. I hope this helps clarify the difference between the two approaches. If you have any further questions, please let me know. Good luck with your code!