|
clc;
clear all;
global f R0 P0 Pa Ki pi rho sigma mu w Pv v AL c r
f=0.256*10^5;
time=1./f;
ttime=2.*time;
R0=10*10^(-6);
P0=1.013*10^5;
Pa=3.77*10^5;
Ki=1.33;
pi=3.14;
rho=1000;
Pv=3540;
sigma=7.25*10^(-2);
mu=0.001;
c=1500;
r=3*R0;
w=2*pi*f;
AL=2*10^(-6);
v=2*pi*f*AL*cos(w*time);
[t,y]=ode45(@qipao,[0 ttime],[R0 0]);
R=y(:,1);dR=y(:,2);
figure(1)
plot(t,R/R0,'k-')
xlabel('t')
ylabel('R(t)/R_0')
figure(2)
dR1=dR(1:(end-1),1);
R=R(1:(end-1),1);
Psc=P0*(R0^3)/((R.^2)*r)+sigma*((R0^2)/((R.^2)*r)-3/R+R*(R0^2-(R.^2))/(r^4))+Pa+Pa*((R0^3+4*(R.^3))/(3*r*(R.^2))-((R0^3)*R+(R.^4))/3*(r^4));
t=t(1:(end-1));
plot(t,Psc,'k-')
xlabel('t')
ylabel('Psc(t)')
function dy=qipao(t,y)
global R0 P0 Pa Ki rho sigma mu w Pv v c
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=((P0+2*sigma/R0)*(R0/y(1))^(3*Ki)-2*sigma/y(1)-4*mu*y(2)/y(1)+Pv-(P0-Pa*sin(w*t)*1))/rho/y(1)-(3/2*((y(2)^2)+v^2))/y(1)-3*(P0+2*sigma/R0)*Ki*((R0/y(1))^(3*Ki))*y(2)/(rho*c*y(1));
end |
|