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

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);

Answer 1

[Sebastian K]

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.

[Pirc1 Pirc2]

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