Hello everyone,
I am quite new to Basilisk and CFD in general, so apologies in advance if I am missing something obvious.
I am trying to implement a 2D Herschel–Bulkley pressure-driven channel flow using navier-stokes/centered.h. My final goal is to validate the numerical solution against the analytical Herschel–Bulkley Poiseuille solution before moving to more complex geometries.
I would appreciate any advice on whether my implementation is conceptually correct and what might be causing the issue described below.
/**
Herschel-Bulkley channel flow
Periodic streamwise direction
No-slip walls
*/
#include "navier-stokes/centered.h"
//==================================================
// DOMAIN
//==================================================
double Lx = 0.256; // domain length
double H = 0.05; // channel height
//==================================================
// FLUID PROPERTIES
//==================================================
double rho_val = 930.; // kg/m3
double tauy = 443.5; // Pa
double K = 78.37; // Pa.s^n
double nHB = 0.462;
double mu_max = 5000.; // regularisation limit
// imposed pressure gradient
double dpdx = 35000.; // Pa/m
//==================================================
face vector muv[];
face vector av[];
scalar un[];
//==================================================
int main()
{
L0 = Lx;
origin (-L0/2., -L0/2.);
N = 256;
periodic (right);
stokes = true;
DT = 1e-4;
TOLERANCE = 1e-3;
run();
}
//==================================================
// INITIALIZATION
//==================================================
event init (t = 0)
{
// Channel geometry
mask (fabs(y) > H/2. ? top : none);
// Constant density
const scalar rhoc[] = rho_val;
rho = rhoc;
const face vector alphac[] = {1./rho_val,1./rho_val};
alpha = alphac;
// Variable viscosity
mu = muv;
// Body force
a = av;
foreach() {
u.x[] = 0.;
u.y[] = 0.;
un[] = 0.;
}
}
//==================================================
// PRESSURE GRADIENT
//==================================================
event acceleration (i++)
{
foreach_face(x)
av.x[] = dpdx/rho_val;
foreach_face(y)
av.y[] = 0.;
}
//==================================================
// HERSCHEL-BULKLEY VISCOSITY
//==================================================
event properties (i++)
{
foreach_face() {
double D11 =
(u.x[] - u.x[-1])/Delta;
double D22 =
((u.y[0,1] - u.y[0,-1]) +
(u.y[-1,1] - u.y[-1,-1]))
/(4.*Delta);
double D12 =
0.5*
(
((u.x[0,1] - u.x[0,-1]) +
(u.x[-1,1] - u.x[-1,-1]))
/(4.*Delta)
+
(u.y[] - u.y[-1])/Delta
);
double D2 =
sqrt(sq(D11) + sq(D22) + 2.*sq(D12));
double mu_eff;
if (D2 > 1e-12) {
mu_eff =
tauy/(sqrt(2.)*D2)
+
K*pow(D2/sqrt(2.), nHB - 1.);
mu_eff = min(mu_eff, mu_max);
}
else
mu_eff = mu_max;
muv.x[] = fm.x[]*mu_eff;
}
}
//==================================================
// LOGGING
//==================================================
event logfile (i += 10)
{
double umax = 0.;
foreach()
umax = max(umax, fabs(u.x[]));
double du = change (u.x, un);
fprintf (stderr,
"i=%d t=%g dt=%g umax=%g du=%g mgp=%d mgu=%d\n",
i, t, dt, umax, du,
mgp.i, mgu.i);
}
Expected behaviour:
A steady Herschel–Bulkley Poiseuille flow should develop in the channel.
Observed behaviour:
The simulation remains stable, but the maximum velocity keeps increasing. For example:t t= 3.54 umax = 47.72 m/s
t = 3.56 umax = 47.89 m/s
t = 3.58 umax = 48.03 m/s
The solution does not seem to approach a steady state.
Since I am new to Basilisk, I would like to ask:
- Is my implementation of the Herschel–Bulkley viscosity correct?
- Is the strain-rate invariant computed correctly?
- Is a body force a.x[] = dpdx/rho the correct way to impose a pressure gradient in a periodic channel?
- Is the channel geometry created correctly using the mask() function?
- Is there anything obviously wrong in this implementation that could explain the continuously increasing velocity?
Since I am still learning Basilisk, any comments on the overall implementation are very welcome.
Thank you very much for your time and help.
Best regards,
Valeria
