Velocity field in two vertical density layers(left=3200,right=3300).-Toy模板网

这篇具有很好参考价值的文章主要介绍了Velocity field in two vertical density layers(left=3200,right=3300).。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

%set parameter
Nx = 31; %ordinary
Ny = 21; 
L1 = 1000; 
L2 = 1500;
g_y = 10;
constant_viscosity = 10.^21;%unit:(Pa*s)
density_left = 3200;
density_right = 3300;
deltax = L1 / (Nx-1);
deltay = L2 / (Ny-1);
x_middle = L1/2;
Kcont = 2*constant_viscosity./(deltax+deltay);
%grids has Ny+1 rows and Nx+1 columns
[x,y] = meshgrid(0:deltax:L1,0:deltay:L2);
%unknown_number is (Ny+1)*(Nx+1)*3
unknown_number = (Nx+1)*(Ny+1)*3;
coe = sparse(unknown_number,unknown_number);
R = zeros(unknown_number, 1);
in_node = zeros(Ny+1,Nx+1);
in_vx = zeros(Ny+1,Nx+1);
in_vy = zeros(Ny+1,Nx+1);
in_p = zeros(Ny+1,Nx+1);
%boundary condition : free slip = -1 ; no slip = 1//why
bc_left=-1;
bc_right=-1;
bc_top=-1;
bc_bottom=-1;
density = zeros(Ny,Nx);
%two vertical layers
for j = 1:Nx
    for i = 1:Ny
        if (j-1)*deltax < x_middle
            density(i,j) = density_left;
        else
            density(i,j) = density_right;
        end
    end
end
%index
for i = 1:Ny+1
    for j = 1:Nx+1
        in_node(i,j) = (j-1)*(Ny+1)+i;
        in_vx(i,j) = 3*in_node(i,j)-2;
        in_vy(i,j) = 3*in_node(i,j)-1;
        in_p(i,j) = 3*in_node(i,j);
    end
end

for i = 1:Ny+1
    for j = 1:Nx+1
        %condition of Vx
        if(i==1 || i==Ny1 || j==1 || j==Nx || j==Nx1)
            %boundary condition
            coe(in_vx(i,j),in_vx(i,j)) = 1;% vx(i,j) = 0
            %right hand side
            R(in_vx(i,j),1) = 0;
            %top free slip
            if(i==1 && j>1 && j<Nx) % vx(i,j)-vx(i+1,j) = 0
                coe(in_vx(i,j),in_vx(i+1,j))=bctop; 
            end
            %bottom free slip
            if(i==Ny+1 && j>1 && j<Nx)% vx(i,j)-vx(i-1,j) = 0
                coe(in_vx(i,j),in_vx(i-1,j))=bcbottom; 
            end
        %x-Stokes equation    
        else
            coe(in_vx(i,j),in_vx(i,j)) = -2*constant_viscosity.*(1/(deltax.^2)+1/(deltay.^2));%vx3
            coe(in_vx(i,j),in_vx(i,j-1)) = constant_viscosity./(deltax.^2);%vx1
            coe(in_vx(i,j),in_vx(i,j+1)) = constant_viscosity./(deltax.^2);%vx5
            coe(in_vx(i,j),in_vx(i-1,j)) = constant_viscosity./(deltay.^2);%vx2
            coe(in_vx(i,j),in_vx(i+1,j)) = constant_viscosity./(deltay.^2);%vx4
            coe(in_vx(i,j),in_p(i,j+1)) = -Kcont./deltax;%P2'
            coe(in_vx(i,j),in_p(i,j)) = Kcont./deltax;%P1'
            %right hand side    
            R(in_vx(i,j),1)=0;    
        end
        %condition of Vy
        if(j==1 || j==Nx+1 || i==1 || i==Ny || i==Ny+1)
            %boundary condition
            coe(in_vy(i,j),in_vy(i,j)) = 1;% vy(i,j) = 0
            %right hand side
            R(in_vy(i,j),1) = 0;
            %left free slip
            if(j==1 && i>1 && i<Ny) % vy(i,j)-vy(i,j+1) = 0
                coe(in_vy(i,j),in_vy(i,j+1))=bcleft; 
            end
            %right free slip
            if(j==Nx+1 && i>1 && i<Ny)% vx(i,j)-vx(i,j-1) = 0
                coe(in_vy(i,j),in_vx(i,j-1))=bcright; 
            end
        %y-Stokes equation    
        else
            coe(in_vy(i,j),in_vy(i,j)) = -2*constant_viscosity.*(1/(deltax.^2)+1/(deltay.^2));%vy3
            coe(in_vy(i,j),in_vy(i,j-1)) = constant_viscosity./(deltax.^2);%vy1
            coe(in_vy(i,j),in_vy(i,j+1)) = constant_viscosity./(deltax.^2);%vy5
            coe(in_vy(i,j),in_vy(i-1,j)) = constant_viscosity./(deltay.^2);%vy2
            coe(in_vy(i,j),in_vy(i+1,j)) = constant_viscosity./(deltay.^2);%vy4
            coe(in_vy(i,j),in_p(i+1,j)) = -Kcont./deltay;%P2'
            coe(in_vy(i,j),in_p(i,j)) = Kcont./deltay;%P1'
            %right hand side    
            R(in_vy(i,j),1) = -g_y.*( density(i,j-1)+density(i,j) )./2;     
        end
        %condition of P
        if(i==1 || j==1 || i==Ny+1 || j==Nx+1 || (i==2 && j==2))        
            % Boundary Condition
            coe(in_p(i,j),in_p(i,j))=1; %p(i,j) = 0 
            % right hand side
            R(in_p(i,j),1)=0; 
            if(i==2 && j==2)
                coe(in_p(i,j),in_p(i,j))=1*Kcont;% p(2,2) = 0,p=p'*Kcont
                R(in_p(i,j),1)=0; 
            end
        %continuity equation    
        else
            coe(in_p(i,j),in_vx(i,j-1)) = -Kcont./deltax;%vx1
            coe(in_p(i,j),in_vx(i,j))   =  Kcont./deltax;%vx2
            coe(in_p(i,j),in_vy(i-1,j)) = -Kcont./deltay;%vy1
            coe(in_p(i,j),in_vy(i,j))   =  Kcont./deltay;%vy2
            % Right hand side
            R(in_p(i,j),1)=0;
        end
    end
end

%compute the solution u1:include vx,vy and p
u1 = coe \ R

vx_matrix = zeros(Ny+1,Nx+1);
vy_matrix = zeros(Ny+1,Nx+1);
p_matrix = zeros(Ny+1,Nx+1);
%decompose the solution of global matrix.
for j = 1:Nx+1
    for i = 1:Ny+1
        vx_matrix(i,j) = u1(in_vx(i,j),1);
        vy_matrix(i,j) = u1(in_vy(i,j),1);
        p_matrix(i,j) = u1(in_p(i,j),1);
    end
end
% Velocity field in two vertical density layers(left=3200,right=3300)
figure(1)
scale = 1
quiver(x,y,vx_matrix(1:Ny,1:Nx),vy_matrix(1:Ny,1:Nx),scale,color = 'r')
xlabel('Horizontal(Km)')
ylabel('Vertical(Km)')
title('Velocity field in two vertical density layers(left=3200,right=3300)')
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')

% Velocity_x field in two vertical density layers(left=3200,right=3300)
figure(2)
pcolor(x,y,vx_matrix(1:Ny,1:Nx))
xlabel('Horizontal(Km)')
ylabel('Vertical(Km)')
title('Velocity_x field in two vertical density layers(left=3200,right=3300)')
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
shading interp;
% Velocity_y field in two vertical density layers(left=3200,right=3300)
figure(3)
pcolor(x,y,vy_matrix(1:Ny,1:Nx))
xlabel('Horizontal(Km)')
ylabel('Vertical(Km)')
title('Velocity_y field in two vertical density layers(left=3200,right=3300)')
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
shading interp;
% Pressure field in two vertical density layers(left=3200,right=3300)
figure(4)
pcolor(x(2:Ny,2:Nx),y(2:Ny,2:Nx),p_matrix(2:Ny,2:Nx))
xlabel('Horizontal(Km)')
ylabel('Vertical(Km)')
title('Pressure field in two vertical density layers(left=3200,right=3300)')
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
colorbar;
shading interp;
% Distribution of two vertical density layers(left=3200,right=3300)
figure(5)
pcolor(x,y,density)
xlabel('Horizontal(Km)')
ylabel('Vertical(Km)')
title('Distribution of two vertical density layers(left=3200,right=3300)')
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
colorbar;
shading interp;

Velocity field in two vertical density layers(left=3200,right=3300).,matlab,velocity field,numerical,mantle