|
[Sponsors] |
Computational fluid dynamics problem (Matlab) |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
March 25, 2013, 13:05 |
Computational fluid dynamics problem (Matlab)
|
#1 |
New Member
srikanth
Join Date: Mar 2013
Posts: 3
Rep Power: 13 |
Hi,
While using the MIT18086_NAVIERSTOKES.m code I encountered a problem to which I need some guidance from to come up with a solution. I modified the MIT18086_NAVIERSTOKES.m code to suit my problem. But I'm unable to achieve the correct required result when the code is executed. I have explained my problem below. A rectangular obstacle (a heat source) is placed in a rectangular domain of incompressible fluid flow. The rectangular obstacle is placed where its left boundary lying on the left boundary of the rectangular domain (obstacle moved to the left side). The top or bottom boundaries of the obstacle do not touch the rectangular domain boundary. This can be seen in the image attached. I want to apply finite difference method to this problem. But I can't figure out the boundary conditions for the obstacle for the velocity components, temperature and also the outflow temperature condition. I have attached the modified code as well. Please be kind enough to help me out to sort this issue. Thank you in advance. The code is given below: ------------------------------------------------------------------------ Main code: global Re Pr nx ny hx hy dt out %==============================MAIN PROGRAM================================ Re = 50; % Reynolds number Pr = 7; % Prandtl number dt = 1e-2; % time step tf = 8e-0; % final time lx = 1; % width of box ly = 1; % height of box nx = 40; % number of x-gridpoints ny = 40; % number of y-gridpoints nsteps = 10; % number of steps with graphic output in = 10; % number of cells for inflow out = 12; % number of cells for outflow RT = 30; % Room temperature obl = 1; obr = 24; obb = 19; obt = 22; %defining the obstacle %-------------------------------------------------------------------------- close all nt = ceil(tf/dt); dt = tf/nt; x = linspace(0,lx,nx+1); hx = lx/nx; y = linspace(0,ly,ny+1); hy = ly/ny; ax = avg(x); ay = avg(y); %-------------------------------------------------------------------------- %Initial conditions U0 = zeros(nx-1,ny); V0 = zeros(nx,ny-1); %T = zeros(nx,ny)+RT; %Boundary conditions uN = x*0; vN = ax*0; uS = x*0; vS = ax*0; uW = [ay(1:end-in)*0 ay(end-in+1:end)*0+1]; vW = y*0; uE = ay*0; vE = y*0; %-------------------------------------------------------------------------- fprintf('initialization') Lq = kron(speye(ny-1),K1(nx-1,hx,2))+kron(K1(ny-1,hy,2),speye(nx-1)); perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq'; figure('un', 'normalized', 'pos',[0.01 0.05 0.98 0.85]) fprintf(', time loop\n--20%%--40%%--60%%--80%%-100%%\n') for k=1:nt vS(1:out) = V0(1:out,1)'; %outflow B.C uS(2:out+1) = U0(1:out,1)'; %outflow B.C %calling Energy function %T = Energy(U,uW,uE,V,vN,vS,T,TN,TS,TW,TE); %calling Navierstokes function [U,V] = Navierstokes(U0,uN,uS,uE,uW,V0,vN,vS,vW,vE); V(obl:obr,obt) = V(obl:obr,obt)*0; V(obl:obr,obb-1) = V(obl:obr,obb-1)*0; U(obr,obb:obt) = U(obr,obb:obt)*0; % U(obl-1,obb:obt) = U(obl-1,obb:obt)*0; U(obl:obr-1,obb) = -U(obl:obr-1,obb-1); U(obl:obr-1,obt) = -U(obl:obr-1,obt+1); % V(obl,obb:obt-1) = -V(obl-1,obb:obt-1); V(obr,obb:obt-1) = -V(obr+1,obb:obt-1); U(obl:obr-1,obb:obt) = U(obl:obr-1,obb:obt)*0; V(obl:obr,obb:obt-1) = V(obl:obr,obb:obt-1)*0; %visualization if floor(25*k/nt)>floor(25*(k-1)/nt), fprintf('.'), end if k==1 | floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt) %stream function rhs = reshape(diff(U')'/hy-diff(V)/hx,[],1); q(perq) = Rq\(Rqt\rhs(perq)); Q = zeros(nx+1,ny+1); Q(2:end-1,2:end-1) = reshape(q,nx-1,ny-1); Q(obl:obr+1,obb:obt+1)=Q(obl:obr+1,obb:obt+1)*0; clf, %contourf(avg(x),avg(y),P',20,'w-'), hold on contour(x,y,Q',20,'k-'); Ue = [uS' avg([uW;U;uE]')' uN']; Ve = [vW;avg([vS' V vN']);vE]; Len = sqrt(Ue.^2+Ve.^2+eps); quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'b-') % streamline((Ue./Len)',(Ve./Len)',2,4) hold off, axis equal, axis([0 lx 0 ly]) % p = sort(p); caxis(p([8 end-7])) title(sprintf('Re = %0.1g ,Pr = %0.1g ,t = %0.2g',Re,Pr,k*dt)) drawnow end U0=U; V0=V; end fprintf('\n') ---------------------------------------------------------------------- function 1: function [Unew,Vnew,P] = Navierstokes(U,uN,uS,uE,uW,V,vN,vS,vW,vE) % function NAVIERSTOKES % Solves the incompressible Navier-Stokes equations in a % rectangular domain with prescribed velocities along the % boundary. The solution method is finite differencing on % a staggered grid with impicit diffusion and a Chorin % projection method for the projection. %---------------------------------------------------------------------- global Re nx ny hx hy dt Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hx^2+... [uW;zeros(nx-3,ny);uE]/hy^2); Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hx^2+... [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hy^2); Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye (nx)); Lp(1,1) = 3/2*Lp(1,1); perp = symamd(Lp);Rp = chol(Lp(perp,perp)); Rpt = Rp'; Lu = speye((nx-1)*ny)+dt/Re*(kron(speye(ny),K1(nx-1,hx,2))+... kron(K1(ny,hy,3),speye(nx-1))); peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut =Ru'; Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),K1(nx,hx,3))+... kron(K1(ny-1,hy,2),speye(nx))); perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv'; % treat nonlinear terms gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1); Ue = [uW;U;uE]; Ue = [2*uS'-Ue(:,1) Ue 2*uN'-Ue(:,end)]; Ve = [vS' V vN']; Ve = [2*vW-Ve(1,:);Ve;2*vE-Ve(end,:)]; Ua = avg(Ue')'; Ud = diff(Ue')'/2; Va = avg(Ve); Vd = diff(Ve)/2; UVx = diff(Ua.*Va-gamma*abs(Ua).*Vd)/hx; UVy = diff((Ua.*Va-gamma*Ud.*abs(Va))')'/hy; Ua = avg(Ue(:,2:end-1)); Ud = diff(Ue(:,2:end-1))/2; Va = avg(Ve(2:end-1,:)')'; Vd = diff(Ve(2:end-1,:)')'/2; U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/hx; V2y = diff((Va.^2-gamma*abs(Va).*Vd)')'/hy; U = U-dt*(UVy(2:end-1,:)+U2x); V = V-dt*(UVx(:,2:end-1)+V2y); %implicit viscosity rhs = reshape(U+Ubc,[],1); u(peru) = Ru\(Rut\rhs(peru)); U = reshape(u,nx-1,ny); rhs = reshape(V+Vbc,[],1); v(perv) = Rv\(Rvt\rhs(perv)); V = reshape(v,nx,ny-1); %pressure correction rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1); p(perp) = -Rp\(Rpt\rhs(perp)); P = reshape(p,nx,ny); U = U-diff(P)/hx; V = V-diff(P')'/hy; Unew = U; Vnew = V; ------------------------------------------------------------------- function 2: function A = K1(n,h,a11) %a11: Neumann=1, Dirichlet=2, Dirichlet mid=3 A = spdiags([-1 a11 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)/h^2; -------------------------------------------------------------------- function 3: function B = avg(A) if size(A,1)==1, A = A'; end B = (A(2:end,:)+A(1:end-1,:))/2; if size(A,2)==1, B = B'; end ---------------------------------------------------------------------- function 4: function Tnew = Energy(U,uW,uE,V,vN,vS,T,TN,TS,TW,TE) %computes the temperature at each discretized point according %to the velocity components using Energy equation %-------------------------------------------------------------------------- global Re Pr nx ny hx hy dt LT = speye(nx*ny)-dt/(Re*Pr)*(kron(speye(ny),K1(nx,hx,2))+... kron(K1(ny,hy,2),speye(nx))); gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1); Ue = [uW;U;uE]; Ve =[vS' V vN']; Te = [2*TW-T(1,:); T; 2*TE-T(end,:)]; Ta = avg(Te); Td = diff(Te)/2; UTx = diff(Ue.*Ta-gamma*abs(Ue).*Td)/hx; Te = [2*TS'-T(:,1) T 2*TN'-T(:,end)]; Ta = avg(Te')'; Td = diff(Te')'/2; VTy = diff((Ve.*Ta-gamma*abs(Ve).*Td)')'/hy; T = reshape(T,[],1); Tx = reshape(UTx,[],1); Ty = reshape(VTy,[],1); T1 = LT*T-(Tx+Ty)*dt; Tnew = reshape(T1,nx,ny); |
|
March 26, 2013, 08:49 |
|
#2 |
Senior Member
Jonas T. Holdeman, Jr.
Join Date: Mar 2009
Location: Knoxville, Tennessee
Posts: 128
Rep Power: 18 |
An observation: a lot of simulations use Pr=0.7, I haven't stopped to figure why, I have used it too, but I notice you set Pr=7.
|
|
March 26, 2013, 09:32 |
|
#3 |
New Member
srikanth
Join Date: Mar 2013
Posts: 3
Rep Power: 13 |
I assumed that the fluid is water. And the Pr (prandl number) for water is 7.
|
|
February 3, 2015, 11:56 |
how to create the obstacle geometry in matlab?
|
#4 |
Member
amine
Join Date: Jun 2012
Posts: 65
Rep Power: 14 |
Hello,
I 'd like to know the program that create a square or rectangular obstacle in matlab??? Thank you |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
An error has occurred in cfx5solve: | volo87 | CFX | 5 | June 14, 2013 18:44 |
Terrible Mistake In Fluid Dynamics History | Abhi | Main CFD Forum | 12 | July 8, 2002 10:11 |
computational fluid dynamics | Amy Moseley | Main CFD Forum | 10 | July 1, 1999 09:46 |
Computational Fluid Dynamics software | T Barron | Main CFD Forum | 5 | December 1, 1998 08:43 |
CFD - Trends and Perspectives | Jonas Larsson | Main CFD Forum | 16 | August 7, 1998 17:27 |