- #1
SebastianF
- 3
- 1
Homework Statement
Hi guys,
we are supposed to program the velocity and temperature fields for the 2-dim Navier-Stokes, using the pressure correction method and Boussinesq. approx for the temperature field. The relevant set of equations and B.C. is displayed below:
Homework Equations
<br /> \frac{\partial u}{\partial x}+\frac{\partial v}{\partial x}=0<br />
<br /> \frac{\partial u}{\partial t}+\frac{\partial P}{\partial x}=-\frac{\partial u^2}{\partial x}-\frac{\partial uv}{\partial y}+Pr\Big( \frac{\partial^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2}\Big)<br />
<br /> \frac{\partial v}{\partial t}+\frac{\partial P}{\partial y}=-\frac{\partial v^2}{\partial y}-\frac{\partial uv}{\partial x}+Pr\Big( \frac{\partial^2 v}{\partial x^2} +\frac{\partial^2 v}{\partial y^2}\Big)+f<br />
<br /> \frac{\partial T}{\partial t}=-\frac{\partial vT}{\partial y}-\frac{\partial uT}{\partial x}+Pr\Big( \frac{\partial^2 T}{\partial x^2} +\frac{\partial^2 T}{\partial y^2}\Big)<br />
where f = Ra * Pr * T
with Pr=Prandtl-number and Ra=Rayleigh Number
Boundary Conditions are:
x=0 : u=v=0 dT/dx=0
x=L_x: u=v=0 dT/dx=0
y=0 : u=v=0 T=T_bottom=0
y=L_y: u=v=0 T=T_top=0
Initial conditions are:
u=v=0
T=T_bottom+(y/L_y)*(T_top-T_bottom)
The Attempt at a Solution
I attempted to implement the above in Matlab:
The code runs but does not result in the expected vortices, but gives a strange velocity source at (1,1), I'm unsure where my error lies. Anyone have an idea?
Main Code:
Code:
% Notes for the Rayleigh-Bernard problem
% 1. Approach
% Careful: Omega from the exercise sheet is a temperature difference, not
% the actual temperature
% Try to imitate the behaviour of V and see if that works for omega
lid_driven_cavity=0;
if (lid_driven_cavity==1)
% Parameters for test case I: Lid-driven cavity
% The Richardson number is zero, i.e. passive scalar.
Pr = 0.71; % Prandtl number
Re = 10.; % Reynolds number
Ri = 0.; % Richardson number
dt = 0.00275; % time step
Tf = 20; % final time
Lx = 1; % width of box
Ly = 1; % height of box
Nx = 30; % number of cells in x
Ny = 30; % number of cells in y
ig = 200; % number of iterations between output
% Boundary and initial conditions:
Utop = 1.;
% IF TEMPERATURE: Tbottom = 1.; Ttop = 0.;
% IF TEMPERATURE: namp = 0.;
else
% Parameters for test case II: Rayleigh-Bénard convection
% The DNS will be stable for Ra=1705, and unstable for Ra=1715
% (Theoretical limit for pure sinusoidal waves
% with L=2.01h: Ra=1708)
% Note the alternative scaling for convection problems.
Pr = 0.71; % Prandtl number
Ra = 2715; % Rayleigh number
Re = 1./Pr; % Reynolds number
Ri = Ra*Pr; % Richardson number
dt = 0.0005; % time step
Tf = 20; % final time
Lx = 10.; % width of box
Ly = 1; % height of box
Nx = 200; % number of cells in x
Ny = 20; % number of cells in y
ig = 200; % number of iterations between output
% Boundary and initial conditions:
Utop = 0.;
% IF TEMPERATURE: Tbottom = 1.; Ttop = 0.;
namp = 0.1;
Tbottom = 1.;
Ttop = 0.;
end%------------------------------------------
% Number of iterations
Nit = ceil(Tf/dt);
% Spatial grid: Location of corners -> position of the vel. nodes
x = linspace(0,Lx,Nx+1);
y = linspace(0,Ly,Ny+1);
% Grid spacing
hx = Lx/Nx;
hy = Ly/Ny;
% Boundary conditions:
uN = x*0+Utop; vN = avg(x,2)*0;
uS = x*0; vS = avg(x,2)*0;
uW = avg(y,2)*0; vW = y*0;
uE = avg(y,2)*0; vE = y*0;
% all row vectors for now
% actual computation of Boundary condtions at Ue, Ve% Initial conditions
% Remember that i=1...x values end up in one colum U(i,j)=U(row,colum) and
% all j=1...y values in one row V(i,j)=V(row,colum)
U = zeros(Nx-1,Ny);
V = zeros(Nx,Ny-1);
% linear profile for T with random noise
% IF TEMPERATURE: T = ... + namp*rand(Nx,Ny)
T = ones(Nx,Ny)* Tbottom + (Ttop - Tbottom)*repmat(avg(y,2),Nx,1) + namp*rand(Nx,Ny);
% Omega = diff(T,1,1);
% T needs wider Boundary conditions vector: Nx, Ny
% Note: Although the spacing is incorrect this is ultimately irrelavant. xT
% will be used to save (B.C.) temperatures, not positions!
xT = linspace(0,Lx,Nx);
yT = linspace(0,Ly,Ny);
% temperature boundary conditions:
tN = xT*0+Ttop; tS = xT*0+Tbottom;
% tW = T(1,:); tE = T(end,:);
% Time series
tser = [];
Tser = [];
%-----------------------------------------
% Compute system matrices for pressure
% First set homogeneous Neumann condition all around
% Laplace operator on cell centres: Fxx + Fyy
Lp = kron(speye(Ny), DD(Nx,hx) ) + kron( DD(Ny,hy) ,speye(Nx));
% Set one Dirichlet value to fix pressure in that point
Lp(1,:) = 0 ; Lp(1,1) = 1 ;
% rhs = 0 inside the loop; further down
% For more speed, you could pre-compute the LU decomposition
% [LLp,ULp] = lu(Lp);
%-----------------------------------------
% Progress bar
fprintf(...
'[ | | | | ]\n')
%-----------------------------------------
% Main loop over iterations
for k = 1:Nit
% include all boundary points for u and v (linear extrapolation
% for ghost cells) into extended array (Ue,Ve)
% add Sidevectors
Ue = [uW; U; uE;];
% rightmost column (translates to bottommost row in actual matrix)
% choosen, so that the average value implemented in Ua will yield uN
Ue = [(2*uS')-Ue(:,1) Ue (2*uN')-Ue(:,end)];
Ve = [vS' V vN'];
Ve = [(2*vW)-Ve(1,:); Ve; (2*vE)-Ve(end,:);];
% boundary conditions for the Temperature
% Add column Side vectors
Te = [(2*tS')-T(:,1) T (2*tN')-T(:,end)];
% Add row top & bottom vectors.
Te = [Te(1,:); Te; Te(end,:);];
% averaged (Ua,Va) of u and v on corners
% Dimensionaltity of Ua,Va should be (nx+1)*(ny+1) for both
% (nx+1)*(ny+2) => (nx+1)*(ny+1)
% averaging in y-direction (Up/down)
Ua = avg(Ue,2);
% (nx+2)*(ny+1) => (nx+1)*(ny+1)
% averaging in x-direction (left/right)
Va = avg(Ve,1);
% averaging projects on cell corners
% Example RayleighBenardConvection used average value, so stay parallel
% and see what happens; Note that alone averaging will return:
% (nx+2)*(ny+2) => (nx+1)*(ny+2) => projection on V
Ta = avg(Te,1);
% (nx+1)*(ny+2) => (nx+1)*(ny+1) => projection on U
Ta = avg(Ta,2);
% Note: taking the difference with the lower boundary value
% (y-direction) should happen at some point -> Omega!
% construct individual parts of nonlinear terms
dUVdx = diff(Ua.*Va, 1, 1)/hx;
dUVdy = diff(Ua.*Va, 1, 2)/hy;
% diff in x-direction (up/down) in T-Matrix
% leaving v-distro
dUOdx = diff(Ua.*Ta, 1, 1)/hx;
% diff in y-direction (left/right) in T-Matrix
% leaving u-distro
dVOdy = diff(Va.*Ta, 1, 2)/hy;
dU2dx = diff(avg(Ue(:,2:end-1).^2,1), 1, 1)/hx;
dV2dy = diff(avg(Ve(2:end-1,:).^2,2), 1, 2)/hy;
% treat viscosity explicitly
viscu = diff( Ue(:,2:end-1),2,1 )/hx^2 + diff( Ue(2:end-1,:),2,2 )/hy^2;
viscv = diff( Ve(:,2:end-1),2,1 )/hx^2 + diff( Ve(2:end-1,:),2,2 )/hy^2;
% (nx+2)*(ny+2)=>(nx)*(ny) + (nx+2)*(ny+2)=>(nx)*(ny)
viscT = diff( Te(:,2:end-1),2,1 )/hx^2 + diff( Te(2:end-1,:),2,2 )/hy^2;
% buoyancy term
% IF TEMPERATURE: fy = ...
% (nx)*(ny)
% dUOdx(:,2:end-1) or dUOdx(:,1:end-1) or dUOdx(:,2:end)
T = T + dt*viscT - dt*(dUOdx(:,1:end-1)+dVOdy(1:end-1,:));
% T = T + viscT - (dUOdx(:,2:end)+dVOdy(2:end,:));
% driving force is the temperature DIFFERENCE:
fy = Ri*diff(T,1,2);
% compose final nonlinear term + explicit viscous terms
U = U + dt/Re*viscu - dt*( dU2dx + dUVdy(2:end-1,:));
% V = V + dt/Re*viscv - dt*( dV2dy + dUVdx(:,2:end-1)); % + % IF TEMPERATURE: dt*fy;
V = V + dt/Re*viscv - dt*( dV2dy + dUVdx(:,2:end-1)) + dt*fy;
% pressure correction, Dirichlet P=0 at (1,1)
rhs = (diff([uW; U; uE],1,1)/hx + diff([vS' V vN'],1,2)/hy)/dt;
rhs = reshape(rhs,Nx*Ny,1);
rhs(1) = 0;
P = Lp\rhs;
% alternatively, you can use the pre-computed LU decompositon
% P = ULp\(LLp\rhs);
P = reshape(P,Nx,Ny);
% apply pressure correction
U = U - dt*diff(P,1,1)/hx;
V = V - dt*diff(P,1,2)/hy;
% Temperature equation: Recomputation with pressure corrected velocities!
% is enforcement of B.C. necessary at this step?
% boundary conditions for the Temperature
%Add column Side vectors
Te = [(2*tS')-T(:,1) T (2*tN')-T(:,end)];
% Add row top & bottom vectors.
Te = [Te(1,:); Te; Te(end,:);];
% IF TEMPERATURE: Tu = ...
Ta = avg(Te,1);
Ta = avg(Ta,2);
dUOdx = diff(Ua.*Ta, 1, 1)/hx;
% IF TEMPERATURE: Tv = ...
dVOdy = diff(Va.*Ta, 1, 2)/hy;
% IF TEMPERATURE: H = ...
viscT = diff( Te(:,2:end-1),2,1 )/hx^2 + diff( Te(2:end-1,:),2,2 )/hy^2;
% IF TEMPERATURE: T = T + dt*H;
T = T + dt*viscT - dt*(dUOdx(:,2:end)+dVOdy(2:end,:));
%-----------------------------------------
% progress bar
if floor(51*k/Nit)>floor(51*(k-1)/Nit), fprintf('.'), end
% plot solution if needed
if k==1|floor(k/ig)==k/ig
% compute divergence on cell centres
if (1==1)
div = diff([uW;U;uE])/hx + diff([vS' V vN'],1,2)/hy;
figure(1);clf; hold on;
contourf(avg(x,2),avg(y,2),div');colorbar
axis equal; axis([0 Lx 0 Ly]);
title(sprintf('divergence at t=%g',k*dt))
drawnow
end
% compute velocity on cell corners
% renew Ua, Va matrixes with freshly computed values of U and V
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,2);
Va = avg(Ve,1);
% Ua = ...
% Va = ...
Len = sqrt(Ua.^2+Va.^2+eps);
figure(2);clf;hold on;
%contourf(avg(x,2),avg(y,2),P');colorbar
contourf(x,y,sqrt(Ua.^2+Va.^2)',20,'k-');colorbar
quiver(x,y,(Ua./Len)',(Va./Len)',.4,'k-')
axis equal; axis([0 Lx 0 Ly]);
title(sprintf('u at t=%g',k*dt))
drawnow
% IF TEMPERATURE: % compute temperature on cell corners
% IF TEMPERATURE: Ta = ...
% boundary conditions for the Temperature
% Add column Side vectors
% Te = [(2*tS')-T(:,1) T (2*tN')-T(:,end)];
% Add row top & bottom vectors.
% Te = [Te(1,:); Te; Te(end,:);];
Ta = avg(Te,1);
Ta = avg(Ta,2);
figure(3); clf; hold on;
contourf(x,y,Ta',20,'k-');colorbar
quiver(x,y,(Ua./Len)',(Va./Len)',.4,'k-')
axis equal; axis([0 Lx 0 Ly]);
title(sprintf('T at t=%g',k*dt))
drawnow
% Time history
if (1==1)
figure(4); hold on;
tser = [tser k*dt];
Tser = [Tser Ue(ceil((Nx+1)/2),ceil((Ny+1)/2))];
plot(tser,abs(Tser))
title(sprintf('Probe signal at x=%g, y=%g',...
x(ceil((Nx+1)/2)),y(ceil((Ny+1)/2))))
set(gca,'yscale','log')
xlabel('time t');ylabel('u(t)')
end
end
end
fprintf('\n')Function: avg
Code:
function B = avg(A,idim)
% AVG(A,idim)
%
% Averaging function to go from cell centres (pressure nodes)
% to cell corners (velocity nodes) and vice versa.
% avg acts on index idim; default is idim=1.
%
% This function belongs to SG2212.m
% nargin returns number of function arguments -> idim not entered
if nargin<2, idim = 1; end
if (idim==1)
B = [ A(2:end,:) + A(1:end-1,:) ]/2;
elseif (idim==2)
B = [ A(:,2:end) + A(:,1:end-1) ]/2;
else
error('avg(A,idim): idim must be 1 or 2')
endFunction: DD
Code:
function A = DD(n,h)
% DD(n,h)
%
% One-dimensional finite-difference derivative matrix
% of size n times n for second derivative:
% h^2 * f''(x_j) = -f(x_j-1) + 2*f(x_j) - f(x_j+1)
%
% Homogeneous Neumann boundary conditions on the boundaries
% are imposed, i.e.
% f(x_0) = f(x_1)
% if the wall lies between x_0 and x_1. This gives then
% h^2 * f''(x_j) = + f(x_0) - 2*f(x_1) + f(x_2)
% = + f(x_1) - 2*f(x_1) + f(x_2)
% = f(x_1) + f(x_2)
%
% For n=5 and h=1 the following result is obtained:
%
% A =
%
% -1 1 0 0 0
% 1 -2 1 0 0
% 0 1 -2 1 0
% 0 0 1 -2 1
% 0 0 0 1 -1
%
% This function belongs to SG2212.m
% produce vector, filled with ones, for further use in A:
% Note that while diag will contain n numbers and therefore will
% be larger than required for the n-1 length sidediagonals, spdiags
% will simply cut-off the superflous information.
diag = ones(n,1);
A = spdiags([diag -2*diag diag], -1:1, n, n)/h^2;
A(1,1) = -1/h^2;
A(n,n) = -1/h^2;Attachment contains the MATLAB m-files of the code fragments above and PDF-pictures of the Temp and velocity fields.