关于喷嘴计算的程序问题指教

daiweisc

自己编写的程序，运行后出现下列问题，请高手指教，程序见附件

Optimization terminated: first-order optimality is less than options.TolFun.
f =
         0   -1.4086
Warning: Default trust-region dogleg method of FSOLVE cannot
handle non-square systems; using Gauss-Newton method instead.
> In fsolve at 248
  In swirlnozzle at 35
f =
         0   -1.4086

f =
         0   -0.3986

f =
         0   -0.1218

f =
         0    0.6373

f =
         0   -0.1440

f =
         0   -0.1218

f =
         0   -0.1218

f =
         0    0.1749

f =
         0    0.0159

f =
  1.0e-003 *
         0   -0.8681

f =
  1.0e-003 *
         0   -0.8681

f =
  1.0e-003 *
         0   -0.8680

f =
  1.0e-004 *
         0    0.5751

f =
  1.0e-003 *
         0    0.4645

f =
  1.0e-009 *
         0   -0.1381

f =
  1.0e-009 *
         0   -0.1381

f =
  1.0e-007 *
         0    0.5120
Optimization terminated: directional derivative along
search direction less than TolFun and infinity-norm of
gradient less than 10*(TolFun+TolX).
f =
         0   -0.1066
Warning: Default trust-region dogleg method of FSOLVE cannot
handle non-square systems; using Gauss-Newton method instead.
> In fsolve at 248
  In swirlnozzle at 45
f =
         0   -0.1066

f =
         0    0.0257

f =
         0    0.0993

f =
  1.0e-003 *
         0   -0.3986

f =
  1.0e-003 *
         0   -0.3986

f =
  1.0e-003 *
         0   -0.3986

f =
  1.0e-004 *
         0   -0.6689

f =
  1.0e-004 *
         0    0.8004

f =
  1.0e-011 *
         0   -0.2870

f =
  1.0e-011 *
         0   -0.2870

f =
  1.0e-007 *
         0    0.1555
Optimization terminated: directional derivative along
search direction less than TolFun and infinity-norm of
gradient less than 10*(TolFun+TolX).
             单位：mm
排放口直径=6.2,修正后的入口直径=6.2旋流半径=4.97,喷射角=48.5    0.5616
m选的合适，不需要再调整，m=
    0.8000
参数合适，不需要调整
旋流室直径=16.2>>

[ 本帖最后由 eight 于 2008-3-10 16:39 编辑 ]

sigma665

function swirlnozzle
% xFsolve.m
% 喷嘴计算，采用《简式旋流喷嘴的设计》文章中的方法，作者“姚增权．马进”

clear all
clc
% -------------------------------------------------------------------------
% 初始参数设定
Kv = 0.8667;                                         % 流量系数Kv
Q = Kv/3600;                                         % 流量，单位m3/s
i = 3;                                               % 旋流入口数目                    
alpha = 60;                                          % 雾化角，alpha单位为度
m = 0.8;                                             % 旋流腔半径与排放口直径之比,m取0.5～3
a = 2*pi*alpha/360;                                  % a将alpha单位转为弧度
Pl = 1e5;                                            % 压力为100kPa
density = 1000;                                      % 水的密度
miu = 1.0034e-6;                                     % 水的运动粘度
x0 = 0.5;
% -------------------------------------------------------------------------
theta = fsolve(@NonlinEqs,x0,[],a);                  % 第一步计算填充系数theta
K = sqrt(2)*(1-theta)/(theta*sqrt(theta));           % 计算几何常数K
Cd = theta*sqrt(theta/(2-theta));                    % 计算流量系数Cd
% -------------------------------------------------------------------------
Dout = sqrt(4*Q/(pi*Cd*sqrt(2*Pl/density)));         % 第二步计算排放口直径Dout
% -------------------------------------------------------------------------
R = m*Dout;                                          % 第三步初步选定旋流半径R
Dp = sqrt(2*R*Dout/(i*K));                           % 计算喷嘴入口直径Dp
% -------------------------------------------------------------------------
Re = 4*Q/(pi*miu*sqrt(i)*Dp);                        % 第四步计算雷诺数Re
% -------------------------------------------------------------------------
s = 10^(25.8/(log10(Re))^2.58-2);                    % 第五步计算摩擦系数
Ky = R*(Dout/2)/(i*(Dp/2)^2+(s/2)*R*(R-Dout/2));     % 计算修正后的几何常数
% -------------------------------------------------------------------------
% 迭代计算填充系数theta1、流量系数Cd1、排放口直径Dout1、旋流半径R1、喷嘴入口直径Dp1
theta1 = fsolve(@NonlinEqs1,x0,[],Ky);
Cd1 = theta1^1.5/sqrt(2-theta1);
Dout1 = sqrt(4*Q/(pi*Cd1*sqrt(2*Pl/density)));
R1 = m*Dout1;
K1 = sqrt(2)*(1-theta1)/(theta1^1.5);      
Dp1 = sqrt(2*R1*Dout1/(i*K1));
% -------------------------------------------------------------------------
% 计算修正后的入口直径Dp2
Dp2 = Dp1/sqrt(0.9);                                    % 修正后的入口直径Dp2
% -------------------------------------------------------------------------
S = fsolve(@NonlinEqs2,x0,[],Cd1,K1);                   % 计算空气芯半径
spray = 2*atan((2*Cd1*K1)/sqrt((1+S)^2-4*(Cd1*K1)^2));  % 计算喷射角
% -------------------------------------------------------------------------
disp('             单位：mm')
fprintf('排放口直径=%.1f,修正后的入口直径=%.1f',Dout1*1000,Dp2*1000)
fprintf('旋流半径=%.2f,喷射角=%.1f',R1*1000,spray/(2*pi)*360)
disp(Cd1)
% -------------------------------------------------------------------------
% 最后参数校验
% 旋流腔半径与排放口直径之比m的校验
if(R1>(Dp2/2))
    disp('m选的合适，不需要再调整，m=');
    disp(m);
else disp('m选的不合适，需要再调整')
end
% 旋流室直径Dk的校验
Dk = 2*R1+Dp2;

if(or(Dk<(2.5*Dout1),Dk>(3.3*Dout1)))
    disp('参数需要调整')
else disp('参数合适，不需要调整')
    fprintf('旋流室直径=%.1f',Dk*1000)
end
% -------------------------------------------------------------------------
% 程序结束
% -------------------------------------------------------------------------


% ------------------------------------------------------------------
function f = NonlinEqs(x,a)
f(1) = tan(a/2)-(2^1.5*(1-x)/(sqrt(x)*(1+sqrt(1-x))));
% -------------------------------------------------------------------
function f = NonlinEqs1(x,Ky)
f(2) = Ky-sqrt(2)*(1-x)/(x^1.5)
% -------------------------------------------------------------------
function f = NonlinEqs2(x,Cd1,K1)
f(2) = Cd1-(sqrt(1-(Cd1*K1)^2)-x*sqrt(x^2-(Cd1*K1)^2)-(Cd1*K1)^2*log((1+sqrt(1-(Cd1*K1)^2))/(x+sqrt(x^2-(Cd1*K1)^2))))



还有请具体说一下是出现什么问题,我看只有警告.

[ 本帖最后由 sigma665 于 2008-3-9 08:56 编辑 ]

daiweisc

这个警告会不会影响结果

sigma665

Warning: Default trust-region dogleg method of FSOLVE cannot
handle non-square systems; using Gauss-Newton method instead.

用了高斯牛顿法,方法不一样,可能精度会不一样