Implementing correct heat flux boundary conditions

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?

[Torsten]

I wonder how you can assume that the radius of the bubble remains constant if there is an evaporating flux at the boundary. And I further wonder how there can be a radial convective heat flux in a spherical object.

Look up "Stefan problem".

Best wishes

Torsten.

[Filip Butula]

Thanks for the info! Looked it up and the moving boundary seems to be what I'm looking for.

To be honest I agree that the radial convective heat flux is a strange addition, but it's what the literature model was using in the form of , so I wasn't entirely sure what to do other than include it.

Would you be able to comment on how the flux boundary could be indicated? As in how a boundary condition "q" can be expressed in the above equation, if "q" is not explicitly used?

Really appreciate the help, cheers!

Answer 1

[Torsten]

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

[Filip Butula]

I forgot to include the DuDr terms on f and s, of course.

I can't believe I didn't catch that, I really should have. Thanks again for helping out, I highly appreciate it!