%% 喷管 MOC 特征线设计与迭代收敛脚本（MATLAB）
% 说明：基于特性线法（MOC）实现五类特征点分类（对称轴、转向区点、
% 内点、末线点、壁面点），并采用迭代直到 x 坐标变化绝对误差 < 1e-8。
% 最后绘制三幅图：
% 1) 特征线网格及五类点标记；
% 2) 优化后离散壁面节点（点网格）；
% 3) 优化后喷管外轮廓（节点连线）。
% 并将最终壁面坐标导出到 Excel。
% 注释均为中文。
clc; close all; clearvars;
%% ===================== 1. 输入参数 =====================
p1 = 4e5; % 室内燃烧室总压 p1 (Pa)
T1 = 1000; % 室内燃烧室总温 T1 (K)
M_e = 2.0; % 设计喷管出口马赫数
gamma = 1.38; % 气体比热比 γ
Rgas = 287; % 气体常数 R (J/(kg·K))
ALT = 10000; % 飞行高度 ALT (m)
%% ===================== 2. 大气环境计算（可选） =====================
if ALT < 11000
T_env_C = 15.04 - 0.00649*ALT;
p0 = 1000*(101.29*((T_env_C+273.1)/288.08)^5.256);
elseif ALT < 25000
T_env_C = -56.46;
p0 = 1000*(22.65*exp(1.73-0.000157*ALT));
else
T_env_C = -131.21 + 0.00299*ALT;
p0 = 1000*(2.488*((T_env_C+273.1)/216.6)^(-11.388));
end
T_env_K = T_env_C + 273.15;
if p1 <= p0
error('室内总压 p1 必须大于环境静压 p0，才能形成超音速膨胀！');
end
%% ===================== 3. Prandtl–Meyer 与喉部条件 =====================
M_throat = 1.0; % 假设喉部Ma = 1
A = sqrt((gamma+1)/(gamma-1));
B = (gamma-1)/(gamma+1);
vPM = @(M) A*atan(sqrt(B*(M.^2-1))) - atan(sqrt(M.^2 - 1));
nu_e = vPM(M_e); % 出口Prandtl–Meyer角 (rad)
nu_e_deg = nu_e * (180/pi); % (deg)
T_max = 0.5 * nu_e_deg; % 最大偏转角 (deg)
%% ===================== 4. 特征线网格初始化 =====================
% 离散化偏转角：从 theta_start 到 theta_start+T_max，步长 0.25°
theta_start = mod(90 - T_max, 0.25);
if theta_start == 0
theta_start = 0.25;
end
theta_deg = theta_start : 0.25 : (T_max + theta_start);
N_char = length(theta_deg); % 特征线条数
theta_rad = theta_deg * (pi/180); % 转为弧度
% 预分配
M_char = zeros(1, N_char);
P = zeros(1, N_char);
RR = zeros(1, N_char);
SL = zeros(1, N_char);
TR = 0.0545; % 喉部半径 (m)
for i = 1:N_char
f = @(M) vPM(M) - theta_rad(i);
M_char(i) = fzero(f, [1, M_e + 0.1]); % 求解 Prandtl–Meyer 反函数
P(i) = TR * tan(theta_rad(i)); % R 特征线与轴线交点 x
RR(i) = -TR / P(i); % R 特征线斜率
mu_i = asin(1 / M_char(i));
SL(i) = tan(theta_rad(i) + mu_i); % L 特征线斜率
end
%% ===================== 5. 初步壁面节点计算 =====================
TM_rad = T_max * (pi/180); % 最大偏转角 (rad)
num_wall = N_char + 1; % 壁面点数量
% 临时存储 N_char 部分（不含喉部）
temp_xw = zeros(1, N_char);
temp_yw = zeros(1, N_char);
% 第一段壁面
temp_xw(1) = (TR + SL(1)*P(1)) / (SL(1) - tan(TM_rad));
temp_yw(1) = tan(TM_rad)* temp_xw(1) + TR;
% 中间各段
for i = 2:N_char-1
s_i = tan(TM_rad) - (i-1)*(tan(TM_rad)/(N_char-1));
b_i = temp_yw(i-1) - s_i * temp_xw(i-1);
temp_xw(i) = (b_i + SL(i)*P(i)) / (SL(i) - s_i);
temp_yw(i) = s_i * temp_xw(i) + b_i;
end
% 最后一段（斜率 = 0）
s_last = 0;
b_last = temp_yw(N_char-1);
temp_xw(N_char) = (b_last + SL(N_char)*P(N_char)) / SL(N_char);
temp_yw(N_char) = b_last;
% 合成最终壁面节点，插入喉部 (0,TR)
xw = [0, temp_xw];
yw = [TR, temp_yw];
%% ===================== 6. 构造点集并初始化分类 =====================
num_axis = 1;
num_turn = N_char;
num_wall = num_wall;
num_last = 1;
num_inner = (N_char - 1)^2;
total_pts = num_axis + num_turn + num_inner + num_wall + num_last;
points_all = repmat(struct('x',NaN,'y',NaN,'theta',NaN,'M',NaN,'type',''), 1, total_pts);
idx = 1;
% 1) 对称轴点
points_all(idx).x = 0;
points_all(idx).y = TR;
points_all(idx).theta = 0;
points_all(idx).M = M_throat;
points_all(idx).type = 'axis';
idx = idx + 1;
% 2) 转向区点
for i = 1:num_turn
points_all(idx).x = P(i);
points_all(idx).y = 0;
points_all(idx).theta = theta_rad(i);
points_all(idx).M = M_char(i);
points_all(idx).type = 'turn';
idx = idx + 1;
end
% 3) 内部点 (初始化为0，迭代后更新)
for i = 1:num_inner
points_all(idx).x = 0;
points_all(idx).y = 0;
points_all(idx).theta = NaN;
points_all(idx).M = NaN;
points_all(idx).type = 'inner';
idx = idx + 1;
end
base_idx_inner = num_axis + num_turn;
% 4) 壁面点
for i = 1:num_wall
points_all(idx).x = xw(i);
points_all(idx).y = yw(i);
points_all(idx).theta = NaN;
points_all(idx).M = NaN;
points_all(idx).type = 'wall';
idx = idx + 1;
end
base_idx_wall = num_axis + num_turn + num_inner;
% 5) 最后一条特征线点
points_all(idx).x = xw(end);
points_all(idx).y = yw(end);
points_all(idx).theta = theta_rad(end);
points_all(idx).M = M_char(end);
points_all(idx).type = 'last';
%% ===================== 7. 迭代求解特征线网格 =====================
tolerance = 1e-8;
max_iter = 1000;
conv_flag = false;
for iter = 1:max_iter
dx_max = 0;
% 更新内部点
for m = 1:(N_char-1)
for n = 1:(N_char-1)
idx_inner = base_idx_inner + (m-1)*(N_char-1) + n;
% 上方点索引
if m == 1
idx_up = num_axis + n; % 转向区
else
idx_up = base_idx_inner + (m-2)*(N_char-1) + n;
end
% 左方点索引
if n == 1
idx_left = idx_up; % 若在行首，则与上方相同
else
idx_left = base_idx_inner + (m-1)*(N_char-1) + (n-1);
end
% 获取斜率
theta_up = points_all(idx_up).theta;
M_up = points_all(idx_up).M;
mu_up = asin(1 / M_up);
slope_Rp = tan(theta_up + mu_up);
theta_left = points_all(idx_left).theta;
M_left = points_all(idx_left).M;
mu_left = asin(1 / M_left);
slope_Lp = tan(theta_left - mu_left);
x_up = points_all(idx_up).x;
y_up = points_all(idx_up).y;
x_left = points_all(idx_left).x;
y_left = points_all(idx_left).y;
% 解特征线交点
A_mat = [ -slope_Rp, 1; -slope_Lp, 1 ];
b_vec = [ y_up - slope_Rp * x_up; y_left - slope_Lp * x_left ];
sol = A_mat \ b_vec;
x_new = sol(1);
y_new = sol(2);
dx = abs(x_new - points_all(idx_inner).x);
dx_max = max(dx_max, dx);
points_all(idx_inner).x = x_new;
points_all(idx_inner).y = y_new;
% 更新 θ 与 M
nu_new = theta_up + mu_up;
theta_new = nu_new / 2;
M_new = fzero(@(M) vPM(M) - 2*theta_new, [1, M_e + 0.1]);
points_all(idx_inner).theta = theta_new;
points_all(idx_inner).M = M_new;
end
end
% 更新壁面点
for k = 1:num_wall
idx_wall = base_idx_wall + k;
if (k == 1) || (k == num_wall)
continue; % 喉部和末点不变
else
i = k - 1;
s_i = tan(TM_rad) - (i-1)*(tan(TM_rad)/(num_wall-2));
b_i = yw(i-1) - s_i * xw(i-1);
xL0 = P(i);
yL0 = 0;
slope_Lp = SL(i);
A2 = [ -slope_Lp, 1; -s_i, 1 ];
b2 = [ -slope_Lp*xL0; b_i ];
sol2 = A2 \ b2;
x_new = sol2(1);
y_new = sol2(2);
dx = abs(x_new - points_all(idx_wall).x);
dx_max = max(dx_max, dx);
points_all(idx_wall).x = x_new;
points_all(idx_wall).y = y_new;
points_all(idx_wall).theta = 0; % 壁面流动角假设为 0
points_all(idx_wall).M = M_char(i);
end
end
% 更新末线点
idx_last = num_axis + num_turn + num_inner + num_wall;
idx_ref = base_idx_wall + (num_wall - 1);
points_all(idx_last).x = points_all(idx_ref).x;
points_all(idx_last).y = points_all(idx_ref).y;
points_all(idx_last).theta = points_all(idx_ref).theta;
points_all(idx_last).M = points_all(idx_ref).M;
% 收敛检查
if dx_max < tolerance
fprintf('迭代收敛: 最大 x 变化 = %.2e, 迭代次数 = %d\n', dx_max, iter);
conv_flag = true;
break;
end
end
if ~conv_flag
warning('达到最大迭代次数(%d)但未收敛，dx_max = %.2e', max_iter, dx_max);
end
%% ===================== 8. 分类收集最终点 =====================
idx = 1;
% 对称轴
points.axis(1,:) = [points_all(idx).x, points_all(idx).y, ...
points_all(idx).theta, points_all(idx).M];
idx = idx + 1;
% 转向区
for i = 1:num_turn
points.turn(i,:) = [points_all(idx).x, points_all(idx).y, ...
points_all(idx).theta, points_all(idx).M];
idx = idx + 1;
end
% 内部
for i = 1:num_inner
r = base_idx_inner + i;
points.inner(i,:) = [points_all(r).x, points_all(r).y, ...
points_all(r).theta, points_all(r).M];
end
% 壁面
for i = 1:num_wall
r = base_idx_wall + i;
points.wall(i,:) = [points_all(r).x, points_all(r).y, ...
points_all(r).theta, points_all(r).M];
end
% 末线点
points.last(1,:) = [points_all(idx).x, points_all(idx).y, ...
points_all(idx).theta, points_all(idx).M];
%% ===================== 9. 绘制特征线网格及五类点 =====================
figure('Name','MOC 特征线网格与点分类','NumberTitle','off');
hold on; box on; grid on;
xlabel('轴向 x (m)'); ylabel('径向 r (m)');
title('MOC 特征线分布及五类点');
% 对称轴 (虚线)
plot([0,0], [0, TR], '--k', 'LineWidth', 1);
% R 特征线 (灰色)
for i = 1:N_char
xR = [0, points.turn(i,1)];
yR = [TR, points.turn(i,2)];
plot(xR, yR, 'Color', [0.5,0.5,0.5], 'LineWidth', 1);
end
% L 特征线 (蓝色)
for i = 1:N_char
xL = [points.turn(i,1), points.wall(i+1,1)];
yL = [points.turn(i,2), points.wall(i+1,2)];
plot(xL, yL, 'b', 'LineWidth', 1);
end
% 标记五类点
scatter(points.axis(:,1), points.axis(:,2), 50, 'k', 'filled'); % 对称轴点
scatter(points.turn(:,1), points.turn(:,2), 50, 'm', 'filled'); % 转向区点
scatter(points.inner(:,1), points.inner(:,2), 8, 'b'); % 内部点
scatter(points.last(:,1), points.last(:,2), 50, 'g', 'filled'); % 末线点
scatter(points.wall(:,1), points.wall(:,2), 50, 'r', 'filled'); % 壁面点
legend({'对称轴','R 特征线','L 特征线','轴线点','转向区点','内部点','末线点','壁面点'}, ...
'Location','BestOutside');
axis equal;
xlim([0, xw(end)*1.02]);
ylim([0, max(yw)*1.02]);
%% ===================== 10. 绘制离散壁面节点（点网格） =====================
figure('Name','离散壁面节点','NumberTitle','off');
plot(xw, yw, 'ro', 'MarkerSize', 5, 'MarkerFaceColor','r');
xlabel('轴向 x (m)');
ylabel('径向 r (m)');
title('优化后离散壁面节点（点网格）');
axis equal;
xlim([0, xw(end)*1.02]);
ylim([0, max(yw)*1.02]);
%% ===================== 11. 绘制连接后的喷管外轮廓 =====================
figure('Name','喷管外轮廓','NumberTitle','off');
plot(xw, yw, 'k-', 'LineWidth', 2);
xlabel('轴向 x (m)');
ylabel('径向 r (m)');
title('优化后喷管外轮廓');
axis equal;
xlim([0, xw(end)*1.02]);
ylim([0, max(yw)*1.02]);
%% ===================== 12. 导出壁面坐标 =====================
n_pts = length(xw);
range_x = sprintf('A1:A%d', n_pts);
range_y = sprintf('B1:B%d', n_pts);
xlswrite('PARAMS.xlsx', xw.', 'PTS', range_x);
xlswrite('PARAMS.xlsx', yw.', 'PTS', range_y);
fprintf('\n最终优化壁面坐标已写入 PARAMS.xlsx 的 \"PTS\" 工作表，共 %d 点。\n', n_pts);
本文详细介绍了CodeIgniter框架更新的主要变化，包括移除旧组件、优化系统类库、新增缓存类及验证码辅助函数等，并增加了对JavaScript的支持。
扫一扫
28
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
