sandbox/Antoonvh/filaments.h
Vector induction by line
This is a solver for the equation:
\displaystyle \mathbf F(\mathbf{x}) = \int_C \mathrm{d}\mathbf{F} = \int_C f(r)\mathrm{d}\mathbf{l} .
Where, \mathbf{F} is a vector field, C is a line or loop, and f is a function that that takes as argument the distance of the point \mathbf{x} and the line segment with length and direction \mathrm{d}\mathbf{l}.
An example
Consider an infinite straight line from \mathbf{x_1} = \{ 0,0,- \infty \} to \mathbf{x_2} = \{ 0,0,\infty \} and f(r) = \frac{\Gamma}{\pi^{3/2}a^3} e^{-\frac{r^2}{a^2}}. Then, using the relevant coordinates \rho = \sqrt{x^2 + y^2} \text{\ and\ } z
\displaystyle \mathbf{F}(\mathbf{x}) = \mathbf{F}(\rho) = \int_{-\infty}^{\infty} \Gamma e^{-\frac{\rho^2 - z^2}{a^2}} d\mathrm{z} = \displaystyle \Gamma e^{-\frac{\rho^2}{a^2}}\int_{-\infty}^{\infty}e^{-\frac{z^2}{a^2}} \mathrm{d}z \mathbf{e_z} = \frac{\sqrt{\pi}a}{\pi^{3/2}a^3}\Gamma e^{-\frac{\rho^2}{a^2}}\mathbf{e_z} = \\ \frac{\Gamma}{\pi a^2}e^{-\frac{\rho^2}{a^2}} \mathbf{e_z}.
Which is closely related to the Lamb-Oseen vortex or Gaussian vortex of size a and strength \Gamma.
Implementation
We follow A. Castillo’s filaments.h/c functions, and compute F from a discrete parameteric curve \mathbf{C}(t) and its parametric direction vector \mathrm{d}\mathbf{l}(t). Noting that the magnitude of the direction vector can be used to have a varying strength along the line.
trace
coord get_vorticity (coord X, int n_seg, coord * c, coord * dl, double (* fun)(double r)) {
coord F = {0, 0, 0};
for (int i = 0; i < n_seg; i++) {
double r = 0;
foreach_dimension()
r += sq(X.x - c[i].x);
double f = fun(r);
foreach_dimension()
F.x += dl[i].x*f;
}
return F;
}
trace
void get_vor_vector (vector omg, coord (* C)(double t), double t0,
double te, int n_seg, double (* fun)(double r)) {
coord c[n_seg], dl[n_seg];
double Tl = te - t0, Dt = Tl/(double)n_seg;
int i = 0;
double sqrt3 = sqrt(3);
for (double ti = t0 + Dt/2; ti < te; ti += Dt) {
c[i] = C(ti);
coord d1 = C(ti - Dt/2), d2 = C(ti + Dt/2);
foreach_dimension()
dl[i].x = d2.x - d1.x;
i++;
}
foreach() {
foreach_dimension()
omg.x[] = 0.;
coord cc = {x, y, z}, vort = {0, 0, 0};
foreach_child() {
coord ccc;
foreach_dimension()
ccc.x = cc.x + child.x*Delta/sqrt3;
coord vor = get_vorticity (ccc, n_seg, c, dl, fun);
foreach_dimension()
vort.x += vor.x/(double)(1 << dimension);
}
foreach_dimension()
omg.x[] = vort.x;
}
boundary ((scalar*){omg});
}
double Gamma_LO = 1, a_LO = 1;
double Lamb_Oseen (double sqr) {
return Gamma_LO/(pow(pi, 3./2.)*cube(a_LO))*exp(-(sqr/sq(a_LO)));
}