-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathget_S0_SLS.m
125 lines (92 loc) · 3.25 KB
/
get_S0_SLS.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
% Initial stress calculation in linear Maxwell material
%
% Zhiren, Dec. 2024
% =========================================================================
% This is my implementation of Sawyer's strategy to estimate initial
% stress integral in a linear Maxwell material during bubble growth stage
% of LIC experiment.
% =========================================================================
function S_peak = get_S0_SLS(Lam,Re,De,Ca)
%% Parameter definition
% Lam = Amplification factor of bubble
% Ca = Cauchy number
% Re = Reynolds number
% De = Deborah number
params = struct;
params.Lam = Lam;
params.Ca = Ca;
params.Re = Re;
params.De = De;
%% Iterate for initial velocity
R0 = 1/Lam;
V0_guess = sqrt( (2/3) * (R0^(-3) - 1) ); % This uses Rayleigh's solution for void collapse
Rmax_err = @(X) abs(findRmax(X,params) - 1);
% V0_opt = fminsearch(Rmax_err,V0_guess);
% fm_options = optimset('Display','iter','PlotFcns',@optimplotfval);
[V0_opt, fval, exitflag, output] = fminsearch(Rmax_err, V0_guess); %, fm_options); % Get access to iteration info
%% Run simulation again to get stress integral at end:
% Set up other parameters
Rb0 = 1/(params.Lam); % Initial radius
pb0 = 1;
SI0 = 0; % Stress integral
tspan_star = 2.0; % Run long enough to reach end of growth
options = odeset('RelTol', 1E-7);
% Run simulation
X0 = [Rb0,V0_opt,pb0,SI0]';
[t_sol,X_sol] = ode23tb(@(t,X) KM_bubble(t,X,params), [0,tspan_star], X0, options);
R_sol = X_sol(:,1);
[Rmax,indx] = max(R_sol);
tol = 1E-3;
if abs(Rmax - 1) > tol
S_peak = nan;
else
S_peak = X_sol(indx,4);
% Deduct the NH part
S_NH = -(5 - Lam^(-4) -4/Lam)/(2*Ca);
S_peak = S_peak - S_NH;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SUBROUTINES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (A) Find Rmax for given initial velocity
function Rmax = findRmax(Vb0,params)
% Set up other parameters
Rb0 = 1/(params.Lam); % Initial radius
pb0 = 1;
SI0 = 0; % Stress integral
tspan_star = 2.0; % Run long enough to reach end of growth
options = odeset('RelTol', 1E-7);
% Run simulation
X0 = [Rb0,Vb0,pb0,SI0]';
[t_sol,X_sol] = ode23tb(@(t,X) KM_bubble(t,X,params), [0,tspan_star], X0, options);
% Get output
R_sol = X_sol(:,1);
Rmax = max(R_sol);
end
%% (B) Watered-down Keller-Miksis simulation
function dxdt = KM_bubble(t,X,params)
Ca = params.Ca;
Re = params.Re;
De = params.De;
Lmax = params.Lam;
% Parameters we won't touch for now. Define here directly:
kap = 1.4;
c_star = 1484/10.0751;
% Current values:
Rb = X(1); % Bubble radius
Vb = X(2); % Bubble wall velocity
pb = X(3); % Isobaric pressure
SI = X(4); % Stress integral
Lb = Lmax*Rb; % Current stretch
pbdot = -3*kap*pb*Vb/Rb;
SIdot_NH = -(2/Ca)*(Vb/Rb)*(Lb^(-4) + 1/Lb);
SI_NH = -(5 - Lb^(-4) - 4/Lb)/(2*Ca);
SIdot = -((SI - SI_NH) + (4/Re)*Vb/Rb)/De + SIdot_NH;
Rdot = Vb;
KM_RHS = (1 + Vb/c_star) * (pb + SI - 1) + (Rb/c_star)*(pbdot + SIdot);
KM_LHS = (3/2)*(1 - Vb/(3*c_star))*(Vb^2);
KM_denom = (1 - Vb/c_star)*Rb;
Vdot = (KM_RHS - KM_LHS)/KM_denom;
dxdt = [Rdot; Vdot; pbdot; SIdot];
end