回复 #15 donkeyxu 的帖子
正如donkeyxu所说的，解析解没有数值解好，而且对于非线性问题很难得到准确的解析解。比如我这个例子，只要将程序改为：syms va c
%P为螺距;
n=400;r=0.006;P=0.010;ha=0.002;
theta=atan(P/(pi*(r+r-ha)));
beita=0.2;
gamma=ha/c;
h1=c+ha;
h2=c;
w=2*pi*n/60;
B1=(1-beita)*P*cos(theta);
B2=beita*P*cos(theta);
u1=w*(r-c*gamma)*cos(theta)-va*sin(theta);
w1=w*(r-c*gamma)*sin(theta)-va*cos(theta);%槽底速度
u2=w*r*cos(theta)-va*sin(theta);
w2=w*r*sin(theta)-va*cos(theta);%槽面速度
SB1=(r-c*gamma)*(1-beita);%与2*pi*L无关，去掉2*pi*L
SB2=r*beita;%与2*pi*L无关，去掉2*pi*L
S=((r/2*sqrt(r^2+P^2/(4*pi^2))+(P^2/(4*pi^2)/2)*log(r+sqrt(r^2+P^2/(4*pi^2))))-((r-c*gamma)/2*sqrt((r-c*gamma)^2+P^2/(4*pi^2))+(P^2/(4*pi^2)/2)*log((r-c*gamma)+sqrt((r-c*gamma)^2+P^2/(4*pi^2)))))/P;%圈数等于长度除以螺距，去掉长度后剩下/P
dp1=(6*w1*w2/c^2)*(beita*gamma/((1-beita)*w1+(1+gamma)^3*w2));
dp2=-(6*w1*w2/c^2)*((1-beita)*gamma/((1-beita)*w1+(1+gamma)^3*w2));
Fx1=(-u1/h1)*SB1;
Fx2=(-u2/h2)*SB2;
Fz1=(-(h1/2)*dp1-w1/h1)*SB1;
Fz2=(-(h2/2)*dp2-w2/h2)*SB2;
Fzz=-dp1*B1*S;%此处用的是p的导数，所以应该有B1.此处正确不用改动，采用的是最大压力值。
Fa=(Fx1+Fx2)*sin(theta)-(Fz1+Fz2)*cos(theta)-Fzz*cos(theta);%加上正压力
%Ff=3*va*r/L;%考虑Ff就需要将Ff除以2*pi*miu*L；Ff=6*pi*va*miu*r;而考虑Ff就引入了长度L，因此此处不考虑Ff
F=Fa;
v=solve(F,'va');
ff=inline(v(2));
gg=fsolve(ff,0.0005)
即可求得解析解，但是不收敛，结果太差了。
只要将其转化为内联函数即可。
页:
1
[2]