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

15 views (last 30 days)
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);

Answers (1)

Sebastian Kammer
Sebastian Kammer on 23 Nov 2016
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 Comment
Pirc1 Pirc2
Pirc1 Pirc2 on 23 Nov 2016
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

Sign in to comment.

Categories

Find more on Simscape Fluids in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!