www.gusucode.com > 动态二维液体模拟仿真源码程序 > 动态二维液体模拟仿真源码程序/2dliquidsimulation/liquid_2d_rectangle.m

    a=3;
b=4;
n=50;

R=sqrt(2*a*b/(n*3*sqrt(3)));

tm=100;
dt=0.02;

m=1; % mass

vis=0.2; % viscosity
visw=0.2; % viscosity from wall

cmp=1; % compressability

rectangle('position',[0 0 a b]);
hold on;

xst=sqrt(3)*R;
yst=3*R/2;

r=[];
ylc=1;
for y=yst/4:yst:(b-yst/8)
    if mod(ylc,2)==0
        x1=xst/2;
    else
        x1=xst;
    end
    x1=x1-xst/4;
    for x=x1:xst:(a-xst/8)
        r=[r, [x;y]];
    end
    ylc=ylc+1;
end
% r=[a*rand(1,n);
%    b*rand(1,n)];
szr=size(r);
n=szr(2);

v=zeros(size(r)); % velocities
%v(2,1)=1;

% rotation
vr0=0.5;
dr=zeros(size(r));
dr(1,:)=r(1,:)-a/2;
dr(2,:)=r(2,:)-b/2;
for nc=1:n
    dr2=dr(:,nc)'*dr(:,nc);
    drl=sqrt(dr2);
    drn=dr(:,nc)/drl;
    drp=[-drn(2);drn(1)];
    v(:,nc)=vr0*drl*drp;
end

pres=ones(1,n); % pressures

rr=r;
rr=[rr [2*a-r(1,:); r(2,:)]]; % left
rr=[rr [-r(1,:); r(2,:)]]; % right
rr=[rr [r(1,:); 2*b-r(2,:) ]]; % up
rr=[rr [r(1,:); -r(2,:) ]]; % down

hp=plot(r(1,:),r(2,:),'.k');
%voronoi(rr(1,:),rr(2,:));
axis equal;

[V,C] = voronoin(rr');
sV=size(V);
VL=sV(1);

As0=zeros(1,n); % cell areas
for nc=1:n
    CC=C{nc};
    As0(nc)=polyarea(V(CC,1),V(CC,2));
end
As=zeros(1,n); % current cell areas

% r(1,1)=r(1,1)+a/20;

% closes point to 0.75*a 0.5*b
drdr=zeros(size(r));
drdr(1,:)=r(1,:)-0.75*a;
drdr(2,:)=r(2,:)-0.5*b;
drdrl=zeros(1,n);
for nc=1:n
    drdrl(nc)=sqrt(drdr(:,nc)'*drdr(:,nc));
end

[tmp mi]=min(drdrl);


for t=0:dt:tm
    
    
    
    rr=r;
    rr=[rr [2*a-r(1,:); r(2,:)]]; % left
    rr=[rr [-r(1,:); r(2,:)]]; % right
    rr=[rr [r(1,:); 2*b-r(2,:) ]]; % up
    rr=[rr [r(1,:); -r(2,:) ]]; % down
    
    [V,C] = voronoin(rr');
    
    cla;
    hp=plot(r(1,:),r(2,:),'.k');
    hold on;
    hold on;
    voronoi(rr(1,:),rr(2,:));
    fill([a 2*a 2*a a],[0 0 b b],'w','Edgecolor','none');
    fill([0 -a -a 0],[0 0 b b],'w','Edgecolor','none');
    rectangle('position',[0 0 a b]);
    xlim([0 a]);
    ylim([0 b]);
    
    axis equal;
    drawnow;
    
    E=zeros(5,10*n); % edges info
    EL=0; % number of edges
    % E(1,i) small number of vertex
    % E(2,i) big number of vertex
    % E(3,i) number of domain
    % E(4,i) number of another domain
    % E(5,i) number of domains
    for nc=1:n
        CC=C{nc};
        As(nc)=polyarea(V(CC,1),V(CC,2));
        LCC=length(CC);
        for p=1:LCC % for each edge
            if p~=LCC
                p1=p;
                p2=p+1;
            else
                p1=p;
                p2=1;
            end
            % 
            % sort CC(p1)  CC(p2):
            if CC(p1)<CC(p2)
                CC1=CC(p1);
                CC2=CC(p2);
            else
                CC2=CC(p1);
                CC1=CC(p2);
            end
            % check all previouse edges:
            ind=find((E(1,1:EL)==CC1)&(E(2,1:EL)==CC2));
            if isempty(ind)
                % add new adge:
                EL=EL+1;
                E(1,EL)=CC1;
                E(2,EL)=CC2;
                E(3,EL)=nc;
                E(5,EL)=1;
            else
                % add info to exist edge:
                E(4,ind(1))=nc;
                E(5,ind(1))=2;
            end
        end
    end
    
    % find pressures
    for nc=1:n
        dAr=(As(nc)-As0(nc))/As0(nc);
        pres(nc)=1-cmp*dAr;
    end
        
    % find forcies:
    F=zeros(size(r));
        
    for nc=1:n
        CC=C{nc};
        LCC=length(CC);
        P=0; % cell perimeter
        F1s=zeros(2,1); % wighted sum forces from neighbours
        F2=zeros(2,1); % friction force from wall if any
        F3=zeros(2,1); % pressures
        
        
        for p=1:LCC % for each edge
            if p~=LCC
                p1=p;
                p2=p+1;
            else
                p1=p;
                p2=1;
            end
            
                        
            % 
            % sort CC(p1)  CC(p2):
            if CC(p1)<CC(p2)
                CC1=CC(p1);
                CC2=CC(p2);
            else
                CC2=CC(p1);
                CC1=CC(p2);
            end
            
            % find edge length:
            dv=V(CC1,:)-V(CC2,:);
            dv2=dv*dv';
            dvl=sqrt(dv2);
            nn=dv'/dvl; % along
            pp=[-nn(2);nn(1)]; % perpendicular
            dvm=(V(CC1,:)+V(CC2,:))'/2; % midddle of edge
            insd=r(:,nc)-dvm; % vector directed from edge to center
            % orient pp strictlly inside cell:
            if (pp'*insd)<0
                pp=-pp;
            end
            P=P+dvl; % add to perimeter
            
            % find this edge:
            ind=find((E(1,1:EL)==CC1)&(E(2,1:EL)==CC2));
            if E(5,ind)==1
                % contact with wall
                F2=-visw*dvl*v(:,nc);
                F3=F3+pres(nc)*dvl*pp; % wall reaction
            else
                % E(5,ind)==2 contact with neighbour
                nnc=find(E(3:4,ind)~=nc);
                nc1=E(2+nnc,ind); % connected domain
                dvv=v(:,nc1)-v(:,nc);  % speed difference
                F1s=F1s+dvl*dvv;
                
                % pressure
                F3=F3+pres(nc1)*dvl*pp;
            end
        end
        %F1=vis*F1s/P; % friction force from neighbours
        F1=vis*F1s; % friction force from neighbours
        
        F(:,nc)=F1+F2+F3; % neighbours+walls+pressures
    end
    
    % integration:
    aa=F/m;
    v=v+aa*dt;
    r=r+v*dt;
    
    i1=find(r(1,:)<0);
    r(1,i1)=0;
    
    i1=find(r(1,:)>a);
    r(1,i1)=a;
    
    i1=find(r(2,:)<0);
    r(2,i1)=0;
    
    i1=find(r(2,:)>b);
    r(2,i1)=b;
    
%     w=2*pi*1*t;
%     r(1,mi)=a/2+(a/4)*cos(w);
%     r(2,mi)=b/2+(b/4)*sin(w);
   
    %set(hp,'Xdata',r(1,:),'YData',r(2,:));
    %drawnow;
end