global r0 dr0 pAtm;
pAtm = 1.01325e5; % standard atmospheric pressure, Pa.
r0 = 4.500e-6; % initial radius of the cavity, m.
dr0 = 0.000; % , m/s.
t0 = 0.000; % start time, s.
t1 = 40.000e-6; % end time, s.
dt = 0.010e-6; % time step, s.
%% ---------- Solve ---------- %%
[t, u] = ...
ode15s(@RayleighPlesset, ...
(t0 : dt : t1).', ...
[r0; dr0], ...
odeset('RelTol', 1.000e-6, ...
'AbsTol', 1.000e-12, ...
... 'NormControl', 'on', ...
'InitialStep', 0.010 * dt, ...
'MaxStep', 0.200 * dt));
r = u(:, 1);
dr = u(:, 2);
%% ========== Plot Result ========== %%
axesWidth = 0.820;
axesHeight = 0.800;
fontName = 'Times New Roman';
fontSize = 16;
textInterpreter = 'latex';
hFigure = ...
figure( ...
'WindowState', 'maximized' ...
);
%% ---------- Pressure Axis ---------- %%
hAxesPressure = ...
axes(hFigure, ...
'Position', [0.666667 * (1.000 - axesWidth), 0.100, axesWidth, axesHeight], ...
'FontName', fontName, ...
'FontSize', fontSize, ...
'Color', 'none', ...
'NextPlot', 'add', ...
'YAxisLocation', 'right');
hAxesPressure.XAxis.Visible = 'off';
hAxesPressure.YAxis.Color = 'red'; % <-- axis color, RED.
hAxesPressure.YLabel.Interpreter = textInterpreter;
hAxesPressure.YLabel.String = '$p$ (kPa)';
hLinePressureScratch = plot(hAxesPressure, t * 1.000e6, pAtm * ones(length(t), 1) * 0.001, ...
t * 1.000e6, DrivingPressure(t) * 0.001);
hLinePressureScratch(1).Color = 'none';
hLinePressureScratch(2).Color = 'none';
%% ---------- Velocity Axis ---------- %%
hAxesVelocity = ...
axes(hFigure, ...
'Position', [0.333333 * (1.000 - axesWidth), 0.100, axesWidth, axesHeight], ...
'FontName', fontName, ...
'FontSize', fontSize, ...
'Color', 'none', ...
'NextPlot', 'add', ...
'YAxisLocation', 'right');
hAxesVelocity.XAxis.Visible = 'off';
hAxesVelocity.YAxis.Color = 'blue'; % <-- axis color, BLUE.
hAxesVelocity.YLabel.Interpreter = textInterpreter;
hAxesVelocity.YLabel.String = '$\dot{R}$ (m/s)';
hLineVelocityScratch = plot(hAxesVelocity, t * 1.000e6, dr);
hLineVelocityScratch.Color = 'none';
%% ---------- Radius Axis ---------- %%
hAxesRadius = ...
axes(hFigure, ...
'Position', [0.333333 * (1.000 - axesWidth), 0.100, axesWidth, axesHeight], ...
'FontName', fontName, ...
'FontSize', fontSize, ...
'Box', 'on', ...
'Color', 'white', ... % <-- background color, WHITE.
'NextPlot', 'add', ...
'YAxisLocation', 'origin');
hAxesradius.XLabel.Interpreter = textInterpreter;
hAxesRadius.XLabel.String = '\it{t} \rm{({\mu}s)}';
hAxesRadius.YLabel.Interpreter = textInterpreter;
hAxesRadius.YLabel.String = '$R$ ($\mu$m)';
hLineRadiusScratch = plot(hAxesRadius, t * 1.000e6, r * 1.000e6);
hLineRadiusScratch.Color = 'none';
%%Plot Radius,Velocity,Pressure%%
hAxesPlot=...
axes(hFigure,...
'Position',hAxesRadius.Position,...
'FontName',fontName,...
'FontSize',fontSize,...
'Color','none',...
'NextPlot','add',...
'YAxisLocation','origin',...
'YLimMode','manual',...
'YLim',[0.000,1.000],...
'Visible','off');
hLineRadiusVisible=...
plot(hAxesPlot, hLineRadiusScratch.XData, (hLineRadiusScratch.YData - hAxesRadius.YLim(1))/(hAxesRadius.YLim(2) - hAxesRadius.YLim(1)));
hLineRadiusVisible.Color = 'black';
hLineRadiusVisible.LineStyle = '-';
hLineVelocityVisible = ...
plot(hAxesPlot, hLineVelocityScratch.XData, (hLineVelocityScratch.YData - hAxesVelocity.YLim(1))/(hAxesVelocity.YLim(2) - hAxesVelocity.YLim(1)));
hLineVelocityVisible.Color = 'blue';
hLineVelocityVisible.LineStyle = '-';
hLineAtmosphericPressureVisible = ...
plot(hAxesPlot, hLinePressureScratch(1).XData, (hLinePressureScratch(1).YData - hAxesPressure.YLim(1))/(hAxesPressure.YLim(2) - hAxesPressure.YLim(1)));
hLineAtmosphericPressureVisible.Color = 'red';
hLineAtmosphericPressureVisible.LineStyle = '--';
hLineDrivingPressureVisible = ...
plot(hAxesPlot, hLinePressureScratch(2).XData, (hLinePressureScratch(2).YData - hAxesPressure.YLim(1))/(hAxesPressure.YLim(2) - hAxesPressure.YLim(1)));
hLineDrivingPressureVisible.Color = 'red';
hLineDrivingPressureVisible.LineStyle = '-.';
legend(hAxesPlot, 'cavity radius', 'rate of cavity radius', 'farfield pressure', 'driving pressure');
hAxesPlot.Legend.Location = 'best';
hAxesPlot.Legend.FontSize = fontSize;
% 仿真参数
t = linspace(0, 0.005, 100); % 时间序列 (s)
R = r0*(1 + 0.5*sin(2*pi*200*t)); % 示例气泡半径变化 (正弦波动)
% 可视化参数
cMap = jet(256); % 颜色映射
pressureRange = [-2e5 2e5]; % 压力显示范围 (Pa)
%% 创建图形窗口
fig = figure('Color','white','Position',[100 100 1200 600]);
% 创建子图布局
ax1 = subplot(1,2,1); % 气泡形态显示
ax2 = subplot(1,2,2); % 压力场显示
%% 初始化气泡形态显示
axes(ax1)
bubble = animatedline('Color','k','LineWidth',2);
axis equal tight
xlabel('X (mm)'); ylabel('Y (mm)');
title('Bubble Contour')
grid on
%% 初始化压力场显示
axes(ax2)
[X,Y] = meshgrid(linspace(-2,2,50), linspace(-2,2,50));
pressureField = peaks(50)*1e4; % 示例压力场数据
hImage = imagesc(ax2, X(1,:), Y(:,1), pressureField);
axis xy equal tight
colormap(ax2, cMap);
caxis(pressureRange)
colorbar
xlabel('X (mm)'); ylabel('Y (mm)');
title('Pressure Field (Pa)')
%% 创建时间轴控件
timeSlider = uicontrol('Style','slider',...
'Position',[100 20 800 20],...
'Min',1,'Max',length(t),...
'Value',1,...
'SliderStep',[1 10]/(length(t)-1),...
'Callback',@updateFrame);
%% 动画更新函数
function updateFrame(src,~)
persistent hBubble hPressure
% 获取当前时间步
currentTimeStep = round(src.Value);
% 更新气泡轮廓
axes(ax1)
theta = linspace(0, 2*pi, 100);
x = R(currentTimeStep)*cos(theta);
y = R(currentTimeStep)*sin(theta);
clearpoints(bubble)
addpoints(bubble,x,y)
% 更新压力场(示例数据)
axes(ax2)
pressureField = (peaks(50) + 0.1*currentTimeStep)*1e4;
set(hImage,'CData',pressureField)
% 更新标题
title(ax1,sprintf('Bubble Contour (t = %.1f μs)', t(currentTimeStep)*1e6))
title(ax2,sprintf('Pressure Field (t = %.1f μs)', t(currentTimeStep)*1e6))
drawnow
end
%% RayleighPlesset方程
function du = RayleighPlesset(t, u)
try % FOR DEBUGGING
persistent gasConst mu b theta0;
persistent pVapor rFarfield pFarfield rhoLiquid eta sigma gamma;
persistent r0 dr0;
persistent pDryAir0 rhoDryAir0;
global pAtm;
if isempty(gasConst)
% All parameters are at temperature of 20℃.
gasConst = 8.314; % gas constant, J/(mol·K).
mu = 0.029; % mocular weight, kg/mol.
% Take nytrogen as an approximation for dry air.
% b = 0.03913e-3; % covolume, m3/mol. Nobel-Abel EOS is adopted.
b = 0.000; % ideal gas.
theta0 = 273.15 + 20.000; % initial temperature, K.
pVapor = 2338.8; % vapor pressure of water, Pa.
% pVapor = 0.000;
rFarfield = Inf; % radius at which properties of farfield are evaluated, m.
% It is usually set as Inf.
% rFarfield = 200.0e-6;
pFarfield = pAtm; % pressure at farfield, Pa.
rhoLiquid = 998.2; % density of water at, kg/m3.
eta = 1.002e-3; % dynamic viscosity of water, Pa·s.
% eta = 0.000; % inviscid.
sigma = 7.28e-2; % interfacial tension between the liquid and vapor phases of pure water, N/m.
% sigma = 0.000;
gamma = 1.000; % specific heat ratio of gas.
% gamma = 1.000 for isothermal process;
% gamma = 1.400 for adiabatic process.
r0 = 4.500e-6; % initial radius, m.
dr0 = 0.000; % initial rate of radius, m/s.
pDryAir0 = pFarfield ...
+ DrivingPressure(0.000) ...
+ 2 * sigma / r0 + 4 * eta * dr0 / r0 - pVapor; % pressre of dry air, Pa.
rhoDryAir0 = 1.000 / (gasConst / mu * theta0 / pDryAir0 + b / mu);
end % /* if isempty(pVapor) */
r = u(1); % radius of the cavity.
dr = u(2); % rate of radius of the cavity.
pDriving = DrivingPressure(t);
rhoDryAir = rhoDryAir0 * (r0 / r)^3;
pDryAir = pDryAir0 * ((1.000 / rhoDryAir0 - b / mu) / (1.000 / rhoDryAir - b / mu))^gamma;
ddr = ( ...
( ...
( ...
((pDryAir + pVapor) - (pFarfield + pDriving) ...
- (4.000 * eta * dr + 2.000 * sigma) / r) ...
/ rhoLiquid) ...
- 0.500 * ((r / rFarfield)^4 - 1.000) * dr * dr) ...
/ (1.000 / r - 1.000 / rFarfield) - 2 * r * dr * dr) / r / r;
du = [dr; ddr];
catch exceptionDebug
end % BREAKPOINT HERE.
end % /* RayleighPlesset */
%驱动压
function pd = DrivingPressure(t)
persistent pAmplitude period omega;
if isempty(pAmplitude)
pAmplitude = 1.000 * 1.01325e5; % amplitude of driving pressure, Pa.
period = 40.000e-6; % period of driving pressure, s.
omega = 2 * pi / period; % circular frequency of driving pressure, rad/s.
end % /* if isempty(pAmplitude) */
if length(t) == 1
% if t <= period
pd = pAmplitude * cos(omega * t);
% else
% pd = pAmplitude;
% end % /* if t <= 40.000e-6 */
else
dataSize = size(t);
pd = zeros(dataSize);
t = t(:);
pd = pd(:);
for ii = 1 : length(t)
pd(ii) = DrivingPressure(t(ii));
end % /* for ii */
pd = reshape(pd, dataSize);
end
end % /* DrivingPressure */
请给我详细解析每行代码的含义
最新发布