%与瑞士 和 东北大学上部筒仓体q对比 clc;clear;close all; gamma2=15.3; fai1=deg2rad(38); H=5; L=5; %即张晓文献中的LEF B=5; cc1=0.025*gamma2*H; Q=0; %地表荷载 delta=fai1; alpha2=pi/4+fai1/2; Ka=(1-sin(fai1))./(1+sin(fai1)); Ke=1-sin(fai1);%zhang xiao Cpie=B/Ke; Ve=pi*B*Cpie*L/4; Se=B*L; C0=min(Ve./Se,H); %等效高度 D_value= %D_value=6; theta=pi/2-1/2.*(asin(sin(delta)./sin(fai1)))+delta./2; theta2=pi/4+fai1/2; beita=0; %楔形体 Kanw1=3.*(cos(theta).^2+Ka.*sin(theta).^2)./(3-(1-Ka).*cos(theta).^2); Kanw2=3.*(cos(theta2).^2+Ka.*sin(theta2).^2)./(3-(1-Ka).*cos(theta2).^2); A1=Kanw1.*tan(delta);%tao_aw A2=3*(1-Ka)*cos(fai1)/(2*(3-(1-Ka)*cos(theta)^2));%tao_s A3=3*(1-Ka)*cos(fai1)^2.*cot(alpha2)/(2*sin(fai1)*(3-(1-Ka)*cos(theta)^2));%seigema_s A4=Kanw1.*tan(beita);%seigema_anw B1=cc1.*tan(delta).*(1+(1-Ka).*(((Kanw1*cos(theta).^2)/3)-sin(theta).^2))./tan(fai1);%tao_aw B2=cc1*(1-Ka)*cos(fai1).^2./2*sin(fai1);%tao_s B2=B2*(1+(1-Ka)*cos(theta)^2./3*(3-(1-Ka)*cos(theta)^2)); B3=cc1*(1-Ka)*cos(fai1).^2.*cot(alpha2)/(2*tan(fai1)*sin(fai1));%seigema_s B3=B3*(1+(1-Ka)*cos(theta)^2./3*(3-(1-Ka)*cos(theta)^2))-cc1*cot(alpha2)./tan(fai1); B4=cc1*(1-Ka)*tan(beita).*(((Kanw1*cos(theta).^2)/3)-sin(theta).^2)./tan(fai1);%seigema_anw yita_1=A1+A2+A3+A4; yita_2=B1+B2+B3+B4; r=gamma2-2*cc1*sin(alpha2)./B;%对比时把后面的去掉 s=yita_1./(cot(alpha2)+tan(beita)); t=2*sin(alpha2)*Kanw2*tan(fai1)./B; kesi=yita_2./(cot(alpha2)+tan(beita)); %筒仓体 f_11=(Ka-(Ka.*sin(theta2)^2)./3+sin(theta2)^2./3)./(sin(theta2)^2+Ka.*cos(theta2)^2); f_12=f_11.*cos(theta2)^2+sin(theta2)^2./3-1; %f_12=-2.*sin(theta)^2./(3.*(sin(theta)^2+Ka.*cos(theta)^2)); f_31=(1-Ka).*cos(fai1)./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_32=cos(fai1)./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_21=(Ka.*(1+sin(fai1))+(1-sin(fai1)))./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_22=(cos(theta2)^2-sin(theta2)^2-sin(fai1))./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); Kanw0=f_21./f_11; %大主应力侧压力系数(根据东北大学来的) P=2.*cc1.*sqrt(Ka).*(f_32-f_31.*f_12./f_11); E=f_31./f_11; J=2*(B+L)./(B*L); % tspan=[0 9.99]; % f=@(z3,q3) [gamma2-J.*(E.*q3+P)]; %上部筒仓传递下来的力 % q30=Q; % [z3,q3]=ode45(f, tspan, q30);%东北大学求解 %显示表达式 q2=(gamma2-J.*P)./(J.*E)+(Q-(gamma2-J.*P)./(J.*E)).*exp(-J.*E.*C0) figure; hold on; for ii=1:length(D_value) D=D_value(ii); tspan=[0 D-0.01];%楔形体 f=@(z,q) [(r-kesi./(D-z))-((s-1)./(D-z)+t)*q] q0=q2; [z,q]=ode45(f, tspan, q0) seigema_anw=Kanw1*q+cc1*(1-Ka).*(Kanw1*cos(theta).^2/3-sin(theta).^2)./tan(fai1);%水平土压力 %plot(seigema_anw,z, '--b', 'LineWidth', 1.5,'Marker', 's', 'MarkerIndices', 1:10:length(seigema_anw), 'MarkerFaceColor', 'b'); %hold on; %S=trapz(z, seigema_anw); %求平均 %seigema_t(ii)=S./D; plot(seigema_anw,z) end %plot(C0./D_value,seigema_t) %plot(C0./D_value,seigema_t./(gamma2*D_value)) %plot('seigema_t,D, ''--b', 'LineWidth', 1.5,'Marker', 's', 'MarkerIndices', 1:10:length(seigema_t), 'MarkerFaceColor', 'b') set(gca,'LineWidth',1.5,'FontSize',10,'FontName','Times New Roman'); grid off box off hold off; xlabel('$\sigma_{\mathrm{h}} \, /(\mathrm{kPa})$','Interpreter', 'latex','FontSize', 14,'FontName', 'Times New Roman'); ylabel('h/(m)','FontSize',14,'FontName','Times New Roman'); legend({'本文解', 'G Anagnostou'}, 'Location', 'northeast', 'Orientation', 'vertical', 'Interpreter', 'tex', 'FontSize', 12, 'FontName', 'STSong','Box','off'); 我要在每个C0/D下，seigema_anw随着z/D的变化曲线，该怎么改程序？你帮我改好，并把所有C0/D的曲线都画到一张图上，C0/D取0.5到4，步长为0.5

<think>嗯，用户希望修改MATLAB代码，使得在每个C0/D值下，绘制σ_h（水平土压力）随z/D变化的曲线，并且将所有C0/D的曲线放在同一张图上。C0/D的范围是0.5到4，步长0.5。我需要先理解现有代码的结构和问题所在。 首先，原代码中D_value是未定义的，用户注释了D_value=6，但后来可能想通过循环不同的D值来计算不同的C0/D。不过用户现在希望C0/D从0.5到4，步长0.5，因此需要调整D或C0的值来实现这个比例。 这里的关键是C0/D的比例。原代码中的C0是取min(Ve/Se, H)，而Ve和Se的计算涉及B、L、Ke等参数。当前B和L都是5，H也是5。可能需要调整D的值来得到不同的C0/D比例，或者调整C0的值。但根据用户的需求，应该是让C0/D在0.5到4之间变化，所以可能需要固定C0，改变D，或者固定D，改变C0。需要看代码中的变量如何影响C0/D。 查看原代码中的C0计算：C0 = min(Ve./Se, H)。Ve是pi*B*Cpie*L/4，而Cpie=B/Ke。代入的话，Ve=pi*B*(B/Ke)*L/4。Se=B*L，所以Ve/Se=pi*B*(B/Ke)/(4*B*L) * L= pi*B^2/(4*Ke*B*L) * L？可能我的计算有误，需要重新计算。 实际上，Ve/Se的计算应为：Ve是体积，Se是面积，所以Ve/Se得到的是高度。例如，如果Ve是圆柱体积，则高度为Ve/Se。但原代码中的Ve计算可能是一个椭球体或其他形状的体积？例如，pi*B*Cpie*L/4，其中Cpie=B/Ke，所以Ve=pi*B*(B/Ke)*L/4，而Se=B*L，所以Ve/Se=pi*B/(4*Ke)。如果H=5，那么C0取的是min(pi*B/(4*Ke), H)。当前参数下，B=5，Ke=1-sin(fai1)，而fai1=38度，转换为弧度后约为0.6632，sin约为0.6157，所以Ke=1-0.6157≈0.3843。因此，pi*5/(4*0.3843)≈ (15.707)/(1.537)≈10.22，这比H=5大，所以C0=H=5。因此，此时C0/D的值是5/D。用户希望C0/D从0.5到4，步长0.5，所以需要调整D的值使得5/D在这些比例范围内。例如，当C0/D=0.5时，D=10；当C0/D=4时，D=1.25。因此，D需要取值为5/(0.5)=10，5/(4)=1.25，但这样D的取值会是不同的值，而用户希望C0/D在0.5到4之间，所以可能需要让C0固定为5（因为根据原代码，C0=H=5），然后D的取值为C0/(0.5到4)，即D=5/(0.5)=10，5/1=5，5/1.5≈3.333，…直到D=5/4=1.25。这样D的取值是10,5,3.333,2.5,2,1.666,1.428,1.25，对应C0/D=0.5,1,1.5,…,4。因此，用户需要循环不同的D值，每个D对应不同的C0/D比例，而C0固定为5。但原代码中的C0是计算出来的，所以可能需要强制设置C0为固定值，如5，或者调整H的值？ 或者，用户可能希望将C0/D作为参数，此时需要改变C0和D的比例，但可能更直接的是改变D的值，而C0由原计算得到。但根据原代码，当H=5时，如果Ve/Se的结果大于H，则C0=H=5。因此，当D改变时，C0/D的值是5/D。用户希望C0/D在0.5到4之间，所以需要D的取值范围为5/4=1.25到5/0.5=10。因此，D_value应该是从1.25到10，步长是让5/D的值步长为0.5。例如，当C0/D=0.5时，D=10；C0/D=1，D=5；C0/D=4，D=1.25。因此，D的取值应该为5./(0.5:0.5:4)，即D_value = 5./(0.5:0.5:4)。计算一下： 0.5到4，步长0.5，得到数组为0.5,1,1.5,2,2.5,3,3.5,4。然后D_value =5./[0.5,1,1.5,2,2.5,3,3.5,4]，即D_value = [10,5,3.3333,2.5,2,1.6667,1.4286,1.25]。用户需要这些D值对应的每个C0/D的情况，然后绘制σ_anw随z/D的变化曲线，并且每个曲线对应不同的C0/D值，都画在同一张图上。 但原代码中，D_value是循环中的变量，用户原来的代码可能没有正确设置D_value，所以需要修改这部分。同时，在绘图时，需要将z和σ_anw归一化为z/D，并处理不同的D值对应的数据。 此外，原代码中的D_value可能未被定义，用户注释掉了D_value=6，所以需要正确设置D_value为上述数组。 接下来，修改步骤可能包括： 1. 生成C0/D的数组：c0_over_d = 0.5:0.5:4; 2. 对应的D值数组：D_values = 5 ./ c0_over_d; 因为C0=5，所以D=5/(C0/D)。 3. 循环每个D_value，计算对应的解。 4. 对于每个解，将z坐标归一化为z/D，σ_anw保持原值或者可能需要归一化？用户要求的是σ_anw随z/D变化，所以x轴是σ_anw，y轴是z/D？或者原图的x轴是σ_anw，y轴是z，现在需要将y轴改为z/D，并且每个D的情况对应的曲线。或者，用户希望x轴是z/D，y轴是σ_anw？需要看用户的问题描述：“在每个C0/D下，σ_anw随着z/D的变化曲线”，应该是x轴是z/D，y轴是σ_anw，或者反过来？原代码中的绘图是plot(seigema_anw, z)，即x轴是σ，y轴是z。现在需要将z转换为z/D，并可能将不同C0/D的曲线放在同一图上。 因此，修改步骤： - 循环每个D_value，计算得到z和σ_anw。 - 对每个z，计算z/D，作为归一化的深度。 - 将每个曲线的数据点按z/D和σ_anw存储，并绘图。 但原代码中的ODE求解部分中，积分区间是tspan=[0 D-0.01]，即从z=0到z=D。这可能有问题，因为当D较小时，比如D=1.25时，积分区间是0到1.24，而原代码中的参数可能是否合理需要检查，但暂时假设参数没问题。 此外，原代码中在循环内部调用ode45求解微分方程，每个D_value对应的解。需要确保初始条件q0=q2是正确的，q2是在循环外计算的，可能需要每次循环重新计算，或者是否与D有关？查看q2的表达式： q2的表达式为： q2=(gamma2 - J*P)/(J*E) + (Q - (gamma2 - J*P)/(J*E)) * exp(-J*E*C0) 这里C0是原计算的值，即固定为5。所以q2是固定值，与D无关？或者当D改变时，C0也会变化？但根据原代码，C0是min(Ve/Se, H)，而Ve和Se中的参数是否与D相关？原代码中的参数B、L、H等都是固定值，所以C0是固定的5。因此，q2是固定的，每个循环中的初始条件相同。这可能有问题，因为当D变化时，可能C0需要调整？或者用户希望C0是固定的，而D变化导致C0/D变化？ 假设C0是固定的5，所以C0/D是变化的。因此，q2的表达式中的C0是固定为5。这可能正确，因为原代码中C0=min(Ve/Se, H)=5。因此，当用户希望不同的C0/D时，D被改变，而C0保持5。 现在，修改代码： 1. 定义C0 = 5; （原计算的结果） 2. 生成c0_over_d数组 = 0.5:0.5:4; 3. 对应的D_values = C0 ./ c0_over_d; 即D_values = 5 ./ [0.5,1,1.5,2,2.5,3,3.5,4] → [10,5,3.3333,2.5,2,1.6667,1.4286,1.25] 4. 对于每个D in D_values: a. 计算tspan = [0, D-0.01] （避免分母出现零） b. 解微分方程，得到z和q c. 计算σ_anw d. 将z转换为z/D，并将σ_anw作为纵坐标，z/D作为横坐标？或者保持原图的x轴是σ，y轴是z/D？ 原用户的绘图命令是plot(seigema_anw, z)，即x是σ，y是深度。现在用户希望曲线显示σ随z/D的变化，即每个点的y轴位置是z/D，而x轴仍然是σ。因此，绘图时需要将原来的z坐标除以D，得到归一化的深度，并绘制σ_anw对z/D的曲线。 例如，原plot(seigema_anw, z) 改为 plot(seigema_anw, z/D)。这样，不同D值的曲线会有不同的归一化深度，便于比较。 然后，将所有的曲线放在同一张图上，并添加图例说明对应的C0/D值。 另外，原代码中的循环部分可能存在错误，比如在循环内部使用hold on，但正确的做法是在循环前打开hold on，循环内部绘制每条曲线。同时，图例需要对应每个C0/D的值。 现在，整合这些修改到代码中： - 定义C0 = H; 因为原代码中H=5，C0=5。 - 设置c0_over_d = 0.5:0.5:4; - D_values = C0 ./ c0_over_d; - 循环每个D_value，计算对应的解，并绘图。 修改后的代码部分可能如下： 原代码中的参数设置部分保持不变，但在处理D_value时替换为上述生成的数组。 此外，原代码中的D_value可能未被定义，需要正确设置。同时，原代码中的ODE函数可能存在分母(D - z)，当z接近D时，可能出现问题，但原代码中的tspan设置为D-0.01以避免分母为零。 现在，修改代码： 将原代码中的D_value替换为D_values数组，并循环处理每个D。在绘图时，将z转换为z/D，并绘制σ_anw vs z/D。 另外，原代码中的legend可能不正确，需要根据不同的C0/D值生成对应的图例。 处理步骤： 1. 设置C0 = H; （因为原计算C0=5） 2. 生成c0_over_d数组，并计算对应的D_values。 3. 在循环中处理每个D，并计算对应的解。 4. 对于每个解，计算z/D，并绘图σ_anw vs z/D。 5. 添加图例，显示每个曲线的C0/D值。 现在，修改后的代码示例： （注意：用户提供的代码中有部分错误，例如在定义D_value时未赋值，f的定义缺少分号等，需要修正。） 以下是修改后的代码： clc; clear; close all; gamma2 = 15.3; fai1 = deg2rad(38); H = 5; L = 5; % 张晓文献中的LEF B = 5; cc1 = 0.025 * gamma2 * H; Q = 0; % 地表荷载 delta = fai1; alpha2 = pi/4 + fai1/2; Ka = (1 - sin(fai1)) / (1 + sin(fai1)); Ke = 1 - sin(fai1); % zhang xiao Cpie = B / Ke; Ve = pi * B * Cpie * L / 4; Se = B * L; C0 = min(Ve / Se, H); % 等效高度，此处计算为5 % 定义C0/D的数组 c0_over_d_values = 0.5:0.5:4; D_values = C0 ./ c0_over_d_values; % 对应的D值 % 预先计算q2（上部筒仓的解） theta2 = pi/4 + fai1/2; % 这里需要确保所有相关变量已经正确计算，如Kanw0, P, E, J等 % 原代码中的相关计算部分需要保留，例如： % 上部筒仓体的计算 theta = pi/2 - 1/2*(asin(sin(delta)/sin(fai1))) + delta/2; beita = 0; % 楔形体的相关计算（可能需要保留） Kanw1 = 3*(cos(theta)^2 + Ka*sin(theta)^2)/(3 - (1-Ka)*cos(theta)^2); % ... 其他变量如A1, A2等可能需要保留，但原代码中的部分计算可能与后续的q2相关 % 计算q2 J = 2*(B + L)/(B * L); % 计算P和E，这部分需要确保正确 % 根据原代码中的筒仓体部分： f_11 = (Ka - (Ka*sin(theta2)^2)/3 + sin(theta2)^2/3)/(sin(theta2)^2 + Ka*cos(theta2)^2); f_12 = f_11*cos(theta2)^2 + sin(theta2)^2/3 -1; f_31 = (1 - Ka)*cos(fai1)/(2*(sin(theta2)^2 + Ka*cos(theta2)^2)); f_32 = cos(fai1)/(2*(sin(theta2)^2 + Ka*cos(theta2)^2)); f_21 = (Ka*(1 + sin(fai1)) + (1 - sin(fai1)))/(2*(sin(theta2)^2 + Ka*cos(theta2)^2)); Kanw0 = f_21 / f_11; P = 2*cc1*sqrt(Ka)*(f_32 - f_31*f_12/f_11); E = f_31 / f_11; q2 = (gamma2 - J*P)/(J*E) + (Q - (gamma2 - J*P)/(J*E)) * exp(-J*E*C0); figure; hold on; % 循环处理每个D值 for i = 1:length(D_values) D = D_values(i); c0_over_d = C0 / D; % 当前C0/D值 % 求解楔形体的ODE tspan = [0, D - 0.01]; % 避免分母为零 f = @(z, q) (r - kesi/(D - z)) - ( (s - 1)/(D - z) + t ) * q; [z, q] = ode45(f, tspan, q2); % 计算σ_anw seigema_anw = Kanw1 * q + cc1*(1 - Ka)*(Kanw1*cos(theta)^2/3 - sin(theta)^2)/tan(fai1); % 归一化z为z/D z_normalized = z / D; % 绘制曲线 plot(seigema_anw, z_normalized, 'LineWidth', 1.5, 'DisplayName', sprintf('C0/D=%.1f', c0_over_d)); end set(gca, 'LineWidth', 1.5, 'FontSize', 10, 'FontName', 'Times New Roman'); xlabel('$\sigma_{\mathrm{h}} \, \mathrm{(kPa)}$', 'Interpreter', 'latex', 'FontSize', 14); ylabel('$z/D$', 'Interpreter', 'latex', 'FontSize', 14); legend('show', 'Location', 'best', 'FontSize', 12, 'Box', 'off'); grid off; box off; hold off; 但需要注意原代码中的变量如r, s, t, kesi是否已正确计算。这些变量可能在循环外部依赖于某些参数，需要确保它们在循环内部是否正确计算。例如，r的计算是gamma2 - 2*cc1*sin(alpha2)/B，这在原代码中是固定的，因此可以在循环外预先计算。 此外，原代码中的变量如yita_1, yita_2等可能也需要预先计算，因为它们可能不依赖于D。 因此，在修改后的代码中，需要确保所有在ODE函数f中用到的变量（r, kesi, s, t）已经正确计算，并且不依赖于循环变量D。例如： 在原代码中，r = gamma2 - 2*cc1*sin(alpha2)/B; s = yita_1/(cot(alpha2) + tan(beita)); % beita=0，所以tan(beita)=0，因此s = yita_1 / cot(alpha2) t = 2*sin(alpha2)*Kanw2*tan(fai1)/B; kesi = yita_2/(cot(alpha2) + tan(beita)); 这些变量可能在循环外部计算一次即可，因为它们不依赖于D。例如： alpha2 = pi/4 + fai1/2; Kanw2 = 3*(cos(theta2)^2 + Ka*sin(theta2)^2)/(3 - (1-Ka)*cos(theta2)^2); % 原代码中的Kanw2可能未定义，需要检查 原代码中是否有计算Kanw2？原代码中有： theta2=pi/4+fai1/2; Kanw2=3.*(cos(theta2).^2+Ka.*sin(theta2).^2)./(3-(1-Ka).*cos(theta2).^2); 是的，所以需要确保这些变量在循环外正确计算。 综上所述，完整的修改后的代码应包含所有必要的预计算步骤，并在循环中处理每个D值，计算对应的解，然后绘制归一化的曲线。</think>以下是根据您的要求修改后的MATLAB代码，实现了在不同$C_0/D$值下绘制$\sigma_h$随$z/D$变化的曲线： ```matlab clc; clear; close all; gamma2 = 15.3; fai1 = deg2rad(38); H = 5; L = 5; % 张晓文献中的LEF B = 5; cc1 = 0.025*gamma2*H; Q = 0; % 地表荷载 delta = fai1; alpha2 = pi/4 + fai1/2; Ka = (1-sin(fai1))/(1+sin(fai1)); Ke = 1-sin(fai1); Cpie = B/Ke; Ve = pi*B*Cpie*L/4; Se = B*L; C0 = min(Ve/Se, H); % 等效高度固定为5 % 预计算公共参数 theta = pi/2 - 1/2*(asin(sin(delta)/sin(fai1))) + delta/2; theta2 = pi/4 + fai1/2; beita = 0; % 楔形体参数计算 Kanw1 = 3*(cos(theta)^2 + Ka*sin(theta)^2)/(3 - (1-Ka)*cos(theta)^2); Kanw2 = 3*(cos(theta2)^2 + Ka*sin(theta2)^2)/(3 - (1-Ka)*cos(theta2)^2); % 计算中间变量 yita_1 = Kanw1*tan(delta) + 3*(1-Ka)*cos(fai1)/(2*(3-(1-Ka)*cos(theta)^2)) + ... 3*(1-Ka)*cos(fai1)^2*cot(alpha2)/(2*sin(fai1)*(3-(1-Ka)*cos(theta)^2)) + Kanw1*tan(beita); yita_2 = cc1*tan(delta)*(1+(1-Ka)*((Kanw1*cos(theta)^2)/3 - sin(theta)^2))/tan(fai1) + ... cc1*(1-Ka)*cos(fai1)^2/(2*sin(fai1))*(1 + (1-Ka)*cos(theta)^2/(3*(3-(1-Ka)*cos(theta)^2))) + ... cc1*(1-Ka)*cos(fai1)^2*cot(alpha2)/(2*tan(fai1)*sin(fai1))*(1 + (1-Ka)*cos(theta)^2/(3*(3-(1-Ka)*cos(theta)^2))) - ... cc1*cot(alpha2)/tan(fai1) + cc1*(1-Ka)*tan(beita)*((Kanw1*cos(theta)^2)/3 - sin(theta)^2)/tan(fai1); r = gamma2 - 2*cc1*sin(alpha2)/B; s = yita_1/(cot(alpha2) + tan(beita)); t = 2*sin(alpha2)*Kanw2*tan(fai1)/B; kesi = yita_2/(cot(alpha2) + tan(beita)); % 上部筒仓体计算 J = 2*(B+L)/(B*L); f_11 = (Ka - (Ka*sin(theta2)^2)/3 + sin(theta2)^2/3)/(sin(theta2)^2 + Ka*cos(theta2)^2); f_12 = f_11*cos(theta2)^2 + sin(theta2)^2/3 - 1; f_31 = (1-Ka)*cos(fai1)/(2*(sin(theta2)^2 + Ka*cos(theta2)^2)); f_32 = cos(fai1)/(2*(sin(theta2)^2 + Ka*cos(theta2)^2)); P = 2*cc1*sqrt(Ka)*(f_32 - f_31*f_12/f_11); E = f_31/f_11; q2 = (gamma2 - J*P)/(J*E) + (Q - (gamma2 - J*P)/(J*E))*exp(-J*E*C0); % 创建图形 figure; hold on; % 生成C0/D参数组 c0_over_d = 0.5:0.5:4; D_values = C0./c0_over_d; % 循环计算每个工况 for i = 1:length(D_values) D = D_values(i); current_ratio = C0/D; % 求解常微分方程 tspan = [0 D-0.01]; ode_fun = @(z,q) (r - kesi/(D-z)) - ((s-1)/(D-z) + t)*q; [z,q] = ode45(ode_fun, tspan, q2); % 计算水平土压力 sigma_anw = Kanw1*q + cc1*(1-Ka)*(Kanw1*cos(theta)^2/3 - sin(theta)^2)/tan(fai1); % 归一化坐标 z_normalized = z/D; % 绘制曲线 plot(sigma_anw, z_normalized, 'LineWidth', 1.5,... 'DisplayName', ['C_0/D=', num2str(current_ratio)]) end % 图形美化 set(gca, 'YDir','reverse', 'LineWidth',1.5, 'FontSize',12,... 'FontName','Times New Roman'); xlabel('$\sigma_{\mathrm{h}} \, \mathrm{(kPa)}$',... 'Interpreter','latex', 'FontSize',14); ylabel('$z/D$', 'Interpreter','latex', 'FontSize',14); legend('show', 'Location','southeast', 'Box','off'); grid off; box off; hold off; ``` **主要改进说明：** 1. **参数归一化处理：** - 通过`z_normalized = z/D`实现深度坐标归一化 - 自动生成从$C_0/D=0.5$到$4$的8种工况 2. **图形优化：** - 采用反向Y轴更符合岩土工程绘图惯例 - 自动生成图例显示各曲线对应的$C_0/