MATLAB: Implementing correct heat flux boundary conditions

heatheat flux

Hi! I'm having quite a bit of trouble implementing the following boundary condition (BC) into a PDEPE solver for temperature change with respect to time inside a bubble:
The PDE is: dT/dt = alpha*d2T/dr2 – vbubble*dT/dr
BC: q = h*(T – Tinf) – m*Hv
Where:
r = radius
alpha = coefficient (constant)
vbubble = coefficient (constant)
q = heat flux
h = local heat transfer coefficient (constant)
T = temperature in bubble (variable)
Tinf = ambient temperature around bubble (constant)
m = evaporation flux into bubble (we can assume it to be constant for now)
Hv = enthalpy of vaporisation of water (constant)
My problem really arises when I run the code, as T should equal Tinf at the edge of the bubble (where radius = r = 1), but for some reason I'm getting the opposite. The code is as follows:
m = 2; % specifies spherical symmetry for PDE solver
r = linspace(0,r,100);
t = linspace(0,t,100);
sol_temp = pdepe(m,@mainfunction_temp,@initialconds_temp,@boundaries_temp,r,t);
T = sol_temp(:,:,1);
function [a,f,s] = mainfunction_temp(r,t,T,DuDr)
a = 1;
f = 0.001;
s = -vbubble;
end
function T0 = initialconds_temp(t)
T0 = Tgas_in;
end
function [pl,ql,pr,qr,tt] = boundaries_temp(rl,Tl,rr,Tr,t)
pl = 0;
ql = 1;
pr = h.*(Tr(1) – T_column) – m(i).*Hv;
qr = 1;
end
Another problem is that T is present in pr, and I'm just not sure how to incorporate the variable that I'm trying to solve for into the boundary condition, unless it's just for the first step where T = Tgas_in. I've therefore replaced it with Tr(1) for now instead of just T (as inputting T says it's a cleared variable), but I think that might be the source of the problem when it comes to getting the correct BC. Any suggestions?

Best Answer

m = 2; % specifies spherical symmetry for PDE solver
r = linspace(0,r,100); 
t = linspace(0,t,100); 
sol_temp = pdepe(m,@mainfunction_temp,@initialconds_temp,@boundaries_temp,r,t);
T = sol_temp(:,:,1);
 
    function [a,f,s] = mainfunction_temp(r,t,T,DuDr)
        a = 1;
        f = 0.001*DuDr;
        s = -vbubble*DuDr;
    end
 
    function T0 = initialconds_temp(t)
        T0 = Tgas_in;
    end
 
    function [pl,ql,pr,qr,tt] = boundaries_temp(rl,Tl,rr,Tr,t)
        pl = 0;
        ql = 1;
        pr = h*(Tr - T_column) - mdot*Hv;
        qr = 1;
    end
Related Question