# 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
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);
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);
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!