Could I get help with modelling the Rayleigh-Plesset eqn for bubble growth

10 次查看(过去 30 天)
Can someone help in telling how to get the radius vs tim e graph?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%Yang and Church Validation %%%%%%%%%%%%%%%%%%%%%%
%Definition of Constants
syms R(t) pdel(t)
R0=1e-6; %Initial bubble radius
f=1e06; %Frequency of pressure wave
Pam=3e06;%Pressure amplitude
fc=1e06; %Frequency of pressure
t=0:1/100/fc:6/fc; %Time of presure application
P=Pam*sin(2*pi*fc*t);%Pressure
Pinf=101e03; %Fluid pressure
rho=1060; %Density of Fluid
c=1540; %Speed of sound in fluid
S=0.056; %Surface Tension
k=1.4; %Polytropic Index
u=0.015; %Viscosity
G=[0;1e06]; %Shear Modulus
Pg0=Pinf+(2*S/R0);
DR=diff(R);
pdel=Pg0*((R0/R)^(3*k))-(4*u*DR/R)-Pinf-P-(2*S/R);
test=(1-DR/c)*R*diff(DR,2)+1.5*(1-DR/(3*c))*(DR^2)-(1+DR/c)*((Pg0*((R0/R)^(3*k))-(4*u*DR/R)-Pinf-P-(2*S/R))/rho)-(R/(rho*c))*diff(Pg0*((R0/R)^(3*k))-(4*u*DR/R)-Pinf-P-(2*S/R))==0;
Sol=dsolve(test,DR(0)==0,R(0)==1e-06);

回答(1 个)

Sebastian K
Sebastian K 2016-11-23
Seems like "dsolve" is not able to solve this differential equation because the order of the differential equation does not match the number of initial conditions.
You have two symbolic function variables: "R" and "DR", where the latter is the first derivative of the former.
When you define your differential equation in the variable "test", notice that you also include a term "diff(DR,2)". This means that the differential equation contains the second derivative of "DR", which is also the third derivative of "R".
If this is not a typo and your equation actually contains the third derivative, then you would need to provide three initial conditions in order to solve this differential equation:
- One for the radius, R(0),
- One for the first derivative of radius, DR(0),
- And one for the second derivative of the radius [e.g., DDR(0), where DDR = diff(DR)].
Currently you are only providing the initial conditions for the radius and its first derivative, hence "dsolve" is not able to solve the differential equation.
  1 个评论
Pirc1 Pirc2
Pirc1 Pirc2 2016-11-23
Thanks for the suggestion. I've since modified/simplified the code to implement it using the ode 15s/45/23 solvers. However, I get the 'unable to meet tolerances' error. Could you please tell what's the problem with the code? Here's the new code. Also, there are only 2 derivatives
function dzdt = osc(t,z)
dzdt=zeros(2,1);
Radiusinit=(10^-6); %Initial bubble radius
Pam=1*(10^6);%Pressure amplitude
fc=3*(10^6); %Frequency of pressure
P=Pam*sin(2*3.14*fc*t);%Pressure
Pinf=101*(10^3); %Fluid pressure
density=1060; %Density of Fluid
c=1540; %Speed of sound in fluid
S=0.056; %Surface Tension
k=1.4; %Polytropic Index
u=0.015; %Viscosity
G=10^6; %Shear Modulus
Pg0=Pinf+(2*S/Radiusinit);
RHS1 = (1/ density)*(1 + (z(2))/c)*(Pg0*((Radiusinit/z(1))^(3*k)) - Pinf - P - (2*S/z(1)) - (4*u*z(2)/z(1))- (4*G/(3*((z(1))^3)))*((z(1))^3 -(Radiusinit^3)));
RHS2 = -1.5*(1 - z(2)/(3*c))*((z(2))^2);
RHS3 = (2*S/(density*c*z(1)))*z(2);
RHS4 = -(4*G*(Radiusinit^3)*((z(1))^-3)*z(2))/(density*c);
RHS5 = (Pg0/(density*c))*((z(1))^(-3*k))*(Radiusinit^3)* z(2);
RHS6 = -(Pam/(density*c))*2*3.14*fc*cos(2*3.14*fc*t)*z(1);
RHS7 = (((z(2))^2)*4*u/(density*c*z(1)));
RHS = (RHS1+RHS2+RHS3+RHS4+RHS5+RHS6+RHS7)/(((4*u/(density*c)) + (1-(z(2))/(density*c))*z(1)));
dzdt(1)=z(2); % get z(1)
dzdt(2)=RHS; % get z(2)
end
and this the caling function:
function [T,Y] = call_osc()
tspan=[0 10^-6]; % set time interval
[T,Y] = ode113(@osc,tspan,[10^-6, 0]');
plot(T,Y(:,1))
end

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Fluid Mechanics 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by