1D Heat equation in Matlab with variable heat Flux at one side

I have a similar problem to the one formulated in the following question:

https://de.mathworks.com/matlabcentral/answers/1743940-1d-heat-equation-in-matlab-with-heat-flux-at-one-side

In my case, however, the heat flux density is not constant and is added to the geometry.

I am trying to implement Beck's nonlinear estimation method. For this I have to evaluate the equation in each iteration step once with the heat flux density q and once with q(1+e), where e = 10^-3, to calculate a sensitivity factor. As I understand it, the temperatures with the heat flow density q(1+e) should always be greater than those of q, but this is not always the case for me. I don't know if I have the wrong idea or if there is a bug in my code. I would be very grateful for any help!

This is my current code:

function pdepeSolverTest

xmesh = 0:.0001:0.014;
tmesh = -1:.05:20;
m = 0;
e = 10^-3;
lamb = 26.7;
rho = 7600;
cp = 496;
T0 = 300;
Tr = 300;


currqdt = 10000;
Tq = pdepe(m,@heatConduction,@heatConductionIni,@heatConductionBc,xmesh, tmesh);

currqdt = 10000 *(e+1);
Tqe = pdepe(m,@heatConduction,@heatConductionIni,@heatConductionBc,xmesh, tmesh);

function [pl,ql,pr,qr] = heatConductionBc(xl,ul,xr,ur,t)
    pl = currqdt;
    ql = 1;
    pr = ur - Tr;
    qr = 0;
end

function u0 = heatConductionIni(x)
u0 = T0;
end

function [c,f,s] = heatConduction(x,t,u,dudx)
c = rho*cp;
f = lamb*dudx;
s = 0;
end

end

In the current code, Tq is larger than Tqe in a few places.

Answer 1

[Torsten]

You are doing numerics here. The maximum difference is 1e-13, I think this is acceptable:

pdepeSolverTest()

ans = logical
   1

ans = 22

ans = 1.1369e-13

function pdepeSolverTest

xmesh = 0:.0001:0.014;
tmesh = -1:.05:20;
m = 0;
e = 10^-3;
lamb = 26.7;
rho = 7600;
cp = 496;
T0 = 300;
Tr = 300;


currqdt = 10000;
Tq = pdepe(m,@heatConduction,@heatConductionIni,@heatConductionBc,xmesh, tmesh,odeset('RelTol',1e-12,'AbsTol',1e-12));

currqdt = 10000 *(e+1);
Tqe = pdepe(m,@heatConduction,@heatConductionIni,@heatConductionBc,xmesh, tmesh,odeset('RelTol',1e-12,'AbsTol',1e-12));

any(Tq(:)-Tqe(:)>0)
nnz(Tq(:)-Tqe(:)>0)
max(Tq(:)-Tqe(:))

function [pl,ql,pr,qr] = heatConductionBc(xl,ul,xr,ur,t)
    pl = currqdt;
    ql = 1;
    pr = ur - Tr;
    qr = 0;
end

function u0 = heatConductionIni(x)
u0 = T0;
end

function [c,f,s] = heatConduction(x,t,u,dudx)
c = rho*cp;
f = lamb*dudx;
s = 0;
end

end

[Elias B.]

Thank you very much for the quick reply. The difference is 1e-4 in another place, but if the code is correct, it's probably not because of that. Then I'll probably have to look for my error in my code of the Beck's method...

[Torsten]

For which location and at which time do you want to compute the sensitivity of temperature on heat flux density ?

[Elias B.]

I would like to use Beck's non-linear estimation method to calculate the surface temperature of a die. For technical reasons, this is measured using 3 thermocouples at different positions below the surface. For this purpose, the heat flux densities are adjusted iteratively for each time step until the calculated temperatures at the measuring points match the real temperature. For this purpose, the sensitivity is calculated in each time step and a delta q by which the heat flux density at the current time is adjusted for the next iteration step until the error is small enough, then q is adjusted for the next time step.

Where n is the current time, k is the number of real measurement points, m is the number of future measurement points (3-4) to be included in the calculation and l is the current iteration step. T denotes the calculated temperatures and Y the measured temperatures.

If you are interested, I can share my complete code, but I don't think this is necessarily a Matlab problem :)