Labfans是一个针对大学生、工程师和科研工作者的技术社区。 论坛首页 | 联系我们(Contact Us)
MATLAB爱好者论坛-LabFans.com
返回   MATLAB爱好者论坛-LabFans.com > 工程数学软件 > MATLAB论坛
MATLAB论坛 一切MATLAB相关问题在此讨论。
回复
 
主题工具 显示模式
旧 2012-09-11, 11:24   #1
sas4444
初级会员
 
注册日期: 2011-01-26
帖子: 1
声望力: 0
sas4444 正向着好的方向发展
默认 关于Rayleigh-Plesset方程中的负数开你次方的问题

function Rx=RP(t,R)
P0=100000; %初始压强P0=0.1MPa
Pv0=2330; %饱和蒸汽压强,
T0=298.16; %初始温度T0=273.16+25
Hv=2257.2; %焓Hv=2257.2KJ/Kg
R0=20e-6; %初始半径R0=0.02m
rho=1000; %液体密度ρ=1000Kg/cm2
nu=1*10^-6; %运动粘度ν=1.0*10^-6m2/s
sigma=0.0725; %表面张力γ=0.0725N/m
m=0.018; %摩尔质量m=0.018Kg/mol
Rg=8.3144; %通用气体常数Rg=8.3144KJ/(Kg·mol·K)
alpha=R0/8.86; %常数α=R0/8.86
r=0.008; %距离r=8*10^-3m
k=1.4; %绝热系数κ=1.4
PG0=P0-Pv0+2*sigma/R0;
Rx=zeros(2,1);
Rx(1)=R(2);
Rx(2)=-((P0*exp(Hv*m/Rg*(1/T0-1/(T0*R0^(3*(k-1))/(Rx(1)^3-alpha^3)^(k-1))))+PG0*(T0*R0^(3*(k-1))/(Rx(1)^3-alpha^3)^(k-1)/T0)*(R0/R(1))^(3*k)-(P0-rho/2*(R(1)^2/r^2*R(2))^2))/rho-4*nu*R(2)/R(1)-2*sigma/(rho*R(1))-1.5*R(2)^2)/R(1);
end
然后用ODE45求解
[t,R]=ode45('RP',[0,200e-6],[10e-6,0]);
其中在运算过程中,(Rx(1)^3-alpha^3)^(k-1) 这一项中会出现负数,从而产生复数解,或者出现无穷大,求大侠们给以指点,如何避免这种情况
sas4444 当前离线   回复时引用此帖
回复


发帖规则
不可以发表新主题
不可以发表回复
不可以上传附件
不可以编辑自己的帖子

启用 BB 代码
论坛启用 表情符号
论坛启用 [IMG] 代码
论坛禁用 HTML 代码



所有时间均为北京时间。现在的时间是 21:20


Powered by vBulletin
版权所有 ©2000 - 2025,Jelsoft Enterprises Ltd.