second-order discretization of viscous stress tensor

selig5576 · August 30, 2017, 12:11

I am discretizing the viscous stress tensor using FDM (as I am going to be performing RANS modelling) and was curious if my formulation is correct.

FMDenaro · August 30, 2017, 12:18

First, the stress tensor has onlt first order derivatives, the outer you wrote is the divergence that applies on the stress component.
Then, it is difficult to see the discretization without a figure of the collocation of the velocity components.
Did you already extracted out the isotropic part of the tensor S?

selig5576 · August 30, 2017, 12:32

The viscous stress tensor I am working with attached. In terms of the velocity-pressure arrangement, I am working on a collocated grid. I will upload a figure of my stencil. By the way, this is for incompressible Navier-Stokes.

FMDenaro · August 30, 2017, 12:51

If you extract out the isotropic part of the modelled stress tensor, it is included in a modified pressure term that you solve directly by the Poisson equation enforcing the divergence-free constraint.
Therefore, you work with Div (mu_t *Grad0 v_bar), being now the tensor tau_xy at zero trace.
Now, assuming you have u,v and mu_t on the same nodes, you start with

d (tau_xy)/dy -> (mu_t(i,j+1/2)*tau_xy(i,j+1/2)-mu_t(i,j-1/2)*tau_xy(i,j-1/2))/dy

you can use a linear reconstruction for mu_t:
mu_t(i,j+1/2)=0.5*(mu_t(i,j+1)+mu_t(i,j))
mu_t(i,j-1/2)=0.5*(mu_t(i,j)+mu_t(i,j-1))

and, accordingly, you discretize the tensor tau_xy and so on

selig5576 · August 30, 2017, 13:33

Thank you for the reply. So if I am understanding you correctly I should get for $\tau_{xy}$

$\frac{\partial}{\partial y} \left(\mu_{T} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right) = \frac{\mu_{i,j+1/2} \tau^{xy}_{i,j+1/2} - \mu_{i,j-1/2} \tau^{xy}_{i,j-1/2}}{dy}$
$\tau^{xy}_{i,j+1/2} = \frac{u_{i,j+1} - u_{i,j}}{dy} + \frac{v_{i+1,j+1} - v_{i,j+1}}{dx}$
$\tau^{xy}_{i,j-1/2} = \frac{u_{i,j} - u_{i,j-1}}{dy} + \frac{v_{i,j-1} - v_{i-1,j-1}}{dx}$

FMDenaro · August 30, 2017, 19:16

selig5576 · August 31, 2017, 10:36

Pardon my stupid question, but if we average the derivatives between [i,i+1] and [i,i-1] then I get
$\frac{\partial}{\partial y} \left(\mu_{T} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right) = \frac{\mu_{i,j+1/2} \tau^{xy}_{i,j+1/2} - \mu_{i,j-1/2} \tau^{xy}_{i,j-1/2}}{dy}$
$\tau^{xy}_{i,j+1/2} = \left(\frac{du}{dy}\right)_{j+1/2} + \left(\frac{dv}{dx}\right)_{i+1/2,j+1/2} = \frac{u_{i,j+1} - u_{i,j}}{dy} + \frac{0.5 \left(v_{i+1,j+1} + v_{i+1,j}\right)}{dx}$
$\tau^{xy}_{i,j-1/2} =\left(\frac{du}{dy}\right)_{j-1/2} + \left(\frac{dv}{dx}\right)_{i-1/2,j-1/2} = \frac{u_{i,j} - u_{i,j-1}}{dy} + \frac{0.5 \left(v_{i-1,j} + v_{i-1,j-1} \right)}{dx}$

FMDenaro · August 31, 2017, 11:51

why are you evaluating dv/dx at half node instead of using i,j+1/2 and i,j-1/2 as for du/dy?

However, I suggest to work in a different manner much more similar to the FV formulation.
Define the flux_ovest(i,j), flux_sud(i,j) matrices (for 2d, in 3D you need a further matrix). Use the staggered notation, for example flux_ovest(i,j) is at the face i-1/2,j and flux_s(i,j) is at i,j-1/2. The compute all the flux once over all the half nodes.
Finally simply sum the fluxes according to

flux_ovest(i+1,j)-flux_ovest(i,j)
flux_sud(i,j+1)-flux_sud(i,j)

selig5576 · August 31, 2017, 12:22

I looked in Ferziger's book on FVM and found something related to discretizing tau_xy. Despite his assumption on constant mu, could I adopt this discretization procedure? I realize that since my viscosity is not constant I will have to interpolate accordingly. Given he is working in a FVM formulation I would have (tau_yx)_east - (tau_yx)_west

EDIT: To be more specific I am currently trying

$(\tau_{yx})_{e} = \left( \frac{v_{I+1,j,k} - v_{i,j,k}}{dx} +\frac{u_{i+1,j+1} - u_{i+1,j-1}}{dy}\right)$
$(\tau_{yx})_{w} = \left(\frac{v_{i,j,k} - v_{i-1,j,k}}{dx} + \frac{u_{i-1,j+1} - u_{i-1,j-1}}{dy} \right)$
Therefore we have
$\frac{ (\tau_{yx})_{e} - (\tau_{yx})_{w}}{dx}$

FMDenaro · August 31, 2017, 12:47

Quote:

Originally Posted by selig5576

I looked in Ferziger's book on FVM and found something related to discretizing tau_xy. Despite his assumption on constant mu, could I adopt this discretization procedure? I realize that since my viscosity is not constant I will have to interpolate accordingly. Given he is working in a FVM formulation I would have (tau_yx)_east - (tau_yx)_west

EDIT: To be more specific I am currently trying

$(\tau_{yx})_{e} = \left( \frac{v_{I+1,j,k} - v_{i,j,k}}{dx} +\frac{u_{i+1,j+1} - u_{i+1,j-1}}{dy}\right)$
$(\tau_{yx})_{w} = \left(\frac{v_{i,j,k} - v_{i-1,j,k}}{dx} + \frac{u_{i-1,j+1} - u_{i-1,j-1}}{dy} \right)$
Therefore we have
$\frac{ (\tau_{yx})_{e} - (\tau_{yx})_{w}}{dx}$

But Eq.(7.90) is presented for a staggered arrangement while you are using a colocated one.

Compute only either tau_xy_est or tau_xy_west an memorize the value in a matrix. Owing to the conservative properties you need nothing else to do the update.

selig5576 · August 31, 2017, 13:30

You are 100% right, I was too hasty and overlooked that.

Quote:

Not to go down a rabbit hole. I have a question: Was the formulation below correct except I did not average du/dy + dv/dx by 1/2? Sorry for my clear confusion

$\frac{\partial}{\partial y} \left(\mu_{T} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right) = \frac{\mu_{i,j+1/2} \tau^{xy}_{i,j+1/2} - \mu_{i,j-1/2} \tau^{xy}_{i,j-1/2}}{\mathbf{2} dy}$
$\tau^{xy}_{i,j+1/2} = \frac{u_{i,j+1} - u_{i,j}}{dy} + \frac{v_{i+1,j+1} - v_{i,j+1}}{dx}$
$\tau^{xy}_{i,j-1/2} = \frac{u_{i,j} - u_{i,j-1}}{dy} + \frac{v_{i,j-1} - v_{i-1,j-1}}{dx}$

In terms of a conservation law, when taking the flux face it is composed of half points. Thank you for your patience...

FMDenaro · August 31, 2017, 14:35

yes, just correct the dv/dx derivative and the discrete formula is right... remember that you do not need to compute twice the same term on a face

selig5576 · August 31, 2017, 15:16

After doing some paper work I finally get

$\frac{\partial}{\partial y} \left(\mu_{T} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right) = \frac{\mu_{i,j+1/2} \tau^{xy}_{i,j+1/2} - \mu_{i,j-1/2} \tau^{xy}_{i,j-1/2}}{dy}$
$\tau^{xy}_{i,j+1/2} = \frac{u_{i,j+1} - u_{i,j}}{dy} + \frac{v_{i+1,j+1} - v_{i,j+1}}{2 dx}$
$\tau^{xy}_{i,j-1/2} = \frac{u_{i,j} - u_{i,j-1}}{dy} + \frac{v_{i,j-1} - v_{i-1,j-1}}{2 dx}$

FMDenaro · August 31, 2017, 15:45

I don't agree, in my mind I have this discretization

dv/dx|i,j+1/2 => 0.25*(v(i+1,j+1)-v(i-1,j+1)+v(i+1,j)-v(i-1,j))/h

August 30, 2017, 13:33	Discretization	#5
selig5576 Senior Member Selig Join Date: Jul 2016 Posts: 213 Rep Power: 10	Thank you for the reply. So if I am understanding you correctly I should get for $\tau_{xy}$ $\frac{\partial}{\partial y} \left(\mu_{T} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right) = \frac{\mu_{i,j+1/2} \tau^{xy}_{i,j+1/2} - \mu_{i,j-1/2} \tau^{xy}_{i,j-1/2}}{dy}$ $\tau^{xy}_{i,j+1/2} = \frac{u_{i,j+1} - u_{i,j}}{dy} + \frac{v_{i+1,j+1} - v_{i,j+1}}{dx}$ $\tau^{xy}_{i,j-1/2} = \frac{u_{i,j} - u_{i,j-1}}{dy} + \frac{v_{i,j-1} - v_{i-1,j-1}}{dx}$ Last edited by selig5576; August 31, 2017 at 13:09.

August 31, 2017, 10:36	Viscous stress tensor	#7
selig5576 Senior Member Selig Join Date: Jul 2016 Posts: 213 Rep Power: 10	Pardon my stupid question, but if we average the derivatives between [i,i+1] and [i,i-1] then I get $\frac{\partial}{\partial y} \left(\mu_{T} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right) = \frac{\mu_{i,j+1/2} \tau^{xy}_{i,j+1/2} - \mu_{i,j-1/2} \tau^{xy}_{i,j-1/2}}{dy}$ $\tau^{xy}_{i,j+1/2} = \left(\frac{du}{dy}\right)_{j+1/2} + \left(\frac{dv}{dx}\right)_{i+1/2,j+1/2} = \frac{u_{i,j+1} - u_{i,j}}{dy} + \frac{0.5 \left(v_{i+1,j+1} + v_{i+1,j}\right)}{dx}$ $\tau^{xy}_{i,j-1/2} =\left(\frac{du}{dy}\right)_{j-1/2} + \left(\frac{dv}{dx}\right)_{i-1/2,j-1/2} = \frac{u_{i,j} - u_{i,j-1}}{dy} + \frac{0.5 \left(v_{i-1,j} + v_{i-1,j-1} \right)}{dx}$

August 31, 2017, 11:51		#8
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,785 Rep Power: 71	why are you evaluating dv/dx at half node instead of using i,j+1/2 and i,j-1/2 as for du/dy? However, I suggest to work in a different manner much more similar to the FV formulation. Define the flux_ovest(i,j), flux_sud(i,j) matrices (for 2d, in 3D you need a further matrix). Use the staggered notation, for example flux_ovest(i,j) is at the face i-1/2,j and flux_s(i,j) is at i,j-1/2. The compute all the flux once over all the half nodes. Finally simply sum the fluxes according to flux_ovest(i+1,j)-flux_ovest(i,j) flux_sud(i,j+1)-flux_sud(i,j) selig5576 likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
calculating the viscous stress tensor tau for viscous dissipation heating	dgadensg	OpenFOAM Programming & Development	5	January 17, 2023 15:20
post-processing utility R, calculating Reynolds stress tensor	volker	OpenFOAM Post-Processing	17	May 21, 2020 17:36
Reynolds Stress Tensor	Sergio13	OpenFOAM Post-Processing	12	May 11, 2020 05:40
Convergence of Reynolds Stress Tensor / Fluctuations in OpenFOAM	Anne Lincke	OpenFOAM Running, Solving & CFD	0	July 17, 2012 10:50
Higher order discretization on staggered grid	Chandra Shekhar	Main CFD Forum	9	January 27, 2005 16:31

August 30, 2017, 12:18		#2
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,785 Rep Power: 71	First, the stress tensor has onlt first order derivatives, the outer you wrote is the divergence that applies on the stress component. Then, it is difficult to see the discretization without a figure of the collocation of the velocity components. Did you already extracted out the isotropic part of the tensor S?

August 30, 2017, 12:51		#4
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,785 Rep Power: 71	If you extract out the isotropic part of the modelled stress tensor, it is included in a modified pressure term that you solve directly by the Poisson equation enforcing the divergence-free constraint. Therefore, you work with Div (mu_t Grad0 v_bar), being now the tensor tau_xy at zero trace. Now, assuming you have u,v and mu_t on the same nodes, you start with d (tau_xy)/dy -> (mu_t(i,j+1/2)tau_xy(i,j+1/2)-mu_t(i,j-1/2)tau_xy(i,j-1/2))/dy you can use a linear reconstruction for mu_t: mu_t(i,j+1/2)=0.5(mu_t(i,j+1)+mu_t(i,j)) mu_t(i,j-1/2)=0.5*(mu_t(i,j)+mu_t(i,j-1)) and, accordingly, you discretize the tensor tau_xy and so on

August 30, 2017, 19:16		#6
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,785 Rep Power: 71	be careful, the derivative dv/dx cannot computed without the averaging between two derivatives. For example dv/dx\|j+1/2 = 0.5(dv/dx\|j+1 + dv/dx\|j) dv/dx\|j-1/2 = 0.5(dv/dx\|j-1 + dv/dx\|j)

August 31, 2017, 14:35		#12
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,785 Rep Power: 71	yes, just correct the dv/dx derivative and the discrete formula is right... remember that you do not need to compute twice the same term on a face

August 31, 2017, 15:16		#13
selig5576 Senior Member Selig Join Date: Jul 2016 Posts: 213 Rep Power: 10	After doing some paper work I finally get $\frac{\partial}{\partial y} \left(\mu_{T} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)\right) = \frac{\mu_{i,j+1/2} \tau^{xy}_{i,j+1/2} - \mu_{i,j-1/2} \tau^{xy}_{i,j-1/2}}{dy}$ $\tau^{xy}_{i,j+1/2} = \frac{u_{i,j+1} - u_{i,j}}{dy} + \frac{v_{i+1,j+1} - v_{i,j+1}}{2 dx}$ $\tau^{xy}_{i,j-1/2} = \frac{u_{i,j} - u_{i,j-1}}{dy} + \frac{v_{i,j-1} - v_{i-1,j-1}}{2 dx}$

August 31, 2017, 15:45		#14
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,785 Rep Power: 71	I don't agree, in my mind I have this discretization dv/dx\|i,j+1/2 => 0.25*(v(i+1,j+1)-v(i-1,j+1)+v(i+1,j)-v(i-1,j))/h