Efficient design optimization of a miniaturized thermoelectric generator for electrically active implants based on parametric model order reduction

This research focuses on the design of a miniaturized thermoelectric generator (TEG) for electrically active implants. Its design optimization is performed using the finite element method. A simplified TEG model is obtained by replacing the thermocouple array with a single representative thermopile, which considers the number and fill factor of the thermocouples as parameters. Instead of rebuilding the geometry of a detailed model with multiple thermocouples, the simplified model adapts the material properties of its representative thermopile, facilitating design optimization. We extend the model by integrating the simplified TEG together with a housing inside a human tissue model for thermoelectric analysis. For computation efficiency and applicability of model order reduction (MOR), a thermal model is derived from the thermoelectric one, with the Peltier effect being considered through an effective thermal conductivity. Through parametric MOR, two parametric reduced‐order models are generated from the full‐scale thermoelectric and thermal model, respectively. Furthermore, we demonstrate the design optimization of TEG both in full‐scale and reduced‐order model for maximal power output and sufficient voltage output.

149 million by 2050. 1 This trend will impact the healthcare system and the availability of specialist healthcare personnel and resources. Hence, implantable medical devices (IMDs) will have enormous potential to deal with this issue as they can be utilized not only for disease treatment but also for real-time health monitoring 2 to improve the quality of medical care and relieve the pressure on the health care system. Nevertheless, the service life of IMDs is limited by their integrated batteries. Specifically, most pacemakers operating with lithium batteries have a longevity of about 8 years. 3 When the energy is depleted, patients have to undergo surgery for battery or implant replacement, suggesting higher risk and longer recovery time for older people.
Fortunately, due to low power requirements of common IMDs, which ranges from tens of μW to several mW, 4 the energy harvesting technologies 5 become an attractive alternative to achieve self-sufficient power supply or prolong the life-span of the batteries by harvesting energy from surroundings or human body. Considerable progress has hitherto been made in harvesting technologies for self-powered IMDs. 6,7,8,9,10,11 Researchers have tried to excavate as many energy sources in the human body as possible to take advantage of them in various implantable energy harvesters. Such energy harvesters include biofuel cells 12,13,14 extracting electrochemical energy, piezoelectric 15,16,17 and triboelectric 18,19 energy harvesters that harness mechanical energy, as well as thermoelectric generators (TEGs) 20,21 that transform thermal energy of the body into electricity through the Seebeck effect. TEGs are particularly suitable, as they do not rely on movable components.
TEGs play an important role in many fields such as automotive, power plants and aerospace. In most cases, the thermal resistance of the environment between the heat source and heat sink can be reduced to a certain extent, which is much smaller than the thermal resistance of TEGs. Hence, for the sake of higher power output, researchers concentrate more on the design of the heat sink than the TEG itself. But when we turn to the design of the TEG for electrically active implants, in which case the thermal resistance of the human tissue is comparable to the TEG and not adaptable, the thermal matching between the TEG and its implantation environment is significant. In such a case, the power output largely depends on the array configuration of thermocouple legs in the TEG. Thinner legs have lower electrical resistance, indicating less power loss. But at the same time, the smaller temperature difference across the legs can be maintained due to their lower thermal resistance, leading to lower power output. Therefore, an optimized array configuration exists to balance these two opposite effects and provide maximum power.
Numerical simulations techniques can enable time and cost savings during the design process. Yang et al. 22 proposed a three-layer tissue model with embedded TEG and performed a numerical simulation. However, such a finite element method (FEM) results in high-dimensional ordinary differential equations (ODEs), causing the high computational cost. To speed up the simulations, Jadhav et al. 23 applied model order reduction (MOR) method 24 on the finite element model of TEG to obtain a highly accurate but low-dimensional model. Afterwards, to investigate the impact of the environmental conditions on TEG and to avoid the need to regenerate the new reducedorder model with each change of the parameters, Yuan et al. 25,26 constructed a TEG reduced model preserving the parameters in convection boundary conditions through the multi-moment matching parametric model order reduction (pMOR) method. 27,28,29,30 Furthermore, Yuan et al. 31,32 extracted the fill factor of thermopiles as an optimization parameter by parameterizing the equivalent thermal conductivity of a thermal dummy model. It was found that the fill factor is a crucial parameter for the performance of the implanted TEG. But the Peltier effect was not considered in the model, which affects the overall thermal field distribution and cause significant errors in simulation results.
In this contribution, a simplified TEG model was used for the convenience of parameterization of the array configuration. Instead of changing the geometry parameters of the thermocouple array, its fill factor and number of thermocouple legs were used to specify the material properties of a representative thermopile. This model was validated through comparison against a detailed TEG model. The Seebeck coefficient and electrical resistivity of the thermoelectric material were obtained from experimental results of a commercially available TEG. An extended model integrates a packaged TEG inside human tissue. Furthermore, a solely thermal model was developed using effective thermal conductivity to approximate the Peltier heat in TEG. Moreover, reduced TEG models were obtained through the pMOR method for efficient parameterization studies and system-level co-simulation. This paper is organized as follows. In Section 2, the simplified and detailed TEG models are presented. Furthermore, the simplified thermoelectric and thermal models are introduced and compared using parameter sweep on fill factor and load resistance. Section 3 focuses on reducing the full-scale model introduced in Section 2 using pMOR for efficient simulations. The emphasis of this section is placed on the design optimization of the implanted TEG. Section 4 concludes the paper and provides an outlook on future work.

| FULL ORDER MODEL
In this section, three-dimensional thermoelectric and thermal models were simulated with ANSYS (Version 2021 R1, Ansys Inc.). The results of stationary analysis are compared with experimental data.

| Initial design
The initial design of this work is depicted in Figure 1. The three-layer tissue has a cross section of 100 mm Â 100 mm, consisting of 2 mm-thickness skin layer, 8 mm-thickness fat layer and 35 mm-thickness muscle layer. The blood perfusion Q b and metabolic heat generation Q m are considered in skin, fat and muscle layers. The temperature on the bottom surface of the muscle layer is constant at 37 C to simulate the human body core. The effects of convection, radiation and evaporation (sweating) are applied on the surface of the skin layer as heat losses. The specific bioheat modeling will be elaborated in Section 2.2.
A TEG is implanted under the skin layer, where the maximum power output is expected. 22,23,33 To fully utilize the height of fat layer, the implanted TEG has a size of 20 mm Â 20 mm Â 8 mm, whose surfaces directly contact the skin and muscle layers. This TEG is comprised of two 2 mm-height ceramic plates and several 4 mm-height thermocouples. The array configuration of thermopiles will be optimized in section 3.3.
To assure secure implantation and reliable operation, the encapsulation of implants is indispensable. Therefore, a housing made of polyetheretherketone (PEEK) was designed, which is a bio-compatible material and has high thermal resistivity. The housing can also contain electronics such as a voltage converter, signal processing circuitry as well as a capacitor or battery for energy storage. Due to its high thermal isolation the housing also raises the temperature gradient across the implanted TEG and hence boosts the power output. The housing is 60 mm wide and 8 mm high, whose four edges were formed as a cambered surface with a radius of 6.8 mm to prevent damage to the human body after implantation. The thickness of the housing wall was set as 1 mm.

| Bioheat human tissue model
To provide realistic human thermal environments, the bioheat human tissue model illustrated in Section 2.1 was developed based on the research work of Yuan et al. 26 Tissue properties from the IT'IS database 34 were applied. Pennes' bioheat equation, 35 which accounts for heat conduction, metabolic heat generation, and perfusion effects, was applied to describe the time transient heat transfer within human tissue as follows: where T is the temperature distribution in human tissue. ρ, c and κ are the density, specific heat and thermal conductivity of different tissues, respectively. Q b and Q m are the rates of blood perfusion and metabolic heat generation, which were applied in skin, fat and muscle layers. ρ b , c b and T a are the density, specific heat and temperature of arterial blood. ω is the perfusion rate in different tissue layers. The temperature on the bottom surface of the muscle layer is fixed at 37 C. The effects of convection, radiation and evaporation (sweating) are modeled as follows: where q conv is the convective heat flux and q rad and q eva are the radiation, and evaporation heat fluxes applied at the skin surface. T amb is the temperature of the ambient air and h c is the heat transfer coefficient. The thermal process of radiation is described by the Stefan-Boltzmann equation with the Stefan-Boltzmann constant σ and emissivity of human skin ε. Moreover, the evaporation effect at the skin surface is described by q eva , where h e is the evaporation heat transfer coefficient, P skin is the saturated vapor pressure at the skin surface, P sa is the saturated vapor pressure and ϕ is the relative humidity. According to the Lewis relation, 36 which shows that the evaporation coefficient could be presented in terms of convection coefficient and the Antoine's equation, 37 the evaporation heat loss q eva in Equation (2) could be expressed with skin wetness w skin as:

| Detailed TEG model
A detailed TEG model was generated based on a commercially available TEG, as illustrated in Figure 2. In this device, 31 p-type and 31 n-type cubic bismuth telluride thermopiles (1.4 mm Â 1.4 mm Â 1.4 mm) are connected in series through copper interconnectors (3.8 mm Â 1.4 mm Â 0.25 mm). Two ceramic plates made of aluminum oxide (20 mm Â 20 mm Â 0.7 mm) act as substrates and thermal contacts. Its finite element model comprises 77,844 degrees of freedom (DoF) and contains only solid elements. For simulation efficiency and applicability of MOR, all material properties were defined as temperature-independent. This simplification is applicable in our scenario as the temperature variation across the thermoelectric material is around 4 K; a temperature range which causes a negligible change in properties. Furthermore, the thermoelectric properties of bismuth telluride strongly depend on its doping material and doping concentration, which varies with different suppliers. Our model implements a Seebeck coefficient and electrical resistivity, which were experimentally obtained from the open-circuit voltage and internal resistance of the commercially available TEG (see Table A1). Our model also implements a lumped resistor as a resistive load. Beyond heat conduction, the analysis considers the Seebeck and Peltier effects; the Thomson effect has not been implemented as the Seebeck coefficient is considered to be temperature-independent. Generally, commercially available TEGs feature low-height and high-density thermopile arrays, which leads to a low thermal resistance. This prevents the device from maintaining a sufficiently high temperature difference and thus harvesting energy inside the human body, which motivates the design optimization of implanted TEGs to provide maximum electrical power.

| Simplified TEG model
In this section, a simplified TEG model is introduced to avert frequent changes of geometry in the detailed TEG model during parametric analysis. It helps to facilitate the parameterization of the array configuration and thermopile cross section.
F I G U R E 1 Model illustration of the initial design, where a TEG is encapsulated in a housing, which also integrates electronics and energy storage. The packaged TEG is embedded into the fat layer of a simplified human tissue geometry model. The temperature at the bottom surface of the muscle layer is fixed at 37 C to simulate the human body core; the blood perfusion Q b and metabolic heat generation Q m are applied in the skin, fat and muscle layers; the convection, radiation and evaporation are applied on the top surface of the skin layer

| Thermoelectric TEG model
Owing to the high electrical and thermal conductivities compared to bismuth telluride, the effects of copper interconnectors in the detailed TEG model are negligible. Hence, the copper layers were omitted in the simplified TEG model. The veritable thermocouple legs were modeled as a single representative thermopile with the same height as the original thermopiles. The simplified TEG model for thermoelectric analysis contains 63,656 DoF, which is 14,188 DoF less than the detailed model. The equivalent material properties of the representative thermopile including equivalent Seebeck coefficient α 0 , electrical resistivity ρ 0 and thermal conductivity κ 0 were parameterized according to the number and fill factor of thermocouple legs (see Figure 3).
According to the Seebeck effect, the electromotive force (EMF) of TEG can be written as: where n is the number of thermocouple legs, α p and α n is Seebeck coefficient of p-and n-type thermoelectric material respectively, ΔT is the temperature difference between the hot and cold side of the thermocouple legs. Then the equivalent Seebeck coefficient of the representative thermopile α 0 can be expressed as: The thermocouple legs are thermally connected in parallel and electrically in-series, as illustrated in Figure 4. The cross-sectional area of single thermocouple leg a is equal to where η is the fill factor indicating the ratio of the cross-sectional area of total thermocouple legs to ceramic plate, and A is the cross-sectional area of ceramic plate. The thermal conductance and electrical resistance of the total thermocouple legs are where κ p and κ n is thermal conductivity of p-and n-type thermoelectric material, respectively; and ρ p and ρ n denote the electrical resistivity of p-and n-type thermoelectric material, respectively. Then the equivalent thermal conductivity and electrical resistivity of the representative thermopile can be expressed as: The change of array configuration in the detailed TEG model during parameter sweep was reflected in the equivalent material properties of the representative thermopile in the simplified TEG model. The simplified TEG model (thermoelectric TEG model) was verified by comparing it against the detailed TEG model. Neglecting the copper interconnectors, the dimension of the simplified model was kept the same as in the detailed model. The equivalent material properties of the representative thermopile were calculated according to fill factor and number of thermocouple legs in the detailed model, which is 30.4% and 62, respectively. For validation purposes the temperature at the bottom surface of both TEG models was fixed as 34 C; an outwards heat flow of 1 W was applied on the surface of the top ceramic plate and a thermoelectric analysis was performed. The simulation results of the two models match well, where the power output at electrical load matching exhibits an average error of 0.8% (See Figure 5).

| Thermal TEG model
The thermoelectric TEG model involves nonlinearity, in which case MOR is not applicable. Hence, a solely thermal TEG model was developed, which implements the thermoelectric effects by means of an effective thermal conductivity.
The rate of heat absorbed at the hot side of the representative thermopile with temperature T h is given by: The rate of heat liberated at the cold side of the representative thermopile with temperature T c is expressed as: Approximation of the detailed TEG model by the simplified TEG model; the copper interconnectors are neglected, the  veritable thermocouple legs are replaced by a representative thermopile where I is current flowing through the TEG, R is internal resistance of the representative thermopile, K refers to thermal conductance of the thermopile. The first term of both equations describes the heat conducted through the representative thermopile, the second term is the Peltier heat at the hot and cold side respectively, the Joule heat in the third term is shared equally at the both sides. Through the Equations (4), (7) and (8), the rates of heat at the both sides are written as: where Q K denotes the rate of heat conducted through the representative thermopile, Z 0 is the equivalent figure of merit of the thermopile, m is the resistance coefficient between load resistance R L and internal resistance R of the TEG: Here, the temperature difference presented across the implanted TEG is only a few Kelvins; compared to the absolute temperature at hot and cold sides (around 310 K). Thus, the Joule heat in Equations 13 and 14 can be neglected. In addition, the Peltier heat at both sides in the thermal model is considered as conductive passing through the representative thermopile (see Figure 6). The rate of effective conductive heat _ Q h 0 and _ Q c 0 is formed as where T avg is the average temperature at the hot and cold side of the representative thermopile. The effective thermal conductivity of the representative thermopile κ e is expressed as F I G U R E 5 Power output with electrical load matching in the simplified and detailed TEG models The effective thermal conductivity κ e and average temperature T avg depend on each other. Therefore, the model was iteratively solved, starting with κ e evaluated at T avg = 300 K. During the sequence κ e is updated in the model until convergence is obtained.
As the thermal TEG model does not solve for electric potential it has only 30,596 DoF, which is only a half of the thermoelectric TEG model. In spite of that, the thermal TEG model was proven as accurate as the thermoelectric TEG model. The verification was achieved through the comparison in simulation results (see Figure 7). The average error of power output between the thermoelectric and thermal TEG model is 0.6%. The Peltier effect leads to higher effective thermal conductivity and smaller temperature difference across the representative thermopile. It is mandatory to include the effect of Peltier heat in the thermal model, as it decreases the power output by almost one third.

| Simulation results
After integrating the TEG together with a housing inside the human tissue, as depicted in Section 2.1, a thermoelectric and a thermal model were developed based on the thermoelectric and thermal TEG models, respectively (see Figure 8). The thermoelectric model and thermal model has DoF of 92,427 and 84,093, respectively. The ambient temperature was fixed as 21 C. The overall heat transfer coefficient h in the convection boundary condition was calculated as 8.58 W/m 2 /K. The inner space of the housing has no thermal conductivity. Except for the simplified TEG, no F I G U R E 6 Approximation of the Peltier heat in the thermal TEG model using effective thermal conductivity; the Peltier heat at both sides is modeled such that the average value as a part of the effective conductive heat F I G U R E 7 Power output with load matching in the thermoelectric and thermal TEG model components were integrated inside the housing. According to a mesh study, the maximum element size was chosen as 4 mm, which is adequate for calculation precision. The material properties in the models are given in Table A1.
In simulation results of the thermoelectric model (see Figure 9), the temperature difference across the representative thermopile varies with the fill factor and load resistance. An increasing fill factor of the representative thermopile leads to a decrease of its thermal resistance, resulting in a drop of temperature difference. With the increase of load resistance, the effective thermal conductivity of the representative thermopile decreases, causing a higher temperature difference. The relative error between the Peltier heat at the hot and cold side has a similar trend as the temperature difference across the representative thermopile. Its maximum relative error is 2.4%, which we consider as acceptable for the approximation of the Peltier heat in the thermal model. Figure 10 exhibits the thermal field distribution of the thermoelectric model with an optimized TEG. The design optimization will be presented in Section 3.3. A significant temperature difference can be observed across the housing due to the good thermal insulating effect of the housing and cavity inside. It improves the performance of the TEG by increasing its temperature difference, which is observed as 4.3 C.
With the parameter sweep on the fill factor and load resistance in Table A2, the parametric power output results in the thermal model are consistent with the thermoelectric model as shown in Figure 11. An optimum value exists for the fill factor, which yields maximum power output with electrical load matching. With a higher fill factor, the TEG has a smaller internal resistance, which is reflected in the decrease of the optimum load. The relative error in power output between the thermoelectric and thermal model ranges from À2.4% to 8.8%, the average relative error with all parametric points is calculated as 1.5%.

| Experimental validation
The commercially available TEG was experimentally characterized in an experimental setup as presented in Figure 12. The setup is composed of a measurement environment, a Peltier controller, and a display. In the measurement environment, the TEG under test is clamped between two aluminum blocks with integrated temperature sensors. The Peltier modules are attached to heat sinks and fans. A Peltier controller allows to adjust the temperature of the hot and cold side independently. The display indicates the temperature readings of integrated sensors.
The validation of the models was achieved through a comparison between the experimental data of a commercially available TEG and simulation results of its detailed TEG model as introduced in Section 2.3. The hot and cold side temperature of the TEG was fixed to 34 C and 33 C, respectively, which is consistent with physiological conditions. A variable resistive load was attached to the TEG for load matching. The simulation results match with the experimental data, presenting the average error as 2.8% in voltage output and 5.7% in power output (see Figure 13). When load resistance is below 0.4 Ω, an obvious error in power output can be observed, which we attribute to the limited accuracy of the resistor.

| PARAMETRIZED REDUCED ORDER MODEL
In this section, the multivariate moment-matching pMOR method, which has been demonstrated to be successful in linear models, 38,39,25,26,31,32 was applied to generate highly accurate but compact surrogates with preservation of the main features of the original large-scale system to speed up the simulations of the thermoelectric and thermal models introduced in Section 2. In the pROMs of the thermoelectric and thermal models, the effective thermal conductivity of the representative thermopile, the heat transfer coefficient and the ambient temperature in convection boundary conditions were set as the parameters, which can be changed in the reduced model.

| pMOR for thermal model
In the full-scale thermal model in ANSYS, the heat generation rates provided by perfusion and tissue's metabolism were applied as the internal heat sources in different tissue layers. Meanwhile, on the surface of the skin layer, the convection, radiation and evaporation effects were considered to be the external heat transfer effects to balance the heat generated inside. Based on the research of Yuan et al. 26 the temperature-dependent perfusion, radiation, and evaporation effects were transformed and a linear thermal tissue model was obtained for the convenience of pMOR. The temperature dependent effects in Equations (1) and (2) can be modeled as follows: where the blood perfusion heat generation Q b is transformed as a convective heat transfer boundary condition q b . Its negative heat transfer coefficient yields the provision of heat to the tissue. The heat radiation effect is linearized by defining a heat transfer coefficient h rad . Therefore, the overall heat transfer coefficient is defined as: Finally, the evaporation effect is approximated by applying the average heat flux boundary condition at the skin surface, where s is the total number of nodes at the skin surface and q i is the heat flux through each surface node. The finite element based spatial discretization of Equations (1) with boundary conditions (2), that is, (20), leads to a system of ODEs. In addition, by setting the effective thermal conductivity of the representative thermopile, the heat transfer coefficient, and the ambient temperature in the convection boundary condition as the parameters, the governing equations of the full-scale parametric thermal model can be written as follows: where E∈ℝ NÂN is the heat capacity matrix. A 0 ,A 1 , A 2 ∈ℝ NÂN are the parameter-independent heat conductivity matrices. N = 84,093 DoF is the full-scale dimension of the thermal system. κ e is the effective thermal conductivity and h is the overall heat transfer coefficient in the convection boundary condition as defined in Equation (21). T is the unknown temperature state vector. B∈ℝ NÂm is the input distribution matrix. u∈ℝ m is the input vector. Q m is the constant metabolic heat generation rates in different tissue layers. T amb is the ambient temperature. T a is the blood temperature assumed as constant at 37 C. q eva is the average heat flux input of the evaporation effect at the skin surface. C∈ℝ pÂN is the user-defined output matrix. m, p are the number of inputs and user-defined outputs. The basic idea of the pMOR method is to find a parameter-independent projection subspace V ∈ℝ NÂr with r ( N. Then Equation (22) allows easy parameter preservation for projection-based MOR. The reduced-order model of Equation (22) is obtained as: where E r , A r κ e , h ð Þ∈ℝ rÂr are the reduced heat capacity and conductivity matrices. B r ∈ℝ rÂm and C r ∈ℝ pÂr are the reduced input and output matrices. The original large-scale temperature state vector T∈ℝ N is projected onto the low dimensional subspace: where x∈ℝ r is a reduced temperature state vector. In this contribution, static thermal simulations were performed to obtain the temperature difference across the representative thermopile. Therefore, as proved in the research work of Yuan et al. 32 the projection subspace V can be constructed by merging the only two orthonormal bases of the disjoint local Krylov subspaces with respect to parameter κ e and h: where κ 0 , and h 0 are the fixed expansion points with respect to each parameter. The disjoint local Krylov subspaces in Equation (25) are generated for one parameter while keeping the other parameter as constant. V κ e ∈ℝ NÂr 1 , V h ∈ℝ NÂr 2 and r ¼ r 1 þ r 2 . In this work, the pROM was generated using the software "Model Reduction inside ANSYS" 40 with the dimension of the pROM r = 10 DoF.
Based on the pROM, the electrical load circuit was constructed as illustrated in Figure 14. In the reduced model, the averaged temperature results from the top and bottom faces of the representative thermopile was set as the outputs. The EMF of TEG can be calculated according to Equation (4) and used as the voltage source in the electrical load circuit. Furthermore, according to Equations (9), (16), (17), and (19), the average value of the temperature outputs from the reduced model was implemented as the feedback for the calculation of the effective thermal conductivity: where κ ini is the original thermal conductivity of the bismuth telluride and κ peltier is the additional effective thermal conductivity caused by the Peltier heat effect when the current flows through the TEG. R is the fill factor dependent internal resistance of TEG as listed in Equation (8) and R L is the load resistance in the circuit. The parametric simulations of the electrical load circuit were performed by choosing different parameter values of fill factor (η) and load resistance (R L ) as shown in Table A2. In order to study the impact of these two parameters on TEG, all the other inputs were set as constant values. In Equations (22) and (23), q eva = À12 W/m 2 , T a = 37 C, h = 8.58 W/m 2 /K and T amb = 21 C were applied. Q m = 498.5, 279.8, 841.6 W/m 3 were used as the metabolic heat generation rates in muscle, fat, and skin layer. The power outputs of TEG obtained from both full and reduced thermal models are compared in Figure 15. It is observed that the results obtained from the pROM match with the results from the fullscale thermal model excellently. The relative error between full and reduced-order models ranges from À1.4% to 1.5% and the average absolute relative error is 0.9%, which means that the pROM of the thermal TEG model is a good approximation of the original full-scale system.

| pMOR for thermoelectric model
In this work, the finite element equations for the coupled thermoelectric model are defined as follows 41 : where To perform the stationary analysis of Equation (27) with DoF N = 92,427 DoF, it is time-consuming for parametric studies with different parameter values. Therefore, the pMOR method discussed in subsection 3.1 was also employed on Equation (27), which demands the extraction of the full system matrices from ANSYS. During this, there is no load resistance connected with the TEG device, leading to the open circuit and terminating the flow of current. Besides, it is worth noting that, as introduced in Section 2, the Thomson effect in TEG is negligible in this case, indicating that the Seebeck coefficients of the representative thermopile were no longer temperature-dependent. Furthermore, all the other material properties used in TEG were constant. As an outcome, there is no Peltier heat generated in TEG and the thermoelectric model tends to behave linearly in nature, defined as: The pMOR method was implemented on the linear thermoelectric model of Equation (28) Figure 14. The only difference is that, in the thermoelectric pROM, the EMF is set as the output from the reduced model directly because the Seebeck coefficients were integrated into Equation (28) with the coupling matrix [K vt ]. The parameter sweeps on fill factor and load resistance in Table A2 were performed to verify the accuracy of the thermoelectric reduced model. In Figure 16, it is investigated that the simulation results from the pROM of the thermoelectric model are close to the results from the thermoelectric full model. The relative error ranges from À4.0% to 7.1%, and the average absolute value of them is 0.1%. The computational time for the parametric simulations of the full and reduced TEG models are shown in Table 1. Based on the pROMs, the integration of 225 groups of parametric simulations takes only a few seconds. Comparing to the parameterized full-scale TEG models, the pROMs provide a significant speed-up of the computational time, which is thousands of times faster. This indicates that the small dimensional pROMs are not only highly accurate but also efficient for the simulations of electronic circuits at the system-level.

| Design optimization
Design optimization aims to determine the thermopile cross section of the implanted TEG, which provides maximum power output. As the voltage output of a realistic implanted TEG only reaches tens of mV a voltage boost converter has to be used to supply voltage levels required by the IMD. Therefore, the voltage output of TEG is expected to fulfill the minimum requirement of the booster converter. Typical booster converters for thermoelectric energy harvesting require about 20 mV input voltages. 42 During the design optimization (which also affects TEG resistance), the resistance of the load was adapted to maintain maximum power delivery. The full and reduced models both underwent the parameter sweep on fill factor and number of thermocouple legs. The optimization results of reduced models match well with results of full models (see Figure 17). In full thermoelectric model, a small spike occurs in the power output, which  may be caused by the optimum load, calculated through the Equation 8 and is not exactly the same as the internal resistance of the simplified TEG. When the fill factor is 6%, the maximum power output can be obtained, which is about 65 μW. The number of thermocouple legs has no impact on power output but the relationship with voltage output, and the higher voltage output can be observed with more thermocouple legs. With the optimum fill factor as 6%, the voltage output ranges from 6 mV to 100 mV with several thermocouple legs from 16 to 144. The optimized TEG needs at least 8 by 8 thermocouple arrays with 0.6 mm Â 0.6 mm cross-section for each leg, providing an output voltage of 22 mV with a matched resistive load.

| CONCLUSION AND OUTLOOK
In this paper, a detailed TEG model based on FEM was generated by reference to a commercially available TEG. The Seebeck coefficient and electrical resistivity of the thermoelectric material were experimentally obtained. Due to the low thermal resistance, the commercially available TEGs are not applicable for electrically active implants to acquire enough thermal energy inside human body, which motivates the design optimization of array configuration in TEG for maximum power output and sufficient voltage output. To facilitate the parameterization of array configuration of thermocouple legs, a simplified TEG model for thermoelectric analysis (thermoelectic TEG model) was presented and validated through the comparison in simulation results against the detailed TEG model. Instead of modifying the geometry of thermocouple legs in the detailed TEG model, equivalent material properties of a representative thermopile in the simplified TEG model were parameterized according to the fill factor and number of thermocouple legs. For computation efficiency and applicability of the pMOR algorithms, the thermoelectic TEG model was then decoupled as thermal TEG model for solely thermal analysis through the approximation of the Peltier heat as a part of effective conductive heat using effective thermal conductivity. After integrating housing and simplified human tissue, a thermoelectric model and a thermal model were developed based on the thermoelectic and thermal TEG model respectively and validated through the comparison in simulation results. Considering the power output and operability of implantation, the dimension of the representative thermopile was determined as 20 mm Â 20 mm Â 4 mm. In such thermal environment, an optimum fill factor of TEG was observed to obtain maximum power output with a matched resistive load, where the TEG thermally matched its implantation environment. For efficient design optimization, the reduced-order models were obtained from the full-scaled thermoelectric and thermal models through the pMOR method and demonstrated to be reliable against the full models. The reduced models can accelerate the simulation process thousands of times compared to full models. Moreover, the design optimization of array configuration of the implanted TEG was implemented both in full models and reduced models. With 21 C ambient temperature and 8.58 W/m 2 /K overall heat transfer coefficient, the maximum power output of TEG was found with an optimum fill factor of 6%, which is about 65 μW. With the same fill factor, the number of thermocouple legs has no impact on power output but voltage output. Higher voltage output can be obtained with more F I G U R E 1 7 (A) Dependence of power output on fill factor with the optimum load; (B) Dependence of voltage output on several thermocouple legs with the optimum fill factor and load thermocouple legs. To fulfill the minimum requirement of booster converter for input voltage, the final design of array configuration came out to have 64 thermocouple legs with a cross-section of 0.6 mm Â 0.6 mm, and the voltage output of TEG was expected to reach 22 mV.
For further research, several aspects can be extended or improved. In terms of simulation, the geometry of housing can be improved for higher power output of the packaged TEG inside. In the work of Rao et al. 43 metallic plates were used to replace a part of ceramic, allowing more heat flow through the TEG. The influence of housing depth, housing size and metallic plate size on power output was investigated. Additionally, the effect of clothing and body hair of the human body should be considered in the tissue model. In the aspect of the experiment, we are interested in the fabrication of the implantable TEG with such high-aspect-ratio thermocouple legs, where some challenges can be foreseen in the processes such as pick and place and soldering. Furthermore, aspects such as biocompatibility, hermeticity and thermal expansion have to be considered in encapsulation. We also look forward to testing the performance of our design in vivo. The surgical measurement data will help us to validate the model and make it more realistic.