In the power plant industry, the turbine inlet temperature (TIT) plays a key role in the efficiency of the gas turbine and, therefore, the overall —in most cases combined —thermal power cycle efficiency. Gas turbine efficiency increases by increasing TIT. However, an increase of TIT would i ncrease the turbine component temperature which can be critical (e.g., hot gas attack). Thermal barrier coatings (TBCs) —porous media coatings —can avoid this case and protect the surface of the turbine blade. This combination of TBC and film cooling produce s a better cooling performance than conventional cooling processes. The effective thermal conductivity of this composite is highly important in design and other thermal/structural assessments. In this article, the effective thermal conductivity of a simpli fied model of TBC is evaluated. This work details a numerical study on the steady -state thermal response of two -phase porous media in two dimensions using personal finite element analysis (FEA) code. Specifically, the system response quantity (SRQ) under i nvestigation is the dimensionless effective thermal conductivity of the domain. A thermally conductive matrix domain is modeled with a thermally conductive circular pore arranged in a uniform packing configuration. Both the pore size and the pore thermal conductivity are varied over a range of values to investigate the relative effects on the SRQ. In this investigation, an emphasis is placed on using code and solution verification techniques to evaluate the obtained results. The method of manufactured solut ions (MMS) was used to perform code verification for the study, showing the FEA code to be second -order accurate. Solution verification was performed using the grid convergence index (GCI) approach with the global deviation uncertainty estimator on a serie s of five systematically refined meshes for each porosity and thermal conductivity model configuration. A comparison of the SRQs across all domain configurations is made, including uncertainty derived through the GCI analysis. NOMENCLATURE A = area α = porosity ε = error FS = factor of safety G = conductance matrix Γ = domain boundary h = characteristic mesh size H = mesh number i = index k = thermal conductivity K = thermal conductivity matrix L = cell length N = total quantity N = shape function Ω = domain p = order of accuracy P = vertex heat load vector q = heat flux vector Q = heat flow per unit length r = mesh ratio R = void radius ρ = energy balance residual ρ = energy balance residual vector S = energy source T = temperature T = triangle vertex temperature vector u = system response quantity U = uncertainty W = cell half-width x = x-coordinate y = y-coordinate Subscripts C = cold 2 Copyright © 20 19 by ASME eff = effective f = formal H = mesh number , hot i = index j = index L∞ = L∞ norm MMS = manufactured solution n = normal num = numerical O = observed t = triangle , transcendental v = vertex x = x-direction y = y-direction Superscripts * = dimensionless , global 1. INTRODUCTION The future of increased efficiency in gas turbine -driven power plant energy generation is heavily dependent on the increased temperature of the turbine inlet temperature (TIT) [1- 4]. In spite of the pursuit of higher TIT levels, virgin metal material in the turbine components are susceptible to a variety of aggr essive and high -consequence degradation or failure modes . Such issues include accelerated thermal creep, material degradation due to oxidation, and cycle fatigue [ 5,6]. Thermal - related corrosion issues arise at various temperature levels anywhere in the ra nge of 650 -1700 °C [2,7-11]. To mitigate thermal issues, thermal barrier coatings (TBCs) are used to protect the virgin material. Such coating s, including spray technologies are typically comprised of some form of ceramic composite material [12-14]. These composite materials are porous in nature and provide significant protection to the metal turbine components to allow for better overall performance and efficiency. Understanding the thermophysical properties of TBC materials is critical to the design and a dvancement of gas turbines for plant power generation for current and future installations. As systems continue to produce higher TIT levels, the thermal performance of advanced TBCs becomes more critical. This investigation evaluates the effective thermal transport in simplified models of TBC porous structures, specifically two -phase porous media, as an initial study into TBC material nano -structure effective thermal response. 2. SYSTEM DESCRIPTION In order to analyze the thermal transport phenomena within the TBC as a part the authors’ efforts , a two -dimensional system is constructed, consistent with the illustration given in Figure 1. For this study, a unit half -cell is given with height, L, and width, W, and is oriented in the Cartesian x-y plane with unit thickness out -of-plane . The bulk matrix material is defined as th e first material phase with thermal conductivity k1, where the pore phase is defined to have thermal conductivity k2. The pore pattern is circular, centered in the unit cell, with radius, R. Both the x=0 and x=W boundaries are given periodic boundary condi tions for system symmetry, and the y=0 and y=L boundaries have enforced hot and cold boundary temperatures, TH and TC, respectively. By enforcing Dirichlet boundary conditions at the hot and cold boundaries, heat flow —per unit length —is induced in the posi tive y-direction, perpendicular to each boundary, indicated by QH and QC, respectively. Figure 1. Unit half -cell nano -porous structure Of interest in this work is th e effect of relative pore size and relative pore thermal conductivity on the overall effective thermal conductivity, keff, of the TBC. Thus, the dimensionless parameter, porosity, α, is used to define the relative pore size, where 𝛼=(𝜋𝑅2)/(2𝑊𝐿) (1) The keff value is defined by overall heat flow through the system and the enforced boundary temperatures such that 𝑘𝑒𝑓𝑓=[(𝑄𝐻+𝑄𝐶)𝐿]/[2𝑊(𝑇𝐻−𝑇𝐶)]. (2) where the domain is assumed to have unit thickness out-of- plane. In theory, QH=QC, but because this study is numerical in nature, the average of the two heat flow values is used in Eq. (2) to define the overall induced heat flow through the domain. To account for variable material properties of the matrix and porous phases in different TBCs, the ultimate system response quantity (SRQ) of interest in this analysis is the dimensionless effective thermal conductivity, k*, which is merely the ratio of THTC k1 k2 QHQC L WR xy 3 Copyright © 20 19 by ASME the effective thermal conductivity with the matrix material thermal conductivity, simply e valuated as 𝑘∗=𝑘𝑒𝑓𝑓/𝑘1. (3) This study will show the effects on the SRQ of varying the k2-to-k1 ratio along with varying the α, using numerical verification techniques to quantify uncertainty in that value. Ratios in the two-phase thermal conducti vities are evaluated from 0 to 2, and α values are varied from 5% to 65%. 3. NUMERICAL APPROACH 3.1 DISCRETIZATION In order to determine the SRQ, the temperature distribution , T, and resultant heat source terms , S, must be solved to satisfy the governing partial differential equation (PDE), the heat equation, as expressed in Eq. (4) , where k is the thermal conductivity. 𝑘(𝜕 𝑇 𝜕𝑥 +𝜕 𝑇 𝜕𝑦 )+𝑆=0 (4) To numerically discretize the heat equation, unstructured triangular meshes were generated over the domain within each material phase . To accommodate the verification approaches used in this study, systematically refined meshes were generated for each porosity level, an example of which is shown in Figure 2a through Figure 2e, for mesh number, H, 5 through 1, respectively. Likewise, the fine st meshes for a selection of porosity models of 5%, 25%, 45%, and 65% are given in Figure 3a through Figure 3d, respectively. Each tth triangle is configured with three vertices such that the ith vertex is located at (xt,i,yt,i) with temperature Tt,I and area At, is shown in Figure 4. Each Hth mesh for a given porosity level is then described with a characteristic mesh size hH. If NH,t is the total number of triangles in mesh H, then ℎ𝐻=√1 𝑁 Σ𝑡=1𝑁 𝐴𝑡 (5) Thus, systematic mesh refinement here approximately halves the characteristic mesh size from mesh H+1 to mesh H. The domain discretization is used with a second order accurate Galerkin Finite Element Method (FEM) approach which is employed to solve Eq. (4) over the entire domain [15]. In short, Eq. (4) can be solved over each element , individually, using the discretized PDE matrix equation ∫(𝜵𝑵)𝑇𝑲𝜵𝑵𝑑𝛺𝑻𝛺=∫𝑵𝑇𝑆𝑑𝛺𝛺−∫𝑵𝑇𝒒𝒏𝑑𝛤𝛤 (6) where Ν is the linear shape function row vector used to interpo late field variables across a single element, Κ is the 2x2 thermal conductivity matrix, T is the vector of triangle vertex temperature values, qn, is the boundary normal heat flux on a triangle boundary, Ω is the element domain, and Γ is the element boundary. Furthermore, the different parts of Eq. (6) are condensed and defined as the conductance matrix, G, and the heat load matrix P, where a) b) c) d) e) 4 Copyright © 20 19 by ASME Figure 2. Systematic mesh refinement for α=25% Figure 3. Finest mesh for varying porosity levels Figure 4. Notional linear triangle element 𝑮=∫(𝜵𝑵)𝑇𝑲𝜵𝑵𝑑𝛺𝛺 (7) and 𝑷=∫𝑵𝑇𝑆𝑑𝛺𝛺−∫𝑵𝑇𝒒𝒏𝑑𝛤𝛤 (8) such that 𝑮𝑻=𝑷. (9) 3.2 SOLVER The meshes employed in this study were generated using open -source Gmsh software [16]. The custom code used for this study was written in Fortran as a generic two -dimensional FE M heat transfer solver. A basic successive over -relaxation method was used to iteratively update solu tions to Eq. (9) until some residual level exit criteria was met. The residual vector, ρ, for the system is computed as 𝝆=𝑮𝑻−𝑷 (10) where the L ∞ norm residual, ρL∞, defined as 𝜌𝐿∞=max|𝝆|. (11) A ρL∞ value of 10-8 was used as exit criteria for sol ution completion. 4. VERIFICATION APPROACH 4.1 CODE VERIFICATION In order to verify the formally second order accurate solution, the Method of Manufactured Solutions (MMS) was employed [17,18]. The MMS approach allows for a user - defined temperature distribution to be selected, TMMS, where the operator —in this case defined by the PDE in Eq. (4) —acts on the manufactured solution. The resulting source term, SMMS, falls out of the operator, where 𝑆𝑀𝑀𝑆=−𝑘(𝜕 𝑇 𝜕𝑥 +𝜕 𝑇 𝜕𝑦 ). (12) Thus, by applying SMMS and the boundary conditions that satisfy TMMS to the domain, the solved system of equations should converge to TMMS with mesh refinement. Here, TMMS is defined as 𝑇𝑀𝑀𝑆(𝑥 𝑦)=cos(2𝜋𝑥)sin(𝜋𝑦+0.75) (13) and is applied as a boundary condition at y=0 and y=L as Dirichlet conditions. It follows that the first partial derivatives are given to be a) b) c) d) Tt,3(xt,3,yt,3) AtTt,1(xt,1,yt,1) Tt,2(xt,2,yt,2) y x 5 Copyright © 20 19 by ASME 𝜕𝑇 𝜕𝑥=−2𝜋sin(2𝜋𝑥)sin(𝜋𝑦+0.75) 𝜕𝑇 𝜕𝑦=𝜋cos(2𝜋𝑥)cos(𝜋𝑦+0.75) (14) and are enforced at the boundaries where x=0 and x=W as Neumann boundary conditions. The resulting source term distribution is then 𝑆𝑀𝑀𝑆=5𝜋2𝑘cos(2𝜋𝑥)sin(𝜋𝑦+0.75). (15) For this study, the L ∞ norm metric is used on temperature error, εMMS ,H, over all Nv vertices in mesh H, where 𝜖𝑀𝑀𝑆 𝐻=max𝑖|𝑇𝑖−𝑇𝑀𝑀𝑆 𝑖| (16) based on the solution temperature, Ti, of the ith vertex. The observed order of accuracy for mesh H, pO,H, is determined by 𝑝𝑂 𝐻=ln(𝜖𝑀𝑀𝑆 𝐻/𝜖𝑀𝑀𝑆 𝐻+1)/ln(ℎ𝐻/ℎ𝐻+1). (17) 4.2 SOLUTION VERIFICATION The numerical error, Unum, for th e SRQ, k*, is approximated in this method by performing solution verification using a modern version of a more traditional grid convergence index (GCI) [18,19] approach that incorporates a global deviation estimator for the solution [20]. This approach uses an empirical trend to scale a factor of safety, FS, for Unum based on how closely the convergent observed order of accuracy is to the formal order of accurac y, thus acknowledging how close the solution is to the asymptotic solution region. The modified transcendental order of accuracy, pt, used for this study is given by 𝑝𝑡=ln*(𝑟1 2𝑝 −1)(|𝑢 −𝑢 𝑢 −𝑢 |)+𝑟1 2𝑝 +/ln(𝑟1 2𝑟2 3) (18) where uH represents the SR Q of interest as computed on mesh H, and 𝑟𝑖 𝑗=ℎ𝑗/ℎ𝑖. (19) For this problem, where the SRQ is k*, the global deviation parameter, Δp, is computed as 𝛥𝑝=min(|𝑝𝑓−𝑝𝑡| 4𝑝𝑓 0.95𝑝𝑓) (20) with pf being the formal order of accuracy of the solution method which, in this study, is 2. Lastly, the global order of accuracy, p*, used to compute the global deviation GCI FS is 𝑝∗=𝑝𝑓−𝛥𝑝 (21) such that 𝐹𝑆=3−1.9(𝑝∗/𝑝𝑓)8 (22) This FS value is then used to estimate the numerical uncertainty on the fine mesh SRQ solution, where 𝑈𝑛𝑢𝑚=𝐹𝑆|(𝑢2−𝑢1)/(𝑟1 2𝑝∗−1)| (23) 5. RESULTS AND DISCUSSI ON 5.1 CODE VERIFICATION Using the MMS solution and conditions described in Eq. (13) through Eq. (15) to perform the MMS verification , the trends of logarithmic error against logarithmic mesh size are shown in Figure 5, where the trend of a formal second order accurate solution is also shown in qualitative comparison. Likewise, Figure 6 shows the convergence of the observed order of accuracy with respect to mesh number. Figure 5. MMS error trend with mesh refinement Figure 6. MMS observed order of accuracy convergence Table 1 through Table 7 detail the specific mesh refinement and MMS observed order of accuracy parameters, showing an approximately second order accurate scheme. -8.0-7.0-6.0-5.0-4.0-3.0-2.0-1.00.0 -5.5 -3.5 -1.5ln(εMMS,H ) ln(h)5% 15% 25% 35% 45% 55% 65% Formal 1.01.52.02.53.03.5 0 1 2 3 4 5pO,H H5% 15% 25% 35% 45% 55% 65% Formal 6 Copyright © 20 19 by ASME Table 1. Systematic mesh refinement and MMS observed order of accuracy parameters for α=5% H rH,H+ 1 ln(hH) ln(εMMS,H ) pO,H 5 --- -2.24 -2.15 --- 4 1.46 -2.62 -2.98 2.18 3 2.05 -3.34 -4.22 1.73 2 1.93 -3.99 -5.46 1.89 1 1.93 -4.65 -6.71 1.89 Table 2. Systematic mesh refinement and MMS observed order of accuracy parameters for α=15% H rH,H+1 ln(hH) ln(εMMS,H ) pO,H 5 --- -2.11 -1.53 --- 4 1.66 -2.62 -2.60 2.11 3 1.99 -3.31 -3.91 1.90 2 1.90 -3.95 -5.46 2.42 1 2.02 -4.65 -6.62 1.65 Table 3. Systematic mesh refinement and MMS observed order of accuracy parameters for α=25% H rH,H+1 ln(hH) ln(εMMS,H ) pO,H 5 --- -2.11 -1.87 --- 4 1.70 -2.64 -2.59 1.35 3 2.01 -3.34 -4.12 2.19 2 1.92 -3.99 -5.54 2.18 1 1.96 -4.67 -6.69 1.71 Table 4. Systematic mesh refinement and MMS observed order of accuracy parameters for α=35% H rH,H+1 ln(hH) ln(εMMS,H ) pO,H 5 --- -2.19 -2.15 --- 4 1.61 -2.67 -3.02 1.80 3 1.88 -3.30 -3.84 1.31 2 1.93 -3.96 -5.32 2.25 1 2.01 -4.65 -6.63 1.88 Table 5. Systematic mesh refinement and MMS observed order of accuracy parameters for α=45% H rH,H+1 ln(hH) ln(εMMS,H ) pO,H 5 --- -2.19 -1.89 --- 4 1.64 -2.69 -3.15 2.52 3 1.81 -3.28 -4.11 1.63 2 2.03 -3.99 -5.52 1.99 1 1.97 -4.67 -6.65 1.66 Table 6. Systematic mesh refinement and MMS observed order of accuracy parameters for α=55% H rH,H+1 ln(hH) ln(εMMS,H ) pO,H 5 --- -2.19 -1.54 --- 4 1.67 -2.71 -2.89 2.63 3 1.80 -3.29 -4.32 2.44 2 2.01 -3.99 -5.33 1.44 1 1.95 -4.66 -6.65 1.98 Table 7. Systematic mesh refinement and MMS observed order of accuracy parameters for α=65% H rH,H+1 ln(hH) ln(εMMS,H ) pO,H 5 --- -2.28 -1.39 --- 4 1.53 -2.71 -2.70 3.08 3 1.86 -3.32 -4.13 2.31 2 1.91 -3.97 -5.45 2.06 1 1.98 -4.65 -6.68 1.79 5.2 SOLUTION VERIFICATIO N As mentioned previously, k* is computed using the average numerically -produced values of QH and QC, as described by Eq. (2) and Eq. (3). These heat flow values are computed by summing the residual values, implied from Eq. (10), of all vertices on the QH and QC boundaries where TH and TC are enforced, respectively. After performing the solution verific ation procedures described with Eq. (18) through Eq. (23), the following table s show the computed global deviation parameters and ultimate Unum approximation for the finest mesh in each material configuration model. Note that the Unum values given in these tables are expressed as percent values relative to the presented fine mesh k* values. All Unum values are obtained to be less than 0.2% of the fine mesh k* values. 7 Copyright © 20 19 by ASME Table 8. Global deviation GCI parameters for α=5% k2/k1 k* r1,2 pt Δp p* FS Unum 0.00 0.91 1.93 1.68 0.32 1.68 2.52 0.147% 0.25 0.94 1.93 1.67 0.33 1.67 2.54 0.072% 0.50 0.97 1.93 1.66 0.34 1.66 2.57 0.034% 0.75 0.99 1.93 1.64 0.36 1.64 2.61 0.013% 1.00 1.00 1.93 0.00 1.90 0.10 3.00 0.000% 1.25 1.01 1.93 1.60 0.40 1.60 2.68 0.008% 1.50 1.02 1.93 1.57 0.43 1.57 2.72 0.014% 1.75 1.03 1.93 1.55 0.45 1.55 2.75 0.018% 2.00 1.03 1.93 1.52 0.48 1.52 2.79 0.021% Table 9. Global deviation GCI parameters for α=15% k2/k1 k* r1,2 pt Δp p* FS Unum 0.00 0.74 2.02 1.89 0.11 1.89 1.79 0.097% 0.25 0.84 2.02 1.88 0.12 1.88 1.82 0.047% 0.50 0.90 2.02 1.88 0.12 1.88 1.85 0.023% 0.75 0.96 2.02 1.87 0.13 1.87 1.89 0.009% 1.00 1.00 2.02 0.00 1.90 0.10 3.00 0.000% 1.25 1.03 2.02 1.86 0.14 1.86 1.96 0.006% 1.50 1.06 2.02 1.85 0.15 1.85 1.99 0.009% 1.75 1.09 2.02 1.84 0.16 1.84 2.02 0.012% 2.00 1.11 2.02 1.83 0.17 1.83 2.06 0.014% Table 10. Global deviation GCI parameters for α=25% k2/k1 k* r1,2 pt Δp p* FS Unum 0.00 0.60 1.96 1.91 0.09 1.91 1.70 0.095% 0.25 0.74 1.96 1.91 0.09 1.91 1.70 0.044% 0.50 0.85 1.96 1.91 0.09 1.91 1.71 0.020% 0.75 0.93 1.96 1.90 0.10 1.90 1.72 0.007% 1.00 1.00 1.96 0.06 1.90 0.10 3.00 0.000% 1.25 1.06 1.96 1.90 0.10 1.90 1.74 0.004% 1.50 1.11 1.96 1.90 0.10 1.90 1.75 0.007% 1.75 1.15 1.96 1.89 0.11 1.89 1.77 0.009% 2.00 1.18 1.96 1.89 0.11 1.89 1.78 0.010% Table 11. Global deviation GCI parameters for α=35% k2/k1 k* r1,2 pt Δp p* FS Unum 0.00 0.48 2.01 1.95 0.05 1.95 1.45 0.090% 0.25 0.65 2.01 1.94 0.06 1.94 1.51 0.040% 0.50 0.79 2.01 1.93 0.07 1.93 1.57 0.019% 0.75 0.90 2.01 1.92 0.08 1.92 1.64 0.007% 1.00 1.00 2.01 0.09 1.90 0.10 3.00 0.000% 1.25 1.08 2.01 1.89 0.11 1.89 1.78 0.005% 1.50 1.15 2.01 1.88 0.12 1.88 1.84 0.008% 1.75 1.21 2.01 1.87 0.13 1.87 1.91 0.010% 2.00 1.26 2.01 1.85 0.15 1.85 1.98 0.012% Table 12. Global deviation GCI parameters for α=45% k2/k1 k* r1,2 pt Δp p* FS Unum 0.00 0.37 1.97 2.09 0.09 1.91 1.67 0.115% 0.25 0.57 1.97 2.05 0.05 1.95 1.45 0.038% 0.50 0.74 1.97 2.02 0.02 1.98 1.22 0.013% 0.75 0.88 1.97 1.98 0.02 1.98 1.21 0.005% 1.00 1.00 1.97 0.00 1.90 0.10 3.00 0.000% 1.25 1.11 1.97 1.92 0.08 1.92 1.61 0.004% 1.50 1.20 1.97 1.89 0.11 1.89 1.79 0.008% 1.75 1.28 1.97 1.86 0.14 1.86 1.94 0.011% 2.00 1.35 1.97 1.83 0.17 1.83 2.08 0.013% Table 13. Global deviation GCI parameters for α=55% k2/k1 k* r1,2 pt Δp p* FS Unum 0.00 0.28 1.95 1.95 0.05 1.95 1.44 0.130% 0.25 0.50 1.95 1.96 0.04 1.96 1.35 0.040% 0.50 0.69 1.95 1.97 0.03 1.97 1.30 0.016% 0.75 0.85 1.95 1.98 0.02 1.98 1.26 0.005% 1.00 1.00 1.95 0.00 1.90 0.10 3.00 0.000% 1.25 1.13 1.95 1.99 0.01 1.99 1.16 0.003% 1.50 1.25 1.95 2.00 0.00 2.00 1.11 0.004% 1.75 1.35 1.95 2.01 0.01 1.99 1.15 0.006% 2.00 1.45 1.95 2.01 0.01 1.99 1.20 0.007% 8 Copyright © 20 19 by ASME Table 14. Global deviation GCI parameters for α=65% k2/k1 k* r1,2 pt Δp p* FS Unum 0.00 0.18 1.98 2.02 0.02 1.98 1.26 0.168% 0.25 0.43 1.98 2.03 0.03 1.97 1.34 0.044% 0.50 0.64 1.98 2.03 0.03 1.97 1.29 0.016% 0.75 0.83 1.98 2.02 0.02 1.98 1.24 0.005% 1.00 1.00 1.98 0.10 1.90 0.10 3.00 0.000% 1.25 1.16 1.98 2.02 0.02 1.98 1.21 0.003% 1.50 1.30 1.98 2.02 0.02 1.98 1.25 0.005% 1.75 1.43 1.98 2.03 0.03 1.97 1.33 0.007% 2.00 1.56 1.98 2.05 0.05 1.95 1.42 0.009% 5.3 CHARACTERISTIC TREND S Compiling all of the k* from Table 8 through Table 14 shows a family of similar curves —roughly quadratic in form — relating k* to both α and k2/k1, as illustrated in Figure 7. Figure 7. Dimensionless effective thermal conductivity characteristic curves In Figure 7, it is clear that an increase in the pore material thermal conductivity leads to an increase in overall effective thermal conductivity of the composite material. Also, it is apparen t from Figure 7 that an increase in pore size amplifies the relative effect of k2 on k*. All curves intersect at (k2/k1=1.0,k*=1.0) , due to the fact that the condition of k2/k1=1.0 implies that the composite material behaves as a single homogeneous material from the perspective of steady - state thermal transport. N ote that error bars on the k* values are indistinguishable due to their relative ly low values. 5.4 FUTURE WORKS The study performed in this work has the potential to be extended to describe the effective thermal behavior in an analytical solution with prescribed uncertainty bounds. Likewise, the effects of pore topology and geom etry on effective thermal responses of composite thermal material could be investigated as well as obtaining optimal thermal conductivity for microporous insulations . Similarly , such a study could have implications on the uncertainty incurred by mesh geome try. Further investigations could be made into analysis of transient effects or the implications of more than one pore material type in a given matrix. 6. CONCLUSION This paper has presented a numerical approach to capturing the effective thermal transport behavior of a two -phase porous material with respect to pore phase thermal conductivity and pore size. It was shown that the employed FEM code was roughly second order accurate using the MMS code verification approach, and numerical uncertainty was estimat ed using the global deviation estimator GCI method. A trend was illustrated, tying the effective thermal conductivity of the overall medium to both relative pore size and relative pore thermal conductivity.