- Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
- Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
- Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.
Trouble with array index in a cfd solver
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
Hello,
I’m trying to implement a simple cfd solver for the incompressible naiver stokes equation on rectangle domain (Lid-driven cavity)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1431818/image.png)
Boundary conditions are needed for the velocity field and the pressure Poisson equation. The velocity boundary conditions are a bit tricky since some of the velocity components are not defined on the boundary. For example, to specify u = utop at the top of the domain is not straightforward because u is defined dx/2 away from the top boundary. One solution to this problem is to create a fictitious velocity outside the domain such that the velocity on the domain boundary satisfies the boundary condition.
Unfortunately, I made a mistake, which is puzzling me.
Here is the script so far:
clc,clear all,close all
% Number of cells
nx=5;
ny=5;
% lengthof the domain
Lx=1;
Ly=1;
% Material data
rho=1;
nu=1;
% Time
dt=0.001;
tg=0.01;
kkk=tg/dt;
%CFL=u*dt/dx;
% BC
u_bot=0;
u_top=1;
v_lef=0;
v_rig=0;
% Index extents
imin=2;
imax=imin+nx-1;
jmin=2;
jmax=jmin+ny-1;
% Create mesh
x(imin:imax+1)=linspace(0,Lx,nx+1);
y(jmin:jmax+1)=linspace(0,Ly,ny+1);
xm(imin:imax)=0.5*(x(imin:imax)+x(imin+1:imax+1));
ym(jmin:jmax)=0.5*(y(jmin:jmax)+y(jmin+1:jmax+1));
% Create mesh sizes
dx=x(imin+1)-x(imin);
dy=y(jmin+1)-y(jmin);
dxi=1/dx;
dyi=1/dy;
u=zeros(imax,jmax);
v=zeros(imax,jmax);
for ijk=1:1:kkk
t=ijk*dt;
% Boundary conditions
u(:,jmin-1)=2*u_bot-u(:,jmin)
u(:,jmax+1)=2*u_top-u(:,jmax)
v(imin-1,:)=2*v_lef-v(imin,:)
v(imax+1,:)=2*v_rig-v(imax,:)
% u momentum discretization
for j=jmin:jmax
for i=imin+1:imax
v_here=0.25*(v(i-1,j)+v(i-1,j+1)+v(i,j)+v(i,j+1));
us(i,j)=u(i,j)+dt*...
(nu*(u(i-1,j)-2*u(i,j)+u(i+1,j))*dxi^2 ...
+nu*(u(i,j-1)-2*u(i,j)+u(i,j+1))*dyi^2 ...
-u(i,j)*(u(i+1,j)-u(i-1,j))*0.5*dxi...
-v_here*(u(i,j+1)-u(i,j-1))*0.5*dyi) ;
end
end
% v momentum discretization
for j=jmin+1: jmax
for i=imin : imax
u_here=0.25*(u(i,j-1)+u(i,j)+u(i+1,j-1)+u(i+1,j));
vs(i,j)=v(i,j)+dt* ...
(nu*(v(i-1,j)-2*v(i,j)+v(i+1,j))*dxi ^2 ...
+nu*(v(i,j-1)-2*v(i,j)+v(i,j+1))*dyi ^2 ...
-u_here*(v(i+1,j)-v(i-1,j))*0.5*dxi ...
-v(i,j)*(v(i,j+1)-v(i,j-1))*0.5*dyi) ;
end
end
% Create Laplacian operator for solving pressure Poisson equation
L=zeros(nx*ny,nx*ny);
for j=1:ny
for i=1:nx
L(i+(j-1)*nx,i+(j-1)*nx)=2*dxi^2+2*dyi^2;
for ii=i-1:2:i+1
if(ii>0 && ii<=nx) % Interior point
L(i+(j-1)*nx,ii+(j-1)*nx)=-dxi^2;
else % Neuman conditions on boundary
L(i+(j-1)*nx,i+(j-1)*nx)= ...
L(i+(j-1)*nx,i+(j-1)*nx)-dxi^2;
end
end
for jj=j-1:2:j+1
if(jj >0 && jj<=ny) % Interiorpoint
L(i+(j-1)*nx,i+(jj-1)*nx)=-dyi^2;
else % Neuman conditions on boundary
L(i+(j-1)*nx,i+(j-1)*nx)= ...
L(i+(j-1)*nx,i+(j-1)*nx)-dyi^2;
end
end
end
end
% Set pressure in first cell (all other pressures w.r.t to this one)
L(1,:)=0;L(1,1)=1;
% MATLAB code to compute the right-hand-side R
n=0;
for j=jmin:jmax
for i=imin:imax
n=n+1;
R(n)=-rho/dt* ...
((us(i+1,j)-us(i,j))*dxi ...
+(vs(i,j+1)-vs(i,j))*dyi);
end
end
% The pressure is found by solving Lpn+1=R
pv=LnR;
% Finally, pv is converted to the mesh representation p(i,j)
n=0;
p=zeros(imax,jmax);
for j=jmin:jmax
for i=imin:imax
n=n+1;
p(i,j)=pv(n);
end
end
% Corrector step
for j=jmin:jmax
for i=imin+1:imax
u(i,j)=us(i,j)-dt/rho*(p(i,j)-p(i-1,j))*dxi;
end
end
for j=jmin+1: jmax
for i=imin : imax
v(i,j)=vs(i,j)-dt/rho*(p(i,j)-p(i,j-1))*dyi;
end
end
end
Maybe someone has a clue how to solve this
Greetings
Steffen
2 Commenti
Steven Lord
il 11 Lug 2023
What is the mistake you've made that's puzzling you?
Risposte (0)
Vedere anche
Categorie
Scopri di più su Mathematics and Optimization in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!