function v = vertical_stress(p_i, h, K) % p_i: 剩余滑坡推力合力，单位为MPa % h: 隧道拱顶至水平地面的高度5，单位为m % K: 集中力作用下的应力分布系数 % q_vf : 滑坡推力对拱顶的附加力，单位为MPa % 定义常量 r = 24; % 围岩重度，单位为KN/m^3 u = 30/180*pi; % 顶板土柱两侧摩擦角，单位为弧度 B = 4; % 隧道上方条块的宽度，单位为m Y = 0.5; % 内外侧的侧压力系数 p_i = 500; %剩余滑坡推力合力，单位为MPa % 计算q_vf K = 3 * p_i / (2 * pi * h^2 * sin(u));; q_pv = r * h / (1 - Y) / (B + h * tan(u)); q_vf = K * p_i * sin(deg2rad(30)); % 计算q_pv q_pv = r * h / (1 - Y) / (B + h * tan(deg2rad(30))); % 计算合力q_v q_v = q_vf + q_pv; v = q_v; end

这个代码仍然有问题。在计算q_vf时，应该使用输入参数h而不是固定数值5。因此，可以将q_vf的计算公式修改为：

q_vf = K * p_i * sin(u);

另外，代码中有一个多余的分号符号。可以将第6行代码中的两个分号符号之一删除，即可正常运行。最终的代码如下所示：

function v = vertical_stress(p_i, h, K)
% p_i: 剩余滑坡推力合力，单位为MPa
% h: 隧道拱顶至水平地面的高度，单位为m
% K: 集中力作用下的应力分布系数
% q_vf : 滑坡推力对拱顶的附加力，单位为MPa

% 定义常量
r = 24; % 围岩重度，单位为KN/m^3
u = 30/180*pi; % 顶板土柱两侧摩擦角，单位为弧度
B = 4; % 隧道上方条块的宽度，单位为m
Y = 0.5; % 内外侧的侧压力系数

% 计算q_vf
K = 3 * p_i / (2 * pi * h^2 * sin(u));
q_vf = K * p_i * sin(u);

% 计算q_pv
q_pv = r * h / (1 - Y) / (B + h * tan(u));

% 计算合力q_v
q_v = q_vf + q_pv;

v = q_v;
end