CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Thermophoresis source term causes instabiltiy in population balance model

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 28, 2022, 20:22
Default Thermophoresis source term causes instabiltiy in population balance model
  #1
New Member
 
Join Date: Jan 2015
Posts: 22
Rep Power: 11
adib_mohammad is on a distinguished road
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, D the diffusion coefficient.

T 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.
adib_mohammad is offline   Reply With Quote

Old   March 29, 2022, 16:06
Default
  #2
New Member
 
Join Date: Jan 2015
Posts: 22
Rep Power: 11
adib_mohammad is on a distinguished road
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.
adib_mohammad is offline   Reply With Quote

Old   March 29, 2022, 18:37
Default
  #3
New Member
 
Join Date: Mar 2022
Posts: 4
Rep Power: 4
jines0708 is on a distinguished road
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.
jines0708 is offline   Reply With Quote

Old   March 29, 2022, 18:46
Default
  #4
New Member
 
Join Date: Jan 2015
Posts: 22
Rep Power: 11
adib_mohammad is on a distinguished road
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.
adib_mohammad is offline   Reply With Quote

Old   March 31, 2022, 03:41
Default
  #5
Senior Member
 
Josh Williams
Join Date: Feb 2021
Location: Scotland
Posts: 113
Rep Power: 5
joshwilliams is on a distinguished road
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.
joshwilliams is offline   Reply With Quote

Old   March 31, 2022, 03:42
Default
  #6
Senior Member
 
Josh Williams
Join Date: Feb 2021
Location: Scotland
Posts: 113
Rep Power: 5
joshwilliams is on a distinguished road
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.
joshwilliams is offline   Reply With Quote

Old   March 31, 2022, 19:16
Default
  #7
New Member
 
Join Date: Mar 2022
Posts: 4
Rep Power: 4
jines0708 is on a distinguished road
Hi Josh,
thanks for the reply. I don't see them implementing thermophoretic force. Or am I missing something?
jines0708 is offline   Reply With Quote

Old   April 1, 2022, 10:40
Default
  #8
New Member
 
Join Date: Jan 2015
Posts: 22
Rep Power: 11
adib_mohammad is on a distinguished road
Quote:
Originally Posted by joshwilliams View Post
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.
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?
adib_mohammad is offline   Reply With Quote

Old   May 13, 2022, 08:23
Default
  #9
New Member
 
Join Date: Jan 2015
Posts: 22
Rep Power: 11
adib_mohammad is on a distinguished road
Any suggestions?
adib_mohammad is offline   Reply With Quote

Old   May 21, 2022, 14:35
Default
  #10
New Member
 
Berin 各ta
Join Date: Dec 2014
Location: Tarragona
Posts: 8
Rep Power: 11
Berin 各ta is on a distinguished road
Quote:
Originally Posted by adib_mohammad View Post
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, D the diffusion coefficient.

T 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.
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
Berin 各ta is offline   Reply With Quote

Old   May 24, 2022, 09:05
Default
  #11
New Member
 
Join Date: Jan 2015
Posts: 22
Rep Power: 11
adib_mohammad is on a distinguished road
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).

Last edited by adib_mohammad; May 24, 2022 at 15:02.
adib_mohammad is offline   Reply With Quote

Old   May 24, 2022, 09:16
Default
  #12
New Member
 
Berin 各ta
Join Date: Dec 2014
Location: Tarragona
Posts: 8
Rep Power: 11
Berin 各ta is on a distinguished road
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".
Berin 各ta is offline   Reply With Quote

Old   May 24, 2022, 11:36
Default
  #13
New Member
 
Join Date: Jan 2015
Posts: 22
Rep Power: 11
adib_mohammad is on a distinguished road
Quote:
Originally Posted by Berin Šeta View Post
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.
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 View Post
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".
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")

Last edited by adib_mohammad; May 24, 2022 at 15:05.
adib_mohammad is offline   Reply With Quote

Reply

Tags
source terms

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[swak4Foam] swak4foam for OpenFOAM 4.0 mnikku OpenFOAM Community Contributions 80 May 17, 2022 08:06
[swak4Foam] funkyDoCalc with OF2.3 massflow NiFl OpenFOAM Community Contributions 14 November 25, 2020 03:30
[swak4Foam] Installation Problem with OF 6 version Aurel OpenFOAM Community Contributions 14 November 18, 2020 16:18
centOS 5.6 : paraFoam not working yossi OpenFOAM Installation 2 October 9, 2013 01:41
friction forces icoFoam ofslcm OpenFOAM 3 April 7, 2012 10:57


All times are GMT -4. The time now is 03:30.