Computational fluid dynamics problem (Matlab)

 Register Blogs Members List Search Today's Posts Mark Forums Read

March 25, 2013, 13:05
Computational fluid dynamics problem (Matlab)
#1
New Member

srikanth
Join Date: Mar 2013
Posts: 3
Rep Power: 5
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.

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);
Attached Images
 untitled2.jpg (11.0 KB, 29 views)

 March 26, 2013, 08:49 #2 Senior Member   Jonas T. Holdeman, Jr. Join Date: Mar 2009 Location: Knoxville, Tennessee Posts: 108 Rep Power: 10 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: 5 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: 61 Rep Power: 6 Hello, I 'd like to know the program that create a square or rectangular obstacle in matlab??? Thank you

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post volo87 CFX 5 June 14, 2013 17:44 Abhi Main CFD Forum 12 July 8, 2002 09:11 Amy Moseley Main CFD Forum 10 July 1, 1999 08:46 T Barron Main CFD Forum 5 December 1, 1998 08:43 Jonas Larsson Main CFD Forum 16 August 7, 1998 16:27

All times are GMT -4. The time now is 18:46.