New code used in revised manuscript

clean_circular_turb_jokhul.m

@@ -0,0 +1,183 @@
+% Heat equation in radial coordinates (circular pipe) to understand Nusselt number in turbulent flow
+% Can set up for heated walls or with internal dissipation profile
+
+clear all
+close all
+
+Re_list=[1000000];
+Nu_list=zeros(size(Re_list));
+for Rei=1:length(Re_list)
+
+    r0=0.2;    % Pipe radius (m)
+
+    % Typical 0 deg conditions used in ice sheet subglacial modeling
+    k=0.558; % Thermal conductivity of water (W/m/K)
+    cw=4.22e3; %specific heat capacity of water (J/kg/K)
+    mu=1.787e-3; % Dynamic viscosity of water at 0 degC (Pa s)
+
+    rhow=1000; % Density of water (kg/m3)
+    kappa=k/(rhow*cw); % Thermal diffusivity
+    nu=mu/rhow; % Kinematic viscosity (m2/s)
+
+    Q=1; % Set to 0 to turn off internal heat generation, 1 for internal dissipation
+
+    dr=0.00001*r0;
+    J=ceil(r0/dr+1);
+    z=0:dr:r0;
+
+    % Get velocity profile
+    y=0:dr:r0;
+    Re=Re_list(Rei); % Prescribe bulk Reynolds number
+    vonK=0.4; % Von Karman constant
+    B=5; % Constant to use in log-law
+    Cf=fzero(@(Cf) -sqrt(2/Cf)+1/vonK*log(Re*sqrt(Cf)/2/sqrt(2))-1/vonK+B,0.009); % Friction factor (from log-law, A/A eqn 1.5)
+    u_b=Re*nu/(2*r0); % Bulk mean velocity -- A/A uses 2h, where h is half-width, but Re should depend on hydraulic radius?
+    tau_w=Cf*rhow*u_b^2/2; % Wall shear stress
+    u_tau=(tau_w/rhow)^0.5; % Friction velocity (inner velocity scale)
+    yplus=y.*u_tau./nu; % Non-dimensional y coordinate near wall
+
+    uplus(yplus<=20)=yplus(yplus<=20)-1.2533e-4*yplus(yplus<=20).^4+3.9196e-6.*yplus(yplus<=20).^5; % Velocity smooth profile (transition at yplus=20)
+    uplus(yplus>20)=1/vonK*log(yplus(yplus>20))+B; % Log-layer
+    u=flip(uplus)'.*(u_tau);
+
+    ub=2/r0^2*trapz(u.*z')*dr; % Mean bulk velocity for circular pipe
+
+    % Dissipation profile (based on Abe and Antonia, missing interior portion)
+    % depth-integrated E=<epsilon_bar>+<nu*(du/dz)^2>=(2.54.*log(yplus)+2.41).*u_tau^3;
+    % eps here epsilon_bar, NOT depth-integrated!
+    outer=find((r0-abs(z))/r0>0.2);
+    eps(outer)=(2.45./((r0-abs(z(outer)))./r0)-1.7).*u_tau^3/r0;
+    inner=find((r0-abs(z))/r0<=0.2);
+    eps(inner)=(2.54./((r0-abs(z(inner)))./r0)-2.6).*u_tau^3/r0;
+    % Inner portion of dissipation function (very near wall)
+    hplus=u_tau*r0/nu;
+    [ii,jj]=find(abs(y/r0-30/hplus)==(min(abs(y/r0-30/hplus))));
+    wall=find((r0-abs(z))/r0<30/hplus);
+    eps(wall)=eps(min(wall)-1);
+
+    % % Get eddy diffusivity
+    tau=tau_w.*z./r0; % Laminar shear stress distribution - linear
+    nu_T=-(tau(1:end-1)+tau(2:end))./2./(rhow.*diff(u')./dr)-nu; % Eddy viscosity (at half-nodes)
+    nu_T(nu_T<0)=0;
+    kappa_T=nu_T; % Eddy diffusivity from Reynolds' analogy (at half-nodes)
+
+    dx=r0;
+    nx=5000/dx;        % Number of x steps (iterations)
+
+    z_half=(z(1:end-1)+z(2:end))./2; % Half-node radial coordinates
+    alpha=(kappa+kappa_T').*dx./dr'.^2.*z_half';
+
+    % Define vectors (coefficients for interior nodes)
+    aconst=-alpha./2./z(2:end)';
+    bconst=(u(2:end-1)+alpha(2:end)./2./z(2:end-1)'+alpha(1:end-1)./2./z(2:end-1)');
+    cconst=-alpha./2./z(1:end-1)';
+
+    % Build a,b,c,r vectors
+    a=[0;aconst];
+    b=[0;bconst;0];
+    c=[cconst;0];
+    r=zeros(J,1);
+
+    % Boundary condition at pipe center (symmetry condition)
+    % Type II at z=0 (no-flux)
+    a(1)=0;
+    b(1)=-1;
+    c(1)=1;
+    r(1)=0;
+
+
+    if Q==0 % Type I at wall (heated)
+        a(J)=0;
+        b(J)=1;
+        c(J)=0;
+        r(J)=1;
+    end
+
+    if Q==1 % Type I at wall (T=0)
+        a(J)=0;
+        b(J)=1;
+        c(J)=0;
+        r(J)=0;
+    end
+
+    % % Type II at top (no-flux)
+    % a(J)=-1;
+    % b(J)=1;
+    % c(J)=0;
+    % r(J)=0;
+
+    % Initial condition u(x,0)=0
+    T=2*ones(J,1);
+
+    T_CN=[];    % Pre-allocate matrix to store plot values for each t
+
+    % Loop through distance along x
+    for it=1:nx
+
+        % Interior nodes
+        r(2:J-1)=alpha(1:end-1)./2./z(2:end-1)'.*T(1:J-2)+(u(2:J-1)-alpha(1:end-1)./2./z(2:end-1)'-alpha(2:end)./2./z(2:end-1)').*T(2:J-1)+alpha(2:end)./2./z(2:end-1)'.*T(3:J)+...
+            Q/rhow/cw*dx.*(eps(2:J-1)'.*rhow+mu.*((u(3:J)-u(1:J-2))./dr).^2);
+
+
+        Tnext=thomas(a,b,c,r);
+        T=Tnext';
+        %     Tbar(it)=trapz(T)*dr/r0;
+        Tbar(it)=2/(ub*r0^2)*trapz(u.*T.*z')*dr;
+
+        % Save to plot for selected downstream distances
+        if it*dx==1 | it*dx==5 | it*dx==10 | it*dx==20 | it*dx==50 | it*dx==100
+            T_CN=[T_CN T];
+        end
+
+%         % Brinkman number
+        Br_lam(it)=mu.*ub.^2./(k.*(T(end)-Tbar(it)));
+        Phi=2*pi.*((mu*trapz(((u(2:J)-u(1:J-1))'./dr).^2.*(z(2:J)+z(1:J-1))/2))+trapz(eps.*z)*rhow)*dr;
+        h_coeff=Phi/(2*pi*r0*(Tbar(end)-T(J)));
+        Br(it)=Phi./(2*h_coeff*(2-T(J)));
+
+
+    end
+
+    % Plotting - turn off for timing in problem 3
+    figure(1)
+    hold on
+    plot(T_CN,z./r0,'--','linewidth',1.2)
+    % axis([0 1 0 1])
+    legend('x=0.01','x=0.1','x=0.2','x=0.5','x=1','x=2','x=5','x=10','x=20','x=50','x=100','x=200','x=500','x=1000')
+    xlabel('T')
+    ylabel('r/r_{0}')
+
+    figure(2);hold on;
+    x=(0:(nx-1)).*dx;
+    plot(x./r0,log(Tbar-T(end)),'linewidth',2);
+    xlabel('x/r_0','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+    if Q==0 % For wall-heating case
+        p=polyfit(x,log(T(end)-Tbar),1);
+        Nu=-p(1)*rhow*cw*u_b*r0^2/k;
+    elseif Q==1  % For internal dissipation case -- with constant T=0 walls
+        Phi=2*pi.*((mu*trapz(((u(2:J)-u(1:J-1))'./dr).^2.*(z(2:J)+z(1:J-1))/2))+trapz(eps.*z)*rhow)*dr;
+        h_coeff=Phi/(2*pi*r0*(Tbar(end)-T(J)));
+        Nu=h_coeff*2*r0/k;
+    end
+
+    % Comparison with Siegel and Sparrow for internal heat generation (fig. 2)
+    % -- make sure to set wall b.c. as no-flux for comparison
+    SS=(T(J)-Tbar(end))./(Q*r0^2/k); % Siegel and Sparrow, fig. 2
+
+    % Plot Brinkman number
+    figure(3);hold on;plot(x./r0,Br)
+    xlabel('x/r_0','fontsize',14);ylabel('Brinkman number')
+
+    % Plot non-dimensional Brinkman
+    figure(4);hold on;plot(x./r0,log((Tbar-T(J))./(2-T(J))),'linewidth',2)
+    xlabel('x/r_0','fontsize',14);ylabel('ln((T_b-T_w)/(T_0-T_w))','fontsize',14)
+
+end
+
+
+
+
+
+
+

plot_ReNu_circular_2020revisions_DBonly.m

@@ -0,0 +1,30 @@
+clear all;close all
+
+
+% Circular pipe flow - wall heating - matches well with D-B correlation
+Re=[10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 200000 300000 400000 500000 750000 1000000];
+Nu=[89.2587 165.3526 236.3955 304.3768 370.2091 434.4064 497.2896 559.0866 619.9480 680.0172 1.2511e+03 1.7895e+03 2.3086e+03 2.8093e+03 4.0257e+03 5.1997e+03]; % With finer mesh, from circular_heatedwalls_data
+Pr=13.5146;
+% figure(1);hold on;plot(Re,Nu,'-o','linewidth',2)
+figure(1);semilogx(Re,Nu,'-o','linewidth',2)
+
+% Circular pipe flow with internal dissipation with Abe/Antonia dissipation profile - epsilon*rhow in Phi
+% calculation (compare Phi/rhow with E) - WITH NEAR-WALL EPSILON (CONSTANT from y/h<30/hplus)
+Re=[10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 200000 300000 400000 500000 750000 1000000];
+Nu=[43.3028 85.510 126.5970 166.9143 206.6587 245.9381 284.8349 323.4069 361.6920 399.7326 763.6548 1.1172e+03 1.4616e+03 1.8287e+03 2.6788e+03 3.5119e+03]; % With finer mesh, from circular_dissipation_data
+% figure(1);hold on;plot(Re,Nu,'-o','linewidth',2)
+figure(1);hold on;semilogx(Re,Nu,'-o','linewidth',2)
+
+Pr=13.5146;
+Re2=10000:10000:1000000;
+DB=0.024.*Re2.^0.8.*Pr^0.4;
+Re2=10000:10000:1000000;
+figure(1);hold on;semilogx(Re2,DB,'--','linewidth',2)
+
+Gni=0.012.*(Re2.^0.87-280).*Pr^0.4;
+figure(1);semilogx(Re2,Gni,'--','linewidth',2)
+
+figure(1);set(gca,'FontSize',18)
+axis([10000 1000000 0 6000])
+xlabel('Re');ylabel('Nu');
+legend('Heated Walls','Dissipation','Dittus-Boelter','Gnielinski')

plot_ReNu_combined_2020revisions.m

@@ -0,0 +1,35 @@
+clear all;close all
+
+% Circular pipe flow - wall heating - matches well with D-B correlation
+Re=[10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 200000 300000 400000 500000 750000 1000000];
+Nu=[89.2587 165.3526 236.3955 304.3768 370.2091 434.4064 497.2896 559.0866 619.9480 680.0172 1.2511e+03 1.7895e+03 2.3086e+03 2.8093e+03 4.0257e+03 5.1997e+03]; % With finer mesh, from circular_heatedwalls_data
+Pr=13.5146;
+figure(1);semilogx(Re,Nu,'-o','linewidth',2)
+
+% Circular pipe flow with internal dissipation with Abe/Antonia dissipation profile - epsilon*rhow in Phi
+% calculation (compare Phi/rhow with E) - WITH NEAR-WALL EPSILON (CONSTANT from y/h<30/hplus)
+Re=[10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 200000 300000 400000 500000 750000 1000000];
+Nu=[43.3028 85.510 126.5970 166.9143 206.6587 245.9381 284.8349 323.4069 361.6920 399.7326 763.6548 1.1172e+03 1.4616e+03 1.8287e+03 2.6788e+03 3.5119e+03]; % With finer mesh, from circular_dissipation_data
+figure(1);hold on;semilogx(Re,Nu,'-o','linewidth',2)
+
+% Restart color order
+ax = gca;
+ax.ColorOrderIndex = 1;
+
+% Channel flow - wall heating
+Re=[10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 200000 300000 400000 500000 750000 1000000]; 
+Nu=[162.4 304.0 436.5 563.6 686.8 807.1 925.0 1041 1155.3 1268.1 2.3430e+03 3.3585e+03 4.3391e+03 5.2859e+03 7.5885e+03 9.8136e+03]; % With finer mesh, from channel_heatedwalls_data
+
+figure(1);hold on;semilogx(Re,Nu,'-.x','linewidth',2)
+
+% Channel flow - with Abe/Antonia dissipation profile - epsilon*rhow in Phi
+% calculation (compare Phi/rhow with E) - WITH NEAR-WALL EPSILON (CONSTANT from y/h<30/hplus)
+Re=[10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 200000 300000 400000 500000 750000 1000000];
+Nu=[91.1845 178.9022 264.1219 347.7000 430.0565 511.4583 592.0762 672.0347 751.4468 830.0000 1.5861e+03 2.3193e+03 3.0395e+03 3.8071e+03 5.5886e+03 7.3389e+03]; % With finer mesh, from AA_dissipationintegral_data
+E=[0.0026 0.0173 0.0530 0.1176 0.2186 0.3632 0.5580 0.8099 1.1253 1.5105];
+Phi=[4.98 33.29 101.75 224.15 415.74 684.34 1048 1512.7 2097.6 2794.3];
+figure(1);semilogx(Re,Nu,'-.x','linewidth',2)
+figure(1);set(gca,'FontSize',18)
+axis([10000 1000000 0 10000])
+xlabel('Re');ylabel('Nu');
+legend('Circular conduit - heated walls','Circular conduit - dissipation','Sheet - heated walls','Sheet - dissipation')

plot_fig11a.m

@@ -0,0 +1,24 @@
+% Plot Figure 11
+clear all;close all;
+
+load Re1e5_jokhul_r2.mat
+
+figure(2);hold on;
+x=(0:(nx-1)).*dx;
+plot(x(1:nx/2),log(Tbar(1:nx/2)-T(end)),'linewidth',2);
+xlabel('x (m)','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+clear all;load Re1e5_heatedwalls_r2.mat
+figure(2);hold on;
+x=(0:(nx-1)).*dx;
+plot(x(1:nx/2),log(Tbar(1:nx/2)-T(end)),'--','linewidth',2);
+xlabel('x (m)','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+clear all;load Re1e5_dissipation_r2.mat
+figure(2);hold on;
+x=(0:(nx-1)).*dx;
+plot(x(1:nx/2),log(Tbar(1:nx/2)-T(end)),'--','linewidth',2);
+xlabel('x (m)','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+figure(2);legend('T_0=2 case','Heated wall case','Dissipation case','fontsize',14)
+% title('Re=10^6','fontsize',16)

plot_fig11d_smallr0.m

@@ -0,0 +1,39 @@
+% Plot Figure 11
+clear all;close all;
+
+load Re1e4_jokhul_r2.mat
+figure(2);hold on;
+x=(0:(nx-1)).*dx;
+plot(x(1:15000),log(Tbar(1:15000)-T(end)),'linewidth',2);
+xlabel('x (m)','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+clear all
+load Re5e4_jokhul_r2.mat
+figure(2);hold on;
+x=(0:(nx-1)).*dx;
+plot(x(1:15000),log(Tbar(1:15000)-T(end)),'linewidth',2);
+xlabel('x (m)','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+clear all
+load Re1e5_jokhul_r2.mat
+figure(2);hold on;
+x=(0:(nx-1)).*dx;
+plot(x(1:15000),log(Tbar(1:15000)-T(end)),'linewidth',2);
+xlabel('x (m)','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+clear all
+load Re5e5_jokhul_r2.mat
+figure(2);hold on;
+x=(0:(nx-1)).*dx;
+plot(x(1:15000),log(Tbar(1:15000)-T(end)),'linewidth',2);
+xlabel('x (m)','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+clear all
+load Re1e6_jokhul_r2.mat
+figure(2);hold on;
+x=(0:(nx-1)).*dx;
+plot(x(1:15000),log(Tbar(1:15000)-T(end)),'linewidth',2);
+xlabel('x (m)','fontsize',14);ylabel('ln(T_b-T_w)','fontsize',14)
+
+figure(2);legend('Re=10^4','Re=5x10^4','Re=10^5','Re=5x10^5','Re=10^6','fontsize',14)
+% title('Re=10^6','fontsize',16)