MATLAB 2D Fluid Simulation

Incompressible 2D MATLAB FDM
fluid simulation calculated with Jacobi method and
Runge-Kutta 4 integration for advection.

UPDATE 2024:
I’m excited to announce that the MATLAB code for the fluid simulation is now available on GitHub! Thank you all for your patience and interest in this project. You can find the full code, along with a README, at the following link:

Feel free to check it out, try it for yourself, and let me know what you create with it! If you use the code for your own animations or projects, please give credit to me. Happy coding!

[jackellis]

My computer has had nightmares ever since seeing this

[huyphanducnhat]

wow amazing
i have never code that much in matlab, i 99% did all the things with simulink when i was in uni

[cedricliekonwah]

Is it possible to have a hint on how you code this?

[Hzur]

Now I haven’t studied university level physics (high school physics only) but did you use the navier stokes equations and cfd to simulate this? It’s very fascinating to see! I’m excited since I’m planning to study engineering and I’m sure this will be covered

[stivenmass]

Hi Robert, impressive simulation and nice visualisation. Would you mind sharing any detail, paper or the exact equation you’re solving here ?

Thanks

[YourPhysicsSimulator]

Impressive video
May I ask?
I'm having quite a lot of trouble trying to make a simulation like this. The advection part is giving me a lot of problems.
Could you explain how you apply it with MATLAB? How do you implement RK in this advection part?
Thanks!

[Flo]

Could you give some more information about the general setup? Is this a compressible simulation? Otherwise the box would need an outflow somewhere, right? What is the dimension of the domain? 200 by 400 meters? What is the time scale of this simulation from start to end in real time?

[YourPhysicsSimulator]

Okay, now I understand the problem, Robert! I'm sorry for commenting again, but I am too close to making this work out to give up
So interp2 works magically when the velocities point to a region defined with meshgrid. I have seen the fluid work fine for a moment.
However when the velocities point outside the defined rectangle it returns as default NaN, and interp2 collapses.
How did you avoid this problem?
Do you create a bigger rectangle and confine your fluid in a smaller rectangle inside this one?

[erencolak]

I am having hard time to understand this simulation. Your inlet is not an inlet (because it is a closed enclosure, if you inject more fluid it wouldn't be incompressible), and not really placed at the boundary. How does it even work?

[ingmoser]

Hi @Robert S. really beautiful, can you share the code?

[Slyerkid]

I would want to know how, but I never want to touch matlab again in my life.

[Tidepodious]

Is there any particular reason you used 4point Jacobi stencil instead of the 5point? So [0 1 0; 1 -4 1; 0 1 0] instead of [0 1 0; 1 0 1; 0 1 0]. Was this to give the diffsion effect/smoothing of the pressure as opposed to the more accurate 5point?

[karstenkunneman]

>Runge-Kutta 4
oh God I'm having flashbacks

[polloconcoca]

How can I learn this? is there some book? or notes? or course?

[mohammedzainulla]

Can i now how to increse the velocity of the flow?

[JorgeJimenez-desv]

How many time did it take to simulate it?
Can you share us your computer specs?

[electromecanica_automatismos]

Yo solo recuerdo la intro LATAM de ULTRAMAN

[kingarthurtheth]

Anyone else try to follow one the initial swirls as long as possible?...

[chawlz]

You will not make it in modern engineering if you can’t use matlab.


E-E-E-E-E-E-E-E-E-E

[mukeshkalel]

Hello @iaseaaa. I have follow through your comment and tried to perform this simulation. can you please look into my code and doe this look right. I don't know what to do next?

clear all
clc

s = 200; % height of the sim
[X, Y] = meshgrid(1:s*2, 1:s);
[p, vx, vy] = deal(zeros(s, s*2));

vx(95:105, 10:15) = 0.1;
[dx, dy] = gradient(p);
vx(2:end-1, 2:end-1) = vx(2:end-1, 2:end-1)-dx(2:end-1, 2:end-1);
vy(2:end-1, 2:end-1) = vy(2:end-1, 2:end-1)-dy(2:end-1, 2:end-1);

[pvx, pvy] = RK4(X, Y, vx, vy, -1);
vx = interp2(vx, pvx, pvy);
vy = interp2(vy, pvx, pvy);

rhs = -divergence(vx, vy);
J = [0 1 0; 1 0 1; 0 1 0]/4;

for i=0:100
p = conv2(p, J, 'same')+rhs/2;
end

function [x_new, y_new] = RK4(px, py, vx, vy, h)
k1x = interp2(vx, px, py);
k1y = interp2(vy, px, py);
k2x = interp2(vx, px+h/2*k1x, py+h/2*k1y);
k2y = interp2(vy, px+h/2*k1x, py+h/2*k1y);
k3x = interp2(vx, px+h/2*k2x, py+h/2*k2y);
k3y = interp2(vy, px+h/2*k2x, py+h/2*k2y);
k4x = interp2(vx, px+h*k3x, py+h*k3y);
k4y = interp2(vy, px+h*k3x, py+h*k3y);
x_new =
y_new =
end

can you help me on moving forward?