twophase.m

基于流函数法求解流体颗粒二相绕流问题的源码

function varargout=twophase(t,Y,flag,tspan,Y0)
switch flag
 case ''                     
   varargout{1} = f(t,Y);
 case 'init'                
   [varargout{1:3}] = fi(tspan,Y0);
 case 'events'             
   [varargout{1:3}] = fev(t,Y,Y0);
 otherwise
   error(['Unknown flag ''' flag '''.']);
 end
% -------------------------------------------------------------
function Yd= f(t,Y)
% c=2,pi=3.1416,l=2,w=192/pi,u0=14,r0=0.35,y0=0.3,R=1.47,A=300,a=2
A=3.2,Wy=0.8,R=1,k=2.72;
Yd= [Y(3);
 Y(4);
-A*Y(3)+A*(1+k*sin(pi*Y(5)))*sin(pi*Y(1))*cos(pi*Y(2))+R*pi/2*(1+k*sin(pi*Y(5)))*(Y(3)*cos(pi*Y(1))*cos(pi*Y(2))-Y(4)*sin(pi*Y(1))*sin(pi*Y(2)))+R*pi*(1+k*sin(pi*Y(5)))^2*sin(pi*Y(1))*cos(pi*Y(1))+1.5*R*pi*k*cos(pi*Y(5))*sin(pi*Y(1))*cos(pi*Y(2));
-A*Y(4)-A*(1+k*sin(pi*Y(5)))*cos(pi*Y(1))*sin(pi*Y(2))+R*pi/2*(1+k*sin(pi*Y(5)))*(Y(3)*sin(pi*Y(1))*sin(pi*Y(2))-Y(4)*cos(pi*Y(1))*cos(pi*Y(2)))+R*pi*(1+k*sin(pi*Y(5)))^2*sin(pi*Y(2))*cos(pi*Y(2))+1.5*R*pi*k*cos(pi*Y(5))*cos(pi*Y(1))*sin(pi*Y(2))+A*Wy;
1]

 
% -------------------------------------------------------------
function [ts,y0,options] = fi(tspan,Y0)
 ts=tspan;
 y0=Y0;
 options =odeset('Events','on','Reltol',1e-5,'Abstol',1e-5);;
% ------------------------------------------------------------
function [value,isterminal,direction]=fev(t,Y,Y0)
case1=Y(1);case2=1-Y(1);case3=Y(2);case4=1-Y(2);
value=[case1,case2,case3,case4];
direction=[0,0,0,0];
isterminal=[1,1,1,1];