Source term linearization

lost.identity · April 26, 2010, 08:40

Hi,

I'm trying to linearize a source term in my momentum equations written in spherical coordinates. I'm trying to solve a 1-D problem.

The source term is given by

$\frac{\partial}{\partial{r}}(r\Gamma{u})$

where $\Gamma$ is the diffusivity for example,

is the radial coordinate and

is the velocity in the r-direction.

Using Finite Volume discretization what I've tried to do is integrate this over the control volume to get

$\int_w^e \frac{\partial}{\partial{r}}(r\Gamma{u})\, dr =(r\Gamma{u})_e-(r\Gamma{u})_w$

where

and

are the faces of the control volume. The usual source term linearization, for example, by Patankar (1980) is $S=S_c+S_P{u_P}$ .

The problem I have is that the source term linearization above has the term ${u_P}$ , which is the value of

at the centre of the control volume (i.e. at grid point

) but from my control volume integration I get the values ${u_e}$ and

which are the face values. I know that I can evaluate the face values in terms of the centre coordinates but then this makes things a lot more complicated.

Hamidzoka · April 27, 2010, 02:03

hi,
"Up" does not necessarily equal to "U" in position of "P".
it deoends the upper and lower bounds of your integration.
for example if you integrate the momentum equation to find discrete momentum equation for "e", your source term will have the form of S=Sc+Se*Ue.
and similarly for w, you will find source term to be: S=Sc+Sw*Uw.

does it help?

lost.identity · April 27, 2010, 10:19

Thanks for the reply. Actually that's the problem, because the lower and upper bounds of my integration has to be the cell faces, which are e and w. P is the centre point of this cell.

Imagine the control volume is shown as below.

w |----------P----------| e

I integrate the momentum equation over this control volume. So the unsteady term for example will be of the form

$\int_t^{t+\Delta{t}}\,\int_e^w \frac{\partial{\rho{u}}}{\partial{t}}\, dt\,dx = (\rho_P^1{u}_P - \rho_P{u}_P^0)\Delta{x}$

where the subscript 1 is the current time step and 0 is the previous and delta_x is the grid size. Here it is assumed that the grid-point value of u (i.e. the value at P) prevails over the control volume.

My problem is, as I mentioned in my first post when I integrate the particular source term that I have, I get

(because the source term is actually in the form of a spatial gradient). If this derivative wasn't there I'd get

like the unsteady term above.

In the original form by Patankar, the full discretised equation will be like

(a_E + a_W - S_P)u_P = a_Eu_E + a_Wu_W + Sc

so the Sp term needs to be multiplied by the value of u at the same point (whether it may be e, w or P) as

and

. But I don't know how I could get that. I've seen a number of papers where they've done this for cylinderical coordinates and it should be similar but they never show what the form of Sp they have.

Hamidzoka · April 28, 2010, 00:46

Hi;
you problem seems to be quite fundamental!
firstly, you should use staggered contorl when integrating momentum equations.
this method is utilized to prevent an unrealistic pressure field. consider the following domain:
W----w----P----e----E
you should use one of the following CVs for x-momentum discretization:

W|-----w-----|P (backward staggered)

or

P|-----e------|E (forward staggered)

the control volume you discussed before is just used for scalar equations discretization such as mass, energy, K, epsilon, etc.

Secondly, when you discrete x-momentum equation, convection terms(the right side of momentum equation should not be integrated and probably interpolated just like the way you did. when you have convection term you should use "hybrid", "power-low" or similar methods for stability considerations. after doing this you can organize your source terms.
for example if you use forward scheme, all source terms containing "Ue" will be considered as "SeUe" term(just like your SpUp) and the rest (such as "Uw") will be considered as "Sc". (for more information refer to CFD books like the one written by Patankar)

dose it help?
regards

lost.identity · May 4, 2010, 08:12

Hi,

Thanks for the reply. Firstly even though I drew the control volume as above, I'm actually using a staggered CV for the momentum equations.

Secondly I am trying to use the power law for the convection terms, for the Cartesian equations this is straightforward but for spherical equations I have additional geometric source terms that's giving me the problem.

The full equation which I'm trying to discretise is

$r^2\frac{\partial\overline\rho\tilde{u}}{\partial{t}} + \frac{\partial{J}}{\partial{r}} = -r^2\frac{\partial{\overline{p}}}{\partial{r}} + S$

where
$J = r^2\left( \overline\rho\tilde{u}\tilde{u} - \Gamma\frac{\partial\tilde{u}}{\partial{r}} \right)$

$\Gamma = \frac{4}{3}(\mu_t + \mu)$

$S = -\frac{\partial}{\partial{r}}\left( r^2\Gamma\frac{\tilde{u}}{r}- r^2\frac{2}{3}\overline\rho\tilde{k} \right)$

The only problem that I have is the term $-\frac{\partial}{\partial{r}}r^2\Gamma\frac{\tilde{u}}{r}$ in the source S. This term is a bit like the convective term but I don't know exactly how I would treat it.

I was thinking of putting this in the total convective-diffusive flux term J but then I don't know how to obtain the powerlaw from there. Also this is different to what I've seen in some of the papers. They have the convective-diffusive flux and the source terms as shown above.

As I mentioned in my previous post, the problem is when I discretise this source term. Using the forward staggered grid below

w---------P|-----e------|E--------e1

I get,

$S_P = \Gamma_Er_Eu_E - \Gamma_Pr_Pu_P$

i.e. the velocities here are at the cell faces.

However, the discretisation of the convective terms only gives me the velocity

which are the values at the cell centre. So how can I form the source term when this is the case? one of them is at the cell centre and the other at the cell faces?

I could use the interpolation of the cell centres to get the cell face values for example

$u_E = \frac{2u_eu_{e1}}{u_e+u_{e1}}$

but then I still don't know how I could factorise this interms of

.

George · May 7, 2010, 08:46

Hello,

refering to your first post, you can express the Ue, Uw velocities like: Ue=(UP+UE)/2. With this way when you introduce the above equations in the source term, the variable UP appears (together with the UE), which is then inserted in the Sp.

Hope it helps a little bit.

George

lost.identity · May 8, 2010, 10:56

I realise that I'm going about this the wrong way. I just have to collect the source terms together and multiply it by the volume of the cell.

However, I'm still struggling with the linearization. I've seen a lot of papers and the book by Patankar that explains you how to do this using Taylor series expansion of the source term.

For example

$S = S^* + \left(\frac{\partial{S}}{\partial{u}}\right)^*(u-u*)$

where the superscript * refers to the value at the previous iteration.

My problem now is with the differentiation since
$S = -\frac{\partial}{\partial{r}}\left( r^2\Gamma\frac{\tilde{u}}{r}- r^2\frac{2}{3}\overline\rho\tilde{k} \right)$
$\Gamma = \frac{4}{3}(\mu_t + \mu)$

The eddy viscosity $\mu_t$ is calculated using $\mu_t=C_{\mu}\rho\frac{k^2}{\epsilon}$ where the variables $\rho,k,\epsilon$ are in general functions of

, this would create a more complex source term with a lot of terms and I'm not sure whether this is the right way to go about this?

George · May 8, 2010, 12:07

Hi,

on the other hand, you can use a different methodology, where in the source term "Sc" you include all the terms calculated from the previous iteration. This method is not so robust, since it might need some more time in order to converge. Stability issues might also arise. But according to the problem complexity, it can end up to the same calculations.

siddiquesil · October 10, 2019, 11:29

Quote:

Originally Posted by Hamidzoka

hi,
"Up" does not necessarily equal to "U" in position of "P".
it deoends the upper and lower bounds of your integration.
for example if you integrate the momentum equation to find discrete momentum equation for "e", your source term will have the form of S=Sc+Se*Ue.
and similarly for w, you will find source term to be: S=Sc+Sw*Uw.

does it help?

I am also trying to solve the N.S equation in spherical co-ordinate.
I have a similar type of source terms in the r-momentum equation which is in derivate form.
But as per my last forum post, I was informed that there is no extra source term in the r-momentum equation apart from physical forces which are pressure, centrifugal, Coriolis, gravitation, etc.

check the link of my last post
How to discretize derivative Source Term??

But I have picked up the r-momentum equation from the wiki page which has lots of source terms in it apart from physical forces. The link is here.
https://en.wikipedia.org/wiki/Navier–Stokes_equations

My questions are as follows.

1. Is the conservative r-momentum equation in spherical co-ordinate given in the wiki page correct?

2. Can you provide me any other link for the conservative N.S equation in spherical co-ordinate?

3. Are there extra terms in momentum equation that arise apart from physical forces for spherical co-ordinate while solving by using FVM?

April 26, 2010, 08:40	Source term linearization	#1
lost.identity Senior Member Join Date: Apr 2009 Posts: 118 Rep Power: 17	Hi, I'm trying to linearize a source term in my momentum equations written in spherical coordinates. I'm trying to solve a 1-D problem. The source term is given by $\frac{\partial}{\partial{r}}(r\Gamma{u})$ where $\Gamma$ is the diffusivity for example, is the radial coordinate and is the velocity in the r-direction. Using Finite Volume discretization what I've tried to do is integrate this over the control volume to get $\int_w^e \frac{\partial}{\partial{r}}(r\Gamma{u})\, dr =(r\Gamma{u})_e-(r\Gamma{u})_w$ where and are the faces of the control volume. The usual source term linearization, for example, by Patankar (1980) is $S=S_c+S_P{u_P}$ . The problem I have is that the source term linearization above has the term ${u_P}$ , which is the value of at the centre of the control volume (i.e. at grid point ) but from my control volume integration I get the values ${u_e}$ and which are the face values. I know that I can evaluate the face values in terms of the centre coordinates but then this makes things a lot more complicated.

April 27, 2010, 10:19		#3
lost.identity Senior Member Join Date: Apr 2009 Posts: 118 Rep Power: 17	Thanks for the reply. Actually that's the problem, because the lower and upper bounds of my integration has to be the cell faces, which are e and w. P is the centre point of this cell. Imagine the control volume is shown as below. w \|----------P----------\| e I integrate the momentum equation over this control volume. So the unsteady term for example will be of the form $\int_t^{t+\Delta{t}}\,\int_e^w \frac{\partial{\rho{u}}}{\partial{t}}\, dt\,dx = (\rho_P^1{u}_P - \rho_P{u}_P^0)\Delta{x}$ where the subscript 1 is the current time step and 0 is the previous and delta_x is the grid size. Here it is assumed that the grid-point value of u (i.e. the value at P) prevails over the control volume. My problem is, as I mentioned in my first post when I integrate the particular source term that I have, I get (because the source term is actually in the form of a spatial gradient). If this derivative wasn't there I'd get like the unsteady term above. In the original form by Patankar, the full discretised equation will be like so the Sp term needs to be multiplied by the value of u at the same point (whether it may be e, w or P) as and . But I don't know how I could get that. I've seen a number of papers where they've done this for cylinderical coordinates and it should be similar but they never show what the form of Sp they have.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
momentum source term	zwdi	FLUENT	14	June 27, 2017 16:40
DxFoam reader update	hjasak	OpenFOAM Post-Processing	69	April 24, 2008 02:24
DecomposePar links against liblamso0 with OpenMPI	jens_klostermann	OpenFOAM Bugs	11	June 28, 2007 18:51
UDF Scalar Code: HT 1	Greg Perkins	FLUENT	8	October 20, 2000 13:40
UDFs for Scalar Eqn - Fluid/Solid HT	Greg Perkins	FLUENT	0	October 11, 2000 04:43

April 27, 2010, 02:03		#2
Hamidzoka Senior Member Hamid Zoka Join Date: Nov 2009 Posts: 293 Rep Power: 19	hi, "Up" does not necessarily equal to "U" in position of "P". it deoends the upper and lower bounds of your integration. for example if you integrate the momentum equation to find discrete momentum equation for "e", your source term will have the form of S=Sc+SeUe. and similarly for w, you will find source term to be: S=Sc+SwUw. does it help?

April 28, 2010, 00:46		#4
Hamidzoka Senior Member Hamid Zoka Join Date: Nov 2009 Posts: 293 Rep Power: 19	Hi; you problem seems to be quite fundamental! firstly, you should use staggered contorl when integrating momentum equations. this method is utilized to prevent an unrealistic pressure field. consider the following domain: W----w----P----e----E you should use one of the following CVs for x-momentum discretization: W\|-----w-----\|P (backward staggered) or P\|-----e------\|E (forward staggered) the control volume you discussed before is just used for scalar equations discretization such as mass, energy, K, epsilon, etc. Secondly, when you discrete x-momentum equation, convection terms(the right side of momentum equation should not be integrated and probably interpolated just like the way you did. when you have convection term you should use "hybrid", "power-low" or similar methods for stability considerations. after doing this you can organize your source terms. for example if you use forward scheme, all source terms containing "Ue" will be considered as "SeUe" term(just like your SpUp) and the rest (such as "Uw") will be considered as "Sc". (for more information refer to CFD books like the one written by Patankar) dose it help? regards

May 4, 2010, 08:12		#5
lost.identity Senior Member Join Date: Apr 2009 Posts: 118 Rep Power: 17	Hi, Thanks for the reply. Firstly even though I drew the control volume as above, I'm actually using a staggered CV for the momentum equations. Secondly I am trying to use the power law for the convection terms, for the Cartesian equations this is straightforward but for spherical equations I have additional geometric source terms that's giving me the problem. The full equation which I'm trying to discretise is $r^2\frac{\partial\overline\rho\tilde{u}}{\partial{t}} + \frac{\partial{J}}{\partial{r}} = -r^2\frac{\partial{\overline{p}}}{\partial{r}} + S$ where $J = r^2\left( \overline\rho\tilde{u}\tilde{u} - \Gamma\frac{\partial\tilde{u}}{\partial{r}} \right)$ $\Gamma = \frac{4}{3}(\mu_t + \mu)$ $S = -\frac{\partial}{\partial{r}}\left( r^2\Gamma\frac{\tilde{u}}{r}- r^2\frac{2}{3}\overline\rho\tilde{k} \right)$ The only problem that I have is the term $-\frac{\partial}{\partial{r}}r^2\Gamma\frac{\tilde{u}}{r}$ in the source S. This term is a bit like the convective term but I don't know exactly how I would treat it. I was thinking of putting this in the total convective-diffusive flux term J but then I don't know how to obtain the powerlaw from there. Also this is different to what I've seen in some of the papers. They have the convective-diffusive flux and the source terms as shown above. As I mentioned in my previous post, the problem is when I discretise this source term. Using the forward staggered grid below w---------P\|-----e------\|E--------e1 I get, $S_P = \Gamma_Er_Eu_E - \Gamma_Pr_Pu_P$ i.e. the velocities here are at the cell faces. However, the discretisation of the convective terms only gives me the velocity which are the values at the cell centre. So how can I form the source term when this is the case? one of them is at the cell centre and the other at the cell faces? I could use the interpolation of the cell centres to get the cell face values for example $u_E = \frac{2u_eu_{e1}}{u_e+u_{e1}}$ but then I still don't know how I could factorise this interms of .

May 7, 2010, 08:46		#6
George New Member George Kosmadakis Join Date: Apr 2009 Posts: 3 Rep Power: 17	Hello, refering to your first post, you can express the Ue, Uw velocities like: Ue=(UP+UE)/2. With this way when you introduce the above equations in the source term, the variable UP appears (together with the UE), which is then inserted in the Sp. Hope it helps a little bit. George

May 8, 2010, 10:56		#7
lost.identity Senior Member Join Date: Apr 2009 Posts: 118 Rep Power: 17	I realise that I'm going about this the wrong way. I just have to collect the source terms together and multiply it by the volume of the cell. However, I'm still struggling with the linearization. I've seen a lot of papers and the book by Patankar that explains you how to do this using Taylor series expansion of the source term. For example $S = S^* + \left(\frac{\partial{S}}{\partial{u}}\right)^(u-u)$ where the superscript * refers to the value at the previous iteration. My problem now is with the differentiation since $S = -\frac{\partial}{\partial{r}}\left( r^2\Gamma\frac{\tilde{u}}{r}- r^2\frac{2}{3}\overline\rho\tilde{k} \right)$ $\Gamma = \frac{4}{3}(\mu_t + \mu)$ The eddy viscosity $\mu_t$ is calculated using $\mu_t=C_{\mu}\rho\frac{k^2}{\epsilon}$ where the variables $\rho,k,\epsilon$ are in general functions of , this would create a more complex source term with a lot of terms and I'm not sure whether this is the right way to go about this?

May 8, 2010, 12:07		#8
George New Member George Kosmadakis Join Date: Apr 2009 Posts: 3 Rep Power: 17	Hi, on the other hand, you can use a different methodology, where in the source term "Sc" you include all the terms calculated from the previous iteration. This method is not so robust, since it might need some more time in order to converge. Stability issues might also arise. But according to the problem complexity, it can end up to the same calculations.