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

Inner holes

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.

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.

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