CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   Finite Volume Method For Cylindrical Coordinates (http://www.cfd-online.com/Forums/main/77713-finite-volume-method-cylindrical-coordinates.html)

falopsy July 1, 2010 09:32

Finite Volume Method For Cylindrical Coordinates
 
I am developing a CFD code for a stirred tank but have quite a few issues I will appreciate if anyone can help me with:

1. There are very few books on the discretisation of the Navier-Stokes equation in cylindrical coordinates. I am having problem knowing which terms to put at the surface and which to put at the centre due to the appearance of terms like 1/r^2.

2. Also, I have issues with staggered grid and I dont really understand the collocated grid methods either. My problem with staggered grid has to do with how to store the velocity vectors as matrices as most program dont accept frational index i.e. vij+1/2 used in staggered grid.

3. Most books have the discretised equation for rectangular coordinate but none for cylindrical coordinates. Is there any good book that explained the discetisation of the Navier-Stokes equation in cylindrical coordinates?

4. I am also thinking maybe I should work in rectangular coordinates but with a cylindrical coordinate fitted grids and then approximate the length of the mesh as rectangular i.e dx=dr, dy = rdTheta, dz= dz. My only concern here is that I need to now convert the source term vectors from rectangular to cylindrical and also my boundary conditions. Does anyone understand what I'm trying to say or does it not make any sense.

Any help will be greatly appreciated. Thanks

arjun July 1, 2010 17:11

Quote:

Originally Posted by falopsy (Post 265285)
I am developing a CFD code for a stirred tank but have quite a few issues I will appreciate if anyone can help me with:


Is mesh motion involved??
If it is there then how will you handle that in cylindrical coordinates?


PS: I have written a code for cylindrical coordinate but can not comment much on it as lots of patents etc are still pending. At the year end, there will be some publications that time probably you will find out.

falopsy July 1, 2010 18:57

Thanks arjun,for your response. Mesh motion is not involved. I will just be grateful if you can give me the equations for the coefficients (aP, aE, aW, etc) in the discretised equation so I can compare them to the one I obtained to verify the correctness before I start coding. Thanks

arjun July 1, 2010 20:52

Quote:

Originally Posted by falopsy (Post 265347)
Thanks arjun,for your response. Mesh motion is not involved. I will just be grateful if you can give me the equations for the coefficients (aP, aE, aW, etc) in the discretised equation so I can compare them to the one I obtained to verify the correctness before I start coding. Thanks

Nothing great about them. I use collocated method and SIMPLE algorithm as peric describe. And i treat it as a unstructured mesh.
Download the peric's code and work along that line, you will be fine.

My patents are about mesh motion in large meshes and development of solver for pressure equation. (Solves to machine precision in less than 5 iteration for 1-2 billion cells mesh).

bjohn July 2, 2010 00:14

If it is finite-volume, I don't see the point of solving the equations written in the cylindrical coordinate system (maybe I'm missing something). If finite-volume, you have a control volume and You integrate the equations over the control-volume. What you need is then a numerical flux at the face in the face normal direction.

Maybe I'm missing simething, but I don't know why the equations are discretized in the cylindrical coordinate system. Ignore me if it is not relevant.

arjun July 2, 2010 00:25

Quote:

Originally Posted by bjohn (Post 265356)
If it is finite-volume, I don't see the point of solving the equations written in the cylindrical coordinate system (maybe I'm missing something). If finite-volume, you have a control volume and You integrate the equations over the control-volume. What you need is then a numerical flux at the face in the face normal direction.

Maybe I'm missing simething, but I don't know why the equations are discretized in the cylindrical coordinate system. Ignore me if it is not relevant.

you are not missing very much.

In fact i think in his case it would be better for him to just write a mesh from cylindrical format to any format that openFOAM type solver can understand and use them. It would be better than reinventing the wheel.

falopsy July 2, 2010 05:56

Quote:

Originally Posted by bjohn (Post 265356)
If it is finite-volume, I don't see the point of solving the equations written in the cylindrical coordinate system (maybe I'm missing something). If finite-volume, you have a control volume and You integrate the equations over the control-volume. What you need is then a numerical flux at the face in the face normal direction.

Maybe I'm missing simething, but I don't know why the equations are discretized in the cylindrical coordinate system. Ignore me if it is not relevant.


Thanks bjohn, what I mean is integration over the control volume after which you have expressions with aP, aN, aE etc. What I want is a confirmation of my expressions for these constants as they are much more complex in cylindrical coordinates than rectangular coordinates and I have not seen any book with them. The other option is to approximate the control volume with a rectangular box and just use the more common expression in rectangular coordinates but I dont know if that is visible.

@ arjun, like you, my final aim is not just the flow in a stirred tank. I'm working on crystallisation and just need the flow fields at different time steps to predict my crystal properties and can't just use openFOAM. Even if I can, I'll like to learn how to do it. Any help you can offer will be appreciated. Thanks

arjun July 2, 2010 06:17

Quote:

Originally Posted by falopsy (Post 265418)

@ arjun, like you, my final aim is not just the flow in a stirred tank. I'm working on crystallisation and just need the flow fields at different time steps to predict my crystal properties and can't just use openFOAM. Even if I can, I'll like to learn how to do it. Any help you can offer will be appreciated. Thanks


i think i can help you a lot.

you got to tell me this one thing, is your calculation less than 5 to 10 million cells and are you happy running the simulation on shared memory parallel. (preferably on windows).

If this is true then within a week or two i am going to release inavier version 2, that would solve you navier stokes.
Plus you will be able to add models to it with little efforts.

The solver takes care of navier stokes part you take care of model that you want to add.

arjun July 2, 2010 06:18

if you want, i can send you manual in pdf format. By reading it you will know if it is any use to you or not.

in that case give me your email id.

falopsy July 2, 2010 09:39

Quote:

Originally Posted by arjun (Post 265424)
if you want, i can send you manual in pdf format. By reading it you will know if it is any use to you or not.

in that case give me your email id.

Arjun, I have sent u my email address and will be expecting the manual. Thanks a mil

arjun July 2, 2010 16:34

Quote:

Originally Posted by falopsy (Post 265477)
Arjun, I have sent u my email address and will be expecting the manual. Thanks a mil


I sent you email.


i will be out of town for today and tomorrow so might not mail these two days. (not sure if am taking my laptop with me).

falopsy July 2, 2010 18:25

Thanks, I have seen the message. The iNavier will be very helpful for me for validation. I will surely acknowledge u in my dissertation.My question again is: I have already written an explicit time stepping code but my issue then is the CFL restriction because crystallisation typically runs for hourrs, say 2 hrs (7200s), and the time limitation is such that it will take forever to finish the simulation. That is why I'm trying to rewrite it using implicit time stepping. My questions again are these:

1. Can I choose the time step to be around 1 sec for an implicit scheme. I noticed u used a time step of 1e-3 in the manual, why use something that small for an implicit scheme?
2. is it possible to have steady state in an impeller induced flow? I am asking this because the impeller is always moving and I kind of dont see how the flo w can reach steady state? I dont know if you get what I mean.

arjun July 4, 2010 21:33

Quote:

Originally Posted by falopsy (Post 265537)
Thanks, I have seen the message. The iNavier will be very helpful for me for validation. I will surely acknowledge u in my dissertation.My question again is: I have already written an explicit time stepping code but my issue then is the CFL restriction because crystallisation typically runs for hourrs, say 2 hrs (7200s), and the time limitation is such that it will take forever to finish the simulation. That is why I'm trying to rewrite it using implicit time stepping. My questions again are these:

1. Can I choose the time step to be around 1 sec for an implicit scheme. I noticed u used a time step of 1e-3 in the manual, why use something that small for an implicit scheme?
2. is it possible to have steady state in an impeller induced flow? I am asking this because the impeller is always moving and I kind of dont see how the flo w can reach steady state? I dont know if you get what I mean.


Without looking at the whole thing (equations) it is really difficult to say anything. If you have them in any papers or documents or if you could write down then please email me. I could tell you if using finite volume solver inavier/openfoam would be feasible or not.

Here are few notes though:

1. One of the old versions of inavier, i did solve some curing problems that run for minutes (if not hours) and i used large steps. I did compare my results with polyflow software and they matched well.


2. If you are using cylindrical coordinate and mesh is regular that for uniform poisson equation
nabla^2( phi ) = src ,
you can create a direct solver which is very very fast.

3. Unfortunately when you try to use fully implicit method with SIMPLE type method, the pressure equation is no longer uniform and it is of form:

nabla (k . nabla phi ) = src, where k varies in calculation domain so direct solver is no longer valid.

So it all depends on nature of your equations and your mesh.

falopsy July 7, 2010 08:20

Hi arjun,

Is it possible to work with cylindrical grids but with velocity components in cartesian coordinates. My question is, how do I calculate the area of the faces perpendicular to the velocity components in terms of r, theta and z. I think I will prefer to work in cartesian coordinates because I can easily verify my equations for the coefficients. I have already programmed the solver for my program and have the user interface ready. I just need to make sure the equations for the coefficients are correct. I will send my equations in cylindrical coordinates to your before end of day tomorrow. Thanks a lot for your help.

arjun July 7, 2010 20:55

Quote:

Originally Posted by falopsy (Post 266154)
Hi arjun,

Is it possible to work with cylindrical grids but with velocity components in cartesian coordinates.

Yepp, this is what i did.


Quote:

Originally Posted by falopsy (Post 266154)
Hi arjun,

Is it possible to work with cylindrical grids but with velocity components in cartesian coordinates. My question is, how do I calculate the area of the faces perpendicular to the velocity components in terms of r, theta and z. I think I will prefer to work in cartesian coordinates because I can easily verify my equations for the coefficients. I have already programmed the solver for my program and have the user interface ready. I just need to make sure the equations for the coefficients are correct. I will send my equations in cylindrical coordinates to your before end of day tomorrow. Thanks a lot for your help.


Before we go further, let me clearify one thing that so far i am assuming that you are talking about cylindrical grids which has meshes:
Ni x Nj x Nk ,

where Ni is say theta distribution (0 to 360 deg),
Nj is radial distribution 0 to Rmax
and Nk is z distribution zmin to zmax.

Once we fixed this, you should be able to access every thing by (i,j,k).
That means lets say pressure is given by P(i,j,k).

Now in a finite volume solver each cell has 6 neighbors, those could be accessed by (i,j,k) notation, by i+/-1, j+/-1, k+/-1 .

For discitesation in finite volume sense you will need fluxes on these faces plus you will need area and distance between them.

convection term will be calculated based on discretisation scheme (based on flux).

Diffusional term would be something like = visc_eff_at_face * (Area / dist_between_cells)

Area and dist_between_cells could be easily calcuted, based on your grid.

you can send me equations, as you get time.

lost.identity July 14, 2010 19:34

Hello,

Instead of cylindrical coordinates I'm trying to discretise the governing equations written in spherical coordinates using FVM.

However, I'm considering a spherically symmetric flow (retaining only r-dependence), i.e. it is 1-D. Therefore the general governing equation can be written as

\frac{\partial}{\partial{t}}(\rho\Phi) + \frac{1}{r^2}\frac{\partial}{\partial{r}}(r^2J) = S_{\Phi}

where the total flux J = \rho{u_r}\Phi - \Gamma\partial{\Phi}/\partial{r}.

When I use the FVM in 1-D (only retaining r-dependence) for the following control volume, my discretised equation looks exactly the same as it does for the Cartesian coordinates! This doesn't sound right to me as I don't know where the 1/r^2 terms went.

W---------|w----------P----------|e---------E

However, for a 2-D polar coordinates the equations I obtain were different to the 2-D Cartesian as expected.


For a detailed description of what I'm done, here's how I'm trying to discretise it in 1-D spherical coordinates.

The general conservation equation (in 3-D) can be written as

\frac{\partial}{\partial{t}}\rho\Phi + \nabla\cdot\boldsymbol{J} = S_{\Phi}

where \boldsymbol{J}=\rho\boldsymbol{u}\Phi - \Gamma_{\Phi}\nabla\Phi. I then apply FVM to obtain (note using fully implicit scheme)

\left[ (\rho\Phi)_P-(\rho\Phi)_P^0 \right]\frac{\Delta{V}}{\Delta{t}} + \sum_f\boldsymbol{J}_f\cdot\boldsymbol{A}_f = S_{\Phi}\Delta{V}

where f refers to the cell faces and A are the face area vectors. I think the key here is the face area vectors. So for the 1-D control volume shown above the equation just becomes

\left[(\rho\Phi)_P-(\rho\Phi)_P^0\right]\frac{\Delta{V}}{\Delta{t}} + \boldsymbol{J}_e\cdot\boldsymbol{A}_e + \boldsymbol{J}_w\cdot\boldsymbol{A}_w= S_{\Phi}\Delta{V}

Now this is where I have the biggest doubt, for the control volume above, wouldn't the face area vectors just be

A_e = \hat{\boldsymbol{e}}_r
A_w = -\hat{\boldsymbol{e}}_r

For spherical coordinates the gradient operator is given by

\nabla = \frac{\partial}{\partial{r}}\hat{\boldsymbol{e}}_r + \frac{1}{r}\frac{\partial}{\partial{\theta}}\hat{\boldsymbol{e}}_{\theta} + \frac{1}{r\sin\theta}\frac{\partial}{\partial{\phi}}\hat{\boldsymbol{e}}_{\phi}

which therefore for only r-dependence gives

\boldsymbol{J}_e\cdot\boldsymbol{A}_e = (\rho{u_r}\Phi)_e - \Gamma_e\left(\frac{\partial{\Phi}}{\partial{r}}\right)_e

which is the same as that for 1-D Cartesian coordinates with x replaced by r

And since in 1-D \Delta{V}=\Delta{r}, isn't my discretisation the same as that for Cartesian coordinates?

arjun July 14, 2010 22:36

sorry did not reply to your mail, been very busy.

falopsy July 15, 2010 06:34

Thanks arjun,

Concerning the calculation of area, I wont be able to just use the area of the cylindrical control volume if my velocities are in cartesian coordinates as the will no longer be perpendicular to the cylindrical coordinates. Another issue will be implementation of boundary condition say I want the velocity to be a particular value in radial direction. I hope you understand my concerns.

Anyway, I think I have the right equations now. I have been writing them all in my notebook and will type them up and send to you before end of day today so you can verify for me. I am using the general discretisation scheme written in TVD format and using the deferred correction method and gauss-siedel to solve the algebraic equation. My problem now is that my algorithm seems not to be diverging even when I use the upwind scheme. Is this a normal issue and how do I go around it.

Thanks

falopsy July 15, 2010 19:08

Hi Arjun,

Should I still send my derivation of the equations to you cos I noticed you are yet to reply my last mails. Please let me know if I can send to you. Thanks

arjun July 16, 2010 19:21

Quote:

Originally Posted by falopsy (Post 267554)
Hi Arjun,

Should I still send my derivation of the equations to you cos I noticed you are yet to reply my last mails. Please let me know if I can send to you. Thanks

you can definitely send in, i will look as time permits. I keep very very busy usually. Plus computers at home on which i work are not connected to internet. And the one that is connected is usually taken by my 3year old daughter (she spends lot of time with youtube) so sometimes i need to wait a bit to see my emails etc etc.

:-D


All times are GMT -4. The time now is 10:31.