A Finite Difference Algorithm Applied to the Averaged Equations of the Heat Conduction Issue in Biperiodic Composites—Robin Boundary Conditions

This note deals with the heat conduction issue in biperiodic composites made of two different materials. To consider such a nonuniform structure, the equations describing the behavior of the composite under thermal (Robin) boundary conditions were averaged by using tolerance modelling. In this note, the process of creating an algorithm that uses the finite difference method to deal with averaged model equations is shown. This algorithm can be used to solve these equations and find out the temperature field distribution of a biperiodic composite.


Introduction
The primary purpose of this research was to create a computational algorithm that can be used in the analysis of the heat conduction issue with regard to biperiodic composites. Periodic or functionally graded structures made of two or more materials can exhibit advantageous mechanical or strength properties (in view of the adopted criterion). For this reason, it is important to be able to consider and model these structures. There are many methods which can be used in the analysis of heterogeneous structures such as composites. Among them are asymptotic homogenization [1,2], homogenization with microlocal parameters [3,4], the finite element method [5,6], the boundary element method [7], the meshless methods [8], the higher order theory [9,10], and the methods using the Green-Lindsay model [11]. In the analyzed biperiodic structures, the so-called periodicity cells can be mentally separated. In turn, these cells are characterized by a certain dimension, called the microstructure size or microstructure parameter. Unfortunately, most known methods do not allow including this parameter in the considerations. Therefore, to obtain the equations describing the considered issue, tolerance modelling [12][13][14] was used. Then, a computational algorithm was created to analyze the heat conduction phenomenon in biperiodic composites, using the finite difference method on the derived equations. By using this algorithm, it is possible to obtain the temperature field distribution in the analyzed composite. On selected edges of this structure, the third kind boundary conditions were imposed (the Robin boundary conditions) [15].
The theoretical structure we analyze herein is a composite made of two materials with different properties; confer Figure 1. These material properties are changing in a periodic way along both perpendicular directions (x 1 and x 2 ). The dimension along x 1 is denoted by L 1 , and that in direction x 2 -analogously-by L 2 . The mentioned microstructure parameters along directions x 1 and x 2 , are denoted by l 1 and l 2 , respectively. The microstructure size is directly related to the number of composite layers (cells). The volume share of each material in every cell is determined and constant. The product of the share of the first material in a cell along x1 and along x2 is a volume share of the first material in a cell and is denoted by ν1 = ν1(x1, x2). The volume share of the second material in a cell ν2 = ν2(x1, x2) can be calculated in an analogous way.

Averaged Equations
Equation (1) describing the heat conduction phenomenon with reference to biperiodic composite is an equation with noncontinuous coefficients: where θ is an unknown temperature field; c and ρ specify the material properties such as a specific heat and a mass density, respectively; and kij defines the components of the conductivity tensor. To average this equation, tolerance modelling was used [16][17][18][19][20][21][22]. Tolerance modelling, also called the tolerance averaging technique, introduces in a process of modelling a new concepts, definitions, and assumptions. The new concepts primarily relate to functions that comply with certain conditions. Among them we can distinguish the tolerance-periodic function, the highly-oscillating function, and the slowlyvarying function. Among the assumptions of this method, the most important, from the point of view of the present considerations, is the micro-macro decomposition assumption. In conformity with this assumption, the unknown temperature field can be taken according to the following formula: where θ is an averaged part of the temperature field Θ, called the macrotemperature. The remainder of the equation is the sum of the products of the fluctuation shape functions g1 and g2; and fluctuation amplitudes ψ1 and ψ2. The mentioned sum is called the fluctuating part of the temperature field Θ. The fluctuation amplitudes are the new unknowns, and the fluctuation shape functions are assumed a priori. These functions should be assumed individually, taking into account the specifics of the issues under consideration. In this The volume share of each material in every cell is determined and constant. The product of the share of the first material in a cell along x 1 and along x 2 is a volume share of the first material in a cell and is denoted by ν 1 = ν 1 (x 1 , x 2 ). The volume share of the second material in a cell ν 2 = ν 2 (x 1 , x 2 ) can be calculated in an analogous way.

Averaged Equations
Equation (1) describing the heat conduction phenomenon with reference to biperiodic composite is an equation with noncontinuous coefficients: where θ is an unknown temperature field; c and ρ specify the material properties such as a specific heat and a mass density, respectively; and k ij defines the components of the conductivity tensor. To average this equation, tolerance modelling was used [16][17][18][19][20][21][22]. Tolerance modelling, also called the tolerance averaging technique, introduces in a process of modelling a new concepts, definitions, and assumptions. The new concepts primarily relate to functions that comply with certain conditions. Among them we can distinguish the tolerance-periodic function, the highly-oscillating function, and the slowlyvarying function. Among the assumptions of this method, the most important, from the point of view of the present considerations, is the micro-macro decomposition assumption. In conformity with this assumption, the unknown temperature field can be taken according to the following formula: where θ is an averaged part of the temperature field Θ, called the macrotemperature. The remainder of the equation is the sum of the products of the fluctuation shape functions g 1 and g 2; and fluctuation amplitudes ψ 1 and ψ 2 . The mentioned sum is called the fluctuating part of the temperature field Θ. The fluctuation amplitudes are the new unknowns, and the fluctuation shape functions are assumed a priori. These functions should be assumed individually, taking into account the specifics of the issues under consideration. In this note, it is important to adopt functions that guarantee the continuity of the total temperature, both between the layers (cells), and inside the cells on interfaces between materials. Therefore, some combination of the basic functions (saw-type and piecewise parabolic functions) shown in Figure 2 are assumed.
Materials 2021, 14, x FOR PEER REVIEW  3 of 19 note, it is important to adopt functions that guarantee the continuity of the total temperature, both between the layers (cells), and inside the cells on interfaces between materials. Therefore, some combination of the basic functions (saw-type and piecewise parabolic functions) shown in Figure 2 are assumed. By γ1 and γ2, the shares of the first material in a cell in direction x1 and in direction x2, respectively, are denoted. The fluctuation shape functions g1(x1, x2) and g2(x1, x2) are assumed according to the following formulas: f1, f2, h1, and h2 are depicted in Figure 2. By considering the micro-macro decomposition assumption, and the other assumptions and definitions of the tolerance modelling discussed in [23][24][25], Equation (1) describing the heat conduction issue was averaged, leading to the tolerance model equations: where K stands for the thermal conductivity tensor whose components are kij, ∇ is a gradient operator defined as (∂1, ∂2, ∂3), overlined ∇ is a gradient in x3 direction (0, 0, ∂3), and ∂ is a gradient operator defined as (∂1, ∂2, 0). Equations (5)-(7) can be obtained by using the orthogonalization method [26] or by the extended stationary action principle [27,28].
The averaging operator, denoted in Equations (5)-(7) by triangular brackets, is defined according to the following equation: By γ 1 and γ 2 , the shares of the first material in a cell in direction x 1 and in direction x 2 , respectively, are denoted. The fluctuation shape functions g 1 (x 1 , x 2 ) and g 2 (x 1 , x 2 ) are assumed according to the following formulas: f 1 , f 2 , h 1 , and h 2 are depicted in Figure 2. By considering the micro-macro decomposition assumption, and the other assumptions and definitions of the tolerance modelling discussed in [23][24][25], Equation (1) describing the heat conduction issue was averaged, leading to the tolerance model equations: where K stands for the thermal conductivity tensor whose components are k ij , ∇ is a gradient operator defined as (∂ 1 , ∂ 2 , ∂ 3 ), overlined ∇ is a gradient in x 3 direction (0, 0, ∂ 3 ), and ∂ is a gradient operator defined as (∂ 1 , ∂ 2 , 0). Equations (5)-(7) can be obtained by using the orthogonalization method [26] or by the extended stationary action principle [27,28].

Algorithm of Finite Difference Method
First, the equations of the tolerance model were reformulated and written in index notation. The indices take the following values: the superscripts a, b = 1, 2, the subscripts a, b = 1, 2, 3; and the subscripts α, β = 3: To prepare an algorithm of the finite difference method solving equations of the tolerance model, some assumptions had to be made at this stage. It was assumed that a non-stationary, two-dimensional problem would be considered. The character of the fluctuation shape functions, adopted and shown in the Figure 2, was also considered, which led to modified equations: In the next step, Equations (12)- (14) was reformulated using the differential quotient formulas and written for the node (i,j): where ∆x 1 and ∆x 2 define the distances between nodes in directions x 1 (i) and x 2 (j); confer Figure 3. By t the time coordinate is denoted, and by m, the time steps. The nodes were introduced to the center of each subcell and the interfaces between the cells and between the subcells, so the number of the nodes m1 in direction x1 is equal to 4N1 + 1, where N1 is a number of the cells in direction x1, and analogously, the number of the nodes m2 in direction x2 is equal to 4N2 + 1, where N2 is a number of the cells in direction x2.
Next, the boundary conditions have to be specified, as they affect the process of the algorithm (if we know the values of the desired unknown at the selected nodes, we do not write the equation associated with that unknown at those nodes).

Matrix of Coefficients
The creation of the algorithm was started by building a matrix of coefficients found in the equations.   The nodes were introduced to the center of each subcell and the interfaces between the cells and between the subcells, so the number of the nodes m 1 in direction x 1 is equal to 4N 1 + 1, where N 1 is a number of the cells in direction x 1 , and analogously, the number of the nodes m 2 in direction x 2 is equal to 4N 2 + 1, where N 2 is a number of the cells in direction x 2 .
Next, the boundary conditions have to be specified, as they affect the process of the algorithm (if we know the values of the desired unknown at the selected nodes, we do not write the equation associated with that unknown at those nodes).

Matrix of Coefficients
The creation of the algorithm was started by building a matrix of coefficients found in the equations. The nodes were introduced to the center of each subcell and the interfaces between the cells and between the subcells, so the number of the nodes m1 in direction x1 is equal to 4N1 + 1, where N1 is a number of the cells in direction x1, and analogously, the number of the nodes m2 in direction x2 is equal to 4N2 + 1, where N2 is a number of the cells in direction x2.
Next, the boundary conditions have to be specified, as they affect the process of the algorithm (if we know the values of the desired unknown at the selected nodes, we do not write the equation associated with that unknown at those nodes).

Matrix of Coefficients
The creation of the algorithm was started by building a matrix of coefficients found in the equations.   It was assumed that the total temperature will be known at the nodes of area 0 (for i equals 1 and j from 1 to m 1 ), which translates to the fact that the macrotemperature will also be known, and we do not write Equation (15) for these nodes.
The nodes of area 1 (for i from 2 to m 1 − 1 and for j from 2 to m 2 − 1) are internal nodes and are not affected by the boundary conditions, and it is necessary to write Equation (15) for these nodes. By writing an equation for a given node, the coefficients for the unknown at that node (i, j), but also at the node above it (i − 1, j), below it (i + 1, j), to the right (i, j + 1), and to the left (i, j − 1), were grouped. For the purposes of this work, simplified nomenclature was adopted for the coefficient groups. The coefficients for the unknown (macrotemperature θ) in the node (i, j) in this area are written below: the coefficients in the node (i − 1, j): the coefficients in the node (i + 1, j): the coefficients in the node (i, j − 1): the coefficients in the node (i, j + 1): The nodes of area 2 (for i equal to m 1 and for j from 2 to m 2 − 1) are the nodes on the bottom surface of the composite. On this surface, the Robin boundary conditions are assumed. These conditions are related to a heat exchange and are expressed according to the formula: where U 1 is a heat transfer coefficient and Θ e is an external temperature. This condition was written by using the micro-macro decomposition assumption and the differential quotient formulas: where the averaging operation is valid only in direction x 2 . From the above equation, the component θ i+1,j was determined, which allowed the elimination of the unknown in the virtual node (i + 1, j) from Equation (15). The formulas describing the individual coefficients are shown below: The nodes of area 3 (for i from 2 to m 1 − 1 and for j equals m 2 ) represent the nodes on the right surface of the composite. On this surface, the Robin boundary conditions are also assumed. These conditions are expressed according to the formula: where U 2 is a heat transfer coefficient. This condition was written by using the micro-macro decomposition assumption and the differential quotient formulas: where the averaging operation is valid only in direction x 1 . From the above equation, the component θ i,j+1 was determined, which allowed the elimination of the unknown in the virtual node (i, j + 1) from Equation (15). The formulas describing the individual coefficients are shown below: The node of area 4 (for i equals m 1 and for j equal to m 2 ) is the node in the bottom right corner of the composite. In this area, the boundary conditions described by Equations (23) and (29) were considered. Therefore, the formulas describing the coefficients are as follows: The nodes of area 5 (for i from 2 to m 1 − 1 and for j equal to 1) represent the nodes on the left surface of the composite. It was assumed that this edge is thermally insulated: which leads to the condition: and allows to the elimination of the unknown in the virtual nodes (i, j − 1) from the equation. The formulas describing the individual coefficients are shown below: where the averaging operation is valid only along x 1 . The node of area 6 (for i equals m 1 and for j equal to 1) is the node in the bottom left corner of the composite. In this area, the boundary conditions described by Equations (23) and (21) are considered. Therefore, the formulas describing the coefficients are as follows: where the averaging operation is valid only along x 1 .
In each area, in the case of the macrotemperature θ, in the node (i, j), the coefficient related to the time coordinate was added: and it was assumed that in the initial time the values of all unknowns are known. There are also the fluctuation amplitudes of the temperature ψ 1 and ψ 2 as unknowns in Equation (15). The coefficients for the fluctuation amplitudes ψ 1 are as follows: • areas 3, 4 • areas 5, 6  Equation (16) is related to the fluctuation amplitudes of the temperature ψ1. Consid ering this equation, the composite nodes were grouped into areas; confer Figure 5. It is assumed that the fluctuation amplitudes of the temperature ψ1 will be known fo the nodes of area 4 (for i from 1 to m1 and j equal to m2) and area 5 (for i from 1 to m1 and equals 1), which translates to the fact that fluctuation amplitudes of the temperature ψ are known, and we did not write Equation (16) for these nodes. Analogously, as in th It is assumed that the fluctuation amplitudes of the temperature ψ 1 will be known for the nodes of area 4 (for i from 1 to m 1 and j equal to m 2 ) and area 5 (for i from 1 to m 1 and j equals 1), which translates to the fact that fluctuation amplitudes of the temperature ψ 1 are known, and we did not write Equation (16) for these nodes. Analogously, as in the case of Equation (15), the coefficients were grouped and written separately for each unknown. First, the coefficients on the macrotemperature were written as: Then, the coefficients on the fluctuation amplitudes of the temperature were written: • areas 2, • area 3, In each area, in the case of the fluctuation amplitudes of the temperature ψ 1 , in the node (i, j), the coefficient related to the time coordinate was added:

Equation for Fluctuation Amplitude of the Temperature ψ 2
Equation (17) is related to the fluctuation amplitudes of the temperature ψ 2 (x 1 , x 2 ). Considering this equation, the composite nodes were grouped into areas again; confer Figure 6. It is assumed that the fluctuation amplitudes of the temperature ψ2 will be known for the nodes of area 3 (for i equal to 1 and j from 2 to m2), area 4 (for i equals m1 and j from 2 to m2), and area 5 (for i from 1 to m1 and j equals 1), which translates to the fact that fluctuation amplitudes of the temperature ψ2 are known, and we did not write Equation (17) for these nodes. First, the coefficients of the macrotemperature were written: • area 1, • area 2, Then, the coefficients on the fluctuation amplitudes of the temperature were written: • area 1, • area 2, In each area, in the case of the fluctuation amplitude of the temperature ψ2, in the node (i, j), the coefficient related to the time coordinate was added: (85) It is assumed that the fluctuation amplitudes of the temperature ψ 2 will be known for the nodes of area 3 (for i equal to 1 and j from 2 to m 2 ), area 4 (for i equals m 1 and j from 2 to m 2 ), and area 5 (for i from 1 to m 1 and j equals 1), which translates to the fact that fluctuation amplitudes of the temperature ψ 2 are known, and we did not write Equation (17) for these nodes. First, the coefficients of the macrotemperature were written:

Vector of Free Terms
• area 2, Then, the coefficients on the fluctuation amplitudes of the temperature were written: In each area, in the case of the fluctuation amplitude of the temperature ψ 2 , in the node (i, j), the coefficient related to the time coordinate was added:

Vector of Free Terms
In the case of the free terms vector, the procedure is analogous to the matrix of coefficients. The composite nodes were grouped into areas; confer Figure 7.
In the case of the free terms vector, the procedure is analogous to the matrix of coefficients. The composite nodes were grouped into areas; confer Figure 7. The vector of free terms is connected with the boundary conditions; therefore, free terms do not appear in equations written for area 0. The free terms do not concern the nodes in areas 23 and 24, because at these nodes the equations are not used. According to the boundary conditions, described above, the values of the macrotemperature θ, are known for the top surface of the composite, so these values occurred as free terms in the individual equations written for the nodes of areas 1-4 and 16 (i − 1, j) and areas [20][21][22] The values of the fluctuation amplitudes of the temperature ψ1 are known for the left and right surfaces of the composite, so these values occurred as free terms in the individual equations written for the nodes of areas 4, 7, 11, and 16-18 (i, j; i−1, j; i + 1, j), for the nodes of areas 15 and 19 (i, j; i − 1, j), for the nodes of areas 1, 5, 8, and 12 (i, j − 1), and for the nodes of areas 3, 6, 10, and 14 (i, j + 1).
In turn, the values of the fluctuation amplitudes of the temperature ψ2 are known on the left, top, and bottom surfaces of the composite, so these values occur as free terms in the individual equations written for the nodes of areas 12-14, and 20-22 (i, j; i, j − 1; i, j + 1), for the nodes of areas 16-18 (i, j; i + 1, j; i − 1; j), for the node of area 19 (i, j; i − 1, j; i, j + 1), for the node of area 15 (i, j; i, j − 1), for the node of area 1 (i − 1, j; i, j − 1), for the nodes of area 5 (i, j − 1), for the node of area 8 (i, j − 1; i + 1, j), for the nodes of areas 2-4 (i − 1, j), and for the nodes of areas 9-11 (i + 1, j).
Additionally, for Equation (15) and several areas, there are additional components that need to be added to the free terms. These are associated with the equations expressing the boundary conditions: • areas 4, 7, and 11, • areas 12, 13, 14, and 19, The vector of free terms is connected with the boundary conditions; therefore, free terms do not appear in equations written for area 0. The free terms do not concern the nodes in areas 23 and 24, because at these nodes the equations are not used. According to the boundary conditions, described above, the values of the macrotemperature θ, are known for the top surface of the composite, so these values occurred as free terms in the individual equations written for the nodes of areas 1-4 and 16 (i − 1, j) and areas 20-22 ( The values of the fluctuation amplitudes of the temperature ψ 1 are known for the left and right surfaces of the composite, so these values occurred as free terms in the individual equations written for the nodes of areas 4, 7, 11, and 16-18 (i, j; i−1, j; i + 1, j), for the nodes of areas 15 and 19 (i, j; i − 1, j), for the nodes of areas 1, 5, 8, and 12 (i, j − 1), and for the nodes of areas 3, 6, 10, and 14 (i, j + 1).
In turn, the values of the fluctuation amplitudes of the temperature ψ 2 are known on the left, top, and bottom surfaces of the composite, so these values occur as free terms in the individual equations written for the nodes of areas 12-14, and 20-22 (i, j; i, j − 1; i, j + 1), for the nodes of areas 16-18 (i, j; i + 1, j; i − 1; j), for the node of area 19 (i, j; i − 1, j; i, j + 1), for the node of area 15 (i, j; i, j − 1), for the node of area 1 (i − 1, j; i, j − 1), for the nodes of area 5 (i, j − 1), for the node of area 8 (i, j − 1; i + 1, j), for the nodes of areas 2-4 (i − 1, j), and for the nodes of areas 9-11 (i + 1, j).

Solution
The Crank-Nicolson method was used [29]. The equations, written for the nodes, form a system of non-uniform linear equations that can be written as follows: where parameter α is equal to 0.5. The matrix of coefficients is divided into matrix of coefficients dependent on time coordinate K t and independent of time coordinate K. The free terms vector is denoted by Q, and the unknowns' vector is denoted by q. Each equation in Equations (15)- (17) was written in turn for each node. Solving the above system of equations allows us to know the distribution of the sought unknowns-the macrotemperature θ and the fluctuation amplitudes of the temperatures ψ 1 and ψ 2 . Then, it is possible to know the distribution of the total temperature Θ by using the micromacrotemperature assumption as shown in Equation (2).

An Example of an Application
The problem analyzed using the algorithm described above was a non-stationary, two dimensional problem of a heat conduction issue in a biperiodic composite with dimensions . The calculations were carried out for the number of the cells N 1 = N 2 = 10. In turn, the share of the first material in the cell along both directions x 1 and x 2 was equal to 1/5. By solving the system of non-uniform Equation (89), it was possible to obtain the plots and the maps of the distributions of the sought unknowns. The plot of the distribution of discretized macrotemperature θ = θ(x 1 , x 2 ) is shown in Figure 8, and as a map in Figure 9.    The plot of the distribution of discretized fluctuation amplitudes of the temperature ψ1 = ψ1(x1, x2) is shown in Figure 10, and the map in Figure 11. The plot of the distribution of discretized fluctuation amplitudes of the temperature ψ 1 = ψ 1 (x 1 , x 2 ) is shown in Figure 10, and the map in Figure 11.  The plot of the distribution of discretized fluctuation amplitudes of temperature ψ2 = ψ2(x1, x2) is shown in Figure 12, and the map in Figure 13.  The plot of the distribution of discretized fluctuation amplitudes of temperature ψ2 = ψ2(x1, x2) is shown in Figure 12, and the map in Figure 13. The plot of the distribution of discretized fluctuation amplitudes of temperature ψ 2 = ψ 2 (x 1 , x 2 ) is shown in Figure 12, and the map in Figure 13.  From the above figures, it can be observed that the values of the fluctuation amplitudes effected the values of total temperature most significantly near the surface for x1 = 0. Close by this surface, the fluctuation amplitudes represented almost 10% of the values of the macrotemperature. In the remainder of the composite, the values of the fluctuation amplitudes were close to zero and were directly related to the analyzed boundary conditions.

Conclusions
Considerations that were explored and the creation of the finite difference method algorithm allowed us to formulate the following conclusions:  From the above figures, it can be observed that the values of the fluctuation amplitudes effected the values of total temperature most significantly near the surface for x1 = 0. Close by this surface, the fluctuation amplitudes represented almost 10% of the values of the macrotemperature. In the remainder of the composite, the values of the fluctuation amplitudes were close to zero and were directly related to the analyzed boundary conditions.

Conclusions
Considerations that were explored and the creation of the finite difference method algorithm allowed us to formulate the following conclusions: Figure 13. Map of the distribution of fluctuation amplitudes of temperature ψ 2 .
From the above figures, it can be observed that the values of the fluctuation amplitudes effected the values of total temperature most significantly near the surface for x 1 = 0. Close by this surface, the fluctuation amplitudes represented almost 10% of the values of the macrotemperature. In the remainder of the composite, the values of the fluctuation amplitudes were close to zero and were directly related to the analyzed boundary conditions.

Conclusions
Considerations that were explored and the creation of the finite difference method algorithm allowed us to formulate the following conclusions: 1.
The heat conduction equation is an equation with noncontinuous coefficients with reference to the analyzed biperiodic structure.

2.
Tolerance modelling makes it possible to average the equations and consider the impacts of the microstructure size on the issues analyzed.

3.
The resulting equations are equations of many variables, and it was necessary to solve them numerically.

4.
The Crank-Nicolson method was used to solve the obtained system of non-uniform equations, which ensured convergence of the solutions.

5.
The created algorithm is universal and may allow one to analyze a biperiodic structure composing two materials with arbitrary material properties arranged as in Figure 1. 6.
The created algorithm may allow one to analyze an arbitrary value of the external temperature and the temperature of one on the surfaces of the composite. 7.
Changing the boundary conditions involves modifying the algorithm, because in the finite difference method, equations are written only for nodes, where the values of the sought unknowns are not known.