Skip to content

Commit e96eb81

Browse files
mrodrig6mrodrig6
committed
updated the files
1 parent e038782 commit e96eb81

File tree

2 files changed

+123
-102
lines changed

2 files changed

+123
-102
lines changed

‎src/data_assimilation/f_new.m‎

Lines changed: 122 additions & 101 deletions
Original file line numberDiff line numberDiff line change
@@ -1,65 +1,77 @@
11
function [t,X,tau_del] = f_new(ti_star,tf_star,xi,vars,tau_del)
2-
3-
global NT Pext_type Pext_Amp_Freq disptime Tgrad Tmgrad ...
4-
comp t0 neoHook nhzen sls linkv k chi fom foh We Br A_star ...
5-
B_star Rv_star Ra_star L L_heat_star Km_star P_inf T_inf C_star ...
6-
De deltaY yk deltaYm xk yk2 Pv REq D_Matrix_T_C DD_Matrix_T_C ...
7-
D_Matrix_Tm DD_Matrix_Tm Ca Re
8-
9-
NT= vars{1};Pext_type=vars{2};Pext_Amp_Freq=vars{3};disptime=vars{4};
10-
Tgrad=vars{5};Tmgrad=vars{6};Cgrad=vars{7};comp=vars{8};t0=vars{9};
11-
neoHook=vars{10};nhzen=vars{11};sls=vars{12};linkv=vars{13};k=vars{14};
12-
chi=vars{15};fom=vars{16};foh=vars{17};We=vars{18};Br=vars{19};
13-
A_star=vars{20};B_star=vars{21};Rv_star=vars{22};Ra_star=vars{23};
14-
L=vars{24};L_heat_star=vars{25};Km_star=vars{26};P_inf=vars{27};
15-
T_inf=vars{28};C_star=vars{29};De=vars{30};deltaY=vars{31};yk=vars{32};
16-
deltaYm=vars{33};xk=vars{34};yk2=vars{35};Pv=vars{36};REq=vars{37};
17-
D_Matrix_T_C=vars{38};DD_Matrix_T_C=vars{39};D_Matrix_Tm=vars{40};
18-
DD_Matrix_Tm=vars{41}; tspan_star = vars{42}; NTM = vars{43}; rho = vars{44};
19-
R0 = vars{45}; fung = vars{46}; fung2 = vars{47}; fungexp = vars{48};
20-
fungnlvis = vars{49};
21-
22-
23-
%************************************************
24-
% March equations in time
25-
26-
Uc = sqrt(P_inf/rho);
27-
28-
Br = xi(2*NT+NTM+5);
29-
foh = xi(2*NT+NTM+6);
30-
31-
% Ca = P_inf./(xi(2*NT+NTM+7));
32-
% Re = (P_inf*R0)./((xi(2*NT+NTM+8)).*Uc);
33-
34-
Ca = 1./(xi(2*NT+NTM+7));
35-
Re = 1./((xi(2*NT+NTM+8)));
36-
37-
38-
De = xi(2*NT+NTM+9);
39-
alpha = (xi(2*NT+NTM+10));
40-
lambda_nu = xi(2*NT+NTM+11);
41-
42-
43-
% restrictions on dynamics to keep quantities physical:
44-
45-
xi(3) = exp(xi(3));
46-
xi(5+NT:4+(2*NT)) = max(xi(5+NT:4+(2*NT)),0);
47-
48-
options = odeset('RelTol',1e-5,'AbsTol',1e-8);
49-
% [~ ,X] = ode23tb(@bubble, [ti_star tf_star], xi(1:end-2)',options);
50-
51-
[t ,X] = ode23tb(@bubble, [ti_star tf_star], xi(1:2*NT+NTM+4)',options);
52-
53-
xf = [X(end,:)';Br;foh;xi(2*NT+NTM+7:end)];
54-
xf(3) = log(xf(3));
55-
56-
57-
%%
58-
%*************************************************************************
59-
% Nested function; ODE Solver calls to march governing equations in time
60-
% This function has acess to all parameters above
61-
62-
function dxdt = bubble(t,x)
2+
3+
global NT Pext_type Pext_Amp_Freq disptime Tgrad Tmgrad ...
4+
comp t0 neoHook nhzen sls linkv k chi fom foh We Br A_star ...
5+
B_star Rv_star Ra_star L L_heat_star Km_star P_inf T_inf C_star ...
6+
De deltaY yk deltaYm xk yk2 Pv REq D_Matrix_T_C DD_Matrix_T_C ...
7+
D_Matrix_Tm DD_Matrix_Tm Ca Re
8+
9+
NT= vars{1};
10+
Pext_type=vars{2};Pext_Amp_Freq=vars{3};disptime=vars{4};
11+
Tgrad=vars{5};
12+
Tmgrad=vars{6};Cgrad=vars{7};comp=vars{8};t0=vars{9};
13+
neoHook=vars{10};
14+
nhzen=vars{11};sls=vars{12};linkv=vars{13};k=vars{14};
15+
chi=vars{15};
16+
fom=vars{16};foh=vars{17};We=vars{18};Br=vars{19};
17+
A_star=vars{20};
18+
B_star=vars{21};Rv_star=vars{22};Ra_star=vars{23};
19+
L=vars{24};
20+
L_heat_star=vars{25};Km_star=vars{26};P_inf=vars{27};
21+
T_inf=vars{28};
22+
C_star=vars{29};De=vars{30};deltaY=vars{31};yk=vars{32};
23+
deltaYm=vars{33};
24+
xk=vars{34};yk2=vars{35};Pv=vars{36};REq=vars{37};
25+
D_Matrix_T_C=vars{38};
26+
DD_Matrix_T_C=vars{39};D_Matrix_Tm=vars{40};
27+
DD_Matrix_Tm=vars{41};
28+
tspan_star = vars{42}; NTM = vars{43}; rho = vars{44};
29+
R0 = vars{45};
30+
fung = vars{46}; fung2 = vars{47}; fungexp = vars{48};
31+
fungnlvis = vars{49};
32+
33+
34+
%************************************************
35+
% March equations in time
36+
37+
Uc = sqrt(P_inf/rho);
38+
39+
Br = xi(2*NT+NTM+5);
40+
foh = xi(2*NT+NTM+6);
41+
42+
% Ca = P_inf./(xi(2*NT+NTM+7));
43+
% Re = (P_inf*R0)./((xi(2*NT+NTM+8)).*Uc);
44+
45+
Ca = 1./(xi(2*NT+NTM+7));
46+
Re = 1./((xi(2*NT+NTM+8)));
47+
48+
49+
De = xi(2*NT+NTM+9);
50+
alpha = (xi(2*NT+NTM+10));
51+
lambda_nu = xi(2*NT+NTM+11);
52+
53+
54+
% restrictions on dynamics to keep quantities physical:
55+
56+
xi(3) = exp(xi(3));
57+
xi(5+NT:4+(2*NT)) = max(xi(5+NT:4+(2*NT)),0);
58+
59+
options = odeset('RelTol',1e-5,'AbsTol',1e-8);
60+
% [~ ,X] = ode23tb(@bubble, [ti_star tf_star], xi(1:end-2)',options);
61+
62+
[t ,X] = ode23tb(@bubble, [ti_star tf_star], xi(1:2*NT+NTM+4)',options);
63+
64+
xf = [X(end,:)';
65+
Br;foh;xi(2*NT+NTM+7:end)];
66+
xf(3) = log(xf(3));
67+
68+
69+
%%
70+
%*************************************************************************
71+
% Nested function; ODE Solver calls to march governing equations in time
72+
% This function has acess to all parameters above
73+
74+
function dxdt = bubble(t,x)
6375
% Break x vector into indv. values
6476
R = x(1); % Bubble wall Radius
6577
U = x(2); % Bubble wall velocity
@@ -80,14 +92,15 @@
8092
if (Tmgrad == 1)
8193
if t/tspan_star> 0.001
8294
%Might need to tune 0.001 for convergence:
83-
84-
85-
guess= real(-.001 + tau_del(end)); %+tau_del(end);
86-
95+
96+
97+
guess= real(-.001 + tau_del(end));
98+
%+tau_del(end);
99+
87100
if isnan(guess) || isinf(guess)
88101
guess
89102
end
90-
103+
91104
prelim = fzero(@(x) Boundary(x,Tm,Tau,C,P),guess);
92105
else
93106
guess = -.0001;
@@ -142,7 +155,7 @@
142155
%if (Pext_type == 'sn')
143156
% Pext = -Pext_Amp_Freq(1)/P_inf*sin( Pext_Amp_Freq(2)*t*t0) ;
144157
% P_ext_prime = -Pext_Amp_Freq(2)*t0*Pext_Amp_Freq(1)/P_inf...
145-
% *cos( Pext_Amp_Freq(2)*t*t0) ;
158+
% *cos( Pext_Amp_Freq(2)*t*t0) ;
146159
%elseif (Pext_type == 'RC')
147160
% Pext = Pext_Amp_Freq(1)/P_inf ;
148161
% P_ext_prime = 0;
@@ -151,7 +164,7 @@
151164
% P_ext_prime = 0;
152165
%elseif (Pext_type == 'ip')
153166
% Pext = -Pext_Amp_Freq(1)/P_inf*...
154-
% (1-heaviside(t-Pext_Amp_Freq(2)/t0)) ;
167+
% (1-heaviside(t-Pext_Amp_Freq(2)/t0)) ;
155168
% P_ext_prime = 0;
156169
%end
157170
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -177,14 +190,14 @@
177190
DDC = DD_Matrix_T_C*C;
178191

179192
% Temp. field in the material
180-
193+
181194
DTm = D_Matrix_Tm*Tm;
182195
DDTm = DD_Matrix_Tm*Tm;
183196

184197
%***************************************
185198
% Internal pressure equation
186199
%pdot = 3/R*(Tgrad*chi*(k-1)*DTau(end)/R - k*P*U +...
187-
% + Cgrad*k*P*fom*Rv_star*DC(end)/( T(end)*R* Rmix(end)* (1-C(end)) ) );
200+
% + Cgrad*k*P*fom*Rv_star*DC(end)/( T(end)*R* Rmix(end)* (1-C(end)) ) );
188201
pdot = 3/R*(Tgrad*chi*(k-1)*DTau(end)/R - k*P*U +...
189202
+ Cgrad*k*P*fom*Rv_star*DC(end)/( R* Rmix(end)* (1-C(end)) )) ;
190203
% *****************************************
@@ -208,10 +221,10 @@
208221

209222
Tau_prime = first_term + second_term + third_term;
210223
% if Tmgrad == 0
211-
Tau_prime(end) = 0;
224+
Tau_prime(end) = 0;
212225
% else
213-
% Tau_prime(end) = K; %JY???
214-
% Tau_prime(end) = 0; % What is K?
226+
% Tau_prime(end) = K; %JY???
227+
% Tau_prime(end) = 0; % What is K?
215228
% end
216229
Tau_prime = Tau_prime*Tgrad;
217230

@@ -229,7 +242,8 @@
229242
% C_prime = C_prime*Cgrad;
230243

231244
% Vapor concentration equation
232-
U_mix = U_vel; % + fom/R*((Rv_star - Ra_star)./Rmix).*DC;
245+
U_mix = U_vel;
246+
% + fom/R*((Rv_star - Ra_star)./Rmix).*DC;
233247
one = DDC;
234248
% % JY!!! %
235249
%two = DC.*( -((Rv_star - Ra_star)./Rmix).*DC - DTau./(K_star.*T) );
@@ -254,7 +268,8 @@
254268
%third_term = 3*Br./yk2.^6.*(4/(3*Ca).*(1-1/R^3)+4.*U/(Re.*R)).*U./R;
255269
%Tm_prime = first_term+second_term+third_term;
256270
%Tm_prime(end) = 0; % Sets boundary condition on temp
257-
%Tm_prime(1) = 0; % Previously calculated;
271+
%Tm_prime(1) = 0;
272+
% Previously calculated;
258273
%Tm_prime = Tm_prime*Tmgrad; %Tmgrad makes this quantity zero
259274

260275
% Material temperature equations
@@ -266,7 +281,8 @@
266281

267282
Tm_prime = first_term+second_term+third_term;
268283
Tm_prime(end) = 0; % Sets boundary condition on temp
269-
Tm_prime(1) = 0; % Previously calculated;
284+
Tm_prime(1) = 0;
285+
% Previously calculated;
270286
Tm_prime = Tm_prime*Tmgrad; %Tmgrad makes this quantity zero
271287
%*****************************************
272288

@@ -357,25 +373,27 @@
357373
disp('Not finished in non-IC case for Fung model!')
358374
end
359375
% ====== JY!!! First order Fung G + first order mu model approx ======
360-
% % alpha = 0; lambda_nu = 0.001;
376+
% % alpha = 0;
377+
lambda_nu = 0.001;
361378
% Lv = 1;
362379
%
363380
% zeta = linspace(-1,0.99,200);
364381
%
365382
% tempr = R*( (2./(1-zeta)-1)*Lv + 1 );
366383
% tempr0 = (tempr.^3+R0^3-R^3).^(1/3);
367384
% gammadot = -0.5*( 2*(tempr0.^2)./(tempr.^3) + 1./tempr0 ) *R^2./(tempr.^2) * U;
368-
% % figure; plot(zeta,gammadot); pause;
385+
% % figure;
386+
plot(zeta,gammadot); pause;
369387
% % tempmu = 1/Re .* heaviside(-abs(gammadot)+1/lambda_nu) .* (1-lambda_nu^2*(gammadot.^2));
370388
% tempmu = 1/Re .* exp(-lambda_nu^2*(gammadot.^2));
371389
% tempS = -12*2*tempmu*U/R .* (1-zeta).^2 ./ (2+(1-zeta)*(1/Lv-1)).^4 /(Lv^3);
372390
%
373391
% S = -(1-3*alpha)*(5 - 4/Rst - 1/Rst^4)/(2*Ca) - ...
374-
% 2*alpha*(-27/40 - 1/8/Rst^8 - 1/5/Rst^5 -1/Rst^2 + 2*Rst)/(Ca) + ...
375-
% trapz(zeta,tempS);
392+
% 2*alpha*(-27/40 - 1/8/Rst^8 - 1/5/Rst^5 -1/Rst^2 + 2*Rst)/(Ca) + ...
393+
% trapz(zeta,tempS);
376394
%
377395
% Sdot = -2*U/R*(1-3*alpha)*(1/Rst + 1/Rst^4)/Ca - ...
378-
% 2*alpha*U/R*(1/Rst^8 + 1/Rst^5 + 2/Rst^2 + 2*Rst)/(Ca) + 4/Re*U^2/R^2;
396+
% 2*alpha*U/R*(1/Rst^8 + 1/Rst^5 + 2/Rst^2 + 2*Rst)/(Ca) + 4/Re*U^2/R^2;
379397
elseif fung2 == 1
380398
if strcmp(Pext_type,'IC')
381399
Rst = R/REq;
@@ -463,14 +481,15 @@
463481
udot = ((1+U/C_star)...
464482
*(P + abs(1-Cgrad)*Pv -1/(We*R) + S - 1 - Pext) ...
465483
+ R/C_star*(pdot+ U/(We*R^2) + Sdot -P_ext_prime ) ...
466-
- 1.5*(1-U/(3*C_star))*U^2)/((1-U/C_star)*R);%+JdotA/(C_star));
467-
484+
- 1.5*(1-U/(3*C_star))*U^2)/((1-U/C_star)*R);
485+
%+JdotA/(C_star));
486+
468487
%elseif fungexp==1
469488
% part1 = Sdot;
470489
% part2 = ((1+U/C_star)...
471-
% *(P + abs(1-Cgrad)*Pv -1/(We*R) + S - 1 - Pext) ...
472-
% + R/C_star*(pdot+ U/(We*R^2) + 0 -P_ext_prime ) ...
473-
% - 1.5*(1-U/(3*C_star))*U^2)/((1-U/C_star)*R);
490+
% *(P + abs(1-Cgrad)*Pv -1/(We*R) + S - 1 - Pext) ...
491+
% + R/C_star*(pdot+ U/(We*R^2) + 0 -P_ext_prime ) ...
492+
% - 1.5*(1-U/(3*C_star))*U^2)/((1-U/C_star)*R);
474493
% Sdot = (1+4/Re/C_star)^(-1) * (part1-4/Re/R*part2);
475494
% udot = (1+4/Re/C_star)^(-1) * (part2+R/C_star*part1);
476495
%end
@@ -482,28 +501,29 @@
482501
Sdot = Sdot - SdotA*udot/R; % JY!!! Pay attention to here!
483502
end
484503

485-
dxdt = [rdot; udot; pdot; Sdot; Tau_prime; C_prime; Tm_prime];
504+
dxdt = [rdot;
505+
udot; pdot; Sdot; Tau_prime; C_prime; Tm_prime];
486506
if isreal(rdot)==0
487507
pause;
488508
end
489509
end
490-
%*************************************************************************
491-
492-
% Other nested functions used to carry out repetetive calculations
493-
% throughout the code
510+
%*************************************************************************
511+
512+
% Other nested functions used to carry out repetetive calculations
513+
% throughout the code
494514
function Tw = TW(Tauw)
495515
%calculates the temperature at the bubble wall as a fuction of \tau
496516
Tw = (A_star -1 + sqrt(1+2*Tauw*A_star)) / A_star;
497517
end
498-
518+
499519
function Cw = CW(Tw,P)
500520
% Calculates the concentration at the bubble wall
501521

502522
%Function of P and temp at the wall
503523
theta = Rv_star/Ra_star*(P./(Pvsat(Tw*T_inf)/P_inf) -1);
504524
Cw = 1./(1+theta);
505525
end
506-
526+
507527
function Tauw = Boundary(prelim,Tm,Tau,C,P)
508528
% Solves temperature boundary conditions at the bubble wall
509529
% Create finite diff. coeffs.
@@ -527,21 +547,22 @@
527547
chi*(-coeff*[prelim ;T_trans] )/deltaY + Cgrad*...
528548
fom*L_heat_star*P*( (CW(TW(prelim),P)*(Rv_star-Ra_star)+Ra_star))^-1 *...
529549
(TW(prelim) * (1-CW(TW(prelim),P)) ).^(-1).*...
530-
(-coeff*[CW(TW(prelim),P); C_trans] )/deltaY;
531-
%************************************************************************
550+
(-coeff*[CW(TW(prelim),P);
551+
C_trans] )/deltaY;
552+
%************************************************************************
532553
end
533-
534-
% Gaussian pressure functions
535-
% acoustic pressure
554+
555+
% Gaussian pressure functions
556+
% acoustic pressure
536557
function p = pf(t)
537558
if t<(dt_star-5*tw_star) || t>(dt_star+5*tw_star)
538559
p=0;
539560
else
540561
p = -Pext_Amp_Freq(1)*exp(-(t-dt_star).^2/tw_star^2);
541562
end
542563
end
543-
544-
% time derivative of acoustic pressure
564+
565+
% time derivative of acoustic pressure
545566
function pdot = pfdot(t)
546567
if t<(dt_star-5*tw_star) || t>(dt_star+5*tw_star)
547568
pdot=0;
@@ -550,5 +571,5 @@
550571
end
551572

552573
end
553-
574+
554575
end

‎toolchain/s_test_suite_masstrans.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@
7777
if isempty(failed_tests)
7878
% success
7979
fprintf('✅ All tests PASSED.\n');
80-
exit(0);
80+
exit(0);
8181
else
8282
% fail the workflow
8383
fprintf('❌ Tests FAILED at indices: %s\n', sprintf('%d ', failed_tests));

0 commit comments

Comments
 (0)