asidesi

3

Apr 20th, 2021
678
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. clear all;
  2. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  3. % Set:
  4. % (0) decide whether restar (restart=1) or scratch (restart=0)
  5. % (1) the number of points (number of pressure cells)
  6. % (2) the dimension of the domain
  7. % (3) the Reynolds number and other geometrical characteristics
  8. % (4) (for the moment) a fixed time step
  9. %      and the number of time steps
  10. % (5) the boundary values and rkpm yes/no
  11. % (6) check convergence every iconv time steps
  12. %     eventually save the field and
  13. %     decide whether or not open "on the fly" figures      
  14. %                                                                  
  15. %               Ub(3) and Vb(3)                                      
  16. %           ¡---------------------¡                                  
  17. %   Ub(4)   ¡                     ¡ Ub(2) and Vb(2)                    
  18. %   Vb(4)   ¡                     ¡                                    
  19. %           -----------------------                                      
  20. %               Ub(1) and Vb(1)                                          
  21. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  22. % (0)
  23. restart=0 % decide if from scratch (0) or from former solution (1)
  24. % (1)
  25. ncx=794; % number of nodes for the pressure in x
  26. ncy=547; % number of nodes for the pressure in y
  27. % (2)
  28. lx=36.28884826325411;
  29. ly=25.0; % y-wise extension of the domain
  30. % (3)
  31. Re=100; % Reynolds number of the immersed circle (U upstream * Diameter/kinematic visc.)
  32. D=1;   %Diameter of the cylinder;
  33. rc=D/2;% circle radius
  34. xc=5.461608775137110;  % center of circle (coordinates)
  35. yc=12.5;
  36. % (4)
  37. nt=2000; CFL=0.1; % number of time steps and CFL number
  38. [Ub,Vb]=set_velocity_bc(ncx,ncy,lx,ly); % set velocity on the boundaries
  39. % decide if rkpm or not
  40. irkpm=1 % (1) yes rkpm
  41. uhlm=0 % (1) yes uhlmann
  42. iconv=20; % every iconv steps gives some info on the screen
  43. snapshot_it=100; % every snapshot save the current solution
  44. nsteps=[1:snapshot_it:nt];
  45. % assemble the grid and the initial solution
  46. if restart==0
  47.    [xu,xv,xp,yu,yv,yp,u,v,p]=initialize(ncx,ncy,lx,ly,Ub,Vb); %initialize the field and grid if from scratch
  48.    save xp.mat xp;
  49.    save yp.mat yp;
  50.    save xu.mat xu;
  51.    save yu.mat yu;
  52.    save xv.mat xv;
  53.    save yv.mat yv;
  54. else % otherwise load previous solution and grid
  55.    load uu; load vv; load pp;
  56.    load xp; load yp; load xu; load xv; load yu; load yv;
  57. end
  58. % initialize the matrices for Helmholz and Poisson Problems
  59. nxu=length(xu); nyu=length(yu); %no. of nodes for u (x-comp. of velocity) % CAREFUL STAGGERED GRID
  60. nxv=length(xv); nyv=length(yv); %no. of nodes for v (y-comp. of velocity)
  61. nxp=length(xp); nyp=length(yp); %no. of nodes for p
  62. % no uniform mesh with lx assumed not to be equal to ly!!!!
  63. hx=lx/ncx; hy=ly/ncy; %mesh spacing
  64. Deltat=CFL*hx/max(Ub(4,:)); % Allowed time step using entrance velocity as max
  65. disp('Delta t ')
  66. disp(Deltat)
  67. beta=Deltat/hx; % some parameters
  68. alpha=Re/Deltat;
  69. %h=xp(2)-xp(1);
  70. % prefactor some parts of approximate factorization of viscous terms
  71. % We use D'Yakonov scheme (approximate factorization)
  72. tic
  73. [alxu,bexu,cxu,alyu,beyu,cyu]=predyakonovU(nxu,nyu,hx,hy,Deltat,Re);
  74. [alxv,bexv,cxv,alyv,beyv,cyv]=predyakonovV(nxv,nyv,hx,hy,Deltat,Re);
  75. [PX,PY,PMIX]=assmatFD(ncx,ncy,hx,hy); % diagonalize pressure laplacian operator
  76. toc
  77. disp('finished eig problems ')
  78. [xlag,ylag,theta,dS]=initialize_lag(xu,yu,rc,xc,yc,irkpm,0); % discretize circle with uniformly spaced Lagrange nodes
  79. if irkpm ==1
  80. tic
  81. [epsilonx,bcoefx,ix_x,jy_x,area_x,nsup_x,hx_x,hy_x]=precompute_eps(xu(2:length(xu)-1),yu(2:length(yu)-1),xlag,ylag,dS);
  82. [epsilony,bcoefy,ix_y,jy_y,area_y,nsup_y,hx_y,hy_y]=precompute_eps(xv(2:length(xv)-1),yv(2:length(yv)-1),xlag,ylag,dS);
  83. toc
  84. end
  85. if uhlm==1
  86. [epsilonx,bcoefx,ix_x,jy_x,area_x,nsup_x,hx_x,hy_x]=precompute_eps_norkpm(xu(2:length(xu)-1),yu(2:length(yu)-1),xlag,ylag,dS);
  87. [epsilony,bcoefy,ix_y,jy_y,area_y,nsup_y,hx_y,hy_y]=precompute_eps_norkpm(xv(2:length(xv)-1),yv(2:length(yv)-1),xlag,ylag,dS);
  88. end
  89. disp('finished precomputing immersed boundary data')
  90. %%%% start the time loop
  91. disp('starting time advancement')
  92. isnapshot=0;
  93. Qin=sum(u(nxu,2:nyu-1)); % flow rate at inlet
  94. AB1=1; AB2=0; CGP=0; % first step Euler for nonlinear terms, from the second one Adams-Bashfort
  95. laplax=zeros(nxu-2,nyu-2); NLx=zeros(nxu-2,nyu-2);
  96. laplay=zeros(nxv-2,nyv-2); NLy=zeros(nxv-2,nyv-2);
  97. % initialise snapshot matrix
  98. Number_snap=100;
  99. Asnap=zeros((ncx-1)*(ncy-1),Number_snap);
  100. Snapshot_count=0;
  101. for it=1:nt
  102. disp('step number'); disp(it)
  103. tic
  104.     actime=(it-1)*Deltat;
  105.     if mod(it,iconv) == 0
  106.        uold=u; vold=v;
  107.     end
  108.     if mod(it,10*iconv) == 0
  109.       AB1=1; AB2=0;
  110.     end
  111.  
  112. %%%%%%%% X_MOM EQUATION %%%%%%%%%%%%%%
  113. % FROM A to B prediction of u from x momentum equation without immersed object (fully explicit)
  114. % --------- A --------------
  115. % compute the laplacian of u
  116.    lapla=laplacian(u,hx,hy);
  117. % compute NL terms of X-mom equation and dp/dx
  118.    NL=nonlinear(u,v,Ub,Vb,lx,ly,ncx,ncy,1);
  119.    GP=pressgrad(p,hx,hy,1);
  120. % compute the force for the x-mom equation
  121.    rhs=sumupAB(u,lapla,NL,laplax,NLx,CGP*GP,Deltat,Re,AB1,AB2); %(2.10a)
  122.    rhs=immersed_boundary(rhs,xu(2:length(xu)-1),yu(2:length(yu)-1),xlag,ylag,dS,epsilonx,bcoefx,...
  123.        ix_x,jy_x,area_x,nsup_x,hx_x,hy_x);
  124. % --------- B --------------
  125. % rhs contains the spread force field required to impose bc's on the circle for u
  126.    forcex=hx*hy*sum(sum(rhs)); % compute integral force in x direction on the circle
  127.    Cd=2*forcex/D; %...and the drag coeficient
  128. % assemble rhs for the implicit step
  129. % from C to D redo x-momentum with the spreaded body force on the rhs
  130. % Approximate fact. for viscous terms and Adams Bashfort for the non-linear ones
  131. %---------C ------------
  132.    rhs=rhs-AB1*NL+AB2*NLx-CGP*GP;
  133.    [Ubs]=setBCout(u,v,Deltat,hx,Qin,1);
  134.    ustar=DyakonovU(rhs,u,Ub,Ubs,alxu,bexu,cxu,alyu,beyu,cyu,hx,hy,Re,Deltat);
  135. %---------D ------------
  136.    laplax=lapla; NLx=NL;
  137. % DO STAGE A-B and C-D for the y momentum equation
  138. %%%%%%%% Y_MOM EQUATION %%%%%%%%%%%%%%
  139. % compute the laplacian of v
  140.    lapla=laplacian(v,hx,hy);
  141. % compute NL terms of Y-mom equation and dp/dy
  142.    NL=nonlinear(u,v,Ub,Vb,lx,ly,ncx,ncy,2);
  143.    GP=pressgrad(p,hx,hy,2);
  144. % compute the force for the y-mom equation
  145.    rhs=sumupAB(v,lapla,NL,laplay,NLy,CGP*GP,Deltat,Re,AB1,AB2);
  146.    rhs=immersed_boundary(rhs,xv(2:length(xv)-1),yv(2:length(yv)-1),xlag,ylag,dS,epsilony,bcoefy,...
  147.              ix_y,jy_y,area_y,nsup_y,hx_y,hy_y);
  148.     if mod(it,iconv) == 0
  149.        disp('saving force field')
  150.        save rhs.mat rhs;
  151.     end
  152.    forcey=hx*hy*sum(sum(rhs));
  153.    Cf=2*forcey/D;
  154. % assemble rhs for the implicit step
  155.    rhs=rhs-AB1*NL+AB2*NLy-CGP*GP;
  156.    [Vbs]=setBCout(u,v,Deltat,hx,Qin,2);
  157.    vstar=DyakonovV(rhs,v,Vb,Vbs,alxv,bexv,cxv,alyv,beyv,cyv,hx,hy,Re,Deltat);
  158.    laplay=lapla; NLy=NL;
  159. %%%%%% Get ready for the projection step
  160. % set the boundary values on the boundary frame
  161.    [u,v]=setbc(ustar,vstar,Ub,Vb,Ubs,Vbs);
  162. % prepare the rhs for the pressure equation
  163.     div=divergence(u,v,hx,hy,ncx,ncy);
  164. % solve the Poisson problem for pressure
  165.    phi=solveFD(PX,PY,PMIX,div/Deltat);
  166.    [u,v,p]=finalvelocityVK(phi,u,v,p,Deltat,Re,hx,hy);
  167.   if mod(it,iconv) == 0
  168.      disp('saving the field')
  169.      isnapshot=isnapshot+1;
  170.      dummy=screenout(uold,u,vold,v,actime,hx,hy,ncx,ncy,Cd,Cf);
  171.      forcexx(isnapshot)=forcex; forceyy(isnapshot)=forcey;
  172.      Cdd(isnapshot)=Cd; Cll(isnapshot)=Cf;
  173.      time(isnapshot)=actime;
  174.      dummy=savefield(u,v,p,phi,time,forcexx,forceyy,Cdd,Cll);
  175.      [omega,psi,xo,yo]=vort(u,v,hx,hy);
  176.      Snapshot_count=Snapshot_count+1;
  177.      Asnap(:,Snapshot_count)=reshape(omega(2:ncx,2:ncy),(ncx-1)*(ncy-1),1);
  178.   end
  179. toc
  180. disp('time for 1 time step')
  181. AB1=3/2; AB2=1/2; CGP=1; % Adam Bashfort coef.
  182. end % END of time looping
  183. save SnapMatrix.mat Asnap;
  184. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  185. %%% some post process
  186. [omega,psi,xo,yo]=vort(u,v,hx,hy);
  187. [X,Y]=meshgrid(xo(2:ncx),yo(2:ncy));
  188. %[omega,psi,xo,yo]=vort(u,v,hx,hy);
  189. fp=figure;
  190. ap=axes;
  191. maxval=max(max(psi));
  192. minval=min(min(psi));
  193. vals=linspace(minval,maxval,50);
  194. contour(X,Y,psi',vals);
  195. set(ap,'DataAspectRatio',[1 1 1],'XLim',[0 lx],'YLim',[0 ly]);
  196. title(sprintf(' Stream function  (%g...%g); Re_b=%g',...
  197.                 minval,maxval,Re));
  198. colorbar;
  199.  
  200. fo=figure;
  201. ao=axes;
  202. maxval=max(max(omega));
  203. minval=min(min(omega));
  204. vals=linspace(minval,maxval,100);
  205. contour(X,Y,omega(2:ncx,2:ncy)',vals);
  206. set(ao,'DataAspectRatio',[1 1 1],'XLim',[0 lx],'YLim',[0 ly]);
  207. title(sprintf(' Vorticity (%g...%g); Re_b=%g',...
  208.                 minval,maxval,Re));
  209. colorbar;
  210.  
  211. fx=figure;
  212. ax=axes;
  213. maxval=max(max(u));
  214. minval=min(min(u));
  215. vals=linspace(minval,maxval,50);
  216. contour(xu,yu,u',vals);
  217. set(ax,'DataAspectRatio',[1 1 1],'XLim',[0 lx],'YLim',[0 ly]);
  218. title(sprintf(' u-com velocity (%g...%g); Re_b=%g',...
  219.                 minval,maxval,Re));
  220. colorbar;
  221.  
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×