Very small timestep for two-phase flow

[Brandon Aranda]

Hi all,

I am trying to simulate an experiment of air-water two-phase flow in a 0.1953 m diameter, 8 m length pipe. The fluid properties are taken to be those of the inlet pressure of 250 kPa. The issue I am having is that the timestep becomes very small, on the order of 10^-6 s and decreasing up to 10^-9 s. At these conditions, I don't expect bubbles to have a velocity larger than 3 m/s, if even that. I'm using multigrid and with the refinement level I'm currently using, the total mesh size is ~84 million cells, meaning each cell should be ~1.5 mm. With the CFL condition of 0.5 and the given cell size and velocity (which should be a conservative value), I would expect the timestep should be ~10^-4. Does anyone have any ideas why the timestep is so small, and what I could do to fix it? I would imagine it might be from my inlet setup, in which I have 64 tightly packed injection holes of 0.014 m diameter injecting air at ~2 m/s. I've chosen this as the mass flow rate is quite high and using fewer injection holes would require an even larger inlet velocity which I would imagine would exacerbate the problem. The pressure also fluctuates dramatically through the channel, making me believe it may be related to the boundary conditions. I'd appreciate any suggestions on the best set of boundary conditions for this setup and any other suggestions. I've provided a minimal version of the code here.

Thanks,

Brandon Aranda

[Wojciech Aniszewski]

Hi Brandon,
have you dumped the velocity anywhere to see if it really maxes out at the level you "expect"?

And by the way why do you need the system math.h?

regards
v

> --
> You received this message because you are subscribed to the Google Groups "basilisk-fr" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to basilisk-fr...@googlegroups.com.
> To view this discussion visit https://groups.google.com/d/msgid/basilisk-fr/e06ed267-125a-4743-b899-425ce3263d6dn%40googlegroups.com.

> #include "grid/multigrid3D.h"
> #include "embed.h"
> #include "navier-stokes/centered.h"
> #include "navier-stokes/perfs.h"
> #include "vof.h"
> #include "tension.h"
> #include "two-phase.h"
> #include "view.h"
> #include "maxruntime.h"
> #include <math.h>
> #include "common.h"
> #include "reduced.h"
>
> #define mu(f) (1./(clamp(f, 0, 1)*(1./mu1 - 1./mu2) + 1./mu2))
>
> // Simulation properties to change
> const double MAXTIME = 2.0;
> const int LEVEL = 12;
>
> // Geometrical properties of experiment
> const double PI = 3.1415926;
> const double DIAMETER = 0.1953, WIDTH = 0.2 [1], HEIGHT = 8.0;
> const int NX = 40; //this should be equal to height/width
> const int NY = 1, NZ = 1;
>
> const double VLiqIn = 1.028949;
> const double VAirIn = 2.312438;
> const double Vrad = 0.2;
> const double Dinj = 0.014;
> #define NHOLES 32 // makes 2x this number of holes, two rings of them
> const int Nholes = NHOLES;
> double angles[NHOLES];
>
> // Global variables
> bool master;
> scalar f0[], vInjY[], vInjZ[];
>
> // Boundary conditions
> u.n[top] = dirichlet(0.0);
> u.n[bottom] = dirichlet(0.0);
> u.n[front] = dirichlet(0.0);
> u.n[back] = dirichlet(0.0);
> u.n[left] = dirichlet(f0[] > 1e-3 ? VAirIn:VLiqIn);
> u.n[right] = neumann(0.0);
>
> u.t[top] = dirichlet(0.0);
> u.t[bottom] = dirichlet(0.0);
> u.t[front] = dirichlet(0.0);
> u.t[back] = dirichlet(0.0);
> u.t[left] = dirichlet(vInjY[]);
>
> u.r[top] = dirichlet(0.0);
> u.r[bottom] = dirichlet(0.0);
> u.r[front] = dirichlet(0.0);
> u.r[back] = dirichlet(0.0);
> u.r[left] = dirichlet(vInjZ[]);
>
> u.n[embed] = dirichlet(0.0);
> u.t[embed] = dirichlet(0.0);
> u.r[embed] = dirichlet(0.0);
>
> f[left] = (f0[]);
> p[left] = neumann(0);
> pf[left] = neumann(0);
>
> f[right] = neumann(0);
> p[right] = dirichlet(0.);
> pf[right] = dirichlet(0.);
>
> void updateAperatures()
> {
> vertex scalar phi[];
>
> foreach_vertex()
> {
> phi[] = sq(DIAMETER/2.0) - sq(y) - sq(z);
> }
>
> fractions(phi, cs, fs);
> fractions_cleanup(cs, fs);
>
> // Define angles for injection holes
>
> for (int i = 0; i < Nholes; i++) {
> angles[i] = i*PI/(Nholes/2);
> }
> }
>
> int main (int argc, char * argv[]) {
> maxruntime(&argc, argv);
> #if !_MPI
> master = true;
> #else
> master = (mpi_rank == 0);
> #endif
>
> // Air-water properties at 0.25 MPa operating conditions
>
> rho1 = 2.35 [0]; // air density (average) kg/m^3, 2.79 at inlet, 1.91 at outlet
> rho2 = 997.561 [0]; // water density kg/m^3
> mu1 = 1.85508e-5; // air viscocity Pa-s
> mu2 = 8.8871e-4; // water viscocity Pa-s
> f.sigma = 0.072; // surface tension N/m
> G.x = -9.81; // gravity in x-direction m/s^2
>
> // Simulation domain size setup
>
> size(HEIGHT);
> dimensions(nx = NX, ny = NY, nz = NZ);
> origin (0., -WIDTH/2.0, -WIDTH/2.0);
> init_grid(pow(2, LEVEL));
>
> TOLERANCE = 1e-3 [*];
> CFL = 0.5;
>
> run();
> }
>
> event init (t = 0)
> {
> updateAperatures();
>
> // Initialize the injection holes at the inlet with a void of 1
> // and the corresponding radial velocity
>
> scalar f00[], f01[];
>
> for(int i = 0; i < Nholes; i++)
> {
> // Make first set of holes, 32 around circumference right around the edge
> fraction(f00, sq(Dinj/2.0) - sq(y - (DIAMETER-Dinj)/2.0*cos(angles[i]))
> - sq(z - (DIAMETER-Dinj)/2.0*sin(angles[i])));
> foreach(){
> f0[] = f0[] + f00[];
> vInjY[] = vInjY[] - Vrad*cos(angles[i])*f00[];
> vInjZ[] = vInjZ[] - Vrad*sin(angles[i])*f00[];
> }
>
> // Make second set of holes, 32 around circumference right around the first set
> fraction(f01, sq(Dinj/2.0) - sq(y - (DIAMETER-3.0*Dinj)/2.0*cos(angles[i] + PI/(Nholes)))
> - sq(z - (DIAMETER-3.0*Dinj)/2.0*sin(angles[i] + PI/(Nholes))));
> foreach(){
> f0[] = f0[] + f01[];
> vInjY[] = vInjY[] - Vrad*cos(angles[i])*f01[];
> vInjZ[] = vInjZ[] - Vrad*sin(angles[i])*f01[];
> }
> }
>
> boundary((scalar*){u});
>
> foreach()
> {
> // Initialize velocity everywhere
> u.x[] = (cs[] > 0.) ? VLiqIn:0.;
> }
>
> }
>
> event endImage (t = MAXTIME)
> {
> view(camera = "front", fov = 15., width = 800, ty = 0., tx = -0.5,
> height = 600, bg = {1,1,1});
> box();
> squares("u", alpha = 0.0001, min = 0, max = 3);
> colorbar(horizontal = true, label = "V (m/s)", min = 0, max = 3);
> save("u.ppm");
> clear();
> }

--
/^..^\
( (••) ) Wojciech (Voietec [vôitek]) ANISZEWSKI
(|)_._(|)~ Research Assistant @ Univ-Rouen
GPG ID : AC66485E
Twitter : @echo_dancers3
Mastodon: @w...@fediscience.org
BlueSky : @aniszewski.bsky.social
Scholar : https://tinyurl.com/y28b8gfp
OrcId : https://orcid.org/0000-0002-4248-1194
RG : https://www.researchgate.net/profile/Wojciech_Aniszewski