8

About the code: I have a code which simulates concentration from advection-diffusion-reaction PDE in 2D space (X,Y) with time. The solution is obtained using fully implicit finite-difference method and includes the capability to simulate a media with spatially varying permeability and reaction constant (through upwinding by harmonic mean). I have been able to test the code for homogeneous media and it works fine. It's a solution for the following equation: $$\begin{align} \frac{\partial C}{\partial t} + \nabla.\left(v C - D\nabla{C} \right)= \alpha C \end{align}$$

Issue: The code gives realistic solution for a medium with uniform layers of permeability as shown below. However, once the medium starts to get more heterogeneous it starts to throw unrealistic results (i.e. negative concentration and large fluctuations with time). Any guesses why I am getting unrealistic results for heterogeneous media (with values in the similar range as uniform layered media)?

UPDATE: I think I have figured the issue for getting negative values and sharp fluctuations in concentration. I think the problem is not with the numerical model, but with the physical values. Even though I was already taking realistic parameter values for which I was getting unrealistic concentration, I had to further adjust them (decrease the convection velocity) to get realistic results

haresfur
  • 4,419
  • 16
  • 32
  • could be related to which finite difference scheme is being used(forward,backward,difference).could also be due to not using weighted values – shrey Nov 22 '15 at 05:31
  • Forward difference in time and space. May I know which weights do you mean? –  Nov 22 '15 at 06:32
  • Not an expert, but I wonder if the problem is in your solver for the finite difference. IIRC, the transport equations can be subject to a large amount of numerical dispersion. What happens if you use a finer grid? – haresfur Nov 22 '15 at 21:57
  • haresfur: I have tried finer grids but that didn't help either, though it changed the results to some extent but the problem of sharp oscillations still existed even with finer grids. –  Nov 22 '15 at 22:04

1 Answers1

2

It's not clear exactly what is being modelled here, but it seems to me that there are two ways in which the concentration can 'go negative'. Firstly, the rate of change of concentration can be massive, in which case see what happens when modelling with much smaller time steps. Or, the diffusion term substantially exceeds the advection term, which is physically unlikely. However, you are modelling a large (more than 10-fold) anisotropy, in which almost all the flow will be in the x, horizontal direction - it is almost a one-dimensional problem in fact. So how can the advection be less than the diffusion? In the case of a constant concentration source, such as a saline water body, the negative results don't make sense, in which case the calculations have gone unstable, probably from using too large a time step. In the case of a pulsed source, as the concentration peak passes (along the x direction), the advected concentration rises and falls, such that the lateral (y axis) diffusion increases, then decreases, i.e. the concentration gradient changes direction, becoming negative.

As if this wasn't hard enough, in real life applications the 'sink' terms probably far exceed any back-diffusive effects.

Gordon Stanger
  • 14,238
  • 23
  • 44
  • Gordon: I have implemented fully implicit scheme in which the accuracy of results are independent from the effect of time step size. Additionally, even if diffusion is zero (in most cases it's almost negligible compared to advective flux), I still have this problem. –  Nov 22 '15 at 18:19
  • Then I'm a bit puzzled also. Do you have the source code, or is this an off the shelf package? – Gordon Stanger Nov 22 '15 at 19:24
  • I wrote it myself...so I have the source code, Gordon. –  Nov 22 '15 at 20:55
  • In what language? Fortran I can read. C, I can't – Gordon Stanger Nov 23 '15 at 07:04
  • In python, Gordon –  Nov 23 '15 at 13:12
  • Python? Never heard of it, but then I am an old fashioned fossil myself. Sorry I can't help you with the code, but I will ask some of my former colleagues about this problem when I get back to Oz - I'm in Georgia at present. – Gordon Stanger Nov 24 '15 at 15:14
  • @Pupil Please post your code. – Isopycnal Oscillation Nov 24 '15 at 21:14
  • @IsopycnalOscillation: Can I share with you separately? –  Nov 24 '15 at 21:46
  • @Pupil Sure, i made a temp email at alt.b1-cm892rk@yopmail.com you can send it there if you'd like. I cant promise anything but I'll take a look. – Isopycnal Oscillation Nov 25 '15 at 05:21
  • @Pupil Ignore previous messages. I was trying to move the conversation to chat so we dont spam the comment box but it does not work... I created a chat room http://chat.stackexchange.com/rooms/32129/solution-for-advection-diffusion-reaction-pde-with-heterogeneous-porous-media join me there. – Isopycnal Oscillation Nov 25 '15 at 21:08
  • @GordonStanger: I think I have figured the issue...I think problem is not with the numerical model, but with the physical values. Even though I was already taking realistic parameter values for which I was getting unrealistic concentration, I had to further adjust them (decrease the convection velocity) to get realistic results. –  Nov 26 '15 at 18:59
  • @pupil That is what i meant by scales. I am glad you got it to work. – Isopycnal Oscillation Nov 27 '15 at 07:37
  • Well, that's a relief, I was scratching my head over what could be driving negative values. Happy to hear it is sorted. – Gordon Stanger Nov 27 '15 at 19:29