% A script to fit heat transfer coefficients based on a one term % approximation to the convection solution on a sphere % Created by Kristen Fauria in spring 2017 and modified in September 2017 %% clear all close all %%% Clast dimensions %%% %%% Choose clast and comment out other clast; %%% ML01 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% V =11.049/100^3; %particle volume (m^3) r = 0.01; % location of thermocouple (m) mass = 7.699/10^3; %mass (kg) phi = 0.69; %porosity %%% MLO2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% V=50.489/100^3; %particle volume (m^3) r=0.015; % location of thermocouple in (m) mass=38.896/10^3;% mass (kg) phi = 0.77; % porosity) %%% Pum06 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% V=6.816/100^3; %particle volume (m^3) r=0.008; % location of thermocouple in (m),check this mass=5.703/10^3;% mass (kg) phi = 0.63; % porosity) %%% SM16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% V=7.0332/100^3; %particle volume (m^3) r=0.008; % location of thermocouple in (m),check this mass=4.272/10^3;% mass (kg) phi = 0.74; % porosity) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculate clast radius and surface area rc = (3/4/pi*V)^(1/3); %effective particle radius (meters) RS=r/rc; %dimensionless distance (thermocouple position) A = 4*pi*rc^2; % effective particle surface area (m^2) rho = mass/V; %clast density (kg/,^3) %%% Preliminaries %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Heat capacities % cp = 820; %heat capacity of glass phase (J/kg K) cpW = 4184; % heat capacity of water (J/kg K) % Thermal Conductivities % %kr = 2.5; %thermal conductivity of rock (W/mK) (solid rock: 2-7; quartz: 3) kp = 1.1; %thermal conductivity of glass (W/ mK) kw = 0.58; % thermal conductivity of water (W/ mK) % Thermal diffusivities % dr=kp/(cp*2350);%thermal diffusivity of rock phase dw=kw/(cpW*1000);%thermal diffusivity of water phase dw=0.14*10^-6; % Effective diffusivities and conductivities (take porosity into accoutn) % keff=phi*kw+(1-phi)*kp; %effective thermal conductivity; deff=phi*dw+(1-phi)*dr; %effective diffusivity ~2.5*10^-7 m^2/s - this matches observation of foamy melt from Bagdassrov and Dingwell (1994) % Perhaps thermal diffusivity should be weighted by units of volume and % thermal conductivity by mass %MT=phi*1000+(1-phi)*2300;%total density %deff2=phi*1000/MT*dw+(1-phi)*2300/MT*dr;%weight by mass %keff2=phi*1000/MT*kw+(1-phi)*2300/MT*kp %% %%% Load data %%% % Manually change file name % %Data=textread(fullfile(pathname,'20170630_MLO1_500_dry_temp')); clear Data [filename,pathname]=uigetfile('*.txt'); Data=textread(fullfile(pathname,filename)); %Data=textread('20170630_ML01_500_dry');150 %% % Change which channel the temperature data is coming from depending on the % number of channels recorded. Time=Data(:,1); % time Temp=Data(:,2); % temperature (C) Tinf=mean(Data(:,3)); %background temperature degrees (C) Tinf=20; %Plot the data to manually select the "start" point figure; plot(Time,Temp); ylim([0, 700]); [temp1,temp2]=ginput(2);%select points to allow to zoom into fit xlim([temp1(1), temp1(2)]); [X,Y]=ginput(2); %select start and end points for the fit close %initial temperature (C) %%% Create time=zero %%% ind1 =find(Time>=X(1),1);%find first index where time >= X(1) ind2=find(Time>=X(2),1);%find end of the series Time2=Time(ind1:ind2)-X(1);%create time from zero Temp2=Temp(ind1:ind2); Ti=max(Temp2); %% % Compare theta to approximate solution predictions with different Biot % numbers clear Prediction1 Prediction2 Theta=(Temp2-Tinf)/(Ti-Tinf); % Table to iterate and find the right Biot number % first column is Bi, second is lambda and third is C1 Table1 = [1 1.5708 1.2732 2 2.0288 1.4793 3 2.2889 1.6227 4 2.4556 1.7201 5 2.5704 1.787 6 2.6537 1.8338 7 2.7165 1.8621 8 2.7654 1.9106 9 2.8044 1.9249 10 2.8363 1.9781 20 2.9857 1.9781 30 3.0372 1.9898]; Table2=[1.00001731 4.712392654 -0.424414357 1.083308736 4.73 -0.429971901 1.169152631 4.748 -0.435501085 1.26058984 4.767 -0.441173483 1.352940049 4.786 -0.44668074 1.461103694 4.808 -0.452855607 1.580732566 4.832 -0.459350389 1.702209246 4.856 -0.465599642]; Table3=[1.008 7.855 0.255090874 1.087 7.865 0.259421396 1.166 7.875 0.263716914 1.261 7.887 0.268825211 1.356 7.899 0.273882842 1.468 7.913 0.279719223 1.576 7.9265 0.285281557 1.701 7.942 0.291588364]; %deff=8*10^-7; Fo=deff.*Time2./(rc^2);%dimensionless time and fourier number %get rid of values with Fo<0.2 Start=find(Fo>0.2,1); Fo(1:Start-1)=[]; Theta(1:Start-1)=[]; Time3=Time2(Start:end); figure; for i=1:length(Table1) C1=Table1(i,3); L1=Table1(i,2); %C2=Table2(i,3); %L2=Table2(i,2); %C3=Table3(i,3); %L3=Table3(i,2); %Prediction1(i,:)=C1*exp(-L1^2.*Fo)./(L1*RS).*sin(L1*RS)+C2*exp(-L2^2.*Fo)./(L2*RS).*sin(L2*RS)+C3*exp(-L3^2.*Fo)./(L3*RS).*sin(L3*RS); % look at predictions from the Equation 5.50a Prediction2(i,:)=C1*exp(-L1^2.*Fo).*sin(L1*RS)./(L1*RS); hold on; %plot(Prediction1(i,:),Theta);hold on plot(Prediction2(i,:),Theta); end hold on; plot(Theta,Theta,'k'); xlabel('model predictions'); ylabel('observations'); set(gca,'FontSize',20); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Lumped capactitance fit %% % Create fit to data %LogFit=log((Temp(ind1:ind2)-Tinf)/(Ti-Tinf)); % %%plot and get linear fit to determine H %%% % figure;plot(Time2,Coeff*LogFit); % set(gca,'FontSize',20) % xlabel('Time (s)') % ylabel('ln[(T-T_i_n_f)/(T_i-T_i_n_f)]') % % p = polyfit(Time2,Coeff*LogFit, 1); % %p2= polyfit(LogFit,Time2/Coeff,1); % yfit = polyval(p,Time2); % hold on; %p1 is the slope % plot(Time2,yfit); % H=p(1);%Fitted heat transfer coefficient %units of W,(m2K), H also equals p2(1) % legend('observation','lumped capacitance model fit'); % biot=H*rc/kp; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% % Create Figure to show model fit % % Must manually put in results to calculate this one % % BestFit=5; % figure; % hold on % plot(exp(Fo),Theta,'b') % plot(exp(Fo),Prediction2(BestFit,:)); % set(gca,'FontSize',20); % xlabel('exp(Fo)'); % ylabel('(T(r,t)-T_\infty)/(T_i-T_\infty)'); % %ylabel('\theta(r,t)=\frac{T(r,t)-T_{\infty}}{T_i-T_{\infty}}'); % legend('observation','Equation XX model fit'); %% BestFit=9; figure; hold on Obs=log(Theta); Model=log(Prediction2(BestFit,:)); plot(Fo,Obs,'b') plot(Fo,Model,'k'); set(gca,'FontSize',20); xlabel('\tau'); ylabel('ln[(T(r,t)-T_\infty)/(T_i-T_\infty)]'); %ylabel('\theta(r,t)=\frac{T(r,t)-T_{\infty}}{T_i-T_{\infty}}'); legend('observation','Equation XX model fit');