Как построить двухмерную проблему адвекции при использовании сетки для определения значения скорости? - PullRequest
Новая
       21

Как построить двухмерную проблему адвекции при использовании сетки для определения значения скорости?

0 голосов
/ 13 марта 2019

В настоящее время я пытаюсь решить уравнение 2D-адвекции 2-го порядка, используя схему против ветра.Сначала задача состоит в том, чтобы построить график quiver(), а затем переместить его поверх contourf().При использовании данных для скорости u и v в схеме против ветра я получаю прямолинейные выходы, как показано ниже.Однако я использую начальное условие phi0 = cos(x).При проверке значений для X, Y, u, v все они имеют разные значения в (i, j) местах.Я вижу, что мои phi0 и phi остаются постоянными в каждом столбце, но выход адвекции должен быть разным в разное время.Я следовал своему 1D 2-му порядку, который отлично работает, но, похоже, не может получить откорректированный движущийся график.Любой совет по моей настройке или если вы можете указать, где я ошибаюсь с этим сюжетом, очень помог бы! 

clear all
clc

%Problem 1
%Part B
%Creating a quiver plot for the 2D vector profile
L=2*pi;%Length Lx=Ly = 2pi
L0=0;
N=31; % Nx=Ny=31
%get a value of dx=dy to know distance between steps 
dx=L/N;
dy=L/N
 x=L0-dx*2:dx:L+dx*2;
 y=L0-dy*2:dy:L+dy*2;

[X,Y]=meshgrid(x,y);

u=cos(X).*sin(Y);
v=-sin(X).*cos(Y);

figure
hold on


%Part B & C using the courtf plot Phi=cos(x)
phi0=cos(X);
contourf(X,Y,phi0)
colormap autumn
colorbar
xlabel('Length from 0 to 2*pi with dx spacing')
ylabel('Length from 0 to 2*pi with dy spacing')
title('Quiver plot on the phi=cos(x) Contour Plot')
quiver(X,Y,u,v)
hold off

figure(2)
plot(x,phi0)
%%
%Writing a code to sovle 2D advection for  part D
t=0; %initial time
tmax=10; %Maximum time
dt=0.01; %time step

phi=phi0;
phip1=phi0;

%phi(:)=phi; %Initial Condition
nsteps = tmax/dt
%Add periodic boundary conditions for both x & y direction
for n=1 : nsteps
    phi(1,:) = phi(end-2,:);
    phi(2,:) = phi(end-3,:);
    phi(end,:) = phi(3,:);
    phi(end-1,:) = phi(4,:);
    %Y ghost cells
    phi(:,1) = phi(:,end-2);
    phi(:,1) = phi(:,end-3);
    phi(:,end) = phi(:,3);
    phi(:,end-1) = phi(:,4);

    for i=3:N+1
        for j=3:N+1

          if u > 0 & v>0
               phip1(i,j)= phi(i,j) - u(i)*dt/(2*dx) * (3*phi(i,j)-4*phi(i-1,j)+phi(i-2,j))- v(j)*dt/(2*dx) *(3*phi(i,j)-4*phi(i,j-1)+phi(i,j-2))
          elseif u <0 & v<0
               phip1(i,j)= phi(i,j) - u(i)*dt/(2*dx) * (-3*phi(i,j)+4*phi(i+1,j)-phi(i+2,j))- v(j)*dt/(2*dx) *(-3*phi(i,j)+4*phi(i,j+1)-phi(i,j+2))
          elseif u >0 & v <0
               phip1(i,j)= phi(i,j) - u(i)*dt/(2*dx) * (3*phi(i,j)-4*phi(i+1,j)+phi(i+2,j))- v(j)*dt/(2*dx) *(-3*phi(i,j)+4*phi(i,j-1)-phi(i,j-2))
          elseif u <0 & v >0
               phip1(i,j)= phi(i,j) - u(i)*dt/(2*dx) * (-3*phi(i,j)+4*phi(i-1,j)-phi(i-2,j))- v(j)*dt/(2*dx) *(3*phi(i,j)-4*phi(i,j+1)+phi(i,j+2))
          end
        end
    end
    t=t+dt;
    phi=phip1;


    plot(x,phi)
    %pause(0.5)
end

1 Ответ

1 голос
/ 13 марта 2019

Проблема, которую я вижу, заключается в том, что вы строите phi (который в 2D) против x (который является 1D).

Я не уверен на 100% в правом разделе, который вы хотите построитьно что-то вроде этого должно работать: plot(x,phi0(1,:)).Для этого нужно построить первый срез phi в направлении y.

enter image description here

EDIT

Чтобы визуализировать phi0 как функцию X и Y, вы можете использовать surf(X,Y,phi0) или mesh(X,Y,phi0).

...