韩国小清新网站模板/鞍山seo外包

%%%一个兄弟找的，还没来得及调试。先分享出来~有需要研究这个的朋友大家一起讨论啊~

%%%最近发现了个好东西叫MATDEM，还没来得及详细研究。http://www.matdem.com/list/?94_1.html，开源软件，可做后

%%%处理。

%%%还是老样子，有啥问题Feel free to tell me, 我会的一定解答！  QQ群：293267908。

global dt Nb N h rho mu ip im a;
global kp km dtheta K;
initialize
init_a

for clock=1:clockmax
  XX=X+(dt/2)*interp(u,X);
  ff=spread(Force(XX),XX);
  [u,uu]=fluid(u,ff);
  X=X+dt*interp(uu,XX);
  
  %animation:
  vorticity=(u(ip,:,2)-u(im,:,2)-u(:,ip,1)+u(:,im,1))/(2*h);
  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
end

%%%2019/6/25 22:34:20%%%

function [uuu,uu]=fluid(u,ff)
global a dt rho mu;
w=u-(dt/2)*skew(u)+(dt/(2*rho))*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-dt*skew(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));

%%% 2019/6/25 22:41:20%%%

% ib2D.m
% This script is the main program.

global dt Nb N h rho mu ip im a;
global kp km dtheta K;
initialize
init_a

for clock=1:clockmax
  XX=X+(dt/2)*interp(u,X);
  ff=spread(Force(XX),XX);
  [u,uu]=fluid(u,ff);
  X=X+dt*interp(uu,XX);
  
  %animation:
  vorticity=(u(ip,:,2)-u(im,:,2)-u(:,ip,1)+u(:,im,1))/(2*h);
  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
end

%%% 2019/6/25 22:41:33%%%

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

%%% 2019/6/25 22:41:43%%%

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)+(c*F(k,1))*w;
  f(i1,i2,2)=f(i1,i2,2)+(c*F(k,2))*w;
end