Towards efficient design optimization of a miniaturized thermoelectric generator for electrically active implants via model order reduction and submodeling technique

Thermoelectric generators (TEG) convert the thermal energy into electrical energy and are under investigation as a power supply for medical implants. To improve the performance of TEG, the design optimization process through finite element model simulation is preferred by biomedical engineers. This paper aims to provide an efficient method of speeding up the design optimization process of TEG. A three‐dimensional realistic human torso model incorporating the TEG is investigated, where the internal heat transfer in human tissue is characterized by Pennes bioheat equation. In addition, convection, radiation, and evaporation effects at the skin surface are applied to identify the heat transfer effects between the human body and the environment. To speed up finite element simulation of the large‐scale human torso model, projection‐based model order reduction (MOR) is applied for generation of a compact but highly accurate model. Parametric MOR (pMOR) further enables generating a parameter‐independent compact model. For an efficient design optimization of TEG, this compact human torso model is applied within a thermal submodeling approach. Its temperature distribution results are back‐projected and used as boundary conditions for the TEG submodel. The achieved speed‐up in simulation time, demonstrated in this work, clearly indicates that the design optimization process of TEG is more efficient with the combination of MOR and submodeling techniques.


Abstract
Thermoelectric generators (TEG) convert the thermal energy into electrical energy and are under investigation as a power supply for medical implants. To improve the performance of TEG, the design optimization process through finite element model simulation is preferred by biomedical engineers. This paper aims to provide an efficient method of speeding up the design optimization process of TEG. A three-dimensional realistic human torso model incorporating the TEG is investigated, where the internal heat transfer in human tissue is characterized by Pennes bioheat equation. In addition, convection, radiation, and evaporation effects at the skin surface are applied to identify the heat transfer effects between the human body and the environment. To speed up finite element simulation of the large-scale human torso model, projectionbased model order reduction (MOR) is applied for generation of a compact but highly accurate model. Parametric MOR (pMOR) further enables generating a parameter-independent compact model. For an efficient design optimization of TEG, this compact human torso model is applied within a thermal submodeling approach. Its temperature distribution results are back-projected and used as boundary conditions for the TEG submodel. The achieved speed-up in simulation time, demonstrated in this work, clearly indicates that the design optimization process of TEG is more efficient with the combination of MOR and submodeling techniques.

| INTRODUCTION
The growing societal percentage of elderly people (eg, in Europe) will challenge the health care system in the next few decades. One out of three persons in Germany will be over 65 years old by 2050. 1 This trend leads to a large demand for medical implants for therapies such as bone and cartilage regeneration, deep brain stimulation to treat movement disorders, and fixing abnormal heart rates with pacemakers. However, the batteries currently used in these implants are limited in lifetime and necessitate replacement every 5 to 10 years. Moreover, chemical side effects and the large size of the batteries encourage engineers to develop energy harvesting technologies for self-powered implantable medical devices.
Energy harvesters transform energy present in the environment into electrical energy. Various applications exist together with dedicated energy harvesters to address the growing demand of ubiquitous electronic systems. 2 Medical implants, which integrate miniaturized energy harvesters, shall support the battery and thus prolong the implant's lifetime. Several research groups have focused on energy harvesting technologies in medical implants and several review reports have been published indicating the state-of-the-art of energy harvesting technologies for implantable devices. [3][4][5][6] Multiple harvesting solutions have been applied in medical implants, such as harvesting solar energy [7][8][9] and transferring its power wirelessly to the implants. [10][11][12] Alternatively, the kinetic energy from human motion and thermal energy from human body can be collected [13][14][15] and provided as electrical power.
The thermoelectric energy harvester transforms thermal into electrical energy through the Seebeck effect in thermoelectric materials. Biomedical engineers use such harvesters as a power supply in implantable medical devices by collecting the human body's thermal energy. In Yang et al, 16 the authors designed a model of a square-shaped thermoelectric generator (TEG) and integrated it in a human tissue model. Solving the Pennes bioheat equation enable to determine the temperature difference present across the TEG. Yang et al 15 demonstrated an implantable TEG and its potential as a continuous power supply for in vivo experiments. On the basis of Yang's paper, Jadhav et al 17 performed simulations of a finite element (FE) model of the TEG in a simplified model of human tissue. This model considers a constant metabolic heat generation in muscle tissue as the sole heat source. Mathematical model order reduction (MOR) was applied to this linear thermal model. The conventional MOR algorithms, [18][19][20] in which the original full-scale FE model is projected onto a lower subspace constructed by the orthonormalized Krylov-subspace, have been successfully applied to linear thermal models of various microsystems. 21 Further, Yuan et al 22 applied parametric MOR (pMOR) to generate a parametrically reduced simplified human tissue model for investigating the impact of environmental parameters on the temperature distribution inside the human tissue. In accordance with the multivariate moment-matching algorithms, 23-26 a parameterindependent reduced order model of human tissue could be obtained. pMOR has been successfully applied in various case studies. 22,27,28 Through the application of MOR and pMOR algorithms, a highly accurate and compact human tissue model could be generated for efficient simulations together with TEG model. Such simulations study a simplified human tissue model, which consist of three layers for the muscle, fat, and skin section. Heat loss at the skin surface is modeled via a convection boundary condition. Verma et al 29 built a model based on a realistic geometry of a human forearm model. Convection and radiation effects were considered as the external heat transfer effects, though no MOR algorithm was applied. In this paper, we aim to introduce the FE model simulations of a realistic human torso model incorporating the TEG. Besides convection and radiation, also evaporation from the skin surface (sweating) has been considered. The TEG was positioned in the fat layer, where the maximum temperature gradient exists, and thus maximum power output is to be expected. 16,17,29 To investigate the temperature difference across the TEG model, MOR and pMOR algorithms were applied to generate the reduced-order human torso model for efficient thermal simulations. Because of the nonlinear radiation and evaporation heat transfer effects, several linearization strategies were applied to the model. On the basis of the simulation results from the reduced human torso model, a submodeling technique was further implemented, which separated the simulations of the human tissue model and the TEG submodel. The simulations in the submodel enable the efficient design optimization process of TEG.
In Section 2, the models of the TEG and human torso are presented. In Section 3, the numerical simulation methods are illustrated. The submodeling technique is introduced. Linearization strategies are then used to generate a linear thermal human torso model. MOR and pMOR algorithms are further employed for generating a compact model based on this linearized thermal model. The temperature distribution obtained from the reduced-order human torso model are applied to the submodel. In Section 4, the accuracy and applicability to parametric studies of the parameterized reduced-order human torso model is verified. The accuracy of the temperature distribution across the TEG inside the submodel is demonstrated as well. The conclusion of this work and the outlook for further research on this topic are given in Section 5.

| CASE STUDY
In this section, the three-dimensional FE model of a disk-shaped TEG model and a realistic human torso model were developed in ANSYS (Version 2019 R1, Ansys Inc.). The TEG was further placed in the fat layer in the chest region of the human torso model. The aim of the simulations was to find out the temperature distribution result on TEG.

| Thermoelectric generator
Thermoelectric generators enable the conversion of heat into electrical energy. The assembling setup of an electrically active implant is shown in Figure 1A. It contains a TEG, an energy buffer, and an application-specific integrated circuit (ASIC). The thermocouple legs are connected with copper connections. The implant is positioned in the fat layer of the human tissue where the maximum temperature is found. 16,17,29 In this paper, the geometry of the TEG FE model was constructed based on a commercially available TEG (see Figure 1B) in ANSYS. An array of 16 × 16 p-type (hole transporting) and n-type (electron transporting) bismuth telluride thermocouples with a height of 2.27 mm and a cross section of 0.8 × 0.8 mm were located between the upper and lower ceramic square plates made of aluminum oxide (24.6 × 24.6 × 0.565 mm) without considering the copper interconnections. The overall height of the TEG is 3.4 mm. Further, the TEG was surrounded by a disk-shaped housing made of Teflon with the same height. After spatial discretization, the FE model of TEG contains 127 307 nodes in total. The material properties used for each part in the TEG are shown in Table A1.
Through the Seebeck effect on thermoelectric material, the temperature difference between two plates of TEG creates an electrical potential in the thermocouples. A thermocouple is made up of one p-type thermocouple leg and one n-type thermocouple leg. Thermal energy from the human body is then converted into electrical energy. The Seebeck voltage output from the TEG could be calculated using the equation: where the voltage output generated in TEG is proportional to n and ΔT, which presents the number of thermocouple legs in TEG and the temperature difference generated between the top and bottom surfaces of the thermocouple legs. α 1 ,α 2 are the Seebeck coefficients.

| Human torso model
In order to obtain the realistic human tissue temperature distribution results, a half human torso model consisting of realistic geometry of the solid internal organs, skeleton, and main vessels, as well as muscle, fat, and skin layers, was  Figure 2A). Realistic thermal data and physiologically correct material parameters were assigned to the various tissue sections 31 (see Table A2). Subsequently, as shown in Figure 2B, the TEG model was placed in the fat layer of the human torso model in the chest region. The final realistic human torso FE model incorporating the TEG contains 1 045 923 nodes in total.
To characterize the internal heat transfer in human tissue, the Pennes bioheat equation is used 32 : where ρ, c, and κ are the density, specific heat capacity and thermal conductivity properties of different tissues. T is the unknown nodal temperature distribution results on the human tissue and TEG. Furthermore, in Equation (2), the blood perfusion and metabolic heat generation rates, Q b and Q m , are considered in muscle, fat, and skin layers while no heat is considered by blood perfusion and metabolism in other internal organs. T t and T a are the body tissue temperature and arterial blood temperature, respectively. ρ b and c b are the density and specific heat capacity of the blood and ω is the blood perfusion rate in body tissue. In this case study, instead of using a central blood pool model, 33,34 we assume that the temperature of arterial blood T a is constant at 37 C when it travels from the heart pool to other body parts. This assumption neglects the temperature variations in arterial blood and obviously overestimates the effects of blood perfusion in tissue. However, with the constant arterial blood temperature, it is easier to apply the blood perfusion effect in the model. In the future, if this assumption limits the accuracy of the blood perfusion source term, one can still overcome it by adjusting the arterial blood temperature or the blood perfusion rate in the Pennes bioheat equation. 35 The material properties of the tissues and values of Q m in muscle, fat, and skin tissue layers are shown in Table A2. Meanwhile, for the human torso model, the heat generated inside is balanced instantaneously by the various external heat transfer effects. In this study case, the convection, radiation, and evaporation effects at the skin surface were applied as the main heat losses. The contributions of heat losses for each effect were calculated as 29.0%, 38.1%, and 24.2%, respectively. 36 No boundary conditions were imposed at the truncated neck, extremities, and the lower torso, but only at the skin surface. The total heat transfer between the skin surface and the environment could be described as the equation as shown in as follows 37 : where q conv , q rad , and q eva are the convection, radiation, and evaporation heat fluxes normal to the boundary skin surface. T skin is the temperature distributed at the skin surface and T amb is the ambient temperature. Firstly, in the convection effect, h c is the convection heat transfer coefficient, which is defined as an air velocity-dependent parameter 37 in W/m 2 /K: where v air is the air velocity in m/s. Secondly, for the thermal process of radiation, it is described by the Stefan-Boltzmann equation with the Stefan-Boltzmann constant σ and emissivity ϵ. Finally, 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, 38 which shows that the evaporation coefficient could be presented in terms of convection coefficient, and the Antoine equation, 37 the evaporation heat loss q eva in Equation (3) could be expressed with skin wetness w skin as follows:

| NUMERICAL SIMULATION METHODS
The FE model of human torso model incorporating the TEG was built to find the temperature difference generated in TEG. The normal design optimization process of TEG is shown in Figure 3A. The temperature distribution result on TEG is obtained after the solution of the FE model. However, because of the complex geometry of the realistic human torso model and TEG, a large-dimensional FE model was generated (1 045 923 nodes in total) and the simulations of such model were time-consuming. Therefore, in this paper, we presented a new method, which combines the MOR and submodeling technique (see Figure 3B). With this new method, the simulations of TEG could be handled in a relatively small human tissue submodel, whose dimension of the FE model was much more smaller (132 826 nodes in human tissue submodel incorporating the TEG) than the global human torso model.

| Submodeling technique
The submodeling technique has already been developed and applied in ANSYS (Version 2019 R1, Ansys Inc). The idea of the submodeling technique is to separate the simulation of a complex model into two parts. For example, in our study case, a representative TEG was first placed into the global human tissue model (see Figure 4A). The TEG's structure of thermocouple legs was replaced at this time with a block structure, where fewer elements were generated because of spatial discretization (127 307 nodes in detailed TEG and 2694 nodes in representative TEG). Secondly, the detailed TEG was placed in a submodel. It was surrounded in the simulation by only a little human tissue in order to minimize the influence of human tissue on the TEG (see Figure 4B). The temperature distribution results calculated from the human torso model incorporating the representative TEG were used as temperature boundary conditions in the detailed TEG submodel. With the submodeling technique, it was possible to do simulations of the TEG for design optimization in the submodel only. Given that fewer elements were generated in the submodel, correspondingly less computational effort was used during the simulations. Although the submodeling methodology speeds up simulations by means of simulating the TEG in a relatively small submodel-obviating, as it does, the repeat simulations of the whole human torso model should any modification to the TEG be required-a full-scale torso model simulation is still indispensable for obtaining the same temperature boundary conditions as those used in the submodel. Therefore, the MOR algorithm was applied in this study case to generate an accurate and highly compact model of a human torso model incorporating the representative TEG for the rapid simulation of temperature distribution results. Moreover, if the environmental conditions surrounding the human body are changed, one needs to generate different reduced models reflecting different environmental conditions. To solve this, the pMOR method was applied to generate a parameter-independent reduced model. The variables in the environmental boundary conditions could be considered to be the parameters in the reduced model and the MOR process was only required once.

| Linearization strategies
The conventional MOR and pMOR algorithms are designed for the linear systems. However, the human torso model is presented by the bioheat Equation (2), where the blood perfusion heat generation rate is resulting temperaturedependent input. This leads to the nonlinearity in the thermal model. Moreover, the nonlinear radiation and evaporation heat transfer effects at the skin surface were applied to balance the internal heat conduction in the human tissue. Therefore, the linearization strategies were introduced in this paper to construct a linear thermal model of human tissue.
According to the analogy between the blood perfusion heat generation rate Q b and convection boundary condition q conv , the blood perfusion heat generation rate could be treated as a "convection-type" effect. 39 This nonlinear temperature-dependent blood perfusion heat generation input could be transformed as a convective heat transfer boundary condition: The blood perfusion rate coefficient ρ b c b ω in Equation (6) is analogous to the film coefficient in the convection boundary condition. It was set as negative to present the heat flowing into the human tissue as a heat generation input.
In the next, the radiation effect in Equation (3), which describes the continuous energy interchange between the human body and the environment by means of electromagnetic waves, was linearized if the temperature difference ΔT = T skin − T amb is small. 40 The linear relationship of radiation effect could be obtained by expanding q rad as Taylor series around T amb , Comparing Equation (7) to the convection boundary condition q conv , we got a radiation heat transfer coefficient h rad = 4σϵT amb 3 . The radiation effect was then further integrated with the convection boundary condition through the sum of convection and radiation heat transfer coefficients. The last nonlinearity occurs in the evaporation heat transfer effect in Equation (5). No linear relationship could be found in this heat transfer effect. Therefore, we approximated the average heat flux generated from the evaporation heat transfer effect on each node of skin surface through the snapshot method. After the FE spatial discretization, r nodes were obtained at the skin surface of the human body. From Equation (5), the heat fluxes generated on each node at s points in time (t i ) was collected and merged as a snapshot matrix X ∈ R r × s : where the vector q i (i = 1, 2, 3Á Á Ás) presents the vector of heat fluxes on each node of the skin surface at each time point. The average heat flux on each node was obtained by taking the weighted average of these snapshots: where w i denotes the weights applied to generate the average heat flux on each node of the skin surface. They were defined based on the ratio between the heat flux in each time step and the total time. In addition, in the representative TEG, the linear thermal conductivity property of the block material was calculated based on the experimental data, 41 combining with these linearization strategies, the linear thermal FE model of human torso model incorporating the representative TEG was presented as follows: where A, E ∈ R N × N are the global heat conductivity and heat capacity matrices. The blood perfusion and radiation heat transfer effects in Equations (6) and (7) were applied as convection effects and integrated in global heat conductivity matrix A. B ∈ R N × m is the input distribution array and u is the input vector. Here, the input vector contains the constant metabolic heat generation rate Q m and the approximated average heat flux input q eva on each node at the skin surface. C ∈ R p × N is the user-defined output matrix. N is the dimension of the system and m and p are the number of inputs and user-defined outputs.

| Model order reduction
To speed up the simulations of the human tissue model incorporating the representative TEG in this paper, the Krylovsubspace-based MOR method was applied to generate a highly accurate and compact model based on the linearized FE model as shown in Equation (10) (N = 921 336). This method projects the original large-scale system onto a low dimensional subspace V ∈ R N × n with n < < N. It was constructed by the orthonormal bases of the Krylov subspace. The full-scale temperature state vector T in Equation (10) can then be replaced by the projection equation: where x(t) ∈ R n is a reduced temperature state vector. According to the Krylov-subspace-based moment matching approach, 18-20 the subspace is found based on the moments (Taylor coefficients) in the transfer function of Equation (10): The expansion coefficient m i are called the moments of the transfer function around the expansion point s 0 . Then, the projection matrix V is defined as an orthonormal base of the following Krylov subspace: By applying projection Equation (11) on linear thermal system Equation (10), a low-dimensional model of order n is then obtained: where E r , A r , B r , and C r are the reduced system matrices. The moments in the transfer function of the reduced system presented in Equation (14) is defined as: It has been proven that the first n moments of m i and m i (r) in Equations (12) and (15) are matched, 18 which permits the accuracy of the reduced system.

| Parametric model order reduction
The normal MOR algorithm generates an accurate reduced model for efficient simulations of the linear thermal system. However, in our study case, if the environmental conditions around the human body are changed, the MOR has to be repeated for generating a new reduced model with the newly specified boundary conditions. Therefore, it is essential to generate a reduced model which preserves the parameters in the environmental boundary conditions. In the previous research, 22,27,28 the authors have shown that through the multivariate moment-matching pMOR methods, [23][24][25][26] it is possible to extract the convection film coefficient and the ambient temperature as the parameters in the reduced model. In this paper, based on Equation (10), the parametrized full-scale human torso linear thermal model was obtained as follows: According to the linearized radiation heat transfer effect (Equation 7), the radiation coefficient h rad was added with the convection heat transfer effect h c . The parametric studies were carried with only the parameter h defined in Equation (16): The parameter of ambient temperature in convection boundary condition q conv was moved to the right-hand side of Equation (10) and applied as an input in input vector u.
To obtain a parametric reduced order model, a parameter-independent projection subspaceṼ was constructed based on the transfer function of (16), which contains two physically-independent variables (s and h) 22 : The moments with respect to the variables h and s were obtained separately. When we calculated the moments of one variable, the other variable was kept as constant: Equations (19) and (20) shows the moments m i s and m i h around variable s and h fixed at points h = h 0 and s = s 0 , respectively. The Krylov subspaces of each parameter were generated separately as: colspan The final global parameter-independent projection matrixṼ was constructed by merging the orthonormalized bases of V s and V h : Similar to Equation (14), the reduced parametric order model was obtained as: where E r , A 0r , A 1r ∈R n 1 + n 2 ð Þ× n 1 + n 2 ð Þ , B r ∈R n 1 + n 2 ð Þ× m and C r ∈R p × n 1 + n 2 ð Þ . The moments of the reduced parametric model around variable s and h were defined as: It had been proposed and demonstrated that 22,23,27,28 an accurate parameter-independent reduced model would be generated through neglecting the mixed moments and matching the first n 1 and n 2 moments with respect to s and h separately.

| SIMULATION RESULTS
In this section, submodeling and MOR techniques were applied to the realistic human torso model incorporating the TEG. The simulation process was separated into two parts. Firstly, a representative TEG was placed in the human torso model. The pMOR reduction method was applied to generate an accurate parameter-independent reduced model, where the convection heat transfer coefficient at the skin surface and the ambient temperature were applied as the independent parameters. Then, the reduced temperature state given by the reduced model was projected back to its full size through projection Equation (11), and the full-scale temperature distribution results in the human torso model were approximated. Secondly, these approximated nodal temperature distribution results were used as the temperature boundary conditions in the human tissue submodel incorporating the detailed TEG.
In the full-scale human tissue model incorporating the representative TEG, the temperature distribution results were obtained in two steps. Firstly, an initial thermal state was simulated through the steady state thermal simulation with the convection film coefficient h c = 9.49 W/m 2 /K and the ambient temperature T amb = 15 C. On the basis of this initial state, a transient thermal simulation was carried out in the second step, where the convection film coefficient was changed to h c = 4.61 W/m 2 /K. Further, to verify the accuracy of the parameter-independent reduced human torso model, different values of the convection film coefficient and the ambient temperature h c = 6.98 W/m 2 /K and T amb = 25 C were used. The details of these selected parameter values are shown in Table 1.
The values of the convection film coefficient and radiation coefficient were calculated depending on the air velocity and ambient temperature, respectively. We imitated environmental conditions in an indoor office by selecting a pleasant air velocity between 0.2 and 1.5 m/s, and ambient temperature was selected between 15 C and 25 C according to European Pharmacopoeia.
Because of the reason mentioned in Section 3, in order to apply the submodeling and pMOR techniques, a representative TEG was positioned in the human tissue model to replace the detailed TEG model. Thus, it is essential to investigate the influence of the representative TEG on the temperature distribution result in human tissue. We processed full-scale thermal simulations of the human torso model, with details, and the representative TEG separately. The temperature results in the path selected at the skin surface of the submodels were compared (see Figure 5).
From the comparison result shown in Figure 5, the maximum temperature relative error between the nodes on the selected paths in the detailed and representative TEG implanted submodels is 2.83%. This result shows that the representative TEG will not influence the temperature distribution results to excess and they can be used as the temperature boundary conditions in the submodel.
After that, a parametrically reduced order model was generated by applying the pMOR algorithm in the software "Model Reduction inside ANSYS." 42 The accuracy of the parametrically reduced model was verified as shown in   In the parametric reduced model, n 1 = 31 moments around frequency parameter s was enough and n 2 = 5 additional moments around film coefficient parameter h was required to give the accurate results. It was found that the norm of the first five column vectors in projection matrices V s and V h were the same. A deflation was applied when merging these two projection matrices. In other cases, the accuracy of the reduced model can be improved by increasing the values of n 1 and n 2 . The temperature results selected from the three nodes in the muscle, fat, and skin layers were compared. The maximum relative error between the full and reduced models is 0.778%, which shows that the parametrically reduced model is an excellent approximation of the full-scale human torso thermal model.
Subsequently, we used the accurate temperature results as obtained from the parameter-independent-reduced model to approximate the full-scale temperature distribution result on the human torso model by projecting the reduced temperature state vector back to full size according to Equation (11). These full-scale temperature results were later used as the temperature boundary conditions in the detailed TEG implanted human tissue submodel.
In the human torso model incorporating the detailed TEG, the convection, radiation, and evaporation heat transfer effects at the skin surface and the blood perfusion and metabolic heat generation rates in muscle, fat, and skin layers were all still applied. To verify the accuracy of the thermal simulation in the submodel, we compared temperature results between two detailed TEG models: one was simulated in the global human torso model and the other was simulated in the submodel (see Figure 7). The maximum temperature error between the nodes in the selected path in the TEG model was 0.516 C, which was calculated with 1.64% as its relative error. This result shows that through the MOR and submodeling techniques, the accuracy of temperature distribution results in the TEG model is guaranteed.
Finally, we compared the simulation time used for generating the temperature results on the TEG with and without MOR and submodeling techniques (see Table 2). For the full-scale thermal human torso model incorporating the detailed TEG, the total simulation time was 1515.05 seconds. Comparing with the total computational time in the submodel (669.34 s), the computational effort was improved 2.26 times. In addition, there is no need to repeat the steps of pMOR and reprojection if the environmental condition is not changed. The optimization process of the TEG could be handled with the submodel in just 181.96 seconds, which is 8.33 times faster than the simulation in the full-scale torso model.

| CONCLUSION AND OUTLOOK
In this paper, a newly designed TEG was introduced and positioned in a realistic human torso model, which is geometrically constructed based on MRI data. Blood perfusion effect including the bioheat equation was used to model internal heat transfer conditions in human tissue and the convection, radiation, and evaporation heat transfer effects were applied as the external heat transfer effects at the skin surface of the human torso model. Several linearization strategies were then introduced in sequence to construct a linear thermal model of the human torso for applying the MOR and pMOR algorithms. The blood perfusion heat generation rate was treated as a "convection-type" effect where the film coefficient was defined as negative to present the heat flowing into the human tissue. A linear relationship was found in the radiation heat transfer effect and an equivalent radiation coefficient was used, which was later added to the convection heat transfer coefficient in convection effect. The average evaporation heat flux on each node at the skin surface of the human torso FE model was used as input, which was modeled by calculating the summarized weighted heat flux snapshots at each selected time step. On the basis of this linearized thermal model, we introduced the methodology of MOR and pMOR to speed up the thermal simulation in the human torso model. Furthermore, through the submodeling technique, the simulation of the TEG was separated into two parts. One was a simulation of the global human torso model with a representative TEG. We compared the selected temperature results generated in torso models, into which both representative and detailed TEGs were implanted, respectively. The influence of the representative TEG in the torso model was negligible because of the small error observed from the comparison. Secondly, the pMOR algorithm was applied to generate a parameter-independent-reduced model and the fullscale temperature effect on the torso model was then approximated through projection Equation (11). The accuracy of the parametric reduced model was verified by comparing the temperature results from the selected nodes in full and reduced human torso models. This approximated full-scale temperature result was then used for the temperature boundary conditions in a human tissue submodel incorporating the detailed TEG. The submodeling technique enables the simulations of TEG in a relatively small submodel with less computational effort. In future, further external heat transfer effects, such as respiration, would be taken into consideration in a more realistic human body model, where some but not all of these heat transfer effects could be linearized. Therefore, it will be essential to apply nonlinear MOR and pMOR algorithms to generate the reduced models. In the design optimization process of TEG, we are interested in choosing the right geometry for the thermocouple legs and finding the optimal array configuration of them. The thermocouple legs' cross section area, height, and array configuration will have impacts on both the thermal properties and electrical performance of the TEG. To balance these opposite effects, these geometrical parameters can be considered in the optimization application. Further, system-level simulations based on the reduced TEG model are intended to be constructed for providing a suitable power management circuit for this power supply generator.