TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Interfacial Forces

The interfacial forces are source terms added to the momentum equation that represent the mechanical interactions between the phases; they couple the equations of motion of the phase velocities. In this page, the expression given is generally the force acting on the gas phase; the opposite force acts on the liquid phase. Each model fills coefficient tables in the code; the convention is described per force below.

Drag force

The general expression of the drag force is:

\[\overrightarrow{F_{l\rightarrow g}^D}= - \frac{3}{4} C_D \frac{\alpha_g\rho_l}{d_b} \|\overrightarrow{u_g}-\overrightarrow{u_l}\| (\overrightarrow{u_g}-\overrightarrow{u_l})=-f^{D}\|\overrightarrow{u_g}-\overrightarrow{u_l}\|(\overrightarrow{u_g}-\overrightarrow{u_l}) \]

The base interfacial-drag source (Source_Frottement_interfacial_base) reads the parameters a_res, dv_min, exp_res, beta (defaults a_res=-1., dv_min=0.01, beta=1., exp_res=2) and the frottement_interfacial correlation. The drag operator fills the coefficient table so that \(\text{coeff}(k_1,k_2,0)=f^{D}\|\overrightarrow{u_{k_1}}-\overrightarrow{u_{k_2}}\|\) and \(\text{coeff}(k_1,k_2,1)=f^{D}\) (the derivative with respect to the relative velocity), with the table symmetric in \((k_1,k_2)\).

Availability of drag-force models in TrioCFD/CMFD:

Model Used Validated Test case
Constant Yes Yes TrioCFD/Tube analytique & Trust/Tube analytique
Composant Yes No
Ishii Zuber Deformable Yes No
Ishii Zuber Yes No
Tomiyama Yes Yes TrioCFD/CoolProp & TrioCFD/Gabillet
Weber Yes No
Wallis Yes No
Sonnenburg Yes No
Garnier Yes No
Rusche Yes No
Simmonet Yes No
Zenit Yes No

Constant (Frottement_interfacial_bulles_constant, params coeff_derive, rayon_bulle): if no rayon_bulle is specified, it takes d_bulles. \(f^{D}=\frac{3}{4}\frac{C_d\alpha_g\rho_l}{d_b}\).

Composant (Frottement_interfacial_bulles_composant): \(f^{D}_{ij}=\frac{3}{4}\frac{C_d\alpha_i\alpha_j\rho_m}{d_b}\), with \(\rho_m=\sum_k \alpha_k \rho_k\) and \(i\neq j\).

Ishii-Zuber, viscous regime [80] (Frottement_interfacial_Ishii_Zuber_Deformable):

\[f^{D}=\frac{1}{2}\alpha_g\rho_l\sqrt{\frac{(\rho_l-\rho_g)g}{\sigma}}\frac{1}{\sqrt{\max(1-\alpha_g,\ 0.001)}} \]

If \(\alpha_l<10^{-6}\), then \(f^D\times\alpha_l\times10^{-6}\).

Ishii-Zuber, viscous + particle regime [80] (Frottement_interfacial_Ishii_Zuber):

\[f^{D}=\frac{3}{4}\frac{\max\left(\frac{24}{Re_b}(1+0.1Re_b^{0.75}),\ \frac{2}{3}\sqrt{\frac{(\rho_l-\rho_g)g d_b^2}{\sigma}}\right)\beta\alpha_g\rho_l}{d_b},\quad Re_b=\frac{\rho_l d(u_g-u_l)}{\mu_l} \]

Tomiyama (contaminated) [171] (Frottement_interfacial_Tomiyama, params beta, constante_gravitation, contamination):

\[f_D=\frac{3}{4}\frac{\alpha_g\rho_l}{d}\begin{cases} \max(\min(16/Re_b(1+.15Re^{.687}),48/Re_b),8Eo/(3+12)) & \text{no contamination (0)}\\ \max(\min(24/Re_b(1+.15Re_b^{.687}),72/Re_b),8Eo/(3Eo+12)) & \text{slight contamination (1)}\\ \max(24/Re_b(1+.15Re_b^{.687}),8Eo/(3Eo+12)) & \text{high contamination (2)} \end{cases} \]

with \(Eo=\frac{g(\rho_l-\rho_v)d^2}{\sigma}\), \(Re_b=\frac{\rho_l d_b(u_g-u_l)}{\mu_l}\). Chosen because, as shown in [165], it yields similar results to other closures while allowing the contamination level to be adjusted.

Bubble critical diameter (incoming) [99] (Frottement_interfacial_Weber, param Weber_critique, default \(We_c=8\)): \(f^{D}=\frac{6\alpha_g}{\pi d_b^{*3}}\frac{24}{Re_b}(1+0.1Re_b^{0.75})\), with \(d_b^*=\frac{\sigma}{\rho_l(u_g-u_l)^2}We_c\). (Warning in source: not homogeneous.)

Wallis (annular flow) [177] : \(f^{D}=5\times10^{-3}\rho_g\frac{4\sqrt{\alpha_g}}{D_h}\left(1+300\frac{1-\sqrt{1-\alpha_g}}{2}\right)\).

Sonnenburg (drift flux): \(f^{D}=\rho_l\frac{\alpha_l\alpha_g}{D_h}\left(\frac{16}{9}(1-\alpha_g^*(1-\frac{9}{16}\sqrt{\frac{\rho_g}{\rho_l}}))\frac{1-\alpha_g^{*40}}{\tanh(32\alpha_g^*)}\right)^2\), with \(\alpha_g^*=\min(\max(\alpha_g,0.001),0.999)\).

Swarm corrections apply a multiplier to \(f^D\):

  • Garnier [62] : \(f^{D}_{new}=f^{D}\begin{cases}\alpha_l\times114.2 & \alpha_l<0.5\\(1-\alpha_g^{1/3})^{-2}&\text{else}\end{cases}\) (validated for \(\alpha_g<0.35\), \(D_{sm}<5.5\) mm).
  • Rusche [150] : \(f^{D}_{new}=f^{D}(\exp(3.64\alpha_g)+\alpha_g^{0.864})\) (validated for \(\alpha_g<0.5\)).
  • Simonnet [160] : \(f^{D}_{new}=f^{D}\alpha_l(\alpha_l^{25}+(4.8\frac{\alpha_g}{\alpha_l})^{25})^{-2/25}\) (validated for \(\alpha_g<0.3\), \(D_{sm}<10\) mm).
  • Zenit [190] : \(f^{D}_{new}=f^{D}\frac{(1+3\alpha_g)^2}{\alpha_l^2}\) (validated for \(\alpha_g<0.18\)).

Lift force

The general expression of the lift force is:

\[\overrightarrow{F_{l\rightarrow v}^L}=-C_L\rho_l\alpha_g(\overrightarrow{u_g}-\overrightarrow{u_l})\wedge\overrightarrow{u_l}=-f^L(\overrightarrow{u_g}-\overrightarrow{u_l})\wedge\overrightarrow{u_l} \]

The base source (Source_Portance_interfaciale_base) reads beta, g and the Portance_interfaciale correlation, requires a liquid phase, and creates a vorticite field. The lift operator fills out.Cl so that \(Cl(k_1,k_2)=f^L\) (symmetric).

Model Used Validated Test case
Constant Yes Yes TrioCFD/Tube analytique, TrioCFD/CoolProp, Trust/Tube analytique
Sugrue Yes Yes TrioCFD/CoolProp, TrioCFD/Gabillet
Tomiyama Yes No

Constant (Portance_interfaciale_Constante): \(f^{L}=C_L\rho_l\alpha_g\max(\min(\frac{\alpha_l-0.05}{0.25},1),0)\); the void-fraction correction damps the lift at high void fractions.

Sugrue: \(f^{L}=\rho_l\alpha_g\max(1.0155-0.0154\exp(8.0506\alpha_g),0)\times\min(5.0404-5.0781Wo^{0.0108},0.03)\), with the wobbling number \(Wo=\min(\frac{k_l Eo}{\max((u_g-u_l)^2,10^{-8})},6)\) and \(Eo=\frac{g(\rho_l-\rho_g)d_b^2}{\sigma}\).

Tomiyama (sign reversal) [172] : \(f^{L}=\rho_l\alpha_g\begin{cases}\min(0.288\tanh(0.121Re_b),f(Eo)) & Eo<4\\ f(Eo) & 4\leq Eo\leq10.7\end{cases}\), with \(f(Eo)=0.00105Eo^3-0.0159Eo^2-0.0204Eo+0.474\).

Added-mass force

The general expression is \(\overrightarrow{F_{l\rightarrow v}^{AM}}=C_{AM}\alpha_g\rho_l\frac{D(u_g-u_l)}{Dt}=f^{AM}\frac{D(u_g-u_l)}{Dt}\). The added-mass operator updates the a_r table (with \(f^{AM}\) added/subtracted on the diagonal/off-diagonal of the \((k_1,k_2)\) velocity-derivative coefficients).

Model Used Validated Test case
Constant Yes Yes TrioCFD/Tube analytique, Trust/Tube analytique
Wijngaarden Yes No
Zuber Yes No

Constant (Masse_ajoutee_Coef_Constant, defaults beta=0.5, limiter_liquid_=0.5): \(f^{AM}=\min(\beta\rho_l\alpha_g,\ \rho_l\alpha_l\times\texttt{limiter\_liquid\_})\). Direct void-fraction influence is limited at \(\alpha_{gmax}=\frac{\texttt{limiter\_liquid\_}}{\texttt{limiter\_liquid\_}+\beta}\) (default 0.5). The injected mass flux is \(\dot{m}_{inj}=\min(\beta\rho_l,\texttt{limiter\_liquid\_}\frac{\alpha_l}{\alpha_g})\dot{m}\times\{\texttt{inj\_ajoutee\_gaz\_}\text{ (gas)},\ \texttt{inj\_ajoutee\_liquide\_}\text{ (liquid)}\}\); no limiter part if \(\alpha_g<0.0001\).

Wijngaarden (two-bubble interaction) [10] : \(f^{AM}=\min(\beta(1+2.78\alpha_g)\rho_l\alpha_g,\ \rho_l\alpha_l\times\texttt{limiter\_liquid\_})\) (default \(\alpha_{gmax}\approx0.34\)). Warning: the corrected value in [10] is 3.32 instead of 2.78.

Zuber (swarm of compliant bubbles) [193] : \(f^{AM}=\min(\beta\frac{1+2\alpha_g}{\max(1-\alpha_g,0.001)}\rho_l\alpha_g,\ \rho_l\alpha_l\texttt{limiter\_liquid\_})\) (default \(\alpha_{gmax}\approx0.303\)).

Dispersion force

The general expression of the turbulent dispersion force is \(\overrightarrow{F_{l\rightarrow v}^T}=-f^T\nabla\alpha_g\). The base source (Source_Dispersion_bulles_base) reads beta and requires the Dispersion_bulles correlation, filling out.Ctd so that \(Ctd(k_1,k_2)=f^{T}\) (coefficient in front of \(\nabla\alpha_{k_1}\)), symmetric.

Model Used Validated Test case
Constant bubble Yes Yes TrioCFD/Tube analytique, Trust/Tube analytique
Constant turbulent Yes Yes TrioCFD/Tube analytique
Lopez Yes No
Burns Yes Yes TrioCFD/CoolProp, TrioCFD/Gabillet

Constant bubble [115] : \(f^{T}=D_{td}\rho_l(u_g-u_l)^2\).

Constant turbulent [30] (default C_td_=0.1): \(f^{T}=C_{td}\rho_l k_l\).

Lopez de Bertodano (Stokes regime) [31] : \(f^{T}=2\rho_l k_l\frac{1}{(1+St)St}\), with \(St=\frac{\tau^F}{\tau^t}\), \(\tau^t=\frac{\nu_t}{k_l}\), \(\tau^F=\frac{\frac{4}{3}\rho_g d}{C_d\rho_l(u_g-u_l)}\). (In the literature: \(f^{T}=\rho_lk_l\frac{C_{\mu}^{1/4}}{(1+St)St}\) with \(\tau^t=C_{\mu}^{3/4}\frac{k_l}{\varepsilon_l}\).)

Burns (Favre-averaged drag) [12] :

\[F^{T}=\frac{f^D\|\vec{u_g}-\vec{u_l}\|(\nu_t+\nu_{BIA})}{Pr_t}\left(\frac{1}{\alpha_g}\nabla\alpha_g-\frac{1}{\alpha_l}\nabla\alpha_l\right) \]

with \(\nu_t=C_\mu\frac{k^2}{\varepsilon}\), \(\nu_{BIA}=\texttt{coefBIA\_}\frac{k_{WIT}}{\omega_{WIT}}\), and \(\omega_{WIT}=2\mu_l\frac{Re_b}{C_{\lambda}^2 d^2}\max(\min(\frac{16}{Re_b}(1+0.15Re_b^{0.687}),\frac{48}{Re_b}),8\frac{Eo}{3(Eo+4)})\). (In the literature: \(f^{T}=\frac{f^D\|\vec{u_g}-\vec{u_l}\|\nu_t}{Pr_t}(\frac{1}{\alpha_g}+\frac{1}{1-\alpha_g})\).)

Wall force

Antal (wall lubrication) [7] (Correction_Antal_PolyMAC_P0, defaults Cw1_=-0.1, Cw2_=0.147): \(F^{WL}=C_{WL}\alpha_g\rho_l\frac{(\vec{u_g}-\vec{u_l})^2}{d_b}\vec{n}\), with \(C_{WL}=\max(-C_{W1}+C_{W2}\frac{d_b}{2y},0)\). Developed for fully developed laminar bubbly two-phase flows.

Lubchenko (wall-force damping, PolyMAC only) [112] : a mix of the Lubchenko lift damping and BIF near-wall damping. The lift damping close to the wall is

\[C_L\rightarrow\begin{cases}0 & y/d_b<\tfrac{1}{2}\texttt{portee\_lift\_}\\ C_L(3(\tfrac{2y}{d_b}-1)^2-2(\tfrac{2y}{d_b}-1)^3) & \tfrac{1}{2}\texttt{portee\_lift\_}\leq y/d_b<\texttt{portee\_lift\_}\\ C_L & y/d_b\geq\texttt{portee\_lift\_}\end{cases} \]

The dispersion-balanced part replaces \(\underline{\nabla}\alpha_g\) by \(\alpha_g\frac{1}{y}\frac{d_b\,\texttt{portee\_disp\_}-2y}{d_b\,\texttt{portee\_disp\_}-y}\overrightarrow{n}\) in the range \(y/d<\texttt{portee\_disp\_}\), and a last part cancels the BIF contribution \(0.5d\) from the wall.

Tchen force

Warning
This force is a good example of implementation, but its use is discouraged.

The general expression of the Tchen force [17] is \(\overrightarrow{F_{l\rightarrow g}}=\alpha_g\rho_l\frac{\partial\overrightarrow{u_l}}{\partial t}\), discretized on the right-hand side (secmem) as \(\alpha_g^n\rho_l^n\frac{u_l^n-u^{n-1}}{\Delta t}\). The local void fraction \(\alpha_{\text{loc}}\) at a face is computed from the interlaced volumes (domaine.volumes_entrelaces()) of the cells on both sides of the face; the density is computed the same way. The full term is added in secmem (with sign \(+1\) in gas, \(-1\) in liquid) and the velocity matrix mat is updated without the velocity increment.

Note
The source page included a figure ("ADD FIGURE FROM PDF") illustrating the interlaced-volume management at faces; no image was available, so it is omitted here.