clc;
n = 40; % number of elements
finalT = 0.2; % final time
%boundary = 1; % switch between Dirichlet conditions (1) and Absorbing (0)
k = 6; % polynomial degree
gamma=1.4; % 比热比
%Cr = 0.2; % Courant number -> sets time step size in dt = Cr * h / a
CFL=0.2;
%flux = 1; % flux type, 0 = LF, 1 = HDG
nc = k+1; % number of quadrature points
left = -0.5; % left end of the domain
right = 0.5; % right end of the domain
%plot_accurate = 1; % plot only on nodes (0) or with more resolution (1)
%% set quadrature formula and quadrature nodes for integration
[xg,wg] = legendre_gauss(nc);
% set the node points for the Lagrange polynomials
[xunit,w] = legendre_gauss_lobatto(k+1);
%% START THE SIMULATION
% create mesh y and compute mesh size h
kp1 = length(xunit);
y = zeros(2*n,1);
h = zeros(n,1);
for e=1:n
y(2*e-1) = left + (right-left)*(e-1)/n;
y(2*e) = left + (right-left)*e/n;
h(e) = y(2*e)-y(2*e-1);
end
dx=(right-left)/n;
% create the grid of all node points for interpolation
x=zeros(kp1,n);
for e=1:n
for i=1:kp1
x(i,e) = y(2*e-1)+(y(2*e)-y(2*e-1))*(0.5+0.5*xunit(i));
end
end
% compute time step from Cr number, adjust to hit the desired final time
dt=CFL*dx;
% evaluate reference cell polynomials
%[values,derivatives] = evaluate_lagrange_basis(xunit, pg);
%% MASS MATRIX
[psi,dpsi] = lagrange_basis3(k+1,nc,xunit, xg);
Me = psi * diag(wg) * psi';
m1=inv(Me);
% De = dpsi * diag(wg) * psi';
%M1inv = sparse(kp1*n,kp1*n);
%for e=1:n
% M1inv((kp1*e-k):kp1*e,(kp1*e-k):kp1*e) = inv(0.5*dx*Me);
% De((kp1*e-k):kp1*e,(kp1*e-k):kp1*e) = 0.5*dx*De;
%end
%M=M1inv;
% D=De;
%% 初始条件 - Sod激波管问题
rho = zeros(kp1,n);
u = zeros(kp1,n);
p = zeros(kp1,n);
for e=1:n
for j = 1:kp1
x_pos = x(j,e);
if x_pos < 0.0
rho(j,e) = 1.0;
u(j,e) = 0.0;
p(j,e) = 1.0;
else
rho(j,e) = 0.125;
u(j,e) = 0.0;
p(j,e) = 0.1;
end
end
end
%% 计算初始条件下的守恒变量
E = p/(gamma-1) + 0.5*rho.*u.^2; % 总能量
U = zeros(3,kp1,n);
U(1,:,:) = rho;
U(2,:,:) = rho.*u;
U(3,:,:) = E;
%% run time loop
t=0;
%step_count=0;
%output_interval=100; %输出时间间隔
while t<finalT
if t+dt>finalT
dt=finalT-t;
end
% Runge-Kutta 3
U1 = U+dt*rhs1(U,dx,gamma,xunit,psi,dpsi,wg);
U2 = 0.75*U + 0.25*(U1 + dt*rhs1(U1,dx,gamma,xunit,psi,dpsi,wg));
U = 1/3*U + 2/3*(U2 + dt*rhs1(U2,dx,gamma,xunit,psi,dpsi,wg));
t=t+dt;
end
%% 结果可视化
figure;
x_plot = reshape(x,kp1*n, 1);
rho_plot = reshape(U(1,:,:), kp1*n, 1);
u_plot = reshape(U(2,:,:)./U(1,:,:), kp1*n, 1);
p_plot = reshape((gamma-1)*(U(3,:,:)-0.5*U(2,:,:).^2./U(1,:,:)), kp1*n, 1);
% 绘制图形
subplot(3,1,1);
plot(x_plot, rho_plot, 'b-', 'LineWidth', 1.5);
title('密度分布');
xlabel('x'); ylabel('\rho');
grid on;
subplot(3,1,2);
plot(x_plot, u_plot, 'r-', 'LineWidth', 1.5);
title('速度分布');
xlabel('x'); ylabel('u');
grid on;
subplot(3,1,3);
plot(x_plot, p_plot, 'g-', 'LineWidth', 1.5);
title('压力分布');
xlabel('x'); ylabel('p');
grid on;
%%
%---------------------------------------------------------------------%
%This code computes the Legendre-Gauss points and weights
%which are the roots of the Legendre Polynomials.
%Written by F.X. Giraldo on 4/2008
% Department of Applied Mathematics
% Naval Postgraduate School
% Monterey, CA 93943-5216
%---------------------------------------------------------------------%
function [xgl,wgl] = legendre_gauss(P)
xgl=zeros(P,1);
wgl=zeros(P,1);
p=P-1; %Order of Polynomials (P is the number of points: p+1)
ph=floor( (p+1)/2 );
for i=1:ph
x=cos( (2*i-1)*pi/(2*p+1) );
for k=1:20
[L0,L0_1,L0_2]=legendre_poly(p+1,x); %Compute the N+1 Legendre Polys
%Get new Newton Iteration
dx=-L0/L0_1;
x=x+dx;
if (abs(dx) < 1.0e-20)
break
end
end
xgl(p+2-i)=x;
wgl(p+2-i)=2/( (1-x^2)*L0_1^2 );
end
%Check for Zero Root
if (p+1 ~= 2*ph)
x=0;
[L0,L0_1,L0_2]=legendre_poly(p+1,x);
xgl(ph+1)=x;
wgl(ph+1)=2/( (1-x^2)*L0_1^2 );
end
%Find remainder of roots via symmetry
for i=1:ph
xgl(i)=-xgl(p+2-i);
wgl(i)=+wgl(p+2-i);
end
end
%%
%---------------------------------------------------------------------%
%This code computes the Legendre-Gauss-Lobatto points and weights
%which are the roots of the Lobatto Polynomials.
%Written by F.X. Giraldo on 4/2000
% Department of Applied Mathematics
% Naval Postgraduate School
% Monterey, CA 93943-5216
%---------------------------------------------------------------------%
function [xgl,wgl] = legendre_gauss_lobatto(P)
xgl=zeros(P,1);
wgl=zeros(P,1);
p=P-1;
ph=floor( (p+1)/2 );
for i=1:ph
x=cos( (2*i-1)*pi/(2*p+1) );
for k=1:20
[L0,L0_1,L0_2]=legendre_poly(p,x);
%Get new Newton Iteration
dx=-(1-x^2)*L0_1/(-2*x*L0_1 + (1-x^2)*L0_2);
x=x+dx;
if (abs(dx) < 1.0e-20)
break
end
end
xgl(p+2-i)=x;
wgl(p+2-i)=2/(p*(p+1)*L0^2);
end
%Check for Zero Root
if (p+1 ~= 2*ph)
x=0;
[L0,L0_1,L0_2]=legendre_poly(p,x);
xgl(ph+1)=x;
wgl(ph+1)=2/(p*(p+1)*L0^2);
end
%Find remainder of roots via symmetry
for i=1:ph
xgl(i)=-xgl(p+2-i);
wgl(i)=+wgl(p+2-i);
end
end
%%
%---------------------------------------------------------------------%
%This code computes the Legendre Polynomials and its 1st and 2nd derivatives
%Written by F.X. Giraldo on 4/2000
% Department of Applied Mathematics
% Naval Postgraduate School
% Monterey, CA 93943-5216
%---------------------------------------------------------------------%
function [L0,L0_1,L0_2] = legendre_poly(p,x)
L1=0;L1_1=0;L1_2=0;
L0=1;L0_1=0;L0_2=0;
for i=1:p
L2=L1;L2_1=L1_1;L2_2=L1_2;
L1=L0;L1_1=L0_1;L1_2=L0_2;
a=(2*i-1)/i;
b=(i-1)/i;
L0=a*x*L1 - b*L2;
L0_1=a*(L1 + x*L1_1) - b*L2_1;
L0_2=a*(2*L1_1 + x*L1_2) - b*L2_2;
end
end
%%
%% 函数:计算右侧项
function rhs11=rhs1(U,dx,gamma,xunit,psi,dpsi,weights)
[~,kp1,n]=size(U);
rhs11=zeros(size(U));
rho=U(1,:,:);
u=U(2,:,:)./rho;
E=U(3,:,:);
p=(gamma-1)*(E-0.5*rho.*u.^2);
rhs11_1=zeros(kp1,n);
rhs11_2=zeros(kp1,n);
rhs11_3=zeros(kp1,n);
rho_1=zeros(kp1,1);
u_1 =zeros(kp1,1);
E_1 =zeros(kp1,1);
%% 计算刚度部分
de=0.5*dx*dpsi*diag(weights)*(psi)';
for e=1:n
rho_1(:,1)=rho(:,e);
u_1(:,1) =u(:,e);
E_1(:,1) =E(:,e);
rhs11_1(:,e)=de*rho_1;
rhs11_2(:,e)=de*u_1;
rhs11_3(:,e)=de*E_1;
rhs11(1,:,e)=rhs11_1(:,e);
rhs11(2,:,e)=rhs11_2(:,e);
rhs11(3,:,e)=rhs11_3(:,e);
end
%% 边界通量计算(数值通量)
for e = 1:n
%% 左边界
if e == 1
% 流入边界条件
rho_L = 1.0;
u_L = 0.0;
p_L = 1.0;
E_L = p_L/(gamma-1) + 0.5*rho_L*u_L^2;
U_L = [rho_L; rho_L*u_L; E_L];
% 单元内左端点状态
rho_R = U(1,1,e);
u_R = U(2,1,e)/rho_R;
E_R = U(3,1,e);
p_R = (gamma-1)*(E_R - 0.5*rho_R*u_R^2);
U_R = [rho_R; rho_R*u_R; E_R];
% 计算Roe平均
U_roe = roe_average(U_L, U_R, gamma);
rho_roe = U_roe(1);
u_roe = U_roe(2)/rho_roe;
E_roe = U_roe(3);
p_roe = (gamma-1)*(E_roe - 0.5*rho_roe*u_roe^2);
a_roe = sqrt(gamma*p_roe/rho_roe);
% 计算特征值
lambda1 = u_roe - a_roe;
lambda2 = u_roe;
lambda3 = u_roe + a_roe;
lambda_max = max(abs([lambda1, lambda2, lambda3]));
% 使用更强的Lax-Friedrichs通量
F_L = euler_flux(U_L, gamma);
F_R = euler_flux(U_R, gamma);
F_num = 0.5*(F_L + F_R - 1.5*lambda_max*(U_R - U_L)); % 增加系数提高耗散
% 应用数值通量
rhs11(:,1,e)=rhs11(:,1,e)-F_num;
else
% 内部边界 - 与左侧单元交互
% 当前单元左端点状态
rho_R = U(1,1,e);
u_R = U(2,1,e)/rho_R;
E_R = U(3,1,e);
p_R = (gamma-1)*(E_R - 0.5*rho_R*u_R^2);
U_R = [rho_R; rho_R*u_R; E_R];
% 左侧单元右端点状态
rho_L = U(1,kp1,e-1);
u_L = U(2,kp1,e-1)/rho_L;
E_L = U(3,kp1,e-1);
p_L = (gamma-1)*(E_L - 0.5*rho_L*u_L^2);
U_L = [rho_L; rho_L*u_L; E_L];
% 计算Roe平均
U_roe = roe_average(U_L, U_R, gamma);
rho_roe = U_roe(1);
u_roe = U_roe(2)/rho_roe;
E_roe = U_roe(3);
p_roe = (gamma-1)*(E_roe - 0.5*rho_roe*u_roe^2);
a_roe = sqrt(gamma*p_roe/rho_roe);
% 计算特征值
lambda1 = u_roe - a_roe;
lambda2 = u_roe;
lambda3 = u_roe + a_roe;
lambda_max = max(abs([lambda1, lambda2, lambda3]));
% 使用更强的Lax-Friedrichs通量
F_L = euler_flux(U_L, gamma);
F_R = euler_flux(U_R, gamma);
F_num = 0.5*(F_L + F_R - 1.5*lambda_max*(U_R - U_L)); % 增加系数提高耗散
% 应用数值通量
rhs11(:,1,e)=rhs11(:,1,e)-F_num;
end
%% 右边界
if e == n
% 改进的流出边界条件
rho_R = 0.125;
u_R = 0.0;
p_R = 0.1;
E_R = p_R/(gamma-1) + 0.5*rho_R*u_R^2;
U_R = [rho_R; rho_R*u_R; E_R];
% 使用前一个单元的状态作为外部状态
rho_L = U(1,kp1,e);
u_L = U(2,kp1,e)/rho_L;
E_L = U(3,kp1,e);
p_L = (gamma-1)*(E_L - 0.5*rho_L*u_L^2);
U_L = [rho_L; rho_L*u_L; E_L];
% 计算Roe平均
U_roe = roe_average(U_L, U_R, gamma);
rho_roe = U_roe(1);
u_roe = U_roe(2)/rho_roe;
E_roe = U_roe(3);
p_roe = (gamma-1)*(E_roe - 0.5*rho_roe*u_roe^2);
a_roe = sqrt(gamma*p_roe/rho_roe);
% 计算特征值
lambda1 = u_roe - a_roe;
lambda2 = u_roe;
lambda3 = u_roe + a_roe;
lambda_max = max(abs([lambda1, lambda2, lambda3]));
% 使用更强的Lax-Friedrichs通量
F_L = euler_flux(U_L, gamma);
F_R = euler_flux(U_R, gamma);
F_num = 0.5*(F_L + F_R - 1.5*lambda_max*(U_R - U_L)); % 增加系数提高耗散
% 应用数值通量
rhs11(:,kp1,e)=rhs11(:,kp1,e)-F_num;
else
% 内部边界 - 与右侧单元交互
% 当前单元右端点状态
rho_L = U(1,kp1,e);
u_L = U(2,kp1,e)/rho_L;
E_L = U(3,kp1,e);
p_L = (gamma-1)*(E_L - 0.5*rho_L*u_L^2);
U_L = [rho_L; rho_L*u_L; E_L];
% 右侧单元左端点状态
rho_R = U(1,1,e+1);
u_R = U(2,1,e+1)/rho_R;
E_R = U(3,1,e+1);
p_R = (gamma-1)*(E_R - 0.5*rho_R*u_R^2);
U_R = [rho_R; rho_R*u_R; E_R];
% 计算Roe平均
U_roe = roe_average(U_L, U_R, gamma);
rho_roe = U_roe(1);
u_roe = U_roe(2)/rho_roe;
E_roe = U_roe(3);
p_roe = (gamma-1)*(E_roe - 0.5*rho_roe*u_roe^2);
a_roe = sqrt(gamma*p_roe/rho_roe);
% 计算特征值
lambda1 = u_roe - a_roe;
lambda2 = u_roe;
lambda3 = u_roe + a_roe;
lambda_max = max(abs([lambda1, lambda2, lambda3]));
% 使用更强的Lax-Friedrichs通量
F_L = euler_flux(U_L, gamma);
F_R = euler_flux(U_R, gamma);
F_num = 0.5*(F_L + F_R - 1.5*lambda_max*(U_R - U_L)); % 增加系数提高耗散
% 应用数值通量
rhs11(:,1,e+1)=rhs11(:,1,e+1)-F_num;
end
end
%% 质量矩阵的逆乘右端项:inv(M)*rhs
mass=0.5*dx*psi*diag(weights)*(psi)';
rhs11_11=zeros(kp1,1);
rhs11_22=zeros(kp1,1);
rhs11_33=zeros(kp1,1);
rho_11=zeros(kp1,1);
u_11 =zeros(kp1,1);
E_11 =zeros(kp1,1);
for e=1:n
rho_11(:,1)=rhs11(1,:,e);
u_11(:,1) =rhs11(2,:,e);
E_11(:,1) =rhs11(3,:,e);
rhs11_11(:,e)=mass\rho_1;
rhs11_22(:,e)=mass\u_1;
rhs11_33(:,e)=mass\E_1;
rhs11(1,:,e)=rhs11_11(:,e);
rhs11(2,:,e)=rhs11_22(:,e);
rhs11(3,:,e)=rhs11_33(:,e);
end
end
%% 函数:计算Roe平均
function U_roe = roe_average(U_L, U_R, gamma)
rho_L = U_L(1);
u_L = U_L(2)/rho_L;
E_L = U_L(3);
p_L = (gamma-1)*(E_L - 0.5*rho_L*u_L^2);
rho_R = U_R(1);
u_R = U_R(2)/rho_R;
E_R = U_R(3);
p_R = (gamma-1)*(E_R - 0.5*rho_R*u_R^2);
% Roe平均
sqrt_rho_L = sqrt(rho_L);
sqrt_rho_R = sqrt(rho_R);
rho_roe = sqrt_rho_L * sqrt_rho_R;
u_roe = (sqrt_rho_L*u_L + sqrt_rho_R*u_R) / (sqrt_rho_L + sqrt_rho_R);
H_L = (E_L + p_L) / rho_L;
H_R = (E_R + p_R) / rho_R;
H_roe = (sqrt_rho_L*H_L + sqrt_rho_R*H_R) / (sqrt_rho_L + sqrt_rho_R);
E_roe = rho_roe * (H_roe - 0.5*u_roe^2);
U_roe = [rho_roe; rho_roe*u_roe; E_roe];
end
%% 函数:计算Euler方程通量
function F = euler_flux(U, gamma)
rho = U(1);
u = U(2)/rho;
E = U(3);
p = (gamma-1)*(E - 0.5*rho*u^2);
F = zeros(3,1);
F(1) = rho*u;
F(2) = rho*u^2 + p;
F(3) = u*(E + p);
end
对上述代码分析,给出出现错误的原因和地方,并给出修正后完整的代码
最新发布