Transient planar Poiseuille flow for a viscoelastic Oldroyd-B or FENE-P fluid

The onset of a Poiseuille flow from a fluid at rest after the sudden application of a pressure gradient is characterized, for a Newtonian fluid, by a monotonic exponential increase of the axial velocity.

In the case of a viscoelastic fluid, an elastic oscillation is superposed on the exponential increase. The problem has an analytical solution due to Waters & King (1970),

U(Y,T)=1.5(1Y2)48k=1sin((1+Y)n/2)n3eαnT/2G(T) with n=(2k1)π, αn=1+βEn2/4 and G(T)=sinh(βnT/2)+γnβncosh(βnT/2) with βn=αn2En2andγn=12β4En2

E is the elastic number given by E=λμo/(ρh2). The time is made dimensionless with λ, T=t/λ, and the velocity with the average steady velocity, u=ΔpΔxh23μo

The first modes of the analytical solution imply complex values of βn, which is why we include the complex library.

#include <complex.h>

#define DT_MAX 0.001
// Total viscoelastic viscosity 
#define MU0 1.          
// Ratio of the solvent viscosity to the total viscosity
#define BETA (1/9.)
// Polymer viscosity
#define MUP ((1. - BETA)*MU0)
// Solvent viscosity 
#define MUS (BETA*MU0)
// Relaxation viscoelastic time
#define LAM 1.
// Average velocity steady flow
#define UAVG (1./(3.*MU0))

Note that we are using as dimensioning scales the width of the gap, the density of the fluid and the gradient of pressure (rather than the average velocity for the steady flow).

#include "grid/multigrid.h"
#include "navier-stokes/centered.h"
#include "log-conform.h"
#if FENE_P
# include "fene-p.h"
# define L2 1.

int lev;

int main()
  periodic (right);
  p[left] = dirichlet(1.);
  p[right] = dirichlet(0.);

  DT = DT_MAX;
  const scalar lam[] = LAM;
  λ = lam;
  const scalar mupc[] = MUP;
  mup = mupc;
  const face vector mus[] = {MUS,MUS};
  μ = mus;
#if FENE_P
  lev = 5;
  init_grid (1 << lev);
  L2 = 10.;   run();
  L2 = 50.;   run();
  L2 = 1000.; run();
#else // Oldroyd-B
  for (lev = 4; lev <= 6; lev++) {
    init_grid (1 << lev);

u.t[top] = dirichlet(0);
u.t[bottom] = neumann(0);

event init (i = 0) 
  scalar s = tau_p.y.y;
  s[top] = dirichlet(0.);
  s = tau_p.x.y;
  s[bottom] = dirichlet(0.);
    u.x[] = 0.;
  boundary ((scalar *){u});

double analytical (double Y, double T, int KF)
  double E = LAM*MU0;
  double U = 0.;
  for (int k = 1; k <= KF; k++) {
    double n = (2*k - 1)*M_PI;
    double α = 1. + 0.25*BETA*E*sq(n);
    complex double β = csqrt(sq(α) - E*sq(n));
    double γ = 1. - 0.25*(2. - BETA)*E*sq(n);
    double G = creal(ccosh(0.5*β*T) + γ/β*csinh(0.5*β*T));
    U += 1./n/sq(n)*sin(0.5*(1. + Y)*n)*exp(-0.5*α*T)*G;
  return 1.5*(1. - sq(Y)) - 48.*U;

event uaxis_evolution (t += 0.2; t <= 10.) {

The velocity is made dimensionless with the average velocity. when the flow is fully developed, uaxis/u = 3/2.

  fprintf (stderr, "%g %g %g %d %g\n", t/LAM, 
	   interpolate(u.x, 0.5, 0.)/UAVG, analytical (0, t/LAM, 8), lev, L2);

#if FENE_P
event uprofile (t = end) {
  int np = 20;
  double hh = 0.999999999/np;
  for (int j = 0; j <= np; j++)
    printf ("%g %g %d %g\n", hh*j, interpolate(u.x, 0.5, hh*j)/UAVG, lev, L2);
#endif // FENE_P

For an Oldroyd-B fluid the solution converges to the theoretical solution with grid resolution.

Temporal evolution of the axial velocity

Temporal evolution of the axial velocity

Temporal evolution of the error for different grids

Temporal evolution of the error for different grids

For a FENE-P fluid, the solution converges toward the Oldroyd-B solution when L2.

Temporal evolution of the axial velocity

Temporal evolution of the axial velocity

Velocity profiles at instant t=10

Velocity profiles at instant t=10



E. Gros. Comparing in-house numerical simulations to the analytical solution of the poiseuille flow for the Oldroyd-B model. Master’s thesis, Aachen University, Germany, 2013. [ .pdf ]


ND Waters and MJ King. Unsteady flow of an elastico-viscous liquid. Rheologica Acta, 9(3):345-355, 1970.