# An Efficient Meshless Numerical Method for Heat Conduction Studies in Particle Aggregates
## 2 Materials and Methods
Upon calculating the conductivity gradient, the nodes that will remain near the interface (\(\
abla K\
eq 0\)) will be the ones from the aforementioned resolution. Thus, the grid is constructed automatically and speedily, without user intervention at any creation stage. This process can be continued for successively higher discretizations, doubling the discretization locally every time. The integration part is of key importance for the success of the MLPG method. The size and shape of the sectors on which the weight function is non-zero may vary. Using the MLS approach and the DC PSE method for the trial function, and based on the concept of the MLPG method, the weight function in each local integral can be selected in several ways [11]. Each weight function offers a specific expression of the MLPG method. It has been reported [10; 11] that the most promising expression is the MLPG which uses the step function as a weight function for the integration, which is defined in the space \(\Omega_{s}\) of the integration. In this case, the spatial integrals in the \(\Omega_{s}\) are avoided on account of the Gauss theorem. The overall matrix includes surface integrals only, which makes this method more efficient and attractive, with the solution being faster and more accurate than most other MLPG implementations, due to the simplicity of the calculations and the stability of the resulting solutions. Apart from the form of the weight function, the shape of the integration domain is also critical to the success of the method. Notable increase of the stability and robustness of the method has been reported using rectangular type domains in 2D problems [18]. Usually, equidistant domains around each node are used in MLPG implementations, with a radius around 60% of the local grid distance. Cubic integration domains around each node are used in the present 3D cases and, thus, the whole domain is covered without overlapping. During refinement in regions close to the interfaces, overlapping cannot be avoided without significant increase in the computational time. However, the overlapping is effectively restricted in small areas and does not affect the stability of the method. Figure 2 shows a domain comprising of a test sphere, Figure 2a, and the resulting computational grid, with integration sites around each node in the inlet, Figure 2b, consisting of cubes equal to the local grid spacing of each node. During numerical integration, one can have increased accuracy of the results either by increasing the order of the Gaussian quadrature integration, or by dividing the field of integration into smaller parts while keeping the integration order low. It was reported that the algorithm which computes integrals using subdivisions was more efficient than the one that integrates across the support domain with higher quadrature order [10; 11]. Thus, the surface integrals for the cubic domains are calculated here as the sum of the integrals on the cube sides, and 3-point Gaussian quadrature is used in each side integral. For heat conduction problems, it has been proven [18] that first-order interpolation is computationally efficient, fast and accurate. The optimal size of the support domain for the interpolation is ca 1.7 times the local grid spacing. The process of the grid generation, the construction of the shape functions, and the calculation of the integrals are performed here using Matlab kernels. The biconjugate gradient stabilized method (BiCGSTAB) [25] was used for the numerical solution of the resulting linear system. Figure 1: Nodes contained in the support domain (red squares) of three representative nodes, located at the boundary of the refinement (red circles). For validation purposes, the thermal conductivity was calculated for various particulate systems using Comsol Multiphysics(r) and in-house developed automation Matlab(r) scripts for the generation of the CAD geometry.