N = 192; dt = 0.25; p = zeros(N, N);
u = zeros(N, N); v = zeros(N, N); ustar = zeros(N, N); vstar = zeros(N, N);
k = pi*[0:(N/2-1) (-N/2):(-1)]; [x y] = meshgrid(k, k); ds = -(x.^2 + y.^2);
ds(1,1) = 1; px = 0.0; py = 0.0; gx = 0.0; gy = 0.0;
for t = 0:1000
% constant velocity zone
u(N/2-10:N/2+10, N/2-10:N/2+10) = 0.2;
% solve intermediate velocity U*
[ux uy] = gradient(u);
[vx vy] = gradient(v);
ustar = u + (u.*ux + v.*uy + del2(u)/20 - px/20 + gy)*dt;
vstar = v + (u.*vx + v.*vy + del2(v)/20 - py/20 + gx)*dt;
% solve Poisson equation: del2(p) = f
fhat = fft2(divergence(ustar, vstar));
p = N^2*real(ifft2(fhat./ds));
% enforce del(U) = 0
[px py] = gradient(p);
u = ustar - px*dt;
v = vstar - py*dt;
cla; surface(sqrt(u.^2 + v.^2), 'edgecolor', 'none'); drawnow;
end
Espoo is full of furry animals. |
Just a basic wave equation...
u = zeros(192, 192);
du = zeros(192, 192);
[x y] = meshgrid(linspace(-1, 1, 192), linspace(-1, 1, 192));
for t = 0:0.002:0.4
u = u + exp(-1000*((x-0.7*sin(7*t)).^2+(y-0.7*sin(11*t)-0.2).^2));
u = u + 0.5*exp(-1000*((x+0.7*sin(5*t)).^2+(y+0.7).^2))*sin(3*31*t);
du = du + del2(u);
u = u + du;
cla; imagesc(u); drawnow;
end
u = zeros(192, 192);
du = zeros(192, 192);
[x y] = meshgrid(linspace(-1, 1, 192), linspace(-1, 1, 192));
for t = 0:0.002:0.4
u = u + exp(-1000*((x-0.7*sin(7*t)).^2+(y-0.7*sin(11*t)-0.2).^2));
u = u + 0.5*exp(-1000*((x+0.7*sin(5*t)).^2+(y+0.7).^2))*sin(3*31*t);
du = du + del2(u);
u = u + du;
cla; imagesc(u); drawnow;
end
...
No comments:
Post a Comment