shear_flow.c

Shear flow

Shear flow

In this section, I present the code used to simulate the shear flow. For a step-by-step explanation, feel free to watch my video tutorial: Shear flow.

We start by including the necessary libraries:

#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
#include "vtknew_cell.h"

Here, we work in non-dimension numbers. Thus, we define all the parameters:

Density ratio \rho_r = 1;
Viscosity ratio \mu_r = 1;
Reynolds number Re = 1;
Capillary number Ca = 0.5;
Dimensionless radius R_0 = 1;
Dimensionless length L = 4;
Dimensionless velocity U = 1.

double RHOR = 1.;
double MUR = 1.;
double Re = 1.;
double Ca = 0.5;
double L = 4.;
double R0 = 1.;
double U = 1.;

To close the problem, appropriate boundary conditions are imposed on the velocity field. A shear flow is generated by prescribing opposite tangential velocities at the horizontal boundaries: the top wall moves with a constant velocity U, while the bottom wall moves with velocity −U. This configuration creates a symmetric shear across the domain.

On the vertical boundaries, no normal velocity gradient is imposed, corresponding to homogeneous Neumann conditions for the normal component of the velocity. These conditions ensure that the flow can freely adjust at the inlet and outlet without artificially constraining the normal momentum.

The boundary conditions are implemented as follows:

u.t[top] = dirichlet(U);
u.t[bottom] = dirichlet(-U);
u.n[left] = neumann(0.);
u.n[right] = neumann(0.);

In the main() function, we initialize the domain size (a square domain of side SIZE_BOX), the uniform initial grid resolution 2^8 = 256, and the physical parameters: densities, viscosities, and surface tension.
We then call run() to launch the simulation.

int main() {
    size(L);
    origin(-L/2., -L/2.);
    init_grid(256);

    rho1 = 1.;
    rho2 = rho1 / RHOR;
    mu1 = 1. / Re;
    mu2 = mu1 / MUR;
    f.sigma = 1. / Ca;

    system("mkdir -p shear");
    run();
}

Initial condition

The initial configuration consists of a circular interface of radius R_0, centered in the computational domain. The interface is defined implicitly by the level-set function

x^2 + y^2 = R_0^2.

In the Volume-of-Fluid framework, the interface is initialized through the scalar function

F(x,y,t=0) = x^2 + y^2 - R_0^2,

event init(t=0) {
    fraction(f, sq(x) + sq(y) - sq(R0));

In addition to the interface initialization, an initial velocity field is prescribed throughout the computational domain. The flow is initialized as a simple shear profile, where the horizontal velocity varies linearly with the vertical coordinate, while the vertical velocity component is initially zero.

This corresponds to the velocity field

u_x(x,y,0) = y, \qquad u_y(x,y,0) = 0.

The initialization is implemented as follows:

    foreach() {
        u.x[] = y;
        u.y[] = 0.;
    }
}

To efficiently capture the interface dynamics and velocity gradients while limiting the overall computational cost, an adaptive mesh refinement strategy is employed. The grid adaptation is performed using the wavelet-based refinement algorithm provided by Basilisk.

The refinement criterion is applied to both the volume fraction field f, which tracks the interface, and the velocity field \mathbf{u}. Cells are refined or coarsened according to the estimated discretization error compared to user- defined thresholds.

The maximum level of refinement is set to LEVEL, while the minimum level is restricted to LEVEL-4 in order to maintain a sufficiently coarse background grid.

The adaptation procedure is implemented as follows:

int LEVEL = 8;
event adapt(i++) {
    double femax = 0.001;
    double uxemax = 0.001;
    double uyemax = 0.001;

    adapt_wavelet({f, u}, (double[]){femax, uxemax, uyemax},
                  maxlevel = LEVEL, minlevel = LEVEL - 4);
}

Finally, this event handles data output using output_vtk(),
which saves the volume fraction f, the pressure p, and the velocity field \mathbf{u} in .vtk format, every 0.01 time units. We also save the shape of the interface (f=0.5) using output_facets(f, ptr).

int count = 0;
event vtk(t=0.; t+=0.01; t<=1.) {
    char filename[100];
    sprintf(filename, "vtk/snap-%d.vtk", count);
    FILE * fileptr = fopen(filename,"w");

    scalar * list = {f, u.x, u.y, p};
    output_vtk(list, fileptr);
    fclose(fileptr);

    char name[80];
    sprintf(name, "interface/snap-%d.dat",count);
    FILE *ptr = fopen(name,"w");
    output_facets(f, ptr);
    fclose(ptr);
    
    if (pid() == 0) {
        fprintf(stderr, "t: %-10g\t dt: %-10g\n", t, dt);
    }

    count++;
}