Desktop.zip_时间离散

function LU=Lu_Roe_cylindrical(U,N,dx,dy) %Fjordholm的ERoe格式求解二维潜水波方程，计算数值通量 U=boundary_cylindrical(U,N); g=1; for i=1:N+2 for j=1:N+2 h_L=U{i,j}(1); u_L=U{i,j}(2)/U{i,j}(1); v_L=U{i,j}(3)/U{i,j}(1); h_R=U{i+1,j+1}(1); u_R=U{i+1,j+1}(2)/U{i+1,j+1}(1); v_R=U{i+1,j+1}(3)/U{i+1,j+1}(1); %平均 h_ave=(h_L+h_R)/2; u_ave=(u_R*sqrt(h_R)+u_L*sqrt(h_L))/(sqrt(h_R)+sqrt(h_L)); v_ave=(v_R*sqrt(h_R)+v_L*sqrt(h_L))/(sqrt(h_R)+sqrt(h_L)); c_ave=sqrt((g*(h_R+h_L))/2); %关于通量导的特征向量 r_cap{1}=[1,u_ave-0.5*c_ave,v_ave-0.5*3^(1/2)*c_ave]'; r_cap{2}=[0,-0.5*3^(1/2)*c_ave,0.5*c_ave]'; r_cap{3}=[1,u_ave+0.5*c_ave,v_ave+0.5*3^(1/2)*c_ave]'; %关于通量导的特征向量的逆 l_cap{1}=[(v_ave*3^(1/2)+2*c_ave+u_ave)/(4*c_ave),(-v_ave+u_ave*3^(1/2))/(2*c_ave),(v_ave*3^(1/2)-2*c_ave+u_ave)/(4*c_ave)]'; l_cap{2}=[ -1/(4*c_ave),-3^(1/2)/(2*c_ave),1/(4*c_ave)]'; l_cap{3}=[-3^(1/2)/(4*c_ave),1/(2*c_ave),3^(1/2)/(4*c_ave)]'; %写成特征向量组成的矩阵 R=[r_cap{1},r_cap{2},r_cap{3}]; L=[l_cap{1},l_cap{2},l_cap{3}]; %对角矩阵 A=[abs(u_ave*cos(pi/3)+v_ave*sin(pi/3)-c_ave),0,0;0,abs(u_ave*cos(pi/3)+v_ave*sin(pi/3)),0;0,0,abs(u_ave*cos(pi/3)+v_ave*sin(pi/3)+c_ave)]; H1=[U{i,j}(2)*cos(pi/3)+U{i,j}(3)*sin(pi/3),(U{i,j}(2)^2/U{i,j}(1)+g/2*U{i,j}(1)^2)*cos(pi/3)+(U{i,j}(2)*U{i,j}(3)/U{i,j}(1))*sin(pi/3),(U{i,j}(2)*U{i,j}(3)/U{i,j}(1))*cos(pi/3)+(U{i,j}(3)^2/U{i,j}(1)+g/2*U{i,j}(1)^2)*sin(pi/3)]'; H2=[U{i+1,j+1}(2)*cos(pi/3)+U{i+1,j+1}(3)*sin(pi/3),(U{i+1,j+1}(2)^2/U{i+1,j+1}(1)+g/2*U{i+1,j+1}(1)^2)*cos(pi/3)+(U{i+1,j+1}(2)*U{i+1,j+1}(3)/U{i+1,j+1}(1))*sin(pi/3),(U{i+1,j+1}(2)*U{i+1,j+1}(3)/U{i+1,j+1}(1))*cos(pi/3)+(U{i+1,j+1}(3)^2/U{i+1,j+1}(1)+g/2*U{i+1,j+1}(1)^2)*sin(pi/3)]'; H{i,j}=0.5.*(H2-H1)-0.5*R*A*L*(U{i+1,j+1}-U{i,j}); end end for j=1:N+1 for i=1:N+1 LU{i,j}=-(1/4)*(H{i+1,j+1}-H{i,j}); end end