MATLAB code: Two-phase incompressible and immiscible fluid solver based on IBM
Date:
%Main code
global dt Nb N h rho mu ip im a;
global kp km dtheta K;
initialize
init_a
for clock=1:400 %clockmax
XX=X+(dt/2)*interp(u,X);
ff=spread(Force(XX),XX);
[u,uu]=fluid(u,ff);
X=X+dt*interp(uu,XX);
figure(1);
%animation:
vorticity=(u(ip,:,2)-u(im,:,2)-u(:,ip,1)+u(:,im,1))/(2*h);
contour(xgrid,ygrid,vorticity);
colormap jet;
hold on
plot(X(:,1),X(:,2),’k.-‘,’linewidth’,1,’markersize’,2)
axis([0,L,0,L])
caxis(valminmax)
axis equal
axis manual
drawnow;
axis([0 1 0 1]);
hold off
end
%fluid
function [uuu,uu]=fluid(u,ff)
global a dt rho mu;
w=u-(dt/2)skew(u)+(dt/(2rho))*ff;
w=fft(w,[],1);
w=fft(w,[],2);
uu(:,:,1)=a(:,:,1,1).w(:,:,1)+a(:,:,1,2).w(:,:,2);
uu(:,:,2)=a(:,:,2,1).w(:,:,1)+a(:,:,2,2).w(:,:,2);
uu=ifft(uu,[],2);
uu=real(ifft(uu,[],1));
w=u-dtskew(uu)+(dt/rho)ff+(dt/2)(mu/rho)laplacian(u);
w=fft(w,[],1);
w=fft(w,[],2);
uuu(:,:,1)=a(:,:,1,1).w(:,:,1)+a(:,:,1,2).w(:,:,2);
uuu(:,:,2)=a(:,:,2,1).w(:,:,1)+a(:,:,2,2).w(:,:,2);
uuu=ifft(uuu,[],2);
uuu=real(ifft(uuu,[],1));
%Force
function F=Force(X)
global kp km dtheta K;
F=K(X(kp,:)+X(km,:)-2X)/(dtheta*dtheta);
%init_a
global a;
a=zeros(N,N,2,2);
for m1=0:(N-1)
for m2=0:(N-1)
a(m1+1,m2+1,1,1)=1;
a(m1+1,m2+1,2,2)=1;
end
end
for m1=0:(N-1)
for m2=0:(N-1)
if~(((m1==0)|(m1==N/2))&((m2==0)|(m2==N/2)))
t=(2*pi/N)*[m1;m2];
s=sin(t);
ss=(s*s')/(s'*s);
% a(m1+1,m2+1,:,:)=a(m1+1,m2+1,:,:)-(ss’)/(s’s);
a(m1+1,m2+1,1,1)=a(m1+1,m2+1,1,1)-ss(1,1);
a(m1+1,m2+1,1,2)=a(m1+1,m2+1,1,2)-ss(1,2);
a(m1+1,m2+1,2,1)=a(m1+1,m2+1,2,1)-ss(2,1);
a(m1+1,m2+1,2,2)=a(m1+1,m2+1,2,2)-ss(2,2);
end
end
end
for m1=0:(N-1)
for m2=0:(N-1)
t=(pi/N)*[m1;m2];
s=sin(t);
a(m1+1,m2+1,:,:)=a(m1+1,m2+1,:,:)...
/(1+(dt/2)*(mu/rho)*(4/(h*h))*(s'*s));
end
end
%initialize
L=1.0
N=64
h=L/N
ip=[(2:N),1]
im=[N,(1:(N-1))]
Nb=ceil(pi*(L/2)/(h/2))
dtheta=2*pi/Nb
kp=[(2:Nb),1]
km=[Nb,(1:(Nb-1))]
K=1
rho=1
mu=0.01
tmax=4
dt=0.01
clockmax=ceil(tmax/dt)
for k=0:(Nb-1)
theta=k*dtheta;
X(k+1,1)=(L/2)+(L/4)*cos(theta);
X(k+1,2)=(L/2)+(L/4)*sin(theta);
end
u=zeros(N,N,2);
for j1=0:(N-1)
x=j1*h;
u(j1+1,:,2)=sin(2pix/L);
end
vorticity=(u(ip,:,2)-u(im,:,2)-u(:,ip,1)+u(:,im,1))/(2*h);
dvorticity=(max(max(vorticity))-min(min(vorticity)))/5;
values= (-10dvorticity):dvorticity:(10dvorticity);
valminmax=[min(values),max(values)];
xgrid=zeros(N,N);
ygrid=zeros(N,N);
for j=0:(N-1)
xgrid(j+1,:)=j*h;
ygrid(:,j+1)=j*h;
end
set(gcf,’double’,’on’)
contour(xgrid,ygrid,vorticity,values)
hold on
plot(X(:,1),X(:,2),’ko’)
axis([0,L,0,L])
caxis(valminmax)
axis equal
axis manual
drawnow
hold off
%interp
function U=interp(u,X)
global Nb h;
global N;
U=zeros(Nb,2);
for k=1:Nb
s=X(k,:)/h;
i=floor(s);
r=s-i;
i1=mod((i(1)-1):(i(1)+2),N)+1;
i2=mod((i(2)-1):(i(2)+2),N)+1;
w=phi1(r(1)).*phi2(r(2));
U(k,1)=sum(sum(w.*u(i1,i2,1)));
U(k,2)=sum(sum(w.*u(i1,i2,2)));
end
%laplacian
function w=laplacian(u)
global im ip h;
w=(u(ip,:,:)+u(im,:,:)+u(:,ip,:)+u(:,im,:)-4u)/(hh);
%phi1
function w=phi1(r)
w=zeros(4,4);
q=sqrt(1+4r(1-r));
w(4,:)=(1+2*r-q)/8;
w(3,:)=(1+2*r+q)/8;
w(2,:)=(3-2*r+q)/8;
w(1,:)=(3-2*r-q)/8;
%phi2
function w=phi2(r)
w=zeros(4,4);
q=sqrt(1+4r(1-r));
w(:,4)=(1+2*r-q)/8;
w(:,3)=(1+2*r+q)/8;
w(:,2)=(3-2*r+q)/8;
w(:,1)=(3-2*r-q)/8;
%sk
function f=sk(u,g)
global ip im h;
f=((u(ip,:,1)+u(:,:,1)).*g(ip,:)…
-(u(im,:,1)+u(:,:,1)).*g(im,:)…
+(u(:,ip,2)+u(:,:,2)).*g(:,ip)…
-(u(:,im,2)+u(:,:,2)).g(:,im))/(4h);
%skew
function w=skew(u)
w=u;
w(:,:,1)=sk(u,u(:,:,1));
w(:,:,2)=sk(u,u(:,:,2));
%spread
function f=spread(F,X)
global h N dtheta Nb;
c=dtheta/(h*h);
f=zeros(N,N,2);
for k=1:Nb
s=X(k,:)/h;
i=floor(s);
r=s-i;
i1=mod((i(1)-1):(i(1)+2),N)+1;
i2=mod((i(2)-1):(i(2)+2),N)+1;
w=phi1(r(1)).*phi2(r(2));
f(i1,i2,1)=f(i1,i2,1)+(cF(k,1))w;
f(i1,i2,2)=f(i1,i2,2)+(cF(k,2))w;
end