Heat equation given constant surface heat flux

[user4417]

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).

 

[Hesch]

How would I go about finding temperature distribution in a thin square plate

Suggestion: Divide the plate into a huge number of cubes, all with a volume = (dx)³.
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 ).

 

[user4417]

Suggestion

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.

 

[user4417]

I finally found an explanation that made sense. Thanks for your help.

 

[Chestermiller]

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.

Chet

 

[user4417]

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.

Chet

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).

 

[Chestermiller]

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).

I think you are thinking of forward Euler. Backward Euler will always converge.

Chet

 

[user4417]

I think you are thinking of forward Euler. Backward Euler will always converge.

Chet

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?

 

[Chestermiller]

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?

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:
$$T(-Δx)=T(+Δx)+Q/k$$
This is then substituted into the finite difference equation for x = 0.

Chet

 

[user4417]

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:
$$T(-Δx)=T(+Δx)+Q/k$$
This is then substituted into the finite difference equation for x = 0.

Chet

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.

 

[Chestermiller]

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.

Coding error? Please write out your finite difference equation for general x, and for x = 0.

Error in tridiagonal solution?

Chet

 

[user4417]

Coding error? Please write out your finite difference equation for general x, and for x = 0.

Error in tridiagonal solution?

Chet

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.

 

[user4417]

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?

 

[Chestermiller]

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?

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.

chet

 

[user4417]

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.

chet

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_0\\ T^n_1\\ ...\\ T^{n}_{m+1}\end{array} \right)$​

 

[Chestermiller]

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:

$$\frac{dT_m}{dt}=a\left[\frac{T_{m-1}-2T_{m}+T_{m+1}}{(Δx)^2}\right]\tag{m≠0}$$
$$\frac{dT_0}{dt}=2a\left[\frac{T_{1}-T_{0}}{(Δx)^2}\right]+2a\frac{q}{kΔx}\tag{m=0}$$

Chet

 

[user4417]

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.

 

[user4417]

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.

 

[user4417]

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)$​

 

[Chestermiller]

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)$​

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.

Chet

 

[user4417]

You are right, I substituted the discrete Neumann equation for Tx, i.e. the $T_{-1}^n = T_1^n - 2 \ delta \ x \times \frac{Qleft}{k}$, 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.

 

[user4417]

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 $\frac{y_{j+1} - y_{j-1}}{ 2 \ delta \ x}$. Should I have used $\frac{y_{j} - y_{j-1}}{ delta \ x}$ too?

 

[Chestermiller]

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$​

This difference equation is incorrect, not only because you are using forward euler in time. Please show us your derivation.

Chet

 

[Chestermiller]

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 $\frac{y_{j+1} - y_{j-1}}{ 2 \ delta \ x}$. Should I have used $\frac{y_{j} - y_{j-1}}{ delta \ x}$ too?

Using central differences on the spatial derivatives is fine.

Chet

 

[user4417]

This difference equation is incorrect, not only because you are using forward euler in time. Please show us your derivation.

Chet

$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.$

 

[user4417]

Using central differences on the spatial derivatives is fine.

Chet

Oh, ok, thanks a lot. I'll have to redo the whole thing. Will ask if have further problems.

 

[Chestermiller]

$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$

This is incorrect. Where is the thermal conductivity k in these equations? The sign in front of the 2 is wrong also.

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.$

This is also incorrect. The finite difference equation for forward euler should be:

$$\frac{T_0^{n+1}-T_0^n}{Δt}=a\frac{(T_{-1}^n-2T_0^n+T_1^n)}{(Δx)^2}=a\frac{(-2T_0^n+2T_1^n+2ΔxQ/k)}{(Δx)^2}$$

Voila.

I leave it up to you to get the algebra right in terms of λ.

By the way, I still recommend using backward euler in time for a problem like this.

Chet

 

[user4417]

Lol. Thanks for corrections. My brain's a bit rusty after so many years after school. Sorry!

 

[user4417]

A question. How will my top row of the matrix to be inversed will look like, (1+λ, -λ, 0, 0, ...), or the usual (1+2*λ, -λ, 0, 0, ...)?
When I substitute $U_{j-1}=U_j + Δ\ x \ \frac{Qleft}{k}$ into $-λ \ U_1^n + (1 + 2 \ λ) U_0^n - λ \ U_{-1}^n = U_0^{n-1}$, I get $-λ \ U_1^n + (1 + λ) U_0^n + λ \times Δ\ x \times \frac{Qleft}{k} = U_0^{n-1}$, from which it follows that the top row of matrix b (Ax=b) will be equal to $U_0^{n-1}+λ \times Δ\ x \times \frac{Qleft}{k}$, while the top row of matrix A (Ax=b) will be (1 + λ, -λ, 0, 0, 0, ...), instead of the usual (1 + 2*λ, -λ, 0, 0, 0, ...). But running this simulation, my room temperature in Kelvin of the piece goes to +5 in 10 seconds from 297, and I don't know how to fix it. Could you please check my code too?

        int t_steps=100; // number of timesteps
        double dt = 0.1; // time step dt = 0.2s
        double dx = 5E-4; // space step dx, 0.1mm
        int m=20; // number of nodes on the x
        double L=dx*(m-1); // width of ptfe, 1cm

        int q1=16500; // q1 constant heat flux value at left boundary
        int q2=0;
        double k_ptfe=0.25; // conductivity of ptfe plate
        int alpha_ptfe=8800000; // alpha value for ptfe plate
        double lambda=alpha_ptfe*dt/Math.pow(dx,2); // lambda value for matrix A in x = b * inverse(A)

        double [][] values_A = new double[m+2][m+2]; // matrix A with ghost points at both ends
        double [][] values_b = new double[m+2][1]; // matrix b with temperatures at time = n
        double [][] values_x = new double[m+2][1]; // matrix x with temperatures at time = n + 1

        // preparing matrix A as follows

        for (int m_index_col=0; m_index_col<=m+1; m_index_col++)
            for (int m_index_row=0; m_index_row<=m+1; m_index_row++)
            {
                values_A[m_index_row][m_index_col]=0;
            }

        values_A[0][0]=1+lambda;
        values_A[0][1]=-lambda;
        values_A[m+1][m+1]=1+lambda;
        values_A[m+1][m-1]=-lambda;
       
        for (int m_diag_index=1; m_diag_index<=m; m_diag_index++)
            {
                values_A[m_diag_index][m_diag_index]=1+2*lambda;
                values_A[m_diag_index][m_diag_index-1]=-lambda;        
                values_A[m_diag_index][m_diag_index+1]=-lambda;
            }

        // preparing matrixes b and x as follows

        for (int m_index_row=0; m_index_row<=m+1; m_index_row++)
        {
            values_b[m_index_row][0]=297;
            values_x[m_index_row][0]=0;
        }

        Matrix A = new Matrix(values_A);
        Matrix b = new Matrix(values_b);
        Matrix x = new Matrix(values_x);
 
        // calculation loop

        for (double t=0; t<t_steps*dt; t=t+dt)
        {
            double temp1 = values_b[0][0];
            double temp2 = values_b[m+1][0];

            x = A.solve(b); // solve for matrix x
            values_x = x.getArray();

            values_b = values_x; // increment time

            values_b[0][0]=temp1-lambda*dx*q1/k_ptfe; // set Neumann at left boundary using central difference
            values_b[m+1][0]=temp2-lambda*dx*q2/k_ptfe; // set Neumann at right boundary
   
            Matrix transitory = new Matrix(values_b);
            b=transitory;
            transitory=null;
            System.gc();    
        }

 

[Chestermiller]

It should be (1+2λ, -2λ, 0, ...,0)

Chet

 

[user4417]

It should be (1+2λ, -2λ, 0, ...,0)

Chet

Is it not $\begin{pmatrix} 1+2λ & -λ & 0 & 0 & ... & 0\\ \end{pmatrix}$? Also, could you please check this part of my code, please?

        int t_steps=600;
        double dt = 0.1;
        double dx = 0.01;
        int m=10;
        double L=dx*m;
      
        int q1=16500;
        int q2=0;
        double k=0.25;
        int alpha=8800000;
        double lambda=alpha*dt/Math.pow(dx,2);

        values_A[0][0]=1+2*lambda;
        values_A[0][1]=-lambda;       
        values_A[m+1][m+1]=1+2*lambda;
        values_A[m+1][m]=-lambda;
              
        for (int m_diag_index=1; m_diag_index<=m; m_diag_index++)
            {
                values_A[m_diag_index][m_diag_index]=1+2*lambda;
                values_A[m_diag_index][m_diag_index-1]=-lambda;               
                values_A[m_diag_index][m_diag_index+1]=-lambda;
            }

        for (double t=0; t<t_steps*dt; t=t+dt)
        {     
            x = A.solve(b);
            values_x = x.getArray();
          
            values_b = values_x;

            values_b[0][0]=values_b[1][0]+dx*q1/k;
            values_b[m+1][0]=values_b[m][0]+dx*q2/k

I'm asking because I get

x
T

0.0
660.0000000634616

0.01
6.346153846108459E-8

0.02
5.769230769195468E-8

0.03
5.1923076922824766E-8

0.04
4.615384615369485E-8

0.05
4.038461538456494E-8

0.06
3.461538461543504E-8

0.07
2.8846153846305132E-8

0.08
2.3076923077175224E-8

0.09
1.7307692308045316E-8

0.1
1.1538461538915409E-8

 

[Chestermiller]

Is it not $\begin{pmatrix} 1+2λ & -λ & 0 & 0 & ... & 0\\ \end{pmatrix}$?

If you can't get the algebra right, I can't help you.

Also, could you please check this part of my code, please?

I don't check code.

Chet

 

[Chestermiller]

The elements in each row of the matrix must add up to 1. Try running the code with Q = 0 and see what you get if you use -λ instead of -2λ for the second element in the first row. With Q = 0, your calculation should yield no change in the temperature distribution if the initial temperature distribution is uniform.

Chet

 

[user4417]

The elements in each row of the matrix must add up to 1. Try running the code with Q = 0 and see what you get if you use -λ instead of -2λ for the second element in the first row. With Q = 0, your calculation should yield no change in the temperature distribution if the initial temperature distribution is uniform.

Chet

I don't understand why 2 lambdas needed instead of one. I mean it's all based on the same formulae as other rows, which use single lambda for n and n+1 (after n-1). In any case, results finally stabilized with this solution I don't understand. It seems incredible though, e.g. this is what I get for:

        int t_steps=2*6*10; // number of time steps
        double dt = 0.1; seconds
        double dx = 0.0001; // meters
        int m=30; // number of dx steps

(the thing seems to heat up much more than that, enough for melt-caused flaking of ptfe to occur)

x [m]
T [K]

0.0
316.3741934284825

1.0E-4
309.7741934284825

2.0E-4
309.7741934284825

3.0000000000000003E-4
309.7741934284825

4.0E-4
309.7741934284825

5.0E-4
309.7741934284825

6.000000000000001E-4
309.7741934284825

7.0E-4
309.7741934284825

8.0E-4
309.7741934284825

9.000000000000001E-4
309.7741934284825

0.001
309.7741934284825

0.0011
309.7741934284825

0.0012000000000000001
309.7741934284825

0.0013000000000000002
309.7741934284824

0.0014
309.7741934284824

0.0015
309.7741934284824

0.0016
309.7741934284824

0.0017000000000000001
309.7741934284824

0.0018000000000000002
309.7741934284824

0.0019
309.7741934284824

0.002
309.7741934284824

0.0021000000000000003
309.7741934284824

0.0022
309.7741934284824

0.0023
309.7741934284824

0.0024000000000000002
309.7741934284824

0.0025
309.7741934284824

0.0026000000000000003
309.7741934284824

0.0027
309.7741934284824

0.0028
309.7741934284824

0.0029000000000000002
309.7741934284824

0.003
309.7741934284824

 

[Chestermiller]

I don't understand why 2 lambdas needed instead of one. I mean it's all based on the same formulae as other rows, which use single lambda for n and n+1 (after n-1). In any case, results finally stabilized with this solution I don't understand. It seems incredible though, e.g. this is what I get for:

        int t_steps=2*6*10; // number of time steps
        double dt = 0.1; seconds
        double dx = 0.0001; // meters
        int m=30; // number of dx steps

(the thing seems to heat up much more than that, enough for melt-caused flaking of ptfe to occur)

The finite difference equation for the flux that you wrote in post #29 is a one-sided difference approximation, and is only first order accurate in Δx. The central difference equation that I had written in an earlier post is a central difference approximation, and is second order accurate in Δx (i.e., much more accurate). For the central difference approximation, the elements are (1+2λ) and -2λ.

For the record, the correct difference equation at x = 0 for the central difference approximation to the heat flux should be:

$$(1+2λ)T_0^n-2λT_1^n=T_0^{n-1}+2λΔx\frac{Q}{k}$$

The numerical solution you presented in your previous post at the end of the first time step looks to me like the results of applying forward euler in time, not backward euler. I know this because all the temperatures should have changed significantly with backward euler, not just the temperature at the boundary.

Chet