Computing the deviatoric stress tensor

Lattice Boltzmann Method LBM: Theory
#1

Hello,

I would like to compute the deviatoric stress tensor PI with LB. For BGK, through chapman-enskog expansion if f[sub]i[/sub] = f[sub]i[/sub][sup]eq[/sup] + epsilonf[sub]i[/sub]sup[/sup] + … is possible to show that PI[sub]xy[/sub] = e[sub]ix[/sub]e[sub]iy[/sub](1 - 1/(2tau)f[sub]i[/sub]sup[/sup] (summation over i), right? I’m wondering if computing e[sub]ix[/sub]e[sub]iy[/sub](1 - 1/(2tau)*(f[sub]i[/sub] - f[sub]i[/sub][sup]eq[/sup]) with LB is a proper approximation of deviatoric stress tensor.

Any suggestions?

Thanks

#2

Hi,
indeed it is the right way. Of course you’ll do a slight error since you’ll also have contributions from the order 2 and above terms of the f_i, but these terms are considered to be too small and therefore neglected. Another way of computing it would be

Pi^{tot}{alpha beta}=sum e{i alpha} e_{i beta} f_i
Pi^{0}{alpha beta}=sum e{i alpha} e_{i beta} f_i^{eq}=rhou_{alpha}u_{beta} + c_s^2rhodelta_{alpha beta}

PI_{alpha beta} = Pi^{tot}{alpha beta} - Pi^{0}{alpha beta}

Orestis

#3

Hello,

I have another question concerning the deviatoric stress tensor.
I would like to calculate the tensor
$\sigma_{xy} = \eta (d u_x / dy + d u_y / dx)$.
This can be done for example using (Shaq’s post)
$\sigma_{xy} = (1 - 1 / (2 \tau)) * \sum_i c_{ix} * c_{iy} * (f_i - f_i^eq)$.
Or it can be done using the expression (malaspin’s post)
$\sigma_{xy} = \sum_i c_{ix} * c_{iy} * f_i - \rho * u_x * u_y$.
But isn’t in the latter equation the viscosity missing? I have to put the term $1 - 1 / (2 \tau)$ somewhere. Is this one,
$\sigma_{xy} = (1 - 1 / (2 \tau)) * [\sum_i c_{ix} * c_{iy} * f_i - \rho * u_x * u_y]$,
correct?

Thanks for your help,
Timm

#4

Hi,
I don’t think that you have to add the (1-1/(2tau)) factor. Pi is given by :

Pi=pI+rhouu+2mu*S

where p=c_s^2rho, mu is the dynamic viscosity and S=1/2(grad u + (grad u)^T) the strain rate tensor. Therefore if you want the 2muS part you only have to substract the pI+rhou*u part.

I hope it is clear…

Orestis

#5

Thanks Orestis,

your explanation is pretty clear. However, I have some difficulties comparing the theoretical results of my model with the simulation. There is a factor missing somewhere. I will check this and post again if there are problems left…

Timm

#6

Hello,

here is the problem I have:
I have implemented the shear stress \Pi_{xy} in my code as
$\Pi_{xy} = - \sum_i c_{ix} * c_{iy} * f_i + \rho * u_x * u_y$
without any relaxation parameter dependent terms (as you have suggested, Orestis).
The shape of the profile exactly fits the theoretical profile
$\Pi_{xy} = \nu * \rho * d u_x / d y$ (u_y = 0 in my case)
\nu is the kinematic viscosity here.
But: The two profiles differ by a factor which is NOT a constant for different values of $\tau$. I have run three tests
1: \tau = 0.73094; \nu = 0.07698; factor = 3.08
2: \tau = 0.84641; \nu = 0.11547; factor = 2.38
3: \tau = 0.93301; \nu = 0.14430; factor = 2.10
The factor is the simulated values divided by the analytic values of the shear stress. Now it seems that this factor obeys
$factor = \tau / (3 * \nu)$
Since the velocity profiles always match outstandingly, I think that there is an error calculating the shear stress. But I do not know why this is going wrong.
Does anybody have a flash of genius?

Thank you,
Timm

#7

Maybe if you post the code it will be easier to post an useful answer…

#8

This is the code for calculating the shear stress from the simulated f_i values (for explanations see below):


void computeDevStress()
{
for(k = k1; k < k2; k++)
    {
    if(obstacle[k] == 0)
        {
        devStressXY[k] = 0;
        devStressXZ[k] = 0;
        }
    }

for(k = k1; k < k2; k++)
    {
    if(obstacle[k] == 0)
        {
        for(ii = 0; ii < NRVELO; ii++)
            {
            devStressXY[k] -= vl[ii][0] * vl[ii][1] * ff[ii][k];
            devStressXZ[k] -= vl[ii][0] * vl[ii][2] * ff[ii][k];
            }

        devStressXY[k] += nn[k] * uxs[k] * uys[k];
        devStressXZ[k] += nn[k] * uxs[k] * uzs[k];
        }
    }
}

This is the analytic solution


double devStressXY(double y, double z)
{
double stress = 0;

for(n = 1; n < NUMCON1; n += 2)
    {
    for(m = 1; m < NUMCON1; m += 2)
        stress += 32 / SQ3(M_PI) / m / LY / (SQ(n / LY) + SQ(m / LZ))
            * cos(n * M_PI * y / LY) * sin(m * M_PI * z / LZ);
    }

return nuLB * stress * umaxLB / norm;
}

The velocity is given by


double velX(double y, double z)
{
double velocity = 0;

for(n = 1; n < NUMCON1; n += 2)
    {
    for(m = 1; m < NUMCON1; m += 2)
        velocity += 32 / SQ4(M_PI) / n / m / (SQ(n / LY) + SQ(m / LZ))
            * sin(n * M_PI * y / LY) * sin(m * M_PI * z / LZ);
    }
}

return velocity * umaxLB / norm;
}

‘norm’ is a common normalization factor, ‘umaxLB’ is the velocity of the flow at the center of the rectangular channel, and obviously the shear stress is the derivative of the velocity with respect to y. ‘SQ3’ and ‘SQ4’ are just the powers of the argument. ‘LY’ and ‘LZ’ are the y- and z-extensions of the channel (flow in x-direction), NUMCON1 is the truncation limit of the sum.

#9

I should also add how I calculate the LB values of \tau and \nu.


nuLB = nu * dt / h / h;
tau = 3 * nuLB + 0.5;

‘h’ is the lattice spacing, so hLB = 1, and ‘dt’ is the time step, so dtLB = 1. This means that for a given physical viscosity (say water, nu = 10e-6 m^2 / s) and a fixed lattice spacing h and time step dt the LB values of the viscosity and the relaxation parameter are fixed.

#10

Actually I just found out that if I would use


\Pi_{xy} = [- \sum_i c_{ix} * c_{iy} * f_i + \rho * u_x * u_y] * [1 - 1 / (2 * \tau)]

instead of


\Pi_{xy} = - \sum_i c_{ix} * c_{iy} * f_i + \rho * u_x * u_y

for the simulated value, then everything would just be fine…
Orestis, what do you think? Where is my error in reasoning? Why this additional factor $1 - 1 / (2 * \tau)$?

Timm

#11

That’s something that seems strage to me… Let me ask you two other questions.

When do you compute the Pi_{xy}? After propagation (or before collision)?

And what do you mean by “everything is fine”? What is the order of magnitude of the error and how does is scale when enhancing the accuracy?

#12

Hello Orestis,

let me first ask another question: Given the shear stress Pi in SI units (kg / m^3 * m^2 / s * 1 / s = kg / m^3 m^2 / s^2), am I correct that the LB value PiLB then is PiLB = Pi * dt^2 / dx^2? I have set the density to 1 in lattice units. dt and dx are the time step and the lattice spacing, respectively.
One example with numbers: let the value of Pi at some point be 0.02 kg / m^3 m^2 / s^2 and dx = 10e-5 m and dt = 10e-6 s, then PiLB = 0.0002. Am i right?
Maybe this is the source of my error.

Concerning your first question:
collision => BB => propagation => calc. of velocity and density => calc. of shear stress => collision
For the calculation of the shear stress I thus take the same f_i values as for the calculation of the velocity and density.

Timm

PS: I am just running some additional simulations in order to answer your second question.

#13

Hi Orestis,

I have compared the results of the two methods


for(k = k1; k < k2; k++)
    {
    if(obstacle[k] == 0)
        {
        for(ii = 0; ii < NRVELO; ii++)
            {
            devStressXY[k] += vl[ii][0] * vl[ii][1] * ff[ii][k];
            devStressXZ[k] += vl[ii][0] * vl[ii][2] * ff[ii][k];
            }

        devStressXY[k] -= nn[k] * uxs[k] * uys[k];
        devStressXZ[k] -= nn[k] * uxs[k] * uzs[k];
        devStressXY[k] *= (1 - 1 / 2 / tau);
        devStressXZ[k] *= (1 - 1 / 2 / tau);
        }
    }

and


for(k = k1; k < k2; k++)
    {
    if(obstacle[k] == 0)
        {
        equilibrium();

        for(ii = 0; ii < NRVELO; ii++)
            {
            devStressXY[k] += (1 - 1 / 2 / tau) * vl[ii][0] * vl[ii][1] * (ff[ii][k] - fe[ii]);
            devStressXZ[k] += (1 - 1 / 2 / tau) * vl[ii][0] * vl[ii][2] * (ff[ii][k] - fe[ii]);
            }
        }
    }

And guess what the result is! Both methods exactly result in the same values, the deviations are 0. This means, if the second approach is correct, then the first one is also correct. But it differs from your proposal by a factor of (1 - 1 / 2 tau).

#14

Alright.
Today I have solved the problem, and I would like to explain it in some detail here so that other people can use my enlightenment to save a lot of time.

First, about the correct definition of the deviatoric shear stress in physical units and LB units:
I am using Jonas’ PhD thesis as reference. In eq. (2.66) Jonas writes that


Pi^(1) = - 2 * c_s^2 * tau * rho * S

where c_s is the sound speed, tau is the relaxation parameter, rho is the density and S is the strain rate tensor (1.6),


S = 1 / 2 * (nabla u + (nabla u)^T)

On the other hand, we know that the physical deviatoric shear stress tensor is (1.5)


sigma = 2 * rho * nu * S

where nu is the kinematic viscosity. Moreover the total momentum tensor is (1.3)


Pi = rho * u * u + p * 1 - sigma

The tensor Pi^(1) can be calculated using the LBM (3.23),


Pi^(1) = sum_i f_i * c_i * c_i - rho * u * u - c_s^2 * 1

Now, when I compare everything, then I find the important result


sigma_ij = - (1 - 1 / (2 * tau)) * [sum_k f_k * c_ki * c_kj - rho * u_i * u_j - c_s^2 * delta_ij]

since


nu = c_s^2 (tau - 1 / 2)

This means that I need an additional factor of (1 - 1 / (2 * tau)) and a minus sign for the calculation of the shear stress tensor sigma.

Second, I have wasted a lot of time, because I was really stupid. I tried to implement this factor using the code


stress *= (1 - 1 / 2 / tau);

As you know, 1 and 2 are integers, which means that c++ thinks that 1 / 2 is also an integer, i. e. 0 in this case. 0 divided by tau is always zero, so this factor is always 1. Now I have changed the code to


stress *= (1 - 1 / tau / 2);

and the simulated results for the shear stress match the analytical results. So much about explaining a computer how to compute numbers.

Timm

#15

There is also something in this article:
Chen, Doolen: Lattice Boltzmann Method for Fluid Flows (Annual Review of Fluid Mechanics (30) 329-364)