∂w/∂t + w grad w + grad p/ρ = 𝜈 div grad w + f
div v = 0
excitation: ∂w/∂t = f
this is the simple part: multiply the force field f with the time step and add it to the velocity field.
convection: w grad w
is the nonlinear term of fluid flow, which makes it interesting.
the change of velocity is the local time derivative ∂w/∂t and what flows in if there is a gradient: w grad w.
it's also called self-convection because it's the velocity which is transported by itself.
forward integration can be unstable if the step size is too large.
to keep it stable we integrate backwards: update the velocity with the value of the cells where it comes from.
we use bilinear interpolation between the source cell and it's 4 neighbor cells.
diffusion: ∂w/∂t = 𝜈 div grad v
with the laplace operator on the right side. it's the same form as the heat equation.
differences equal out with time.
this part is not present in the euler equations of inviscous flow.
we solve it under the fft together with the projection step.
fourier
fft:{$[-1+n:#x;(x+r),(x:fft x o)-r:(1@(!n)*-180.%n)*fftx 1+o:2*!n%:2;x]} /radix-2 recursive odd/even split
ifft:{(conj fft conj x)%#x} /inverse is fft under conjugation
rfft:{.5*(x+y;1a270*x-y:conj(x:fft imag[x;y])n!n-!n:#x)} /real241 two for the price of one
fft2:+fft'+fft' /two-dim separate and do it again
projection: handle grad p and div v = 0
we can ignore grad p if we keep the velocity field solenoidal.
there are no sources and the conservation of mass demands div v = 0.
in 2d this is 0=du/dx+dv/dy. in fourier space this is 0=(X*U)+Y*V,
e.g. the scalar product with the wave vector (X;Y) is 0.
this means, the velocity field in fouier space is perpendicular to the wave vector.
we have to project it, to keep it that way.
the projection P keeps U and V orthogonal to the vector of wave numbers X,Y
e.g. it substracts the part in direction of X,Y:
P(U) = (1-XX/RR) U - (XY/RR) V
P(V) = (- XY/RR) U + (1-YY/RR) U, with RR:XX+YY
X and Y are the wavenumbers X:0.+(n*n)#(!m),m,1_(-m)+!m:n%2 in the order after the transform.
we can do the 2d-fft on the complex velocity field if we combine w:u+iv.
to recover U and V we had to do:
U:0.5*(W+conj W I) and V:-i*0.5*(W-conj W I) with I:n/n!n-!n,n (n/decode n!modulo odometer!n,n)
and before transforming back, we combine the parts: P(W) = P(U) + i*P(V).
let's get rid of the decomposition and work on W directly: