radial test cases done, debugging the bubble thermal cases for FD

Author: mrodrig6

src/default_case.m

@@ -3,31 +3,30 @@
 bubtherm        = 0;        % 0 : polytropic assumption, 1: thermal PDE in bubble
 medtherm        = 0;        % 0 : cold fluid, 1: warm fluid assumption
 stress          = 0;        % 1 : NHKV, qKV, 2: linear Maxwell, Jeffreys, Zener, 3: UCM or OldB, 4: PTT, 5: Giesekus
-eps3            = -1;        % this value must be (0, 0.5]
+eps3            = 0;        % this value must be (0, 0.5]
 vapor           = 0;        % 0 : ignore vapor pressure, 1 : vapor pressure
 masstrans       = 0;        % mass transfer, default is no mass transfer 
 
 % solver options
-TFin            = 100e-6;     % final time (s)
+TFin            = 20e-6;     % final time (s)
 TVector         = [0 TFin];
 method          = 45;       % ode45 setting for the time stepper
 spectral        = 0;        % force spectral collocation solution
 divisions       = 0;        % minimum number of timesteps
-Nt              = 12;       % number of points in bubble, thermal PDE
-Mt              = 12;       % number of points outside of bubble, thermal PDE
+Nt              = 10;       % number of points in bubble, thermal PDE
+Mt              = 10;       % number of points outside of bubble, thermal PDE
 Nv              = 150;      % number of points outside of bubble, viscoelastic stress PDE     
 Lv              = 3;        % characteristic length for grid stretching, leave at 3
 Lt              = 3;        % characteristic length for grid stretching, leave at 3
 
 % initial conditions
-R0              = 2.447495043190468e-04;     % initial bubble radius
+R0              = 100E-6;   % initial bubble radius
 U0              = 0;        % initial velocity (m/s)
-Req             = 3.008409399929516e-05;%27E-6;   % Equilibrium radius for pre-stress bubble, see Estrada JMPS 2017    
+Req             = 50E-6;    % Equilibrium radius for pre-stress bubble, see Estrada JMPS 2017    
 
 % output options
 dimensionalout  = 0;        % output result in dimensional variables
 progdisplay     = 0;        % display progress while code running
-plotresult      = 0;        % generate figure containing results
 
 % acoustic options
 rho8            = 1064;%997;              % far-field density (kg/m^3)   

src/f_call_params.m

@@ -151,6 +151,10 @@
     C8 = sqrt(nstate*(P8 + GAM)/rho8); 
 end
 
+if ~medtherm && Mt >= 0
+    Mt = -1;
+    disp('RESETTING Mt = -1 for medtherm == 0');
+end
 check = 1-isnumeric(radial);
 if check || radial > 4 || radial <= 0 
     error('INPUT ERROR: radial must be 1, 2, 3, or 4');

src/f_display.m

@@ -28,7 +28,10 @@
         mass = 'no mass transfer in the bubble';
     end
 
-    if stress == 1
+
+    if stress == 0 
+        const = 'no stress applied';
+    elseif stress == 1
         if Ca == Inf
             const = 'Newtonian fluid';
         else

src/m_imrv2_finitediff.m

@@ -132,7 +132,7 @@
 p = X(:,3); 
 if bubtherm
     Tau = X(:,4:(Nt+3)); 
-    T = (A_star -1 + sqrt(1+2*Tau*A_star)) / A_star; % Temp in bubble
+    T = (alpha -1 + sqrt(1+2*Tau*alpha)) / alpha; % Temp in bubble
     if medtherm
         Tm = X(:,(Nt+4):(2*Nt+3)); 
     end
@@ -203,7 +203,7 @@
         Tau(end) = prelim;
         T = TW(Tau);
 
-        K_star = A_star*T+B_star; 
+        K_star = alpha*T+beta; 
         pVap = (f_pvsat(T(1)*T8)/P8); 
     else
         pVap = p0star;
@@ -238,9 +238,9 @@
               /( T(end)*R*Rmix(end)*(1-C(end)) ) );
 
         % temperature inside the bubble
-        U_vel = (chi/R*(kappa-1)*DTau-yk*R*pdot/3)/(kappa*p);
+        U_vel = (chi/R*(kappa-1)*DTau-y*R*pdot/3)/(kappa*p);
         first_term = (DDTau*chi/R^2+pdot)*( K_star*T/p*kapover);
-        second_term = -DTau*((1/R)*(U_vel-yk*U));
+        second_term = -DTau*((1/R)*(U_vel-y*U));
 
         Taudot= first_term+second_term; 
         Taudot(end) = 0;
@@ -262,16 +262,16 @@
         % internal pressure equation
         pdot = 3/R*(chi*(kappa-1)*DTau(end)/R-kappa*p*U);
         % temperature inside the bubble
-        U_vel = (chi/R*(kappa-1).*DTau-yk*R*pdot/3)/(kappa*p);
+        U_vel = (chi/R*(kappa-1).*DTau-y*R*pdot/3)/(kappa*p);
         first_term = (DDTau.*chi./R^2+pdot).*( K_star.*T/p*kapover );
-        second_term = -DTau.*(1/(R).*(U_vel-yk*U));
+        second_term = -DTau.*(1/(R).*(U_vel-y*U));
 
         Taudot= first_term+second_term; 
         Taudot(end) = 0;
 
     else
         % polytropic gas
-        p = p0star*(1/R)^(3*kappa);
+        p = (p0star-Pv_star)*R^(-3*kappa);
         pdot= -3*kappa*U/R*p;
         pVap = Pv_star;
     end
@@ -338,7 +338,7 @@
 
 function Tw= TW(Tauw)
     %calculates the temperature at the bubble wall as a fuction of \tau 
-    Tw = (A_star -1 + sqrt(1+2*Tauw*A_star)) / A_star;
+    Tw = (alpha -1 + sqrt(1+2*Tauw*alpha)) / alpha;
 end
 
 function Cw= CW(Tw,P)

toolchain/f_generate_test.m

@@ -1,8 +1,6 @@
 clc; clear all; close all;
 
 addpath('../src');
-addpath('../src/spectral/');
-addpath('toolchain/');
 addpath('../unit_tests/');
 load('file_ids.mat');
 
@@ -13,26 +11,32 @@
 for id = testb:teste
     filename = strcat('../unit_tests/',ids{id},'.dat'); 
     [ts,Rs,Us] = m_imrv2_spectral('radial',id);
-    [tf,Rf,Uf] = m_imrv2_spectral('radial',id);
+    [tf,Rf,Uf] = m_imrv2_finitediff('radial',id);
     display(filename)
     if (norm(Rs-Rf,2) < 1E-15)
-        disp('success')
+        disp('----> SUCCESS! <------')
         save(filename,"ts","Rs","Us")
     else
         disp('error radial not working')
     end
 end
-%%
+
 testb = testb + 1;
 teste = teste + 1;
 % bubtherm tests
 for id = testb:teste
     filename = strcat('../unit_tests/',ids{id},'.dat'); 
-    [t,R,U] = m_imrv2_spectral('bubtherm',1);
+    [ts,Rs,Us] = m_imrv2_spectral('bubtherm',1);
+    [t,R,U] = m_imrv2_finitediff('bubtherm',1);
     display(filename)
-    save(filename,"t","R","U")
+    if (norm(Rs-Rf,2) < 1E-15)
+        disp('----> SUCCESS! <------')
+        save(filename,"ts","Rs","Us")
+    else
+        disp('error radial not working')
+    end
 end
-
+%%
 testb = testb + 1;
 teste = teste + 1;
 % medtherm tests