Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Heat equation given constant surface heat flux

  1. Aug 16, 2015 #1
    How would I go about finding temperature distribution in a thin square plate during the the first few milliseconds (or actually a fraction of a millisecond) after t=0s. Initial temperature distribution throughout the plate is known, there is heat flux to one side = Qinj, while heat flux from all other sides = 0. I have seen how they do it for fixed temperatures at wide-end boundaries, but I am not sure how to adopt this method if instead constant flux is given. The equation is d2T/dx2 -a*dT/dt=0 (a is known).
    Last edited: Aug 16, 2015
  2. jcsd
  3. Aug 16, 2015 #2


    User Avatar
    Gold Member

    Suggestion: Divide the plate into a huge number of cubes, all with a volume = (dx)3.
    Supply heat to some area and do a numerical calculation of heat flux/temperatures through/in all the cubes. At the surfaces of the plate, the transfer heat flux will be zero as for the direction out of the plate, but the temperature may change. You are calculating with heat flux, thermal resistance and heat capacity. You may use Kirchhoffs current law by calculations ( an analogy ).

    That's how meteorologists do it, calculating a weather forecast.

    Maybe your computer can calculate these few milliseconds within a few hours?

    ( Remember high precision data/calculations and a small dt ).
    Last edited: Aug 16, 2015
  4. Aug 16, 2015 #3
    I'm not sure how to model it with Kirchoff laws. After more searches I'm trying to figure out how to make that Neumann boundary work. They only give examples for the Dirichlet condition (e.g. here). I'm trying to come up with a formula for Ux,t+1 such that Neumann conditions at boundaries are taken into consideration.
  5. Aug 17, 2015 #4
    I finally found an explanation that made sense. Thanks for your help.
  6. Aug 20, 2015 #5
    You can do it analytically by treating the plate as a semi-infinite slab. During the fraction of a millisecond of interest, the heat would only penetrate a small fraction of the thickness of the plate anyway. So you have a semi-infinite slab being heated on one side by a constant heat flux. See Conduction of Heat in Solids by Carslaw and Jaeger.

  7. Aug 20, 2015 #6
    Yeah, thanks for the link. I think you have a good point. I just need a look at first few molecular layers anyway, which could decrease my array size significantly. I mean if I use backward Euler, the time-steps have to be so small or it will be non-convergent. Just want to examine the rate of heat buildup vs rate of heat pass-over to further layers to determine the depth of surface melting during friction (friction is quite high and heat conductivity is very low).
  8. Aug 21, 2015 #7
    I think you are thinking of forward Euler. Backward Euler will always converge.

  9. Aug 21, 2015 #8
    Yeah, I've read that somewhere. You're the second person mentioning that. Hey, maybe you also have an idea of how to encode a non-zero Neumann boundary using the finite difference. What I got is Ux0,t1 = Ux0,t0 + boundary_constant_heat_flux_value * delta_t. And in calculations, before relooping the code, just add this to whatever the backward Euler calculated for Uxn,t1?
  10. Aug 22, 2015 #9
    A better way is to actually solve the differential equation at the boundary also. This is done by introducing a "fictitious" grid point at x = -Δx, and writing:
    This is then substituted into the finite difference equation for x = 0.

  11. Aug 23, 2015 #10
    Thanks a lot. Though it seems more like (T-1,n) = (T1,n) - 2 * delta_x * (q / k), because it seems that (for the implicit, backward scheme) the derivative of (Tn,0) ~ ( T1,n - T-1,n) / (2 * delta_x) = (q / k), please correct me if I'm wrong.
    Last edited: Aug 23, 2015
  12. Aug 23, 2015 #11
    Coding error??? Please write out your finite difference equation for general x, and for x = 0.

    Error in tridiagonal solution?

  13. Aug 23, 2015 #12
    Thanks but I actually deleted that part of my post which was asking for help with a situation whereby changes in delta x gave big changes in Tm,n+1. I must have just played with thickness of the plate, making it narrow, which would mean temperature stabilization would be quicker, because I stuck a Dirichlet boundary on the right hand side. In reality though, it's not Dirichlet, but a changing heat flux, in the form of heat loss to a metallic plate heat sink located on the right side of ptfe plate. I'm currently trying to figure out what the heat flux formula for this boundary is, because an experiment has shown that during friction on the left side of the ptfe plate, ptfe melts a bit (flaking occurs), so temperature at surface must exceed 327 degrees celcius (melting point), but backward euler with mixed boundary conditions showed that temperature at surface is only 86 degrees celcius.
  14. Aug 23, 2015 #13
    A question: how to write matrix b in Ax=b, if the right boundary is at infinity. Would I just write say { (T1,n) - 2*dx*Qin , (T0,n), (T1,n), (T2,n), ... (Tm+1,n) } ? I.e. just add the ghost point at the end, but do not call it Tright, neither dicretize a Neumann at this point. Correct?
  15. Aug 23, 2015 #14
    I don't understand what you are saying here and I don't understand your notation. That's strange, because I've solved a huge number of problems just like this one during my professional career.

  16. Aug 23, 2015 #15
    That's fine. I'll try to explain my notation (had to pick it up from the Internet, from various sources, as I forgot all my school material). A is the tridiagonal matrix, b is matrix of values at time = n, and x is matrix of values at time = n + 1. So the question is, for matrix b with m+2 entires, is b going to look as follows, if the right-hand-side boundary is neither Neumann nor Dirichlet (i.e. at approximated infinity)?

    \left( \begin{array}{ccc}
    T^n_1 - 2 \times dx \times constant \ boundary \ flux \ Qin\\
    T^{n}_{m+1}\end{array} \right)
    Last edited: Aug 23, 2015
  17. Aug 24, 2015 #16
    The right-hand boundary should be approximated at a large finite value of x (or, in your problem, the actual thickness of the slab). The condition should be either zero flux or no-temperature-change. With zero flux, you would do the same trick at the right boundary as you do at the left boundary.

    The ODEs being solved using the method of lines would be:


  18. Aug 24, 2015 #17
    Yes, thanks, Chet. I made my right boundary for ptfe at 1m, when it's only 3mm, and it happens that after one minute of friction (backward euler) the temperature shoots through the roof, being +/- 600 degrees celcius almost half a meter in (which is about twice the ptfe melting temperature.). This coincides with experiment showing that after a good number of open/close cycles (it's a high pressure valve basically, for 40MPa), the metal casing significantly heats up, while the ptfe rod seal looses its use only after a relatively small number of cycles.
  19. Aug 24, 2015 #18
    That said, I think my algorithm is screwed up, need to recheck it, as if I have 30 time steps, dx = 5E-4m, dt = 0.3s and m=2000, I get surface temperature of 2000 after simulation time elapses.
  20. Aug 24, 2015 #19
    I used

    T^{n+1}_0 = T_1^n - 2 \times delta \ x \times Q left + (1-2 \times lambda)\times T_0^n + lambda \times T_1^n, \\ where \ Qleft \ is \ constant \ boundary \ flux \ at \ left \ boundary\\
    T^{n+1}_m = Tright \ (constant \ temperature \ at \ right \ boundary)
  21. Aug 24, 2015 #20
    I don't like your equation for x = 0 for two reasons
    1. It's forward euler
    2. I have no idea how you came up with this equation. The equation is not consistent with the second order accuracy (spatially) equation that I presented in post #16. Please show the derivation.

  22. Aug 24, 2015 #21
    You are right, I substituted the discrete Neumann equation for Tx, i.e. the [itex]T_{-1}^n = T_1^n - 2 \ delta \ x \times \frac{Qleft}{k}[/itex], into a wrong (explicit) scheme, instead of backward Euler (nevermind the missing k, as most formulas are given for U, heat, not temperature, so I got confused). That's why also with some initial conditions I had oscillations.
    Last edited: Aug 24, 2015
  23. Aug 24, 2015 #22
    Plus in relation to discretizing the boundary condition itself, I was using derivative based on central difference, not backward Euler. Not sure how big of an effect it could have had or whether it matters. Used [itex]\frac{y_{j+1} - y_{j-1}}{ 2 \ delta \ x}[/itex]. Should I have used [itex]\frac{y_{j} - y_{j-1}}{ delta \ x}[/itex] too?
  24. Aug 25, 2015 #23
    This difference equation is incorrect, not only because you are using forward euler in time. Please show us your derivation.

  25. Aug 25, 2015 #24
    Using central differences on the spatial derivatives is fine.

  26. Aug 25, 2015 #25
    [itex]U_x(t_n,x_0) = \frac{U_1^n-U_{-1}^n}{2 \times delta \ x} = Qleft, \ therefore\\
    U_{-1}^n=U_1^n - 2 \times delta \ x \times Qleft \ (1)\\
    Substitute \ (1) \ into \ Forward \ Euler \ for \ j = 0, i.e. \ into \ U_0^{n+1} = U_{-1} + (1-2 \times lambda) \ U_0^n + lambda \times U_1^n, \ and \ voila.[/itex]
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook