## 3 Hybrid Boundary-Type Meshless Method for Modeling Heat Conduction Problems  
In the polar coordinate system, the governing equation of the stationary heat conduction problems is shown as follows:  
\[\frac{\partial^{2}u}{\partial\rho^{2}}+\frac{1}{\rho}\frac{\partial u}{ \partial\rho}+\frac{1}{\rho^{2}}\frac{\partial^{2}u}{\partial\theta^{2}}=0 \tag{4}\]  
where \(\rho\) is the radius and \(\theta\) is the polar angle in the polar coordinate system. In the proposed method, the approximation is established via superpositioning a series of particular non-singular basis functions of the heat equation in terms of many source points located within the domain. An approximation solution of two-dimensional stationary heat conduction problems can be expressed as a linear combination of particular non-singular basis function with unknown coefficients, \(\mathbf{a}_{j,k}\). \[u(\mathbf{x})\approx\sum_{j=1}^{O}\sum_{k=1}^{P}\mathbf{a}_{j,k}\mathbf{H}(x, y_{j}) \tag{5}\]  
where \(\mathbf{x}=(\rho,\theta)\) and \(\mathbf{x}\in\partial\Omega\) is the spatial coordinate that is collocated on the boundary, \(y_{j}=(\rho_{j},\theta_{j})\) is the source point, \(O\) is the number of source points, \(P\) is the order of the particular non-singular basis function, \(\mathbf{a}_{j,k}=[ \begin{array}{cc}c_{jk}&d_{jk}\end{array} ]\) is a vector of unknown coefficients to be determined, \(\mathbf{H}(\mathbf{x},\mathbf{y}_{j})\) is the selected combination of the particular solutions of the heat equation in the two-dimensional polar coordinate system. The particular non-singular basis function of the heat equation can be expressed as  
\[\mathbf{H}(\mathbf{x},\mathbf{y}_{j})=\Big{[} \begin{array}{cc}(\rho_{j}/R)^ {k}\cos(k\theta_{j})&(\rho_{j}/R)^{k}\sin(k\theta_{j})\end{array} \Big{]}^{T} \tag{6}\]  
In this study, we adopt the characteristic length [38] as  
\[R=1.5\times\max(\rho) \tag{7}\]  
where \(\max(\rho)\) is the maximum radial distance in the problem domain. Applying the Dirichlet boundary condition, we obtain  
\[u(\mathbf{x}_{l})\approx\sum_{j=1}^{O}\sum_{k=1}^{P}\mathbf{a}_{j,k}\mathbf{H }_{j,k}(\mathbf{x}_{l},\mathbf{y}_{j})=g(\mathbf{x}_{l}),\;l=1,...,Q \tag{8}\]  
where \(g(\mathbf{x}_{l})\) is the Dirichlet boundary value imposed at the boundary collocation point. \(Q\) is the number of boundary collocation points. The Neumann boundary condition can be written as  
\[\frac{\partial u(\mathbf{x}_{l})}{\partial n}=\frac{\partial u(\mathbf{x}_{ l})}{\partial x}n_{x}+\frac{\partial u(\mathbf{x}_{l})}{\partial y}n_{y} \tag{9}\]  
where  
\[\frac{\partial u(\mathbf{x}_{l})}{\partial x}=\frac{\partial u(\mathbf{x}_{ l})}{\partial\rho}\frac{\partial\rho}{\partial x}+\frac{\partial u(\mathbf{x}_{ l})}{\partial\theta}\frac{\partial\theta}{\partial x}\text{ and }\frac{\partial u(\mathbf{x}_{l})}{\partial y}=\frac{\partial u(\mathbf{x}_{l})}{ \partial\rho}\frac{\partial\rho}{\partial y}+\frac{\partial u(\mathbf{x}_{l} )}{\partial\theta}\frac{\partial\theta}{\partial y} \tag{10}\]  
where \(n_{x}\) and \(n_{y}\) are normal vectors in the \(x\) and \(y\) direction, respectively. \(\partial u/\partial\rho\) and \(\partial u/\partial\theta\) can be rewritten as follows:  
\[\frac{\partial u(\mathbf{x}_{l})}{\partial\rho}=\sum_{j=1}^{O}\sum_{k=1}^{P} \mathbf{a}_{j,k}[ \begin{array}{cc}k(1/R)^{k}\rho_{j}^{k-1}\cos(k\theta_{j })&k(1/R)^{k}\rho_{j}^{k-1}\sin(k\theta_{j})\end{array} ]^{T} \tag{11}\]  
\[\frac{\partial u(\mathbf{x}_{l})}{\partial\theta}=\sum_{j=1}^{O}\sum_{k=1}^{P} \mathbf{a}_{j,k}[ \begin{array}{cc}-k(1/R)^{k}\rho_{j}^{k}\sin(k\theta_{j})&k( 1/R)^{k}\rho_{j}^{k}\cos(k\theta_{j})\end{array} ]^{T} \tag{12}\]  
For the modeling of multi-layer composite materials, we adopted the DDM [37]. The DDM splits the domain into smaller subdomains in which subdomains are intersected only at the interface between layers. For each subdomain, there exists an independent layer with its own thermal conductivity. The boundary and source points are collocated in each subdomain.