Thermo-Elastic Homogenization of TPMS-Based Interpenetrating Phase Composites: A Comparison Between Finite Element and Mechanics of Structure Genome Approaches

1. Introduction

Lightweight, multifunctional materials are in growing demand across the aerospace, automotive, and energy sectors. Structural components in these sectors increasingly need high specific stiffness together with tailored thermal transport and dimensional stability. Examples include structural heat exchangers, battery thermal-management housings, and precision optomechanical mounts, all of which must combine load-bearing capacity with controlled heat flow and minimal thermal distortion. Architected cellular materials and metamaterials have emerged as a design paradigm for meeting these coupled demands. Their effective properties depend not only on the constituent materials but also on a freely tunable unit-cell topology. This gives a much larger design space than conventional stochastic foams or fiber-reinforced composites offer [1,2,3]. However, realizing this potential requires predictive, physically grounded models. These models must resolve the coupled elastic and thermal response of such architectures across the topology and density ranges now capable of being manufactured through Additive Manufacturing (AM).

Among architected topologies, those based on TPMS are a particularly promising class. Their smooth, doubly curved, non-self-intersecting interfaces minimize local stress and heat-flux concentrations relative to strut-based lattices. TPMS geometries were first identified mathematically by Schwarz [4] and Neovius [5] in the nineteenth century. Schoen [6] later classified these surfaces, establishing the gyroid, diamond, and further surfaces now in widespread engineering use. The same morphologies recur in biological and self-assembled systems, from butterfly-wing scales to block-copolymer mesophases, suggesting that evolutionary and thermodynamic optimization alike favor them as load- and transport-bearing interfaces [7]. TPMS RUCs are represented implicitly through trigonometric level-set functions. They can be instantiated as either sheet or solid (network) architectures and continuously tuned in relative density, giving a single geometric framework that spans a wide range of stiffness, energy-absorption, and multifunctional design targets [8]. This tunability is not merely empirical. For the equal-volume-fraction Schwarz P and D surfaces, Jung and Torquato [9] showed that the bi-continuous morphology is simultaneously near-optimal for fluid and heat/electrical transport and, by geometric analogy, for elastic stiffness under a compliant second phase. This provides a rigorous physical basis for TPMS interfaces as natural candidates for multifunctional composite design, rather than an arbitrary geometric choice.

This last observation motivates a natural evolution in TPMS design. Single-phase (sheet or solid) TPMS lattices leave one region void. TPMS-based IPCs instead fill both regions defined by the TPMS level set, with distinct, mutually continuous solid phases. In particulate or short-fiber composites, the reinforcing phase is discontinuous and load transfer relies on interfacial shear. An IPC is different: each phase is independently self-supporting and load-bearing, a topological feature already recognized in early micromechanical treatments of interpenetrating multiphase media [10]. When the TPMS surface itself defines the phase interface, this co-continuity yields smooth, low-concentration load paths between phases. It also improves damage tolerance through crack deflection and delocalization at the curved interface. It also gives access to property combinations that are difficult to realize with either constituent alone, for example, simultaneously high stiffness and high toughness, or tailored stiffness-conductivity pairings [11,12]. Recent work has extended this concept toward increasingly application-relevant configurations. These include spinodal and stochastic TPMS-like topologies with tunable strength-toughness trade-offs [13], and centimeter-scale carbon- and polymer-based architected IPCs. In the latter, matrix infiltration was shown to distribute stress effectively and delocalize failure, achieving specific energy-absorption levels competitive with wound fiber tubes [14]. A recent comprehensive review consolidates this rapidly growing body of work on TPMS-based IPCs, spanning structural, thermal, and even piezoelectric functionalities. It also highlights that quantitative, method-verified predictions of their coupled thermo-elastic and transport properties remain comparatively scarce [15].

Reliable exploitation of this design space requires efficient, accurate homogenization frameworks. These frameworks must predict the effective properties of a RUC as functions of topology and relative density. The dominant benchmark approach is Finite-Element (FE) homogenization of the RUC under Periodic Boundary Conditions (PBCs). It has been applied extensively to TPMS lattices to predict effective elastic moduli [16,17,18]. Complementary numerical strategies include voxel-based asymptotic homogenization [19] and Fast-Fourier-Transform-based homogenization, the latter achieving comparable accuracy at substantially reduced computational cost for periodic microstructures [20]. FE-RUC homogenization is accurate, but it requires a distinct system factorization for every independent load case, typically six for the full elastic stiffness tensor alone. Its cost therefore scales unfavorably when multiple physics and many topology/density combinations must be swept, as is the case for the multifunctional design studies motivating this work.

The MSG offers a computationally efficient alternative. It was developed by Yu and co-workers from its predecessor, the Variational Asymptotic Method for Unit Cell Homogenization (VAMUCH) [21], which our group has also applied and extended [22,23]. MSG casts homogenization as a single variational statement, solved once per physics on the Structure Genome (SG). For three-dimensionally periodic media, this recovers the same effective properties as RUC-based FE analysis, while requiring only inexpensive back-substitutions for each additional load case [24,25,26]. The theoretical equivalence and comparative accuracy of the two approaches were established for single-phase TPMS lattices in our companion study [27]. That work built on earlier applications of MSG to the effective elastic and thermal-conductivity properties of particulate and nano-reinforced composites [28,29]. Practical FE-RUC implementations, in turn, rely on established periodic-boundary-condition and finite-element-analysis frameworks for composite materials [30,31].

The elastic homogenization of TPMS architectures is now comparatively well established. Their thermal characterization is not: both thermal expansion and thermal conductivity have received considerably less systematic attention, despite being equally central to the multifunctional applications motivating this work. Existing studies of TPMS thermal conductivity show that effective conductivity is governed primarily by relative density and topology. This has been demonstrated for pure TPMS lattices manufactured by laser powder bed fusion [32], TPMS-structured phase-change materials for battery thermal management [33], ceramic TPMS lattices [34], and TPMS-based porous media more generally [35]. A smaller body of work extends this analysis to two-phase TPMS foams and IPCs, where effective conductivity and elastic moduli depend jointly on topology and phase contrast [36,37]. The underlying MSG-based thermal-conductivity homogenization framework has itself been established for composite materials with general interfacial conditions, including imperfect interfaces [29]. The CTE of TPMS-based IPCs is still markedly less explored. To our knowledge, only a handful of studies address it explicitly. The most notable is a numerical investigation of the combined thermo-elastic response of gyroid-based composites, including the role of phase-modulus contrast and Poisson’s ratio [38]. To the best of our knowledge, no study to date has jointly addressed effective elastic properties, CTE, and thermal conductivity for the same set of TPMS-based IPC topologies. Still less has any study compared FE-RUC and MSG predictions of all three property groups over a broad relative-density range.

This gap is significant. TPMS lattices, TPMS-based IPCs, FE-RUC homogenization, and the MSG/VAMUCH framework have each been well studied individually. Yet no systematic investigation has compared FE-RUC and MSG predictions of effective elastic properties, CTE, and thermal conductivity for multiple TPMS-based IPC topologies over a broad range of relative densities. Closing this gap matters for two reasons. First, MSG offers an intrinsic computational advantage for multi-physics, multi-configuration design studies; validating it against FE-RUC would establish it as a lower-cost substitute for the multifunctional design of TPMS-based IPCs. Second, it is a prerequisite for understanding how topology and relative density govern the coupled stiffness-CTE-conductivity trade-offs relevant to lightweight thermal-management and dimensionally stable structural applications.

To address this gap, we develop two-phase RUC models with conforming meshes for four TPMS-based IPC topologies: gyroid, diamond (Schwarz D), PMY, and F-Rhombic Dodecahedron (F-RD). Relative densities ranging from $0.1$ to $0.5$ are considered. For each configuration, the effective elastic properties, coefficients of thermal expansion, and thermal conductivities are determined using two homogenization approaches applied to the same discretization: FE-RUC in Abaqus, with periodic boundary conditions imposed through the Abaqus Micromechanics Plugin [30,39], and MSG/VAMUCH through the in-house CmbsFE implementation. The MSG predictions are then quantitatively assessed against the FE-RUC results for all three groups of effective properties. In doing so, the present work extends the purely elastic, single-phase comparison reported in our companion study [27] to the thermo-elastic and thermal-transport behaviour of two-phase IPCs. Finally, the combined influence of topology and relative density on stiffness, thermal expansion, and conductivity is examined to identify topology-dependent trade-offs relevant to the design of lightweight structures with coupled mechanical and thermal requirements.

The remainder of the paper is organized as follows. Section 2 describes the generation of the TPMS-based IPC RUCs and the FE- and MSG-based homogenization procedures. Section 3 presents and discusses the resulting effective properties, including validation against the literature and a quantitative comparison of the two homogenization approaches. Section 4 presents concluding remarks.

2. Methodology

This section describes the procedures used to generate the TPMS-based IPC representative unit cells, the constituent phase properties, and the two homogenization strategies compared in this work, namely FE-based RUC homogenization and MSG-based homogenization. For each strategy, the effective elastic properties, CTEs, and thermal conductivities are obtained following analogous, physically consistent procedures, so that the two approaches can be compared on an equal footing.

2.1. Generation of TPMS-Based Interpenetrating Phase Composite Unit Cells

The TPMS-based IPC microstructures are generated with the open-source Python library Microgen [40], following the same modeling workflow adopted in the companion study on single-phase TPMS architectures [27]. Four TPMS topologies are investigated, matching those considered in the companion study: the gyroid, the diamond (Schwarz D) surface, the PMY surface, and the F-Rhombic Dodecahedron (F-RD, Schoen F-RD) surface. Each topology is implicitly defined, within a cubic cell, by a level-set function $\varphi (x_1,x_2,x_3)=0$. In Microgen, this function is expressed in terms of dimensionless coordinates $x_i=2\pi y_i/L$ ($i=1,2,3$), where $y_i$ are the physical Cartesian coordinates of a point $\mathbf{y}$ within the cubic RUC of edge length L. This is the same local position vector $\mathbf{y}$ used to describe the RUC/SG domain throughout Section 2.3 and Section 2.4:

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{gyroid}}=& \sin x_1\cos x_2+\sin x_2\cos x_3+\sin x_3\cos x_1,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{diamond}}=& \sin x_1\sin x_2\sin x_3+\sin x_1\cos x_2\cos x_3+\cos x_1\sin x_2\cos x_3+\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & +\cos x_1\cos x_2\sin x_3,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{PMY}}=& 2\cos x_1\cos x_2\cos x_3+\sin 2x_1\sin x_2+\sin x_1\sin 2x_3+\sin 2x_2\sin x_3,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{F}-\mathrm{RD}}=& 4\cos x_1\cos x_2\cos x_3-\left(\cos 2x_1\cos 2x_2+\cos 2x_2\cos 2x_3+\cos 2x_3\cos 2x_1\right),\hfill \end{array}$$

(5)

each topology being generated, in the RUC, as the zero level set $\varphi =0$ of the corresponding function. A level-set sampling resolution of 25 was used to generate all RUCs, except for the F-RD topology at $\rho =0.1$. That case required a resolution of 40. F-RD has a strongly asymmetric level-set range, unlike the other three topologies. At low relative density, this makes the solid network sensitive to under-resolution: the coarser sampling produced a mesh with a higher actual volume fraction than intended. The finer resolution corrects this.

This same configuration, F-RD at $\rho =0.1$, required a second exception. All analyses in this study were run on a workstation with an Intel Core i7-8665U CPU and 48 GB of RAM. Its quadratic-element system exceeded the memory available to the MUMPS direct solver used by CmbsFE, even on this workstation. The MSG analysis for this configuration therefore used linear (first-order) tetrahedral elements, rather than the quadratic elements used elsewhere in this study. For a fair comparison, the FE analysis for this same configuration was also rerun with linear elements. Every other topology-density combination in this study uses quadratic elements for both FE and MSG. Section 3.3 discusses the effect on FE-MSG agreement.

Unlike the single-phase (sheet or solid) TPMS lattices considered previously [27], the microstructures investigated here are two-phase IPCs. For a given topology, the corresponding level-set function of Eqs. (1)–() partitions the cubic RUC of volume V into two fully interpenetrating, co-continuous solid phases. These two phases jointly fill the entire RUC volume, with no void space. The solid TPMS network (phase 1) is generated as a Microgen Tpms object with a prescribed relative density $\rho =V_1/V$, where $V_1$ is the volume enclosed by the TPMS solid region and V is the volume of the cubic RUC bounding box. This matches the terminology of the companion study [27] and Microgen’s own density parameter; $\rho$ is varied in the range $\rho =0.1$–$0.5$ for each of the four topologies. The complementary phase (phase 2, denoted the “matrix” phase) occupies the remaining volume $V_2=V-V_1$, so that $\rho +(1-\rho )=1$. The two solid, conforming phases are obtained using the Boolean cutting procedure implemented in Microgen’s cut_phases function. This function subtracts the TPMS solid region from the cubic box and returns two topologically complementary, non-overlapping solids, sharing a common internal interface coincident with the TPMS level-set surface.

Figure 1 shows the solid TPMS network of each of the four topologies, extracted at a representative relative density ($\rho =0.3$) directly from the level-set functions of Eqs. (1)–(). The color code (gyroid, diamond, PMY, F-RD) is used consistently for these four topologies throughout the remainder of this paper. In the two-phase IPC, the complementary matrix phase occupies the remaining volume of each unit cell, sharing the smooth, doubly curved interface visible in Figure 1.

A conforming, periodic tetrahedral mesh of the two-phase RUC is subsequently generated with Gmsh [41], through Microgen’s mesh_periodic utility. The two-phase geometry is first imported into Gmsh via the OpenCASCADE kernel. The two phases are meshed as separate solids sharing a common interface, so a Boolean fragment operation is applied next. This operation shares (rather than duplicates) the interface faces, edges, and nodes between the two phases, ensuring a conforming mesh across the phase boundary. Periodicity is then enforced geometrically, prior to meshing. For each of the three RUC axes, boundary faces on the negative side of the cube are paired with their corresponding faces on the positive side, identified by matching bounding boxes offset by the RUC edge length. Each pair is declared periodic to Gmsh via an affine (translational) transformation. This constrains the mesher to generate identical node distributions on corresponding master and slave faces. The volume mesh is then generated for a prescribed target element size, using tetrahedral elements optimized with Netgen [42]. The element order is subsequently elevated to obtain quadratic (second-order) tetrahedral elements, consistent with the meshing strategy validated in the companion elastic study [27]. One exception applies, as noted above: F-RD at $\rho =0.1$ uses linear (first-order) elements for both FE and MSG. As a final safeguard, the periodicity of the generated mesh is verified numerically: for each pair of opposite faces, the in-plane coordinates of corresponding nodes are compared to within a specified tolerance. The same periodic mesh is exported in two formats from a single geometry-and-meshing script: a Gmsh .msh file for the MSG-based analysis with CmbsFE (Section 2.4), and an Abaqus .inp file for the FE-based analysis (Section 2.3). Both homogenization approaches therefore operate on an identical discretization of the RUC, so differences between the two sets of predictions cannot be attributed to differences in geometry or mesh discretization. Figure 2 shows the conforming, periodic two-phase meshes for the four TPMS-based IPC topologies at $\rho =0.3$. The mesh resolution adopted for the FE analyses is justified by the convergence study presented in Section 3.1.

2.2. Constituent Phase Properties

Both constituent phases are modeled as homogeneous, isotropic, linear thermo-elastic conducting solids. Each is fully characterized by five properties: a Young’s modulus E, a Poisson’s ratio $\nu$, a coefficient of thermal expansion $\alpha$, a thermal conductivity $\kappa$, and a (mass-)specific heat capacity c. Both phases are treated as generic engineering materials. They are distinguished only by their elastic-modulus, CTE, and thermal-conductivity contrasts, not as specific, named polymer, metal, or ceramic constituents. The elastic moduli, Poisson’s ratios, and CTEs follow the application example [43], itself based on the gyroid IPC studied by Chawla and Kiran [38]. Chawla and Kiran [38] report an elastic-modulus ratio of 20 between the two phases. On this basis, the COMSOL example [43] assumes representative values of 200 GPa and 10 GPa for the TPMS and matrix phases, respectively, together with a common Poisson’s ratio of $0.30$. The corresponding CTEs are $0.8\times {10}^{-6}$ 1/K (TPMS phase) and $44.0\times {10}^{-6}$ 1/K (matrix phase). Neither [38] nor the COMSOL example [43] reports a thermal conductivity or specific heat capacity, since both are purely elastic/thermo-elastic studies. Representative thermal conductivities of $0.1$ W/(m·K) (TPMS phase) and $0.3$ W/(m·K) (matrix phase) are therefore assigned for the present thermal-conductivity study, together with a common specific heat capacity of 1650 J/(kg·K) for both phases. The specific heat capacity is required by the CmbsFE/MSG solver formulation of Section 2.4.3. However, only steady-state (time-independent) problems are considered here, so its value does not influence any of the effective properties reported in this work. The resulting properties are summarized in Table 1.

The matrix and TPMS phases are assigned different isotropic thermal conductivities, with a contrast ratio ${\kappa }_{\mathrm{matrix}}/{\kappa }_{\mathrm{TPMS}}=3$; the anisotropy of the effective conductivity tensor of the IPC therefore reflects the combined effect of the geometric arrangement of the two phases (topology and relative density) and this conductivity contrast, consistent with standard two-phase conductivity homogenization.

2.3. Finite-Element-Based Homogenization

The FE-based homogenization of the TPMS-based IPC RUCs is performed in Abaqus, using the Abaqus Micromechanics Plugin [30,31,39]. It follows the general RUC homogenization framework validated for single-phase TPMS lattices in [27]. Here, this framework is extended to the coupled thermal-expansion and heat-conduction problems, using the numerical procedures described in [44].

2.3.1. Periodic Boundary Conditions and Macro-Homogeneity

For a two-phase RUC of volume V, let $C_{ijkl}\left(\mathbf{y}\right)$, ${\alpha }_{ij}\left(\mathbf{y}\right)$ and ${\kappa }_{ij}\left(\mathbf{y}\right)$ denote, respectively, the local (position-dependent) stiffness, CTE, and thermal-conductivity tensors at a point $\mathbf{y}$ within the RUC. Each phase is itself homogeneous and isotropic (Section 2.2), so these piecewise-constant fields take the matrix-phase properties of Table 1 when $\mathbf{y}$ lies in the matrix phase, and the TPMS-phase properties when $\mathbf{y}$ lies in the TPMS phase. Throughout this section, a superposed asterisk (e.g., $C_{ijkl}^{*}$) denotes an effective (homogenized) property of the whole IPC. An overbar (e.g., ${\overline{\epsilon }}_{ij}$) denotes a macroscopic (volume-averaged) field. Unmarked symbols (e.g., ${\sigma }_{ij}\left(\mathbf{y}\right)$) denote local fields that vary from point to point within the RUC. The local displacement field $u_i$ at a point $\mathbf{y}$ within the RUC is decomposed into a macroscopic, affine part and a periodic fluctuation ${\chi }_i$ [27]:

$$u_i\left(\mathbf{y}\right)={\overline{\epsilon }}_{ij}{y}_j+{\chi }_i\left(\mathbf{y}\right),$$

(6)

where ${\overline{\epsilon }}_{ij}$ is the prescribed macroscopic strain. Analogously, the local temperature field $\theta \left(\mathbf{y}\right)=T\left(\mathbf{y}\right)-T_0$ (with $T_0$ a stress-free reference temperature) is decomposed as

$$\theta \left(\mathbf{y}\right)=\overline{{\nabla }_j\theta }{y}_j+w\left(\mathbf{y}\right),$$

(7)

where $\overline{{\nabla }_j\theta }$ is the prescribed macroscopic temperature gradient and w is a periodic temperature fluctuation. Periodic boundary conditions require that ${\chi }_i$ and w take equal values at corresponding points on opposite faces of the RUC [27,44], i.e.,

$${\chi }_i\left({\mathbf{y}}^{+}\right)={\chi }_i\left({\mathbf{y}}^{-}\right),w\left({\mathbf{y}}^{+}\right)=w\left({\mathbf{y}}^{-}\right),$$

(8)

together with anti-periodicity of the corresponding tractions and normal heat fluxes on opposite faces. In Abaqus, Eq. (8) is implemented through multi-point constraint equations, linking the degrees of freedom of matching node pairs on opposite RUC faces [27,30,44]. This is enabled by the conforming periodic mesh described in Section 2.1. PBCs are used throughout this study since among the classical boundary-condition choices (kinematic-uniform, static-uniform, and periodic), they are widely recognized as giving the most accurate estimate of the effective properties for a given RUC size. Also, they satisfy the Hill–Mandel macro-homogeneity condition of equivalent strain energy (respectively, thermal dissipation) between the heterogeneous RUC and its homogenized counterpart [27,44].

2.3.2. Effective Elastic Properties

The effective elastic response of the homogenized IPC is described by the generalized Hooke’s law, ${\overline{\sigma }}_{ij}=C_{ijkl}^{*}{\overline{\epsilon }}_{kl}$, where ${\overline{\sigma }}_{ij}$ and ${\overline{\epsilon }}_{ij}$ denote the volume-averaged (macroscopic) stress and strain,

$${\overline{\sigma }}_{ij}=\langle {\sigma }_{ij}\rangle ={\displaystyle \frac{1}{V}}{\int }_V{\sigma }_{ij}\mathrm{d}V,{\overline{\epsilon }}_{ij}=\langle {\epsilon }_{ij}\rangle ={\displaystyle \frac{1}{V}}{\int }_V{\epsilon }_{ij}\mathrm{d}V,$$

(9)

and $C_{ijkl}^{*}$ is the sought effective stiffness tensor of the IPC. Following [27], the full $6\times 6$ effective stiffness matrix is populated column-by-column from six independent linear-elastic analyses. In each analysis, a single component of the macroscopic strain tensor is prescribed to unity through the PBCs of Eq. (8), while all other components are set to zero. This gives three axial load cases (${\overline{\epsilon }}_{11}$, ${\overline{\epsilon }}_{22}$, ${\overline{\epsilon }}_{33}$) and three shear load cases (${\overline{\gamma }}_{12}$, ${\overline{\gamma }}_{13}$, ${\overline{\gamma }}_{23}$). For each load case, the resulting local stress field ${\sigma }_{ij}\left(\mathbf{y}\right)$, obtained from a static analysis of the two-phase RUC, is volume-averaged according to Eq. (9). This yields the corresponding column of $C_{ijkl}^{*}$, expressed in contracted (Voigt) notation as the $6\times 6$ matrix $C_{mn}^{*}$. The standard index pairs $m,n=1,\dots ,6$ map the tensor indices as $11\to 1$, $22\to 2$, $33\to 3$, $23\to 4$, $13\to 5$, $12\to 6$. Only this direct, unit-strain extraction procedure is used in the present study; the least-squares regression alternative of [44], relevant for more general, non-isolating strain states, is not required here. The effective compliance matrix is obtained by inversion, $S_{mn}^{*}={\left[C_{mn}^{*}\right]}^{-1}$. The nine independent engineering constants of the (generally orthotropic) homogenized IPC are then recovered from the components of $S_{mn}^{*}$ [44]: the effective Young’s moduli $E_1^{*},E_2^{*},E_3^{*}$, shear moduli $G_{12}^{*},G_{13}^{*},G_{23}^{*}$, and Poisson’s ratios ${\nu }_{12}^{*},{\nu }_{13}^{*},{\nu }_{23}^{*}$.

2.3.3. Effective Coefficients of Thermal Expansion

The effective CTE tensor is obtained following the numerical procedure of [44], based on the Duhamel–Neumann thermo-elastic constitutive law relating the macroscopic stress to the macroscopic total strain and a uniform temperature change $\Delta T$:

$${\overline{\sigma }}_{ij}=C_{ijkl}^{*}\left[{\overline{\epsilon }}_{kl}-{\alpha }_{kl}^{*}\Delta T\right],$$

(10)

where ${\alpha }_{kl}^{*}$ is the sought effective CTE tensor. A single thermal-expansion analysis is performed. The PBCs of Eq. (8) are applied with a zero prescribed macroscopic strain, ${\overline{\epsilon }}_{kl}=0$, while a uniform temperature change $\Delta T$ is imposed over the entire RUC. Under these conditions, Eq. (10) reduces, in contracted notation, to ${\overline{\sigma }}_m=-C_{nm}^{*}{\alpha }_n^{*}\Delta T$. The effective CTE vector is then recovered from the homogenized thermal stresses and the previously computed effective compliance matrix (Section 2.3.2) as

$${\alpha }_n^{*}=-{\displaystyle \frac{1}{\Delta T}}{S}_{nm}^{*}{\overline{\sigma }}_m.$$

(11)

2.3.4. Effective Thermal Conductivity

The effective thermal conductivity tensor is obtained from a steady-state heat-conduction analysis of the RUC. This is based on Fourier’s law, $q_i=-{\kappa }_{ij}{\nabla }_j\theta$, and the local energy balance ${\nabla }_i\left({\kappa }_{ij}{\nabla }_j\theta \right)=0$ in the absence of internal heat generation [44]. Three thermal load cases are considered. The macroscopic temperature gradient is prescribed through the PBCs of Eq. (8), along each of the three orthogonal directions in turn: ${\overline{{\nabla }_j\theta }}^k={\delta }_{jk}\Theta$ ($k=1,2,3$), where $\Theta$ is the magnitude of the imposed temperature difference across the RUC. For each load case, the resulting local heat flux $q_i\left(\mathbf{y}\right)$ and temperature gradient ${\nabla }_i\theta \left(\mathbf{y}\right)$ fields are volume-averaged,

$${\overline{q}}_i=\langle q_i\rangle ={\displaystyle \frac{1}{V}}{\int }_V{q}_i\mathrm{d}V,\overline{{\nabla }_i\theta }=\langle {\nabla }_i\theta \rangle ={\displaystyle \frac{1}{V}}{\int }_V{\nabla }_i\theta \mathrm{d}V,$$

(12)

The effective conductivity tensor ${\kappa }_{ij}^{*}$ is then obtained from the homogenized Fourier’s law, ${\overline{q}}_i^k=-{\kappa }_{ij}^{*}{\overline{{\nabla }_j\theta }}^k$, applied to the three load cases $k=1,2,3$ simultaneously:

$${\kappa }_{ij}^{*}=-{\overline{Q}}_{ik}{\left[{\overline{\nabla T}}^{-1}\right]}_{jk},{\overline{Q}}_{ik}\equiv {\overline{q}}_i^k,{\overline{\nabla T}}_{jk}\equiv {\overline{{\nabla }_j\theta }}^k.$$

(13)

2.4. Mechanics of Structure Genome-Based Homogenization

The MSG-based homogenization is performed with CmbsFE, the in-house C++ implementation of the MSG framework. CmbsFE is built on GetFEM [45] for finite element assembly, MUMPS [46,47] for direct linear-system solution, and MuParser [48] for symbolic expression parsing. MUMPS is run in serial mode. CmbsFE was previously used for the effective elastic properties of nano-composites and for the effective thermal conductivity of composites with imperfect interfaces [28,29], and was extended to TPMS architectures in the companion elastic study [27]. MSG is a unified, semi-analytical multiscale constitutive theory [24,25], based on the Structure Genome concept: the smallest mathematical building block of a material or structure that contains its complete constitutive information. For a material with three-dimensional heterogeneity, such as the TPMS-based IPCs considered here, the SG coincides geometrically with the RUC used for the FE analysis of Section 2.3. The same periodic mesh (Section 2.1) is used for both approaches [27].

2.4.1. Variational Foundation

MSG constructs the homogenized (effective) model by minimizing, over the SG, the loss of information, expressed, for a linear thermo-elastic conducting material, in terms of strain energy and thermal “energy”, between the original heterogeneous body and the homogenized body [24,26,27]. The kinematic and thermal fields within the SG are decomposed, analogously to Eqs. (6)–(7), as

$$u_i(\mathbf{x},\mathbf{y})={\overline{u}}_i\left(\mathbf{x}\right)+{\chi }_i(\mathbf{x},\mathbf{y}),\theta (\mathbf{x},\mathbf{y})=\overline{\theta }\left(\mathbf{x}\right)+w(\mathbf{x},\mathbf{y}),$$

(14)

Here, $\mathbf{x}$ and $\mathbf{y}$ denote, respectively, the macroscopic (structural-scale) and local (SG-scale) position vectors, distinct from the dimensionless level-set coordinates $x_i$ of Section 2.1; $\mathbf{y}$ coincides with the local RUC position used in Section 2.3. ${\overline{u}}_i=\langle u_i\rangle$ and $\overline{\theta }=\langle \theta \rangle$ are the volume-averaged, homogenized displacement and temperature fields. ${\chi }_i$ and w are unknown, periodic fluctuation functions satisfying $\langle {\chi }_i\rangle =0$ and $\langle w\rangle =0$ [27,49]. The corresponding local strain and temperature-gradient fields are

$${\epsilon }_{ij}={\overline{\epsilon }}_{ij}+{\chi }_{(i,j)},{\nabla }_i\theta =\overline{{\nabla }_i\theta }+w_{,i},$$

(15)

where the comma denotes differentiation with respect to the local (micro) coordinates, and parentheses denote the symmetric part of the gradient.

2.4.2. Effective Elastic Properties

For the purely mechanical problem, the strain-energy functional to be minimized over the SG, subject to the periodicity and zero-average constraints on ${\chi }_i$, is [26,27]

$$\Pi =〈\frac{1}{2}{C}_{ijkl}\left({\overline{\epsilon }}_{ij}+{\chi }_{(i,j)}\right)\left({\overline{\epsilon }}_{kl}+{\chi }_{(k,l)}\right)〉,$$

(16)

whose Euler–Lagrange stationarity condition yields the governing equation for the fluctuation function,

$${\left[C_{ijkl}\left({\overline{\epsilon }}_{kl}+{\chi }_{k,l}\right)\right]}_{,j}=0,$$

(17)

with ${\chi }_i$ periodic on the SG boundary. Equation (17) is solved numerically via finite element discretization of the SG in CmbsFE. The solution is linear in the macroscopic strain, ${\chi }_k=H_k^{mn}{\overline{\epsilon }}_{mn}$, where $H_k^{mn}$ is an influence-function tensor relating the macroscopic strains to the local displacement fluctuations. The effective stiffness tensor is then obtained by volume-averaging the strain-energy density [26,27]:

$$C_{ijkl}^{*}=〈C_{ijkl}+C_{ijmn}{H}_{(m,n)}^{kl}〉.$$

(18)

Here, $H_{(m,n)}^{kl}$ is the symmetric part of the spatial gradient of $H_m^{kl}$ with respect to the local coordinate $y_n$, following the same comma-derivative and symmetrization convention used for the local strain field in Eq. (15). The local strain and stress fields can subsequently be recovered by dehomogenization [26]:

$${\epsilon }_{ij}={\overline{\epsilon }}_{ij}+H_{(i,j)}^{mn}{\overline{\epsilon }}_{mn},{\sigma }_{ij}=\left(C_{ijmn}+C_{ijkl}{H}_{(k,l)}^{mn}\right){\overline{\epsilon }}_{mn}.$$

(19)

2.4.3. Effective Coefficients of Thermal Expansion

For the coupled thermo-elastic problem, the strain-energy functional of Eq. (16) is extended to include the local thermal-stress coupling and thermal-storage terms, following the two-step MSG thermo-elastic homogenization of [49]:

$$\Pi =〈\frac{1}{2}{C}_{ijkl}{\epsilon }_{ij}{\epsilon }_{kl}+b_{ij}{\epsilon }_{ij}\theta +\frac{1}{2}{c}_v{\displaystyle \frac{{\theta }^2}{T_0}}〉,$$

(20)

where $b_{ij}$ is the local thermal-stress coefficient tensor of a constituent phase, related to the local CTE and stiffness by $b_{ij}=-C_{ijkl}{\alpha }_{kl}$. $c_v$ is the local specific heat per unit volume at constant volume: a volumetric quantity, distinct from the mass-specific heat capacity c of Table 1. The two would coincide only if scaled by a mass density, which is not required elsewhere in this study and hence not assigned a value. Discretizing the fluctuation functions ${\chi }_i$ and w over the SG mesh, and minimizing Eq. (20), yields a linear system for the nodal fluctuation unknowns. This system is linear in both ${\overline{\epsilon }}_{ij}$ and $\overline{\theta }$ [49]. Substituting the solution back into Eq. (20) gives the homogenized strain-energy density in the reduced form

$$\Pi =\frac{1}{2}{\overline{\epsilon }}^{\mathrm{T}}{C}^{*}\overline{\epsilon }+{\overline{\epsilon }}^{\mathrm{T}}{b}^{*}\overline{\theta }+\frac{1}{2}{c}_v^{*}{\displaystyle \frac{{\overline{\theta }}^2}{T_0}},$$

(21)

where $C^{*}$ is the effective stiffness matrix of Eq. (18), $b^{*}$ is the effective thermal-stress coefficient vector, and $c_v^{*}$ is the effective specific heat. The effective CTE vector then follows directly as [49]

$${\alpha }^{*}=-{\left(C^{*}\right)}^{-1}{b}^{*}.$$

(22)

2.4.4. Effective Thermal Conductivity

For the heat-conduction problem, MSG minimizes an analogous thermal “energy” functional over the SG [29,50]:

$$\Phi =〈\frac{1}{2}\left(\overline{{\nabla }_i\theta }+w_{,i}\right){\kappa }_{ij}\left(\overline{{\nabla }_j\theta }+w_{,j}\right)〉,$$

(23)

subject to $\langle w\rangle =0$ and periodicity of w on the SG boundary. The temperature fluctuation is discretized over the SG finite element mesh as

$$w=\mathbf{S}\mathbf{V},$$

(24)

where $\mathbf{S}$ is the matrix of standard finite element shape functions defined over the SG mesh, of the same type as used for the FE analysis of Section 2.3. $\mathbf{V}$ is the vector of unknown nodal values of the fluctuation function w. Imposing stationarity of Eq. (23) with respect to $\mathbf{V}$, subject to Eq. (24), yields a linear system of the form $\mathbf{K}\mathbf{V}=-{\mathbf{D}}_{hg}\overline{\nabla \theta }$. Here $\mathbf{K}=\langle {\left(\mathbf{h}\mathbf{S}\right)}^{\mathrm{T}}\mathbf{\kappa }\left(\mathbf{h}\mathbf{S}\right)\rangle$ is the SG conductance matrix. ${\mathbf{D}}_{hg}=\langle {\left(\mathbf{h}\mathbf{S}\right)}^{\mathrm{T}}\mathbf{\kappa }\rangle$ and ${\mathbf{D}}_{gg}=\langle \mathbf{\kappa }\rangle$ are further SG operator matrices, assembled from the local conductivity tensors $\mathbf{\kappa }$ of the two phases and the finite element shape functions of the SG mesh. $\mathbf{h}$ is the gradient (differential) operator matrix relating nodal shape-function values to their spatial derivatives over the SG [50]. The solution of this linear system is linear in the macroscopic temperature gradient, $\mathbf{V}={\mathbf{V}}_0\overline{\nabla \theta }$, obtained from a single factorization of $\mathbf{K}$. Substituting back into Eq. (23) gives the reduced thermal energy $\Phi =\frac{1}{2}{\overline{\nabla \theta }}^{\mathrm{T}}{\kappa }^{*}\overline{\nabla \theta }$, from which the effective thermal conductivity tensor is identified as [50]

$${\kappa }^{*}={\mathbf{V}}_0^{\mathrm{T}}{\mathbf{D}}_{hg}+{\mathbf{D}}_{gg}.$$

(25)

The local temperature gradient and heat-flux fields within the SG can subsequently be recovered by dehomogenization [50]:

$${\nabla }_i\theta =\overline{{\nabla }_i\theta }+{\left(\mathbf{h}\mathbf{S}{\mathbf{V}}_0\right)}_i\overline{{\nabla }_i\theta },q_i=-{\kappa }_{ij}{\nabla }_j\theta .$$

(26)

2.5. Comparison Framework Between the FE and MSG Approaches

Both approaches solve, over the same RUC/SG domain and the same periodic mesh (Section 2.1), governing equations that are formally equivalent for a three-dimensional, periodic heterogeneous medium [26]. The FE displacement- and temperature-based boundary value problems of Section 2.3, and the MSG variational statements of Section 2.4, reduce to the same Euler–Lagrange equations when periodic boundary conditions are used. In the absence of numerical discretization and solver error, the two methods are therefore expected to yield identical effective properties and local fields. They differ, however, in computational implementation. The FE-based approach requires one full system factorization and solve for each of the ten load cases considered in this study: six for the elastic problem, one for the CTE problem, and three for the thermal-conductivity problem. This is because the assembled system differs between mechanical and thermal degrees of freedom, and between the individual load cases under PBCs [26,27]. The MSG-based approach, by contrast, requires only a single factorization of the SG system matrix for each physics (mechanical or thermal), followed by multiple, computationally inexpensive back-substitutions for the different load cases [26,27,50]. This intrinsic computational advantage was already reported for the purely elastic problem in [27]. It is expected to be more pronounced in the present multifunctional study, since three separate sets of effective properties (elastic, CTE, and thermal conductivity) are extracted for every topology and relative density considered.

3. Results and Discussion

This section first establishes the adequacy of the mesh resolution used in the numerical analyses (Section 3.1). The FE and MSG frameworks are then validated against previously published single-phase TPMS results (Section 3.2). Next, the effective elastic properties, coefficients of thermal expansion, and thermal conductivities of the two-phase IPCs are presented property-by-property (Section 3.3, Section 3.4 and Section 3.5), followed by a joint discussion of their multifunctional response (Section 3.6) and a quantitative comparison of the two homogenization approaches (Section 3.7).

3.1. Mesh Convergence Study

Before validating the homogenization frameworks and presenting the two-phase IPC results, a mesh convergence study was performed for the gyroid topology at the lowest and highest relative densities considered, $\rho =0.1$ and $\rho =0.5$. The target element edge length $l_e$ was varied while the level-set sampling resolution used to generate the TPMS surface was kept fixed. For each mesh, the effective axial Young’s modulus $E_1^{*}$, coefficient of thermal expansion ${\alpha }_1^{*}$, and thermal conductivity ${\kappa }_1^{*}$ were evaluated using the FE-based RUC procedure described in Section 2.3.

Table 2 and Table 3 show that the predicted properties converge rapidly at both relative densities. For $\rho =0.1$, reducing the target element edge length from 10 mm to 3 mm changes $E_1^{*}$, ${\alpha }_1^{*}$, and ${\kappa }_1^{*}$ by less than $0.05\%$, $0.04\%$, and $0.01\%$, respectively. The corresponding changes at $\rho =0.5$ are below $0.21\%$, $0.09\%$, and $0.02\%$, respectively.

Based on these results, a target element edge length of $l_e=6$ mm was selected for the analyses reported in this study. For a 100 mm unit cell, this corresponds to approximately $94,300$ elements at $\rho =0.1$ and $88,100$ elements at $\rho =0.5$. The selected mesh agrees with the finest mesh to within $0.14\%$ for all three effective properties at both relative densities, while maintaining a computational cost suitable for the full set of topology–density configurations.

3.2. Validation of the FE and MSG Frameworks

Before presenting the new two-phase IPC results, the FE and MSG homogenization procedures of Section 2 are validated against two previously published reference cases. This uses the same RUC/SG generation, meshing, and solution workflow, applied here to the single-phase (bulk-material) limit of the gyroid and diamond topologies, following the validation strategy adopted in the companion elastic study [27]. The two reference cases are the titanium gyroid lattice of Refai et al. [16] ($\rho =0.1$, bulk properties $E_s=110$ GPa, $G_s=41$ GPa, ${\nu }_s=0.34$), and the diamond lattice of Zhang et al. [18] ($\rho =0.3$, bulk properties $E_s=190$ GPa, $G_s=73$ GPa, ${\nu }_s=0.30$); the latter includes both a numerical (FEM) and two experimental measurements (EXP). Table 4 compares the effective properties obtained in this study with these literature values.

Both the FE and MSG predictions agree with the gyroid reference of Refai et al. [16] to within $1.9\%$, and with the diamond FEM reference of Zhang et al. [18] to within $0.1\%$. The largest deviation, $4.51\%$, is against Zhang et al.’s first experimental measurement; this is consistent with the level of numerical-experimental scatter reported by the original authors. These deviations are consistent with, and in fact tighter than, those reported for the same validation cases in the companion elastic study [27]. This confirms that the FE and MSG frameworks employed here reproduce established literature results, before being applied to the new two-phase IPC configurations discussed below.

3.3. Effective Elastic Properties

Figure 3 presents the effective axial Young’s modulus $E_1^{*}$, shear modulus $G_{12}^{*}$, and Poisson’s ratio ${\nu }_{12}^{*}$ of the four TPMS-based IPC topologies as a function of relative density $\rho$, predicted by the FE and MSG approaches.

Table 5 reports the numerical values underlying Figure 3, for all four topologies and the full relative-density range investigated.

For all four topologies, $E_1^{*}$ and $G_{12}^{*}$ increase monotonically with relative density. At $\rho =0.1$, $E_1^{*}$ ranges from about $17.8$ to $20.1$ GPa across the four topologies. At $\rho =0.5$, it ranges from about 62 to 76 GPa. This increase is consistent with the growing volume fraction of the stiffer TPMS phase ($E=200$ GPa), relative to the compliant matrix phase ($E=10$ GPa, Table 1). The F-RD topology is the stiffest at every relative density investigated. At $\rho =0.1$, F-RD leads by about $13\%$ over the softest topologies, gyroid and PMY. This gap widens steadily with relative density, reaching about $22\%$ at $\rho =0.5$. The gyroid topology is consistently the most compliant. Diamond and PMY fall between gyroid and F-RD, with their relative ranking varying with $\rho$. The effective Poisson’s ratio decreases with increasing relative density for all four topologies, from ${\nu }_{12}^{*}\approx 0.28$–$0.30$ at $\rho =0.1$ towards ${\nu }_{12}^{*}\approx 0.26$–$0.29$ at $\rho =0.5$. This reflects the progressively smaller influence of the compliant, higher-Poisson’s-ratio matrix phase, as the stiffer TPMS network comes to dominate the load path. The same qualitative trend was reported for single-phase TPMS lattices in the companion study [27].

Only the axial components $E_1^{*}$, $G_{12}^{*}$ and ${\nu }_{12}^{*}$ are shown in Figure 3. Across all four topologies, all five relative densities, and both homogenization methods, the remaining tensor components differ from their axial counterparts by at most $1.3$ (Young’s modulus), $1.5$ (shear modulus), and $1.4\%$ (Poisson’s ratio). This confirms that the effective elastic response of all four TPMS-based IPCs remains quasi-isotropic (more precisely, cubic-symmetric, consistent with the cubic point-group symmetry of the underlying TPMS unit cells) over the density range investigated, as already observed for the corresponding single-phase lattices in [27].

3.4. Effective Coefficients of Thermal Expansion

Figure 4 presents the effective axial CTE ${\alpha }_1^{*}$ as a function of relative density.

Table 6 reports the numerical values underlying Figure 4, for all four topologies and the full relative-density range investigated.

The effective CTE decreases monotonically with increasing relative density for all four topologies. At $\rho =0.1$, ${\alpha }_1^{*}$ ranges from about $23.8$ to $24.8\times {10}^{-6}$/K across the four topologies. At $\rho =0.5$, it ranges from about $5.7$ to $6.3\times {10}^{-6}$/K. This decrease occurs because the low-CTE TPMS phase ($\alpha =0.8\times {10}^{-6}$/K) increasingly constrains the thermal expansion of the high-CTE matrix phase ($\alpha =44\times {10}^{-6}$/K). The reduction is markedly stronger than a rule-of-mixtures estimate would suggest. For example, at $\rho =0.1$, a Voigt-type volume average gives $\alpha \approx 0.1\times 0.8+0.9\times 44\approx 39.7\times {10}^{-6}$/K, roughly $1.6$–$1.7\times$ the values in Figure 4. This reflects the elastic-constraint (Turner-type) mechanism inherent to the Duhamel–Neumann formulation of Section 2.3.3: the stiffer TPMS network suppresses the free thermal expansion of the compliant matrix, well beyond what a volume-weighted average would predict. Unlike the elastic properties, no single topology has the lowest CTE across the whole density range. Diamond has the lowest CTE at $\rho =0.1$ and $\rho =0.2$. PMY takes over at $\rho =0.3$ and $\rho =0.4$. F-RD has the lowest CTE only at $\rho =0.5$, where it is also the stiffest topology (Section 3.3). This shifting ranking shows that the CTE-stiffness relationship is not a simple, universal trend across topologies. It depends on the specific network geometry at each relative density.

3.5. Effective Thermal Conductivity

Figure 5 presents the effective axial thermal conductivity ${\kappa }_1^{*}$ as a function of relative density.

Table 7 reports the numerical values underlying Figure 5, for all four topologies and the full relative-density range investigated.

The effective conductivity decreases monotonically with relative density for all four topologies. At $\rho =0.1$, ${\kappa }_1^{*}$ ranges from about $0.2696$ to $0.2701$ W/(m·K) across the four topologies: a narrow spread of less than $0.2\%$. At $\rho =0.5$, the spread widens to about $0.176--0.180$ W/(m·K). This decrease occurs as the low-conductivity TPMS phase ($\kappa =0.1$ W/(m·K)) progressively displaces the higher-conductivity matrix phase ($\kappa =0.3$ W/(m·K)). At $\rho =0.1$, the predicted values for all four topologies fall within the Hashin–Shtrikman bounds for a two-phase isotropic conductor of the same constituent conductivities and volume fractions (${\kappa }_{\mathrm{HS}}^{-}=0.269$, ${\kappa }_{\mathrm{HS}}^{+}=0.275$ W/(m·K)). This provides an additional, physics-based consistency check on the conductivity homogenization procedure, beyond the mesh convergence study of Table 2. As with the CTE, no single topology has the lowest conductivity across the whole density range. Diamond is lowest at $\rho =0.1$; PMY is lowest at $\rho =0.2--0.4$; F-RD is lowest only at $\rho =0.5$, mirroring the CTE ranking of Section 3.4.

3.6. Cross-Property Discussion

Figure 3, Figure 4 and Figure 5 reveal a consistent multifunctional trade-off with density. Increasing relative density raises stiffness. At the same time, it lowers both CTE and thermal conductivity. All three effective properties are governed by the same underlying competition: the stiff, low-CTE, low-conductivity TPMS phase against the compliant, high-CTE, higher-conductivity matrix phase. The four topologies respond to this trade-off differently. Stiffness ranking is stable with density: gyroid is the most compliant topology, and F-RD is the stiffest, at every relative density investigated. CTE and thermal conductivity rankings are not stable. Diamond has the lowest CTE and conductivity at low relative density. PMY takes over at intermediate density. F-RD only becomes the lowest-CTE, lowest-conductivity topology at $\rho =0.5$, where it is also the stiffest. At that density, F-RD combines the best stiffness, CTE, and conductivity of the four topologies simultaneously. At lower densities, the stiffest topology and the best thermal-transport topology are not the same. From a design standpoint, this has a direct implication: at low relative density, a designer should choose diamond for the lowest thermal expansion and conductivity, or F-RD for the highest stiffness; a single topology cannot deliver both simultaneously. At high relative density, F-RD is the clear choice: it dominates all three properties at once, at low mass. The gyroid topology, with its low stiffness and comparatively strong thermal transport at high density, remains preferable where higher thermal conductivity is the priority, such as in heat-exchanger applications.

3.7. Comparison Between the FE and MSG Approaches

Figure 6 presents parity plots comparing the FE- and MSG-predicted effective Young’s modulus, CTE, and thermal conductivity for all four topologies and all five relative densities considered.

All data points lie close to the $y=x$ line across the full range of topologies, relative densities, and properties, confirming the close agreement between the two homogenization approaches anticipated on theoretical grounds in Section 2.5. Table 8 summarizes the relative error between the MSG and FE predictions, $|\mathrm{MSG}-\mathrm{FE}|/\mathrm{FE}\times 100\%$, for the representative axial elastic modulus, shear modulus, Poisson’s ratio, CTE, and thermal conductivity. Values are maximized over the five relative densities for each topology.

The maximum relative error between the two approaches is below $2.6\%$ for the elastic properties, $1.7\%$ for the CTE, and $0.4\%$ for the thermal conductivity, across all four topologies and the full relative-density range investigated. This level of agreement is consistent with that reported for the purely elastic, single-phase TPMS lattices in the companion study [27], and with the theoretical equivalence of the two formulations established in Section 2.5. The largest discrepancies are observed for the elastic properties of the gyroid topology and for the CTE of the PMY topology. These differences are mesh-independent (Table 2) and consistent with the validation-level agreement of Section 3.2. They are attributed to the differing numerical discretizations underlying the FE (Abaqus, six-load-case direct extraction) and MSG (CmbsFE, single-factorization variational) solution procedures, rather than to a methodological inconsistency between the two frameworks. The F-RD topology at $\rho =0.1$ shows the closest agreement of the entire study (Table 8). This is the only configuration in which both FE and MSG use the same linear element order, as noted in Section 2.1. Its tight agreement supports the view that most of the small remaining differences, elsewhere in the study, stem from the different solution procedures, rather than from the geometry or mesh itself.

4. Conclusions

This study has compared FE-based and MSG-based homogenization for predicting the effective elastic properties, coefficients of thermal expansion (CTE), and thermal conductivities of TPMS-based interpenetrating phase composites (IPCs). Four topologies (gyroid, diamond, PMY, and F-Rhombic Dodecahedron (F-RD)) were investigated over relative densities from $0.1$ to $0.5$. Conforming-mesh RUC models were analyzed identically by both approaches.

Both frameworks were validated against published literature results for single-phase TPMS lattices, with predictions agreeing to within $1.9\%$. This confirms that the homogenization procedures developed here reproduce established results before being extended to the new two-phase IPC configurations. All four IPC topologies exhibit a consistent multifunctional trade-off with density. Increasing relative density raises stiffness. It also reduces both CTE and thermal conductivity. All three properties are governed by the same competition between the stiff, low-CTE, low-conductivity TPMS phase and the compliant, high-CTE, higher-conductivity matrix phase. The F-RD topology is the stiffest of the four at every relative density studied. It combines the lowest CTE and lowest conductivity only at the highest relative density, $\rho =0.5$, where it is the clear multifunctional choice. At lower relative densities, diamond and PMY offer lower CTE and conductivity than F-RD, even though F-RD remains stiffer. Topology selection therefore depends on both the target relative density and the property being prioritized. MSG predictions agree with FE-RUC to within $2.6\%$ for elastic properties, $1.7\%$ for CTE, and $0.4\%$ for thermal conductivity, across all topologies and densities. MSG also has an intrinsic computational advantage: a single system factorization per physics, rather than one per load case. Together, this level of agreement and this efficiency support the use of MSG as a validated, lower-cost substitute for FE-RUC homogenization in multifunctional design studies of TPMS-based IPCs.

Future work could extend this framework to additional TPMS topologies and volume-fraction combinations, incorporate nonlinear or temperature-dependent constituent behavior, and experimentally validate the predicted thermo-elastic and transport properties of the two-phase IPC configurations considered here.