simulationCases/dropImpact.c
Drop Impact Simulation
This file contains the simulation code for drop impact on a solid surface using a Volume of Fluid (VOF) method with viscoelastic fluid modeling capabilities.
The simulation uses an adaptive mesh refinement approach to efficiently resolve the interface dynamics during impact. The model can handle both Newtonian and viscoelastic fluids through a log-conformation formulation for the polymeric stress tensor.
Physical Model
The simulation solves the Navier-Stokes equations for incompressible flow coupled with a viscoelastic constitutive equation. The interface between the two fluids is tracked using a VOF method with surface tension.
Key dimensionless parameters:
- Weber number (\(We\)): Ratio of inertia to surface tension
- Ohnesorge number (\(Oh\)): Ratio of viscous to inertial and surface tension forces
- Deborah number (\(De\)): Ratio of relaxation time to characteristic flow time
- Elastocapillary number (\(Ec\)): Ratio of elastic to capillary forces
@file dropImpact.c @author Vatsal Sanjay @version 0.2 @date Oct 18, 2024
// #include "axi.h"
#include "grid/octree.h"
// #include "grid/quadtree.h"
#include "navier-stokes/centered.h"
#define VANILLA 036
#if VANILLA
#include "log-conform-viscoelastic.h"
#define logFile "logAxi-vanilla.dat"39
#else
#if AXI
#include "log-conform-viscoelastic-scalar-2D.h"
#define logFile "logAxi-scalar.dat"43
#else
#include "log-conform-viscoelastic-scalar-3D.h"
#define logFile "log3D-scalar.dat"46
#endif
#endif
#define FILTERED // Smear density and viscosity jumps
#include "two-phaseVE.h"
#include "navier-stokes/conserving.h"
#include "tension.h"
Simulation Parameters
Timing Parameters
- tsnap: Time interval between snapshots (1e-2)
#define tsnap (1e-2)6263
Numerical Error Tolerances
- fErr: Error tolerance in VOF function f1 (1e-3)
- KErr: Error tolerance in VOF curvature calculated using height function method (1e-6)
- VelErr: Error tolerances in velocity (1e-2)
- Use 1e-2 for low Oh and 1e-3 to 5e-3 for high Oh/moderate to high J
#define fErr (1e-3)71
#define KErr (1e-6)72
#define VelErr (1e-2)7374
Domain Configuration
- xDist: Initial horizontal displacement of droplet center from boundary (5e-2)
- R2: Function to calculate squared distance from droplet center
#define xDist (5e-2)80
#define R2(x, y, z) (sq(x - 1. - xDist) + sq(y) + sq(z))8182
Boundary Conditions
Dirichlet boundary condition for volume fraction at left boundary. Note: No-slip boundary conditions are commented out for testing purposes.
// u.t[left] = dirichlet(0.); // todo: later on use no-slip. For testing, free-slip is faster.
// u.r[left] = dirichlet(0.); // todo: later on use no-slip. For testing, free-slip is faster.
f[left] = dirichlet(0.0);
Global Variables
- MAXlevel: Maximum level of grid refinement
- We: Weber number of the drop
- Oh: Solvent Ohnesorge number
- Oha: Air Ohnesorge number
- De: Deborah number
- Ec: Elasto-capillary number
- tmax: Maximum simulation time
- nameOut: Output filename for snapshots
- dumpFile: Filename for restart dumps
int MAXlevel;
double We, Oh, Oha, De, Ec, tmax;
char nameOut[80], dumpFile[80];
Main Function
Initializes the simulation parameters and grid, then starts the simulation.
Parameters
- argc: Command-line argument count
- argv: Command-line argument array
Return Value
Returns 0 on successful completion
int main(int argc, char const *argv[]) {
dtmax = 1e-5;
L0 = 4.0;
// Values taken from the terminal
MAXlevel = 6;
tmax = 3.0;
We = 5.0;
Oh = 1e-2;
Oha = 1e-2 * Oh;
De = 1.0;
Ec = 1.0;
init_grid(1 << 4);
// Create a folder named intermediate where all the simulation snapshots are stored
char comm[80];
sprintf(comm, "mkdir -p intermediate");
system(comm);
// Name of the restart file
sprintf(dumpFile, "restart");
// Set fluid properties for both phases
rho1 = 1.0, rho2 = 1e-3;
mu1 = Oh / sqrt(We), mu2 = Oha / sqrt(We);
G1 = Ec / We, G2 = 0.0;
lambda1 = De * sqrt(We), lambda2 = 0.0;
f.sigma = 1.0 / We;
run();
return 0;
}
Initialization Event
Sets up the initial condition for the droplet and velocity field.
The event initializes a spherical droplet using the VOF method and sets an initial horizontal velocity. If a restart file exists, the simulation state is restored from it instead.
@param t Time (set to 0 for initialization)
event init(t = 0) {
if (!restore(file = dumpFile)) {
refine(R2(x, y, z) < (1.1) && R2(x, y, z) > (0.9) && level < MAXlevel);
fraction(f, (1 - R2(x, y, z)));
foreach() {
u.x[] = -f[] * 1.0;
}
}
}
Adaptive Mesh Refinement
Refines the computational mesh based on wavelet error estimates for the tracked fields to efficiently allocate computational resources.
The refinement criteria track errors in: - Volume fraction field (f) - Velocity components (u.x, u.y, u.z)
@param i Current iteration number
event adapt(i++) {
adapt_wavelet((scalar *){f, u.x, u.y, u.z},
(double[]){fErr, VelErr, VelErr, VelErr},
MAXlevel, 4);
}
Snapshot Writing
Creates periodic dumps of the simulation state for visualization and restarts.
@param t Current simulation time @param tsnap Time interval between snapshots @param tmax Maximum simulation time
event writingFiles(t = 0; t += tsnap; t <= tmax) {
dump(file = dumpFile);
sprintf(nameOut, "intermediate/snapshot-%5.4f", t);
dump(file = nameOut);
}
Simulation Termination
Outputs final simulation parameters to standard error when the simulation ends.
@param t End time
event end(t = end) {
if (pid() == 0)
fprintf(ferr, "Level %d, Oh %2.1e, We %2.1e, Oha %2.1e, De %2.1e, Ec %2.1e\n",
MAXlevel, Oh, We, Oha, De, Ec);
}
Simulation Logging
Records simulation progress and checks for numerical stability.
This event: 1. Calculates the total kinetic energy of the system 2. Writes simulation data to the log file 3. Checks for numerical stability based on kinetic energy values 4. Terminates the simulation if energy values indicate instability
@param i Current iteration number @return Returns 1 (terminating the simulation) if instability is detected
event logWriting(i++) {
// Calculate kinetic energy
double ke = 0.;
foreach(reduction(+:ke)) {
ke += (2 * pi * y) * (0.5 * rho(f[]) * (sq(u.x[]) + sq(u.y[]))) * sq(Delta);
}
static FILE *fp;
if (pid() == 0) {
const char *mode = (i == 0) ? "w" : "a";
fp = fopen(logFile, mode);
if (fp == NULL) {
fprintf(ferr, "Error opening log file\n");
return 1;
}
if (i == 0) {
fprintf(ferr, "Level %d, Oh %2.1e, We %2.1e, Oha %2.1e, De %2.1e, Ec %2.1e\n",
MAXlevel, Oh, We, Oha, De, Ec);
fprintf(ferr, "i dt t ke\n");
fprintf(fp, "Level %d, Oh %2.1e, We %2.1e, Oha %2.1e, De %2.1e, Ec %2.1e\n",
MAXlevel, Oh, We, Oha, De, Ec);
fprintf(fp, "i dt t ke rM\n");
}
fprintf(fp, "%d %g %g %g\n", i, dt, t, ke);
fprintf(ferr, "%d %g %g %g\n", i, dt, t, ke);
fflush(fp);
fclose(fp);
}
assert(ke > -1e-10);
// Stability check based on kinetic energy
if (i > 1e1 && pid() == 0) {
if (ke > 1e2 || ke < 1e-8) {
const char *message = (ke > 1e2) ?
"The kinetic energy blew up. Stopping simulation\n" :
"Kinetic energy too small now! Stopping!\n";
fprintf(ferr, "%s", message);
fp = fopen("log", "a");
fprintf(fp, "%s", message);
fflush(fp);
fclose(fp);
dump(file = dumpFile);
return 1;
}
}
}