# Intro

This is a minimal example for implementing a ‘fan’ after geometry examples that can be found here. The mentioned function is updated to be applicable for general phi and theta. Also this script adds a velocity based on a stated fan power. All functions can be found at the end of this file.

# Power to velocity forcing and expected exit velocity

The implementation is based on a momentum forcing. An amount of work per time is done on the flow with a certain power, P. Using the work-energy theorem, which states that the work done by all acting forces on a particle is equal the change in kinetic energy, this can be expressed as \Delta E_{kin}=P\Delta t. This directly implies that,

\displaystyle \vec{v}^2_{t+\Delta t} - \vec{v}^2_{t} = \frac{2\textrm{P}\Delta t}{m}\hat{n}^2

which in turn implies for every spatial dimension,

\displaystyle v^2_{i,t+\Delta t} - v^2_{i,t} = \frac{2\textrm{P}\Delta t}{m} n_i^2

From these equations the direction of \hat{n} is not yet defined and |\hat{n}|=1. When it is assumed that the force applied to the flow is normal to the plane in which the fan blades rotate then \hat{n} is the normal vector of this plane. Rearranging and taking the square root results in, \displaystyle \vec{v}_{i,t+\Delta t} = \sqrt{\vec{v}^2_{i,t} + \frac{2\textrm{P}\Delta t}{m}\hat{n}^2_i }

To estimate the velocity to which air is accelerated by this forcing we assume that it is applied over a finite time from t_{start} to t_{end}, moving a distance w in the process: \displaystyle w=\int_{t_{start}}^{t_{end}} v\left(t\right)\textrm{dt} Assuming a zero entrance velocity, v is given at every instance as: \displaystyle \textrm{E}_{\textrm{kin}}=\textrm{P}\Delta t=0.5mv^2 \displaystyle v\left(t\right)=\sqrt{2\textrm{P}t/m} P, w and m are set by the user thus this leaves two equations and two unknowns. Solving for the exit velocity gives: \displaystyle v_{end} = \left(\frac{3\textrm{P}w}{m}\right)^{1/3}

# The code

First we include relevant files for 2D navier stokes solving and geometric fun.

//#include "grid/octree.h" 		// For 3D --> to demanding for server
#include "navier-stokes/centered.h"     // Navier stokes
#include "fractions.h"	      		// Needed to compute fan fractions

Required functions and structures

void init_rotor();
void rotor_update();
void rotor_coord();
void rotor_forcing();

struct sRotor {
double rampT;		// Time to start up rotor
double P, Prho;		// Power, powerdensity
double R, W, A, V;	// Diameter, Thickness, Area ,Volume
double diaVol;		// Diagnosed rotor volume
double x0, y0, z0;	// Origin of rotor
double xt, yt, zt;      // Translation of origin per rotate event
double theta, phi;	// Polar and Azimuthal angle
double thetat, phit;    // Addition to angles per rotate event
double Work;            // Work done by rotor
double cu;              // Characteristic velocity
bool rotate;            // Rotation yes or no
coord nf, nr;		// Normal vector fan, rotation
};

Some global vars

int minlevel, maxlevel;         // Grid depths
struct sRotor rot;  		// Rotor details structure
scalar fan[];			// Fan volume fraction
scalar lev[]; 			// for movie diagnostic level
scalar Ekin[];  		// for movie kinetic energy

Starting with the most loved c function

int main() {
minlevel = 3;
maxlevel = 8;

L0 = 10.;
X0 = Y0 = Z0 = 0.;

init_grid(1<<7);

foreach_dimension(){
periodic (left);
u.x.refine = refine_linear; 		// Momentum conserved
}

fan.prolongation = fraction_refine;		// Fan is a volume fraction
p.refine = p.prolongation = refine_linear;

DT = 10E-5;					// For poisson solver
TOLERANCE=10E-6;				// For poisson solver
CFL = 0.5;					// CFL condition

run();						// Start simulation
}

Initialisation

event init(t = 0){
rot.rotate = true; 				// Turn on rotation
init_rotor();					// See init function for details
fan.prolongation = fraction_refine;		// Tell basilisk that fan is volumetric field
refine (fan[] > 0. && level < maxlevel); 	// Refine where fan is
eps = 0.07*rot.cu;				// Wavelet estimator based on flow velocity
}

Reset to defaults for poisson solver

event reset_pois(i = 10){
DT = 1;
TOLERANCE = 10E-3;
}

event adapt(i++) {
}

Forcing by the rotor

event forcing(i = 1; i++) {
rotor_coord();
rotor_forcing();
}

Rotate the rotor

event rotate(t+=0.1) {
if(rot.rotate) {
// Change center
rot.x0 += rot.xt;
rot.y0 += rot.yt;
rot.z0 += rot.zt;

// Change angles
rot.theta += rot.thetat;
rot.phi += rot.phit;

rotor_update();
}
}

Visual output

event movies(t+=1){
fprintf(stderr, "t=%g\n", t);
foreach(){
lev[] = level;
Ekin[] = sqrt(sq(u.y[]) + sq(u.x[]));
}
boundary({lev, Ekin});

output_ppm (lev, file = "ppm2mp4 level.mp4", n = 512, linear = false, min = minlevel, max = maxlevel);
output_ppm (Ekin, file = "ppm2mp4 ekin.mp4", n = 512, linear = true);

}

End the simulation

event end(t=180){
}

# Results

Kinetic energy density

Grid depth

# Functions

Function returning the sRotor structure, includign default properties

void init_rotor() {
rot.Work = 0.;
if(!rot.rampT)
rot.rampT = 1.; // fan at full thrust after 1 second, linear ramp
if(!rot.R)
rot.R = L0/80.;
if(!rot.W)
rot.W = rot.R/5;
if(!rot.Prho)
rot.Prho = L0/10.;
if(!rot.x0)
rot.x0 = L0/2.;
if(!rot.y0)
rot.y0 = L0/2.;
if(!rot.z0){
#if dimension == 2
rot.z0 = 0.;
#elif dimension == 3
rot.z0 = L0/2.;
#endif
}
if(!rot.theta)
rot.theta = 90.*M_PI/180.;      // Polar angle
if(!rot.phi)
rot.phi = 90.*M_PI/180.;        // Azimuthal angle

if(rot.rotate) {
// determines for every rotate event how much needs to be added to x, y, z, tehta, phi.
rot.xt = 0;
rot.yt = 0;
rot.zt = 0;
rot.thetat = 0.;
rot.phit = -0.1*M_PI/180.; // note that this rate is linked to the frequency of the rotate event
} else {
rot.xt = 0;
rot.yt = 0;
rot.zt = 0;
rot.thetat = 0.;
rot.phit = 0.;
}
rotor_update();
}

Updating relevant rotor variables

void rotor_update() {
// Normal vectors to the fan plane, note that z and y are turned around for consistency regarding polar coordinates (basilisk uses y for the vertical)
rot.nf.x = sin(rot.theta)*cos(rot.phi);
rot.nf.y = sin(rot.theta)*sin(rot.phi);
rot.nf.z = cos(rot.theta);

rot.nr.x = sin(rot.theta)*cos(rot.phi);
rot.nr.y = sin(rot.theta)*sin(rot.phi);
rot.nr.z = cos(rot.theta);

#if dimension == 2
rot.A = 2*rot.R*rot.W;
#elif dimension == 3
rot.A = sq(rot.R)*M_PI;
#endif

// Volume is slice of a sphere
#if dimension == 2
rot.V = 1.*rot.A;
#elif dimension == 3
rot.V = 4.*M_PI*pow(rot.R,3.)/3. -
2*M_PI*pow(rot.R-rot.W/2., 2.)/3.*(2*rot.R + rot.W/2.);
#endif

rot.P = rot.V*rot.Prho;
rot.cu = pow(3*rot.Prho*rot.W, 1./3.);
}

Function returning the volume fractions of a fan object

void rotor_coord() {
scalar sph[], plnu[], plnd[];
fraction(sph, -sq((x - rot.x0)) - sq((y - rot.y0)) - sq((z - rot.z0)) + sq(rot.R)); // sphere
fraction(plnu,  rot.nr.x*(x - rot.x0) + rot.nr.y*(y - rot.y0) + rot.nr.z*(z - rot.z0) + rot.W/2.); // upper plane
fraction(plnd, -rot.nr.x*(x - rot.x0) - rot.nr.y*(y - rot.y0) - rot.nr.z*(z - rot.z0) + rot.W/2.); // lower plane

foreach () {
fan[] = sph[]*plnu[]*plnd[]; // estimate for the overlappnig volumes, not that this is not exact at the boundaries of the 'fan', this is corrected for later in the focring by comparing diagnosed to theoretical volume.
}
boundary({fan}); // update boundaries
}

Function returning new velocities based on a rotor forcing. This is a function of powerdensity, width, direction, ramp-time, and diagnosed volume

void rotor_forcing(){
double tempW = 0.;
double w, wsgn, damp, usgn, utemp, corrP;
foreach(reduction(+:tempW)) {
if(fan[] > 0.) {
foreach_dimension() {
wsgn = sign(rot.nf.x*u.x[]) + (sign(rot.nf.x*u.x[]) == 0)*sign(rot.nf.x); // Get the direction of the flow relative to forcing direction
damp = rot.rampT > t ? t/rot.rampT : 1.;                                  // Linear ramp for starting up the forcing
corrP = rot.diaVol > 0. ? rot.V/rot.diaVol : 1.;                         // Correction if diagnosed volume differs from theoretical volume
w = wsgn*fan[]*damp*sq(rot.nf.x)*(2./rho[])*(corrP*rot.P/rot.V)*dt;      // additional kinetic energy
tempW += 0.5*rho[]*w*dv(); // save work done to asses total work in diagnostics

// New kinetic energy
utemp = sq(u.x[]) + w;
// New sign of the velocity
usgn = 	  1.*(u.x[] >= 0)*(utemp > 0) +
-1.*(u.x[] >= 0)*(utemp < 0) +
1.*(u.x[] <  0)*(utemp < 0) +
-1.*(u.x[] <  0)*(utemp > 0);
// New velocity
u.x[] = usgn*sqrt(fabs(utemp));
}
}
}
rot.Work += tempW; // Diagnostics
}