Finite volume method to determine steady state temp in cylinder

Question

I am attempting to model the steady state behavior of a cylinder using the finite volume method (FVM) subjected to a variety of boundary conditions in Matlab. First off, I am treating the cylinder as being axisymmetric so I am only determine the temperature profile in the r-Z plane. I set up a 2D grid so that the r coordinate is on the 'x' axis and the Z coordinate is the 'y' axis. Since the cylinder is axisymmetric then the heat flux at r = 0 is 0 (insulated). I am also applying an insulated boundary condition to the bottom surface at Z=0. The side surface has convection with a constant convection coefficient and air temperature of 310 K. The top of the cylinder, z=Z, there is a heat flux of 1000 W/m^2, convection from the air at 310K, and also surface to ambient radiation with the surroundings being 298 K. The initial temperature of cylinder is 298 K (initial guess value). The problem that I am having is that the temperature distribution matches Comsol's results with a grid of 40x40 and a thermal conductivity of 1. When I increases the number of cells to 50x50, 60x60, etc, or if i change the thermal conductivity to something other than 1 then the temperature profile is still what it should be but all the temperature are significantly lower than expected. My control loop in Matlab is as follows:

1. Use initial values and problem constants to solve for initial surface temperature of top and side.
2. Use surfacte temperatures and problem geomety to determine aE, aN, aS, aW, aP, and b for every cell in the grid.
3. Use line by line TDMA in both directions to determine the temperature profile.
4. Calculate the error between the new temperature profile and the previous one.
5. Update the surface temperatures and aP and b because of the nonlinear radiation term that uses previous iterate values.
6. Perform line by line TDMA again and check for convergence. If the error is not less than 0.001 then repeat the procedure.

For the boundary conditions, I used an energy balance for the top surface and the side surface. The top surface includes surface to ambient radiation and I calculated the surface temperature to be $$T_{b}=\frac{\alpha G+hT_{f}+\frac{kT_{p}}{\delta z}+\epsilon \sigma(T_{sur}^{4}+3T_{b}^{*4})}{h+4 \epsilon \sigma T_{b}^{*3}+\frac{k}{\delta z}}$$

I also do the same thing with the surface temperature on the side the is experiencing heating through convection. Since the top surface temperature depends on its previous iterate value, I use an initial guess value of 298 for the top surface. I then calculate the initial surface temperature based on that and the equation i provided for the first iteration,

The problem lies in the fact that the surface tempreature depends on both the thermal conductivity and the grid spacing. When I use many cells like 80x80, the top surface temperature for the first TDMA iteration is almost 298 which is hard to believe because its being heated by a 1000 W/m^2 heat flux and its also being heated by the 310 K air temperature. This causes the first iteration to have very little error so convergence only produces a max temperature of like 325 when it should be 450. But if I use a grid of like 40x40 then the surface temperature for the first iteration is much much higher, and convergence causes the final temperature profile to look almost exactly like the expected.

The same exact thing happens with the thermal conductivity, if i change the value then the initial surface temperatures change dramatically which causes the convergence values to be much lower than the expected.

It seems like it all has to do with that initial surface temperature but I am not sure how else to treat it. Any help would be amazing as I have spent many nights til 5am trying to figure out how I could fix this. My contact email is in my profile if you would like to see my matlab code of if you would like to talk to me more about my problem.

Answer 1

Verification is the first step with any new code. I like using the method of manufactured solutions (MMS) because it permits the use of rich solutions, but there are also closed-form analytic solutions for similar problems. I would verify the diffusion solver first, before adding the radiation term.

The idea is to pick any smooth function $u(r,z)$ (does not need to be "physical") and feed it to your differential operator to "manufacture" a forcing term for which your chosen $u(r,z)$ is a solution. Then run your numerical code and perform a grid convergence study, evaluating $\lVert u^h - u \rVert$ on each grid $h$. Confirm that the design order of the method (e.g., quadratic) has been reached. If not, simplify further until you see the design order.

Your symptoms suggest that there may be a scaling bug, but whatever it is, MMS should diagnose the problem.