%*********************************
%作用:求各层缠绕系数
%*********************************
clc;
clear;
Es=6*10^10; %定义Es
a=6.67; %定义a
R=156.7; %定义R
rr=60*pi/180; %定义γ伽玛
q=3.7; %定义δ德尔塔
E=210*10^9; %定义E
e=4; %定义e
Vr=0.2; %定义Vr
c=atan(0.16); %定义ζ捷塔
Esq=100*10^6; %定义Esq
d=4; %定义d
S=1000; %定义S
nz=10; %定义nz(nz>=2)
% syms Es a R rr q E e Vr c Esq d S a b
A(1)=1;
for i=1:nz-1
She(i)=0;
end
for n=2:nz
for i=1:n-1
syms (['x',num2str(i)]); %定义n个符号变量,分别是x1,x2...x(n-1),分别代表S1n,S2n...S(n-1)n
end
fh=[];
for i=1:n-1
fh=[fh,sym(['x',num2str(i)])];%把符号放在矩阵中,方便引用
end
C0=R/(e*q*E);
C(1)=(1+sin(rr))*(2-Vr*cot(rr+c))/(2*Esq*(2*R+d));
% C0=a;
% C(1)=b;
for i=2:n-1
C(i)=2*sin(rr)*(2-Vr*cot(rr+c))/(Esq*(2*R+(1+2*(i-1)*sin(rr))*d));
% C(i)=1;
end
temp=S;
for k=1:n-1
temp=temp-fh(k);
end
zz=C0*temp;
if n==2
m=2*Es*a/(2*R+d);
jie.x1=(m*C0*S+0.5*m*C(1)*S)/(1+m*C0);
else
for i=1:n-1 %代表xi中的i
if i==1
temp=S;
for k=2:n-1
temp=temp-fh(k);
end
te(1)=1/2*C(1)*temp;
else
for j=1:i-1 %代表C(i)中的i
temp=S;
for k=j+1:n-1
temp=temp-fh(k);
end
te(j)=C(j)*temp;
end
if i+1<=n-1
temp=S;
for k=i+1:n-1
temp=temp-fh(k);
end
else
temp=S;
end
te(i)=0.5*C(i)*temp;
end
plus_te=zz;
for k=1:i
plus_te=plus_te+te(k);
end
result(i)=(2*Es*a/(2*R+(1+2*(i-1)*sin(rr))*d))*plus_te;
% result(i)=simple(result(i));
% result(i)=collect(result(i));
end
result_last=result-fh; %列出最后的方程
jie=solve(result_last); %接出方程组x1=jie.x1....xi=jie.xi
end
for k=1:n-1
She(k)=She(k)+jie.(['x',num2str(k)]);
end
xhe=0;
for k=1:n-1
xhe=xhe+jie.(['x',num2str(k)])/S;
end
A(n)=A(n-1)+(1-xhe);
end
A(nz) %A(j)对应于公式中的Aj