Skip to content

Commit 99a4d30

Browse files
mrodrig6mrodrig6
committed
loaded Bachir_IMR to the repository
1 parent 19f6bb5 commit 99a4d30

File tree

9 files changed

+1879
-0
lines changed

9 files changed

+1879
-0
lines changed

‎src/bachir_imr/Finite_diff_mat.m‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
function [Diff_Matrix ] = Finite_diff_mat(Nodes,order,Tm_check)
2+
% Creates finite diffrence matrices
3+
% Nodes: Number of nodes
4+
% order: order of differentiation
5+
% Tm_check: 0 not for external temp , 1 used for ext. temp
6+
7+
8+
if Tm_check == 0
9+
10+
N = Nodes-1;
11+
deltaY = 1/N;
12+
K = 1:1:N+1;
13+
yk = (K-1)*deltaY;
14+
15+
elseif Tm_check == 1
16+
17+
N = Nodes-1;
18+
deltaY = -2/N;
19+
K = 1:1:N+1;
20+
yk = 1+(K-1)*deltaY;
21+
22+
end
23+
24+
%Diff_Matrix = zeros(Nodes,1); % Zhiren: pre-allocation incorrect.
25+
Diff_Matrix = zeros(Nodes);
26+
27+
if order == 1
28+
29+
30+
31+
for counter=2:(N) %in between
32+
33+
Diff_Matrix(counter,counter+1) = 1/2 ;
34+
Diff_Matrix(counter,counter-1) = -1/2 ;
35+
36+
end
37+
38+
if Tm_check == 0
39+
Diff_Matrix(end,end) = 3/2;
40+
Diff_Matrix(end,end-1) = -2;
41+
Diff_Matrix(end,end-2) = 1/2;
42+
43+
end
44+
45+
Diff_Matrix = Diff_Matrix / deltaY ;
46+
47+
elseif order == 2
48+
49+
50+
for counter=2:(N) %in between
51+
52+
if Tm_check == 0
53+
54+
Diff_Matrix(counter,counter+1) = 1+deltaY/yk(counter) ;
55+
Diff_Matrix(counter,counter) = -2 ;
56+
Diff_Matrix(counter,counter-1) = 1-deltaY/yk(counter) ;
57+
58+
elseif Tm_check == 1
59+
60+
Diff_Matrix(counter,counter+1) = 1 ;
61+
Diff_Matrix(counter,counter) = -2 ;
62+
Diff_Matrix(counter,counter-1) = 1 ;
63+
end
64+
end
65+
66+
if Tm_check == 0
67+
68+
Diff_Matrix(1,1)= -6; Diff_Matrix(1,2) = 6;
69+
70+
end
71+
72+
73+
Diff_Matrix = Diff_Matrix / (deltaY^2) ;
74+
end
75+
76+
% size(Diff_Matrix) % For debug - confirm previously pre-allocated incorrectly.
77+
78+
end

‎src/bachir_imr/IMRCalc_Req.m‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
function [R_eq,P_eq,C_eq] = IMRCalc_Req(R0,Tgrad,Cgrad,Pa,G,G1,mu)
2+
3+
%***********************************************************************
4+
% NOTE:
5+
% Used to estimate the equilibrium radii
6+
% It assumed that at equilibrium the medium is stress free, so
7+
% R_eq is needed to set up the Elastic force on the bubble f(R_eq/R)
8+
%***********************************************************************
9+
10+
%***************************************
11+
% Load Parameters :
12+
Pmt = IMRcall_parameters(R0,G,G1,mu); % Calls parameters script
13+
k = Pmt(1); We = Pmt(7); Rv_star = Pmt(11); Ra_star = Pmt(12);
14+
P_inf = Pmt(19); T_inf = Pmt(20);
15+
16+
%***************************************
17+
18+
Pv = Pvsat(1*T_inf)/P_inf; Pa = Pa/P_inf; ma0 = Pa/Ra_star;
19+
MTotal0 = Pa/Ra_star + Pv/Rv_star;
20+
21+
22+
if Cgrad == 0 && Tgrad == 0
23+
24+
fun = @(x) Pa *(1/x)^(3*k)-(1)+Pv-1/(We*x);
25+
26+
exp2 = 1;
27+
x = fzero(fun,exp2);
28+
29+
while (isnan(x))
30+
31+
exp2 = exp2*1.01;
32+
x = fzero(fun,exp2);
33+
34+
end
35+
36+
R_eq = x;
37+
P_eq = Pa*(1/R_eq)^(3*k)+Pv;
38+
theta = Rv_star/Ra_star*(P_eq/Pv-1);
39+
C_eq = 1/(1+theta);
40+
41+
elseif Cgrad == 0 && Tgrad == 1
42+
43+
% Note: For very small initial mass content this case has a hard time
44+
% finding roots. One fix is to just set Cgrad to 1. That always
45+
% converges
46+
%fun = @(x) Pa*(1/x)^(3)-(1)+Pv-1/(We*x);
47+
fun = @(x) Pa + (-1+Pv)*x^3 -x^2/(We);
48+
exp2=1;
49+
x = fzero(fun,exp2);
50+
while (isnan(x))
51+
exp2 = exp2*1.01;
52+
x = fzero(fun,exp2);
53+
end
54+
R_eq = x;
55+
P_eq = Pa*(1/ R_eq )^(3)+Pv;
56+
theta = Rv_star/Ra_star*(P_eq/Pv-1);
57+
C_eq = 1/(1+theta);
58+
59+
elseif Cgrad == 1 && Tgrad == 1
60+
61+
% Full model:
62+
63+
fun = @(x) Pv*(1+(ma0/x)*(Ra_star/Rv_star))-1-...
64+
1/We*(Pv/(Rv_star*x))^(1/3) ; % parameterized function
65+
66+
67+
exp2=MTotal0;
68+
x = fzero(fun,exp2,optimset('display','off'));
69+
70+
while (isnan(x))
71+
72+
exp2 = exp2/1.11;
73+
x = fzero(fun,exp2,optimset('display','off'));
74+
75+
end
76+
MVE = x;
77+
R_eq =(Rv_star*MVE/Pv)^(1/3);
78+
P_eq = Pv*(1+ma0/MVE*Ra_star/Rv_star);
79+
theta = Rv_star/Ra_star*(P_eq/Pv-1);
80+
C_eq = 1/(1+theta);
81+
82+
83+
84+
end
85+
86+
end

‎src/bachir_imr/IMRcall_parameters.m‎

Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
function [P] = IMRcall_parameters(R0,G,G1,mu)
2+
3+
% Code to create parameter .mat file for RP_Cav to use
4+
5+
% Zhiren - Modified Feb. 2021, add some KVFD parameters
6+
% The idea is that we will spit out multiple parameters. Some will not be
7+
% used anyway.
8+
9+
% Parameters:
10+
11+
%0 < G < 10^6; % (Pa) Medium Shear Modulus
12+
%0 < mu < 10^3; % (Pa s) Viscocity
13+
alpha = 1-eps; % Power law viscosity model exponent
14+
G2 = 1E1; %G2 for the Yeoh model
15+
16+
A = 5.28e-5; % (W/m-K^2)Thermal Conductivity coeff
17+
B = 1.165e-2; % (W/m-K)Thermal Conductivity coeff
18+
D0 = 24.2e-6; %Diffusion Coeff m^2/s
19+
k = 1.4; % Ratio of Specific Heats
20+
S = 0.072;%56; % (N/m) Liquid Surface Tension
21+
T_inf = 298.15; % (K) Far field temp.
22+
P_inf = 101325; % (Pa) Atmospheric Pressure
23+
rho = 1064;%1060.2; % (Kg/m^3) Liquid Density
24+
Km = 0.55; %(W/m-K)Thermal Conductivity Medium
25+
Cp = 4.181e3; % Specific Heat Medium J/Kg K;
26+
Dm = Km /(rho*Cp) ; % Thermal Diffusivity m^2/s
27+
Ru = 8.3144598; % (J/mol-K) Universal Gas Constant
28+
Rv = Ru/(18.01528e-3); % (J/Kg-K) Gas constant vapor
29+
Ra = Ru/(28.966e-3); % (J/Kg-K)Gas constant air
30+
L = 2; % Strech variable to map domain outside the bubble
31+
L_heat = 2264.76e3; % (J/Kg) Latent heat of evaporation
32+
C = 1484;%1540;%1484; % sound speed (m/s)
33+
34+
% Intermidiate calculated variables
35+
36+
K_infy = A*T_inf+B;
37+
Rnondim = P_inf/(rho*T_inf);
38+
Uc = sqrt(P_inf/rho);
39+
Pv = Pvsat(T_inf );
40+
P0 = P_inf + 2*S/R0 ; % need to add Pv_sat at room temp
41+
theta = Rv/Ra*(P0-Pv)/Pv; % mass air / mass vapor
42+
C0 = 1/(1+theta);
43+
44+
45+
% Final non-dimensional variables
46+
47+
chi = T_inf*K_infy/(P_inf*R0*Uc);
48+
fom = D0/(Uc*R0);
49+
foh = Dm/(Uc*R0);
50+
Ca = P_inf/G;
51+
Re = P_inf*R0/(mu*Uc);
52+
We = P_inf*R0/(2*S);
53+
Br = Uc^2/(Cp*T_inf);
54+
CaY = P_inf/G2;
55+
56+
A_star = A*T_inf / K_infy;
57+
B_star = B / K_infy;
58+
Rv_star = Rv/Rnondim;
59+
Ra_star = Ra/Rnondim;
60+
P0_star = P0/P_inf;
61+
t0 = R0/Uc;
62+
L_heat_star = L_heat/(Uc)^2;
63+
Km_star = Km/K_infy;
64+
C_star = C/Uc;
65+
De = (mu/G1)*Uc/R0;
66+
67+
% ZZ - For KVFD: G is still spring, G1 is alpha, while mu = emu = vmu (for
68+
% now) Intermediate values will also be treated here, just for clarity
69+
% in future revision ...
70+
71+
FDCa = P_inf/mu;
72+
FDRe = Re;
73+
afd = G1; % Because alpha was used for power-law
74+
Ze = (FDCa^(1-afd))*(FDRe^afd);
75+
76+
77+
P = [k chi fom foh Ca Re We Br A_star...
78+
B_star Rv_star Ra_star P0_star t0 C0 L L_heat_star Km_star ...
79+
P_inf T_inf C_star De alpha CaY ...
80+
Ze afd]; % ZZ: Add FDKV params.
81+
82+
end

0 commit comments

Comments
 (0)