Heat equation given constant surface heat flux

[user4417]

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

Chet

Is it not \begin{pmatrix}<br /> 1+2λ &amp; -λ &amp; 0 &amp; 0 &amp; ... &amp; 0\\<br /> \end{pmatrix}<br />? 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}<br /> 1+2λ &amp; -λ &amp; 0 &amp; 0 &amp; ... &amp; 0\\<br /> \end{pmatrix}<br />?

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

 

[user4417]

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

Thanks a lot for walking through this whole thing with me, I applaud your patience and keenness in helping people out with physics problems. I figured out what my problem was. I was using alpha as Cp * ro / k, when it is k / Cp / ro. I used backward Euler for both Neuman boundaries (and backward Euler for the equation itself) and results were very close, though not in line with experiemntal data, namely temp only rises to 175C degrees at the surface after 10-20s of constant rubbing, while it should be over 326.8C. The experiement conditions were weird though. Second boundary actually had a heat sink (metal case of valve), and the rubbing was periodic and not constant (by-hand closing).