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