CFD Online Discussion Forums - Thermophoresis source term causes instabiltiy in population balance model

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - Thermophoresis source term causes instabiltiy in population balance model (https://www.cfd-online.com/Forums/openfoam-solving/241940-thermophoresis-source-term-causes-instabiltiy-population-balance-model.html)

Thermophoresis source term causes instabiltiy in population balance model

I have created a new solver based on buoyantSimpleFoam to accomodate an additional transport equation that tracks number concentration of particles. The transport equation looks like this:

$\frac{\partial N_{agg}}{\partial t}+\frac{\partial}{\partial x_i} (u_i. N_{agg})= D.\frac{ \partial^2 N_{agg}}{\partial x_i^2 }+ K_{th}.\frac{ \partial }{\partial x_i\ }(\frac{ N_{agg}}{T} \frac{ \partial T}{\partial x_i})$

where $N_{agg}$ is the number concentration of particles, $K_{th}$ the thermophoresis constant,

the diffusion coefficient.

is the gas temperature, and I get it from thermo object.

Code:

volScalarField& T = thermo.T();

I am simulating particle motion in a laminar flow ( $U_{avg} = 0.36 m/s$ ) in a pipe. The inlet flow has a uniform temperature of 1500 K ( $T_{in}=1500 K$ ) and walls the constant temperature of 500 K ( $T_{w}=500 K$ ).
For low diffusion coefficient ( $10^{-11}$ ) that weakens the diffusion term ( $D.\frac{ \partial^2 N_{agg}}{\partial x_i^2 }$ ), the simulation results large negative and positve values for $N_{agg}$ , but the results are meaningful for high diffusion coefficient ( $10^{-5}$ ).

My NAEqn.h file (solves $N_{agg}$ ) looks like this.

Code:



{



    fvScalarMatrix NAEqn

    (

        fvm::ddt(rho, N_agg)

        + fvm::div(phi, N_agg)

        - fvm::laplacian(diff_prefactor*diff_p*rho, N_agg)

    ==

        fvc::laplacian( K_th * (mu * N_agg/T), T)



    );



    NAEqn.relax();

    fvOptions.constrain(NAEqn);

    NAEqn.solve();

    fvOptions.correct(N_agg);

}

Based on my review of literature, the simulation with thermophoresis and low particle diffusion should be possible to perform but I don't know how to achieve this in OpenFOAM. I know that part of the problem is coupling where the source term of this equation depends on the second derivative of temperature, but I don't know how to solve this issue.

I should also add that velocity and tempreature values are physical, and the number concentration equation divegres. The large values of number concentration start to appear near the wall, and they spread to the rest of the computational domain.

I am very new to this and was looking to add thermophoresis into denseparticleFoam with no luck. Do you think even if it's possible? I am using openfoam9.

I am not familiar with the solver you are using, but I think the principles and the procedures are the same. The transport equation of the variable affected by thermophoresis should include the $K_{th}.\frac{ \partial }{\partial x_i\ }(\frac{ \phi}{T} \frac{ \partial T}{\partial x_i})$ with $\phi$ being the variable of interest.

Hi Abib.

You may consider using the OpenQBMM library to solve your PBE. Indeed, it already has a solver called buoyantPbePimpleFoam (tutorial here), which may be easier for you to use than developing your own. Alternatively, maybe looking at their source code could be helpful for you.

Another option is the aerosolved library, which I have not used personally but I have seen several articles using it in literature and at conferences.

Hi Josh,
thanks for the reply. I don't see them implementing thermophoretic force. Or am I missing something?

Quote:

Originally Posted by joshwilliams (Post 825180)

Thanks for introducing the package. Since the code I am diagnosing is part of my thesis project, I cannot shift to another solver. Do you have any suggestions regarding how to solve the divergence issue my this case?

Quote:

Originally Posted by adib_mohammad (Post 824985)

the diffusion coefficient.

is the gas temperature, and I get it from thermo object.

Code:

volScalarField& T = thermo.T();

Code:



{



    fvScalarMatrix NAEqn

    (

        fvm::ddt(rho, N_agg)

        + fvm::div(phi, N_agg)

        - fvm::laplacian(diff_prefactor*diff_p*rho, N_agg)

    ==

        fvc::laplacian( K_th * (mu * N_agg/T), T)



    );



    NAEqn.relax();

    fvOptions.constrain(NAEqn);

    NAEqn.solve();

    fvOptions.correct(N_agg);

}

Hi,

Without completely seeing your case, I cant be sure if the error is what I assume it is.

Since you stated that large values of concentration appear near walls, that is expected because the temperature gradients will be highest there and you will have thermodiffusion/themophoresis stongest at those zones.

Unrealistically high values lead me to conclusion that your boundary conditions for the concentration is wrong. Could you tell me what kind of boundary conditions you used for concentration at walls with fixed temperature?

If you use zeroGradient, that might be source of divergence, because your mass flux will be non-zero (you will have "generation of particles" or "leakage").

You should make fixedGradient boundary condition, related with temperature gradient that you have in the walls, basically:

mass_flux = 0, this leads to grad(C) = - (K_th/D)*grad(T)

The boundary condition above could be implemented simply using groovyBC.

The similar work regarding multicomponent thermodiffusion is implemented here: https://link.springer.com/article/10.../i2019-11818-7 - for binary system (equivalent to your problem)

The implementation of the specie transport in that work was implement in this way (similar to yours obsviously):

Code:

{

    fvScalarMatrix w1Eqn

    (

        fvm::ddt(w1)

      + fvm::div(phi, w1)

      - fvm::laplacian(D1, w1)   

      - fvc::laplacian(DT*w1*w2, T)

    );



    w1Eqn.relax();



    fvOptions.constrain(w1Eqn);



    w1Eqn.solve();



    fvOptions.correct(w1);



}

and it was implement in the boussinesqBouyantPimpleFoam.

If the above mentioned is reason for the "divergence" of particles number, that could explain why your case for the high diffusion coefficient it works. Since your diffusion coefficient is so high, the Soret effect is small (so separation is very limited) and hence the strong temperature gradient in the zones close to the walls does not separate particles in large scale ((K_th/D) - this value becomes very small).

Hope this helps!

Let me know if this was problem,

BR
Berin

Thanks @Berin

I used a fixed value $N_{agg} = 0$ . I assumed that wall is a perfect particle sink adsorbing any particle in the cell adjacent to the wall. I think the boundary condition (bc) causes the problem.

Fixed gradient bc you mentioned ( $\frac{\partial N_{agg}}{\partial y} = - \frac{K_{th}}{D} \frac{ \partial T}{\partial y}$ ) makes sense, but it turns into a zero gradient bc when temperature gradient is zero ( $\frac{ \partial T}{\partial y} = 0 \rightarrow \frac{\partial N_{agg}}{\partial y}=0$ ). Very low temperature gradients (near zero) occur in the regions near the end of pipe, that fluid reaches wall tempeature.

Moreover, this bc eliminates wall deposition due to diffusion, and the solver will not produce physical results for smaller particles (with moderate diffusion coefficients in the order of 1e-5 or 1e-6).

Quote:

Fixed gradient bc you mentioned ( $\frac{\partial N_{agg}}{\partial y} = - \frac{K_{th}}{D} \frac{ \partial T}{\partial y}$ ) makes sense, but it turns into a zero gradient bc when temperature gradient is zero ( $\frac{ \partial T}{\partial y} = 0 \rightarrow \frac{\partial N_{agg}}{\partial y}=0$ ). Very low temperature gradients (near zero) occur in the regions near the end of pipe, that fluid reaches wall tempeature.

But isnt thermophoresis/thermodiffusion all about this? If you dont have temperature gradient, you should not have temperature driven attraction/separation of particles. When your temperature gradient is zero (like in the end of the pipe), thermophoresis should not exist.

Quote:

Moreover, this bc eliminates wall position due to diffusion, and the solver will not produce physical results for smaller particles (with moderate diffusion coefficients in the order of 1e-5 or 1e-6).

I dont quite understand what you mean by elimination of wall position due to diffusion..

However, I dont see why this boundary condition should not work with different D values. It includes D coefficient in it, and it simply states that Np number is conserved at "steady-state".

Quote:

Originally Posted by Berin Šeta (Post 828561)

You are right. I am simulating a pipe flow where hot gas enters a pipe with cold walls. There is a sharp temperature gradient at the beginning, but as flow goes forward along the pipe, the flow temperature gets closer to the wall tempeature due to the heat transfer. Near the end of pipe, the flow has a nearly uniform temperature inidicating zero thermophoresis force. At this point, the bc is equavalent to zero gradient.

Quote:

Originally Posted by Berin Šeta (Post 828561)

The main difference of my work is that the mass flux of particles is not zero on the wall meaning that the total number of particles is not conserved in the pipe. Some particles deposit on the walls or simply being removed from the pipe flow. (I noticed a typo in my previous post. I meant "deposition" not "position")