关于辊系弹性模型的，有大佬知道这个模型哪里出错了嘛，结果一下出错。

[2301_79899011]

% 初始化参数

b = 110; % 带材宽度, mm

m = 11; % 条元数量

E = 70e9; % 弹性模量, Pa

mu = 0.3; % 泊松比

k = 500; % 剪切变形抗力, MPa

tau = 100; % 平均摩擦应力, MPa

sigma1 = 100; % 平均前张应力, MPa

v1 = 10; % 出口横向平均速度, m/s

h0=1;

h1=0.8;

s=rand(1,m)*(b/m);%随机生成宽度

s=s/sum(s)*b;%归一化

% y=[0,cumsum(s)];%从0开始，累加每个条元的宽度得到节线坐标

% x = (y(1:end-1) + y(2:end)) / 2;% 计算条元的中间坐标

y=zeros(1,m+1);

for i=1:m

y(i+1)=y(i)+s(i);

end

% 预分配入口和出口厚度数组

h_0 = zeros(1, m); % 预分配入口厚度数组

h_1 = zeros(1, m); % 预分配出口厚度数组

%计算入口和出口厚度

for i=1:m

h_0(i)=rand()*2+1;% 随机生成厚度，范围1-3mm

h_1(i)= h_0(i)-rand()*0.2;% 出口厚度小于入口厚度

end

R=50;

l=sqrt(R*(h_0(i)-h_1(i))); %接触弧长 mm100;

xn=y(1:end-1)+diff(y)/2;%中性点坐标

% 预分配

h = zeros(1, m); % 变形区厚度

dt_h = zeros(1, m); % 压下量

l0 = zeros(1, m); % 来料长度

h_n = zeros(1, m); % 中性点厚度

h_c = zeros(1, m); % 平均厚度

c = zeros(1, m); % 计算量

y = linspace(0, b, m); % 横向坐标

for i=1:m

h(i)=h_1(i)+(h_0(i)-h_1(i)).*(xn(i)/l).^2;%变形区条元厚度%%%%%%

dt_h(i)=h_0(i)-h_1(i);% 条元的压下量

l0(i)=l/m;%条元来料长度

h_n(i)=(h_0(i)+h(i))/2;%条元中性点厚度

h_c(i)=(h_0(i)+h_1(i))/2;

c(i)=(dt_h(i)*l)/h_n(i);

end

T = zeros(1, m);

K_A = zeros(1, m);

%计算K值

for i=1:m

T(i)=1+(3*h_c(i)*(1-(mu)^2))/(2*E*dt_h(i));

K_A(i)=((8*tau*h_n(i)*(1-(mu)^2))/(E*h_c(i)*dt_h(i)*l*T(i))).*0.5;

end

%条元出口节线上的出口横向位移值的确定

A = zeros(1, m);

B = zeros(1, m);

C_A = zeros(1, m);

N = zeros(1, m);

T = zeros(1, m);

lambda=zeros(1,m);

e_b=E/((1-(mu)^2)*b);

% 构建系数矩阵和常数向量

A_A= zeros(m+1, m+1);%系数矩阵

C_C = zeros(m+1, 1);%常数向量

%填充系数矩阵和常数向量

A_A(1,1)=A(1)-e_b;

A_A(1,2)=B(1);

A_A(1,m)=e_b;

C_C(1)=-k+C_A(1);

for i=1:m

A(i)=N(i)*coth(K_A(i)*s(i))+lambda(i)*((K_A(i)*s(i))/(sinh(K_A(i)*s(i))^2));

B(i)=-lambda(i)*(K_A(i)*s(i)*cosh(K_A(i)*s(i)))/(sinh(K_A(i)*s(i))^2)-N(i)/(sinh(K_A(i)*s(i)));

C_A(i)=(E/(1-(mu)^2))*(1+(h_1(i)/h1)-(h_0(i)/h0)-(l0(i)/l0(i)));%%%%%

N(i)=(k*h_c(i)*K_A(i))/dt_h(i)+(4*tau*h_n(i))/(h_c(i)*dt_h(i)*l*K_A(i))+(E*K_A(i))/(2*(1-(mu)^2));

lambda(i)=(k*h_c(i)*K_A(i))/dt_h(i)+(4*tau*h_n(i))/(h_c(i)*dt_h(i)*l*K_A(i))+(E*K_A(i))/(2*(1-(mu)^2));

end

% (A(1)-e_b)*u(0)+B(1)*u(1)+e_b*u(m)=-k+C_A(1)

% B(j)*u(j-1)+(A(j)+A(j+1))*u(j)+B(j+1)*u(j+1)=C_A(j+1)-C_A(j)

% e_b*u(0)+B(m)u(m-1)+(A(m)-e_b)u(m)=k-C_A(m)

for j=1:m-1

A_A(j+1,j)=B(j);%第 j+1 个条元与第 j 个条元之间的相互作用。

A_A(j+1,j+1)=A(j)+A(j+1);%第 j+1 个条元自身的相互作用，通常涉及到当前条元和相邻条元的系数。

A_A(j+1,j+2)=B(j+1);% j+1 个条元与下一个条元（第 j+2 个）之间的相互作用。

C_C(j+1)=C_A(j+1)-C_A(j);%第 j+1 个条元上，由于物理条件变化（如力、压力等）导致的净效应。

end

A_A(m+1, 1) = e_b;

A_A(m+1, m-1) = B(m);

A_A(m+1, m+1) = A(m) - e_b;

C_C(m+1) = k - C_A(m);

% u = A_A \ C_C;

% disp('出口横向位移 u:');

% disp(u);

% 检查系数矩阵的条件数

rcondA_A = cond(A_A);

if rcondA_A > 1e10 % 如果条件数很大，矩阵可能接近奇异

u = pinv(A_A) * C_C;% 如果条件数较小，则使用伪逆求解

else

A_reg = A_A + 1e-10 * eye(size(A_A)); % 添加正则化项

u = A_reg \ C_C; % 使用正则化矩阵求解

end