# Computational fluid dynamics problem (Matlab)

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

May 5, 2013, 02:53
Computational fluid dynamics problem (Matlab)
#1
New Member

srikanth
Join Date: Mar 2013
Posts: 3
Rep Power: 5
Hi, :)

I was using the MIT18086_NAVIERSTOKES.m code and encountered an issue which I have discussed below. I modified the MIT18086_NAVIERSTOKES.m code to suit my problem. But I'm unable to achieve the correct result.

A rectangular obstacle 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.

The fluid inflow is from the left upper side/boundary of the domain. The fluid outflow is at the bottom left boundary of the domain. Till 0.2 seconds the fluid flow direction is perfect. But after that, the flow direction is incorrect (reverse). How to overcome this problem?

Please be kind enough to help me out to sort this issue.

The code is given below (Main program and 3 functions): :)
--------------------------------------------------------------------------
function [Unew,Vnew] = Navierstokes(U,uN,uS,uE,uW,V,vN,vS,vW,vE,nx,ny)
% 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 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 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 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;

global Re Pr 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 = 15; % number of cells for inflow
out = 15; % number of cells for outflow
obl = 1; obr = 30; 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
U = zeros(nx-1,ny); V = zeros(nx,ny-1);
%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+0.02]; vW = y*0;
uE = ay*0; vE = y*0;
%--------------------------------------------------------------------------
U1=U(:,obt+1:end); V1=V(:,obt+1:ny-1);
U2=U(obr+1:nx-1,obb:obt); V2=V(obr+1:nx,obb:obt-1);
U3=U(:,1:obb-1); V3=V(:,1:obb-2);
nx1=nx; nx2=nx-obr; nx3=nx;
ny1=obb-1; ny2=obt-obb+1; ny3=ny-obt;
%--------------------------------------------------------------------------
fprintf('initialization')
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
%--------------------------------------------------------------------------
uN1=uN;
uS1=uS; uS1(obr+2:end-1)=U(obr+1:end,obt+1)';
uW1=uW(obt+1:end); uE1=uE(obt+1:end);

vN1=vN;
vS1=vS; vS1(obr+1:end)=V(obr+1:end,obt+1)';
vW1=vW(obt+1:end); vE1=vE(obt+1:end);
[U1,V1] = Navierstokes(U1,uN1,uS1,uE1,uW1,V1,vN1,vS1,vW1,vE1 ,nx1,ny1);
%--------------------------------------------------------------------------
uN2=uS1(obr+1:end);
uS2=uS(obr+1:end); uS2(2:end-1)=U(obr+1:end,obb)';
uW2=uW(obb:obt); uE2=uE(obb:obt);

vN2=vS1(obr+1:end);
vS2=V(obr+1:end,obb)';%
vW2=vW(obb:obt+1); vE2=vE(obb:obt+1);
[U2,V2] = Navierstokes(U2,uN2,uS2,uE2,uW2,V2,vN2,vS2,vW2,vE2 ,nx2,ny2);
%--------------------------------------------------------------------------
uN3=uN; uN3(obr+1:end)=uS2;
uS3=uS; uS3(2:out)=U(1:out-1,1)';
uW3=uW(1:obb-1); uE3=uE(1:obb-1);

vN3=vN; vN3(obr+1:end)=vS2;
vS3=vS; vS3(1:out)=V(1:out,1)';
vW3=vW(1:obb); vE3=vE(1:obb);
[U3,V3] = Navierstokes(U3,uN3,uS3,uE3,uW3,V3,vN3,vS3,vW3,vE3 ,nx3,ny3);
%--------------------------------------------------------------------------
U(:,obt+1:ny)=U1; V(:,obt+1:end)=V1;
U(obr+1:end,obb:obt)=U2; V(obr+1:end,obb-1:obt)=[vS2' V2 vN2'];
U(:,1:obb-1)=U3; V(:,1:obb-2)=V3;

%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
clf,hold on
Ue = [uS3' avg([uW;U;uE]')' uN'];
Ve = [vW;avg([vS3' V vN']);vE];
Ue(obl:obr+1,obb:obt+1)=Ue(obl:obr+1,obb:obt+1)*0;
Ve(obl:obr+1,obb:obt+1)=Ve(obl:obr+1,obb:obt+1)*0;
Len = sqrt(Ue.^2+Ve.^2+eps);
quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'k-')
hold off,axis equal, axis([0 lx 0 ly])
title(sprintf('Re = %0.1g ,Pr = %0.1g ,t = %0.2g',Re,Pr,k*dt))
drawnow
end
end
fprintf('\n')
Attached Images
 CS3.jpg (20.9 KB, 9 views)

 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 srik Main CFD Forum 3 February 3, 2015 11:56 volo87 CFX 5 June 14, 2013 17:44 Jacek Stecki Main CFD Forum 0 November 10, 2002 06:49 Abhi Main CFD Forum 12 July 8, 2002 09:11 Jonas Larsson Main CFD Forum 16 August 7, 1998 16:27

All times are GMT -4. The time now is 06:56.