[Gesetzt]

Hello,

I am trying to simulate a reaction around a bubble. I have used Finite-Difference technique to solve the PDE. I am using second-order upwind scheme, because the advection gets very high in comparison to diffusion. Before that I used central difference schemes and there was an instability at the boundary as one would aspect for high Peclet number.

There are no instabilities at the boundary, but for small reaction constants, the equation only works well in the region of low convection. Somehow the concentration ist rising after a certain radius, which makes no sense at all.
The reaction term is of the form k_a*ca(i,j) where ca is the concentration. Is there another way to handle this source term?

Edit: It should be said, that I am simulating a steady-state case.


Best regards
Dan

[FMDenaro]

Originally Posted by Gesetzt

Hello,

I am trying to simulate a reaction around a bubble. I have used Finite-Difference technique to solve the PDE. I am using second-order upwind scheme, because the advection gets very high in comparison to diffusion. Before that I used central difference schemes and there was an instability at the boundary as one would aspect for high Peclet number.

There are no instabilities at the boundary, but for small reaction constants, the equation only works well in the region of low convection. Somehow the concentration ist rising after a certain radius, which makes no sense at all.
The reaction term is of the form k_a*ca(i,j) where ca is the concentration. Is there another way to handle this source term?

Edit: It should be said, that I am simulating a steady-state case.


Best regards
Dan

upwind produces numerical viscosity and the reaction velocity can be strongly affected.. have you tried to make a grid refinement study to check the improvement?
Central scheme can be used if you are able to work with Pe_h=O(1).

[Gesetzt]

Hello,

thanks for your help. I've tried refinement and it gets a little bit better, but I think it is not enough. I thought central scheme is bad for strong convection and I don't unterstand fully what you mean with central based if I use O(1) ( You mean trunctation error with this ? )

I think that the reaction term ist too small in comparison with reaction so I thought there might be some improvememt I can do with this term. I tried to use 1/2*(ca(i-1,j)+ca(i,j))*k_a but it there is no improve. ( i is in radial direction and j is angular direction ).

Best Regards,
Dan

[FMDenaro]

I mean the cell Peclet number O(1)

[Gesetzt]

Hello,

Unfortunatly I think I can't make my grid so fine to achieve this. So I would need another solution.

[FMDenaro]

to reduce the artificial viscosity you can try a third order upwind, for the steady state solution you can implement the QUICK scheme. Note that it is a Finite Volume method.

However, if the reaction is stiff you need further care.

[Gesetzt]

Hello,

Thanks for your help.
I've tried QUICK. But it isn't much better. I think the problem is about the skewness of my streamlines. I have read, that upwind, quick etc. only works well, when the stream goes right into the grid. Because I want to simulate the flow around the sphere this is not the case everywehre. Is there a simple second order scheme which is easy to implement with skewness correction ?

Best regards

Dan

[FMDenaro]

Originally Posted by Gesetzt

Hello,

Thanks for your help.
I've tried QUICK. But it isn't much better. I think the problem is about the skewness of my streamlines. I have read, that upwind, quick etc. only works well, when the stream goes right into the grid. Because I want to simulate the flow around the sphere this is not the case everywehre. Is there a simple second order scheme which is easy to implement with skewness correction ?

Best regards

Dan

I dont think that this is the problem...and QUICK is a multidimensional scheme that works well for steady solution of flows not aligned with the grid

[Gesetzt]

Hello,

I want to go more in detail to find the problem.

I have a problem with the following equation:




I'm trying to solve this with finite difference technique. My discretization :



A stands for c here.

I first tried it with this, but because I used central-difference in doesn't work for Pe(cell)>2 I had strange oscillations at the boundary. So I thought I use a upwind scheme now.

V_r and V_theta are given through polynoms as a function (r,theta).

So for theta<90 ° I Use second order upwind for the convection term V_r > 0.

And for theta >90 ° in the other direction.
V_theta is always <0, so I only use this form for this.



This is my coordinate system :



I also have seen that the problem is strongly dependent on V_theta. If I multiply V_theta with 0.9 for example there is no problem. The concentration is not rising. I've looked at the speed vectors in more detail. The problem is vanishing when the speed vectors near the zero degree line are horizontal or show away from it. I have a Neumann boundary condition there ( also at 180°). Now I also have seen that for the regions where this "pseudo-source" exist the boundary condition isn't met perfectly. There is a gradient, but it should be zero. So I am thinking maybe this acts like a "pseudo-source ". I also have drichilet boundary condition for r=1 and r-->infinite. ( My calculation begins at r=1).

My grid numbering is the following ( for example I have three radial and five angular steps ). My notation is c(i,j), i is the radial direction and j the angular direction. So say for i=1, it is c(1,j) with j=0...4 from phi=0...180°. For zero degress with constants l_i :
c10*l1+...c12*l2+c14*l4=0.

I need no ghoistpoints here. I just sad c12=c14 and rearranged the equation to c10*l1...+c12(l2+l4).
For 180 ° I need ghoistpoints because of second order upwind scheme I have c(i,j-2) which is outside my domain, but at 180 ° there is no problem.

Should I implement Neumann on another way ? Is this a problem because of variable coefficients ?

Best regards

Dan

[FMDenaro]

What about a plot of the solution to see the problem with second order upwind??

First of all, I strongly suggest to check the code with the first order upwind to see if it works. The strong artificial viscosity of the first order upwind would lead to a solution not enough accurate but this way you are sure that any problem you see is due to the second order upwind.

The range of your domain goes from r=1 to r=?

[Gesetzt]

Hello,

I'm now following your suggestion.
I am looking in the code for first-order upwind to find a fault, so far I haven't found one. The outer boundary of r is normally r=2.5, but I also tried it with other values.

Here are the plots :








The second one is a zoom, the concentration there shouldn't rise. The third plot is about the angular direction. As you can see the zero slope criterium isn't fullfilled. I don't understand why because I just rearranged the equation as described above ( here I have used only first order for the Neumann boundary condition ).

I used 1000 steps in radial and 1000 steps in angular direction. Making the grid finer does not change much. I'm using the direct solver in Matlab, because it is fast enough.

[Gesetzt]

Hello FMDenaro,

I have one (maybe very stupid ) question: As for example I want to write my equation for the point c(1,4). I would take all the coefficients only from this point so I would write V_r(1,4)*(c(2,4)-c(0,4)/2*delta_r)+V_phi(1,4)/r....
Or would I have to average the coefficients in some way ?

Maybe I have misunderstand something about the variable coefficients in my equation... Now I'm thinking this would be the right way ?

Best regards

Dan

[FMDenaro]

Originally Posted by Gesetzt

Hello FMDenaro,

I have one (maybe very stupid ) question: As for example I want to write my equation for the point c(1,4). I would take all the coefficients only from this point so I would write V_r(1,4)*(c(2,4)-c(0,4)/2*delta_r)+V_phi(1,4)/r....
Or would I have to average the coefficients in some way ?

Maybe I have misunderstand something about the variable coefficients in my equation... Now I'm thinking this would be the right way ?

Best regards

Dan

sorry, but in your example (1,4) isn't a point on the boundary????
I need a remind and a sketch of your example ...

[Gesetzt]

Hi,

My domain is going from i=0 to i=4 (radial) and j=0 to j=4 (angular) including the boundary. But I only have Dirichlet-condition for i=0 and i=4, so I would have to set up the equation also for i=1 and j=4 as far as I understand, I just have to be careful because I have a Neumann-condition here, so I would have to change the stencil.

I would take c(1,3) for example then. If I now set up the equation for this which value of V_R , V_Phi would I have to take ? They are all function of the radius and the angular direction. Can I take just ( for this example ) V_R(1,3) and V_Phi (1,3) ?

I am sorry, I hope I could explain my problem well enough.

Best regards

Dan

[FMDenaro]

but at any boundary point you just set the Dirichlet value, you do not write any equation...

[Gesetzt]

No not for the points with Drichilet, but for the points with Neumann I have to set up the equation ?

My Matrix looks like this for the example :



So in the first four rows/columns and the last four I have included the dirichilet condition, because of that there are only ones. Let's say I have A*x=b , in my b-vector I have included the dirichelet values for the first four and the last four points.

In the fifth row I have included Neumann Boundary conditon ( also in row 12 ).

[FMDenaro]

Immagine a simple 1D grid of points going from 1 to N. points 1 is Dirichlet and points N is Neumann. Now the equations are written for points going from 2 to N-1 and the matrix system is of order (N-2)x(N-2). The first row is related to equation written in 2 and the Dirichlet boundary condition makes the entry at node 1 a known value. The last row is related to the Neuman BC that alters the entries according to the discretization of the derivative.

In any case, the equation is never discretized at the boundaries.

[Gesetzt]

Hello,

Thanks for your help and patience.
I don't understand fully, I have sticked to this script ( page 116 ):

http://www.scientificpython.net/uplo...k2_lecture.pdf

I want to use centered difference at my boundary, so I would have to introduce the ghost points c15 for the boundary on 180 ° and c1-1 for the boundary on 0°.

My understanding of this was, that I have to write the differential equation now on the boundary ( for example c14) because I have the additional unknown c15 and combine this two equations. If I would use a one-sided difference instead I don't discretizise on the boundary ?

Best regards

Dan

[FMDenaro]

consider the node N to be a boundary node where you have the Neumann BC df/dx=q. As example, consider the diffusion operator written in the node N-1:

[f(N-2)-2*f(N-1)+f(N)]/dx^2

Now, the value f(N) is not known and you do not write an equation for it. You discretize df/dx=q and express f(N) in terms of surrounding nodes. I suggest the second order backward discretization written for the derivative in N and express f(N).


Note that, a different grid colocation on staggered node would produce a different implementation of the second order discretization of the diffusion operator. In this case you can simply substitute the Neumann condition in

d/dx(df/dx)|N-1/2=[df/dx(N) - df/dx(N-1)]/dx = [q - df/dx(N-1)]/dx

[FMDenaro]

Originally Posted by Gesetzt

Hello,

Thanks for your help and patience.
I don't understand fully, I have sticked to this script ( page 116 ):

http://www.scientificpython.net/uplo...k2_lecture.pdf

I want to use centered difference at my boundary, so I would have to introduce the ghost points c15 for the boundary on 180 ° and c1-1 for the boundary on 0°.

My understanding of this was, that I have to write the differential equation now on the boundary ( for example c14) because I have the additional unknown c15 and combine this two equations. If I would use a one-sided difference instead I don't discretizise on the boundary ?

Best regards

Dan

consider the matrix in form 8.19 in your case