Numerical calculation error

This topic has been permanently closed and is being transferred over to MATLAB Answers. Check back later for the link.

Meta Cody
Follow


Numerical calculation error

Nicholas on 12 Sep 2024
Latest activity Reply by Nicholas on 12 Sep 2024

%I encountered the following problem in the calculation: 1. The calculated H is negative, %and I am unsure if the calculation is correct. Some formulas cannot be simplified and still %exist in the form of 5000/51166. 3. Poor overall code fluency
clear all
close
%% 参数定义parameter definition
P = 42;
c = 800;
E = 15000000;
K = 1.8;
P_ya = 18000;
F = 2;
y = 26.8;
R = 20; % radius
syms H B
%H=100;
% 计算破裂角 a Calculate the rupture angle
if K <= 0.5
a = 90;
elseif K <= 1
a = -90 * K + 135;
elseif K <= 3
a = -22.5 * K + 67.5;
else
a = 0;
end
% 显示计算得到的 a 的值
disp(['当 K = ', num2str(K), ' 时，破裂角 a = ', num2str(a), '°']);
%% 求解初始破裂角相关量 Solving the initial rupture angle related quantities
L = H + R * (1 - sind(a));
G_1 = (y * L^2) / (2 * tand(B)); % 三角形块体的自重
p = atan2(tand(P), F); % 折减后的内摩擦角
C = c * L; % 竖直面上粘聚力合力
C_s = c * L / (F * sind(B)); % 破裂面上粘聚力合力
G_0 = 2 * y * H * cosd(a);
z = 0.9 * P; % 按照围岩等级取值，三级围岩取0.9
%% 定义目标函数 E(B) Define the objective function E (B)
%E_func = @(B) (y ./ (2 .* tand(B))) .* sind(B + p) ./ cosd(B + p - z);
E_func=@(B) (cosd(B+p).*sind(B)).*cosd(B).*cosd(B+p-z)+sind(B+p).*sind(B).*(sind(B+p-z).*cosd(B)+cosd(B+p-z).*sind(B));
%% 数值求导函数 Numerical derivative function
% 使用中心差分法计算导数
dE_func = @(B) (E_func(B + 1e-6) - E_func(B - 1e-6)) / (2e-6);
%% 数值寻找导数为零的 B 值
% 只寻找一个接近的 B 值
B_range = [0, 90]; % B 的取值范围
B_init = 45; % 初始猜测值，设置为 45 度
% 使用 fzero 寻找导数为零的 B 值
try
B_zero = fzero(dE_func, B_init);
% 检查找到的 B 值是否满足条件
if abs(dE_func(B_zero)) < 1e-6
disp(['找到满足条件的 B 值为：', num2str(B_zero)]);
else
disp('没有找到导数接近零的 B 值');
end
catch
disp('fzero 计算失败，未找到满足条件的 B 值');
end
B=B_zero
%% 计算埋深 Calculate burial depth
f1 = ((G_1 - C) .* sind(p + B_zero) + C_s .* cosd(p)) ./ cosd(B_zero + p - z);
f2 = (P_ya - G_0 - 2 .* C) / (2 * sind(z));
% 定义控制方程，解出 H
eqn = f1 - f2 == 0;
% 使用 solve 反解出 H
sol_H_sym = solve(eqn, H);
% 将符号解转换为具体的数值
sol_H_num = double(subs(sol_H_sym));
% 显示结果
disp(['解出的 H 的值为：', num2str(sol_H_num)]);

1 Comment

Comments have been disabled.

Nicholas
Nicholas on 12 Sep 2024
help me! thanks you!
1803905353@qq.com
Nicholas

Nicholas

1
Post
1
Reply
0
Followers
Follow
View profile

Posts by this author

Tags

Go to top of page