# Design of Novel Cooling Systems Based on Metal Plates with Channels of Shapes Inspired by Nature
## 2 Theory
### Finite Differences Discretization  
The diffusion equation, Equation (1), was solved using a transient finite difference implicit numerical scheme. For this purpose, the cooling plate was divided in cells of size \(\Delta x\) in the horizontal direction, and \(\Delta y\) in the vertical direction, so that the temperature can be evaluated at \(N\) nodes on the surface of the cooling plate. In view of the square shape of the cooling plate and for the sake of simplicity of the numerical model, \(\Delta x\) and \(\Delta y\) were chosen of the same size, \(\Delta x=\Delta y\). Hence, the total number of nodes is distributed across the plate in a square mesh \(N=N_{r}N_{r}\). As indicated in Section 2.1, the temperature distribution is considered uniform in the direction perpendicular to the cooling plate. Nevertheless, the size \(\Delta z\) of the plate thickness direction shall be accounted for in the numerical model. Therefore, the diffusion equation, Equation (1), is solved using a transient finite difference scheme in a total of \(N\) nodes along the surface of the cooling plate following Equation (11):  
\[\rho\Delta x\Delta y\Delta zc_{p}\frac{T_{m,n}^{(p+1)}-T_{m,n}^{(p)}}{\Delta t }=q_{\rm k}^{(p+1)}+q_{\rm h,a}^{(p+1)}+q_{\rm h,w}^{(p+1)}+q_{\rm in} \tag{11}\]  
where the superscript \(p\) indicates the time instant and the subscripts \(m\) and \(n\) refer, respectively, to the nodes in the horizontal and vertical directions of the cooling plate. Consequently, node (1, 1) corresponds to the upper left corner of the plate and node (\(N_{r}\), \(N_{r}\)) denotes its lower right corner. In Equation (11), the heat transferred by conduction, convection, and the heat entering the system are discretized for each node. For each node in the system, the heat transferred by conduction to the node can be expressed as the conduction contribution from any adjacent node as follows:  
\[q_{k}^{(p+1)}=q_{k_{m-1,n}}^{(p+1)}+q_{k_{m+1,n}}^{(p+1)}+q_{k_{m,n+1}}^{(p+1) }+q_{k_{m,n-1}}^{(p+1)} \tag{12}\]  
where conduction can be evaluated considering Fourier's law and the lateral area through which conduction occurs. For instance, for the node situated in the right-hand-side of node (\(m\), \(n\)), i.e., node (\(m+1\), \(n\)), the contribution to Equation (12) would read:  
\[q_{k_{m+1,n}}^{(p+1)}=k\Delta y\Delta z\frac{T_{m+1,n}^{(p+1)}-T_{m,n}^{(p+1)}} {\Delta x} \tag{13}\]  
Convection through the top and bottom surfaces of the cooling plate can be also expressed in a finite difference form as:  
\[q_{\rm h,a}^{(p+1)}=h_{\rm h,a}\Delta x\Delta y\Big{(}T_{a}-T_{m,n}^{(p+1)} \Big{)} \tag{14}\]  
where \(h_{\rm h,a}=h_{\rm h,at}\) is at the top surface of the plate and \(h_{\rm h,a}=h_{\rm h,ab}\) is at its bottom surface. In cases in which the (\(m\), \(n\)) node is located in the lateral or corner cells of the cooling plate, natural convection through the top and bottom faces would occur only in half or one quarter of the surface of the node, respectively, and lateral vertical convection appears. This additional term can be expressed, for example, for a side lateral node as:  
\[q_{\text{h,a}}^{(p+1)}=h_{\text{h,a}}\Delta y\Delta z\Big{(}T_{\text{a}}-T_{m,n}^ {(p+1)}\Big{)} \tag{15}\]  
where \(\Delta y\Delta z\) is the surface through which heat is dissipated. Finally, nodes in which the water channel is present have an additional forced convection term due to the water cooling, \(q_{\text{h,w}}^{(p+1)}\). In this case, the heat transfer surface between the node and the channel is \(2\Delta y\Delta z\), neglecting the curvature of the channel, since the channel is in the middle plane of the plate.