MECOM 2021
Solid heat simulations are important in industry due to the large number of problems associated with heat buildup in devices. Finding a way to focus heat flows or to isolate certain geometries becomes relevant to solve these problems.
The finite difference method will be used to solve the heat equation that models the transient and stationary regimes when the system is subjected to a temperature gradient. Geometries with different hole shapes are considered.
We explore the effects of different pattern configurations and hole shapes, to improve the controlled direction of heat flow in a continous solid.
The simulations are done with an own code written in C++ and the data visualization was done through Octave.
Crescent hole of cross-sectional area equal to 0.125 of total surface area with 400°C in the blue side and 200°C in the red side.
Crescent hole of cross-sectional area equal to 0.5 of total surface area with 400°C in the blue side and 200°C in the red side.
Crescent hole of cross-sectional area equal to 0.875 of total surface area with 400°C in the blue side and 200°C in the red side.
Triangular hole of cross-sectional area equal to 0.125 of total surface area with 400°C in the blue side and 200°C in the red side.
Triangular hole of cross-sectional area equal to 0.5 of total surface area with 400°C in the blue side and 200°C in the red side.
Triangular hole of cross-sectional area equal to 0.875 of total surface area with 400°C in the blue side and 200°C in the red side.
In principle, the problem starts from a square domain in which different equations were applied to the different sections of the domain. For the center, the equation of heat conduction in solids was applied. For the edges, both zero flow conditions were applied for the y-edges and fixed temperature conditions for the x-edges.
On the other hand, the initial conditions were for all the simulations of a temperature of 200K while for the edges at "x" they varied between 200 and 400 degrees.
The first equation to be considered for this problem was the equation of heat conduction in two-dimensional solids and in time, taking as spatial dimensions x and y and as time dimension t. Where alpha will be the diffusion constant of the solid
To discretize this equation we used the central finite difference formula for the spatial dimensions and the forward formula for the temporal dimension, so that our equation was left with the form
Where the "i", "j" and "k" are finite representations of "x" , "y" and "t" respectively with a step distance of ∆x, ∆y and ∆t. They are also related to the direction of the node in which we are located.
Thus, using the previous expression and applying a little algebra, we arrive at the expression for the value of the node at the next time step, i.e. to predict the value at the node at the next instant.
Then, using fictitious nodes and a bit of algebra we arrive at the expression of the limits of the domain (not including the corners). This is achieved by equating the boundary conditions with the last nodes within the domain, which will present fictitious nodes that will be unknowns to look for in solving the equation.
For the left edge of our domain seen from a front view, i.e. for the edge when x = 0, it will be:
Where "a0" will be 0 if the condition is a heat flow and on the contrary if it is convection, "a0" will be equal to the convection coefficient h of the fluid.
On the other hand "f0" will be equal to the heat flux if the condition is heat flow but if the condition is convective, then "f0" will be equal to the product of the fluid temperature by the coefficient h of the fluid. Finally the coefficient k is the thermal conduction of the material.
For the case presented in which the left edge (x=0) has a fixed temperature, the equation will be much simpler, being simply: Tk+1 = Tk
It is important to note that for the other edges this formula will be identical, changing only the sub indexes depending on which edge we are on, being adapted to each case. For example for the edges at "y" where the condition we imposed is a flux equal to zero, we will have to replace the values of "a0" and "f0" by the appropriate ones, which will be 0 and 0 respectively.
For example for the edge when y = L will be:
For the corners of the domain, the same clearance is performed, except that instead of one fictitious node, there will be two fictitious nodes. This development is not placed here due to the extension of the final formula, however, it should be clarified that for the case considered in this project it was not necessary to use this formula since in the corners a fixed temperature was set because they belong to the edges of x = 0 and x = L.
For more information on the development of the discretized equations at the corners, see [1].
For the discretization of the holes, extreme care must be taken when simulating the conditions for the nodes bordering the hole, especially for nodes with all their neighboring nodes present.
The first and simplest is the case of the boundary with a horizontal or vertical edge of the hole.
For our particular case the hole will have the zero flow condition so the equation for node type 1 will be:
For case 2 in our example with zero flow, the equation will be:
Finally, for the case of node 3 and for flow equal to zero we will have:
We must pay attention that the nodes that are "in contact" with the edge, will have a 2 multiplying, while the nodes with all their neighboring nodes will only have a 1.
In the creation of the program, this point was key to get consistent results and finding this formula was especially difficult because being a node that has the nodes i+1, i-1, j+1 and j-1, we could assume that it behaves like a normal interior node, however by not having a node of the corners appears this formula where it is weighted with a 2 multiplying the node that is near the edge, which causes a mini compensation of errors that in sum will be key to have consistent results.
To validate the program, tests were performed with systems with analytical solutions and with physically known results even though the analytical solution was not known.
At first the program was tested with a domain without holes where we set the end x=0 to oscillate between 280 k and 320 k while for the edges x=L, y=0 and y=L we imposed a zero flux condition and an initial temperature of 300K.
The analytical solution to this problem is:
We found very similar behavior between the finite difference simulation and the analytical solution.
The next validation performed was for the problem of a domain with a hole in its interior. For this problem what was done was to perform simulations with arbitrary shaped holes in the interior of the domain and simulate until the steady state was reached. At this point, given the boundary conditions (fixed temperature at x=0 at 400 k, x=L at 200 k, flux equal to zero at y=0 and y=L, the whole domain at a temperature of 300 k), the flux entering the domain at the edge x=0, would have to be theoretically equal to the outgoing flux x=L and would also have to be equal to the flux at the middle of the domain x=L/2. This could be achieved with an error of 1.4%.
We can observe how for the case of a triangle with an occupied cross section equal to 1/2, the fluxes at x=0, x=L/2 and x=L converge to the same value.
In principle we have to define the rectification calculation used. To do this we define as hot side the edge of the domain subjected to 400k. In this way we have case A and case B, in which we will have at x=0 a temperature of 200k and at x=L we will have 400k, while for case B we will have exactly the opposite (x=0 at a temperature of 400k and x=L at 200k). In this way we will define the rectification of the system as:
First, we present the steady-state flows for the triangle and crescent shapes as a function of the cross-space they occupy in the domain.
We can observe that naturally, as expected, the flow decreases as the cross-sectional area occupied by the hole increases. We can also observe that the crescent shaped hole will hinder the flow more than the triangle shape for the same size, however this effect starts to be noticed from the 3/4 size of the cross-sectional area hindered.
In this image we can observe the rectification as a function of the occupied section in the domain for the triangular shaped hole, thus we can see how the rectification effect increases considerably as a function of the occupied section. In addition to the latter we can also see that the rectification peak takes longer to appear as we have a smaller hole size, this is because the hole being smaller ends up being further away from the hot edge, which causes the heat passage takes longer to reach the asymmetrical section that causes the hole.
According to the above, we plot the peak values of the rectifications as a function of the cross-sectional area occupied for both figures (Triangle and Crescent).
In this way we can observe that the triangle-shaped hole is a better rectifier because its peak is much larger than the crescent-shaped hole, almost doubling the rectification of the latter.
Finally, we plot the time in which the maximum rectification peaks are reached as a function of the cross section occupied by the hole.
We can observe the curious effect that while for the triangle, as we increase the size of the hole, the time in which peak rectification is reached is lower, for the crescent the opposite happens, the time increases as a function of the size occupied by the hole.
[1] M. Necati Özişik , Helcio R. B. Orlande, Marcelo J. Colaço, Renato M. Cotta, "Finite Difference Methods in Heat Transfer, second edition" .NewYork, USA : CRC Press, 2017.
[2] Theodore L. Bergman, Adrienne S. Lavine, Frank P. Incropera and David P. Dewitt, "Heat and Mass Transfer, seventh edition". United States of America: John Wiley & Sons, 2011.
[3] Michel T. Barakoetal, Nanotechnology, volume29, number 15 (2018).
Send an email to neptp.ungs@gmail.com for more information about the project.