% EG20H1 CA3 优化设计 — 学号 50111491
clc; clear;

% ===== 你的设计参数 =====
F   = 90000;        % N
h   = 0.003;        % m
s   = 0.04;         % m
r   = 0.002;        % m
sigma_all = 300e6;  % Pa

% ===== 搜索范围与步长 =====
% d 从 5mm 扫到 80mm，t 从 0.5mm 扫到 15mm
d_vec = 0.005 : 0.001 : 0.080;
t_vec = 0.0005 : 0.0005 : 0.015;

sol = [];  % [d, t, D, sigma_max, Volume]

for d = d_vec
    for t = t_vec
        D = d + 2*h;

        % --- 应力集中系数 Kt（使用模板第2页给出的公式）---
        ratio = 2*h / D;
        Kt = 2.09 - 0.38*ratio - 1.21*ratio^2 + 0.5*ratio^3;

        % --- 名义应力与最大应力 ---
        sigma_nom = F / (d * t);
        sigma_max = Kt * sigma_nom;

        % --- 体积（目标函数）---
        V = s*t*(D + d) + (4 - pi)/2 * r^2 * t;

        % --- 可行性判断 ---
        if sigma_max <= sigma_all
            sol = [sol; d, t, D, sigma_max, V];
        end
    end
end

if isempty(sol)
    error('搜索范围内没有可行解，请增大 d 或 t 的上限。');
end

% ===== 按体积升序排序，取最优 =====
sol_sorted = sortrows(sol, 5);
op = sol_sorted(1, :);

d_opt   = op(1);
t_opt   = op(2);
D_opt   = op(3);
sig_opt = op(4);
V_opt   = op(5);

% ===== 输出结果 =====
fprintf('===== 最优设计结果 =====\n');
fprintf('d  (mm) = %.3f\n', d_opt*1000);
fprintf('t  (mm) = %.3f\n', t_opt*1000);
fprintf('D  (mm) = %.3f\n', D_opt*1000);
fprintf('最大应力 (MPa) = %.2f\n', sig_opt/1e6);
fprintf('体积 (mm^3)   = %.2f\n', V_opt*1e9);

% 校核：最大应力应 <= 300 MPa，且尽量接近
fprintf('\n许用应力 = %.0f MPa，设计利用率 = %.1f%%\n', ...
    sigma_all/1e6, sig_opt/sigma_all*100);