## 2 Materials and Methods
### Finite Element Formulation  
In this section, the adopted numerical method based on the finite element method (FEM) is employed for numerical investigation of characteristics of flow and heat transfer in the diamond-shaped pinned DPH. FEM was introduced to find the solution of structural engineering problems. With the passage of time, the development of this method made it so powerful that it is extensively employed to solve governing partial differential equations of the problems from various fields, such as fluid mechanics, wave propagation, reaction-diffusion and convection-diffusion processes, and in many other fields. The basic steps involved in the computation of elliptic PDE using FEM are: describing the numerical domain; generating the triangular mesh in the computational domain; discretizing the given equation; the corresponding boundary conditions are adjusted to form the system of equations. For the present problem, the \(h\)-adaptive procedure was used to compute more accurate solutions. It may be noted that the governing momentum and energy equations of this present heat transfer problem are elliptic in nature. FEM formulation of an elliptic model PDE is presented as follows. Consider the elliptic model PDE,  
\[\
abla\ (c\
abla u)+au=f \tag{8}\]  
in the domain \(D\), where \(c\), \(a\), \(u\), and \(f\) are functions defined on \(D\). The Dirichlet and generalized Neumann conditions are imposed on the solid surface and the line of symmetry, respectively. These are mathematically expressed by \(hu=r\) and \(\overrightarrow{n}\ (c\
abla u)+qu=g\), respectively, on the edge \(\Gamma\) of the domain \(D\), where \(h\), \(r\), \(q\), and \(g\) are defined on the boundary \(\Gamma\) and \(\overrightarrow{n}\) is the unit normal. The Dirichlet condition can be approximated from the generalized Neumann boundary condition, and thus, the latter is assumed on the whole boundary, \(\Gamma\). Multiplying Equation (8) with a test function \(\psi\) and integrating on \(D\), gives  
\[\int\limits_{D}((-\
abla\ (c\
abla u))\psi+au\psi)dx=\int\limits_{D}f\psi dx \tag{9}\]  
Integrating by parts yields  
\[\int\limits_{D}((c\
abla u)\ \
abla\psi+au\psi)dx-\int\limits_{\Gamma}( \overrightarrow{n}\ (c\
abla u)\Big{)}\psi ds=\int\limits_{D}f\psi dx \tag{10}\]  
Implementing the boundary condition:  
\[\int\limits_{D}((c\
abla u)\ \
abla\psi+au\psi)dx-\int\limits_{\Gamma}(-qu+g) \psi ds=\int\limits_{D}f\psi dx \tag{11}\]  
The weak form of the model equation is as below, after simplification:  
\[\int\limits_{D}((c\
abla u)\ \
abla\psi+au\psi)dx-\int\limits_{D}f\psi dx-\int \limits_{\Gamma}(-qu+g)\psi ds=0,\ \ \ \ \ \forall\psi \tag{12}\]  
where \(u\) and \(\psi\) belong to finite-dimensional subspace \(V_{M}\) of a function space \(V\), which is the collection of those square-integrable continuous functions that disappear on the Dirichlet boundary. The approach of choosing \(u\) and \(\psi\) from the same subspace is called the [PERSON] approach. Hence, \(u\) can be expressed as \(u(x)=\sum_{j=1}^{M}u_{j}\psi_{j}(x)\), where \(u_{j}\) and \(\psi_{j}\) are the scalar coefficients and linear basis functions, respectively. After substituting the solution approximation \(u(x)\), the resulting system of an algebraic equation is given by  
\[\sum\limits_{j=1}^{M}\left(\int\limits_{D}((c\
abla\psi_{j})\ \
abla\psi_{i}+au\psi_{j}\psi_{i})dx+\int\limits_{\Gamma}q\psi_{j}\psi_{i}ds \right)u_{j}=\int\limits_{D}f\psi_{i}dx+\int\limits_{\Gamma}g\psi_{i}ds,\ \ \ \ \ \forall\psi \tag{13}\]  
where the system (13) is represented in the matrix notation as \((K+M+Q)U=F+G\), where \(K_{ij}=\int\limits_{D}(c\
abla\psi_{j})\ \
abla\psi_{i}dx\), \(M_{ij}=\int\limits_{D}a\psi_{j}\psi_{i}dx\), \(Q_{ij}=\int\limits_{\Gamma}q\psi_{j}\psi_{i}ds\), \(F_{ij}=\int\limits_{D}f\psi_{i}dx\), and \(G_{ij}=\int\limits_{\Gamma}g\psi_{i}ds\). To compute the solution of the problem, the formulation (13) is considered with the relevant boundary conditions.