Topology optimized thermoelectric generator: a parametric study

Exploring new structure

The exploration for a new compliant structure involves the removal of material to determine if it can support structural loads while maximizing power generation. The effort in the study is geared toward optimizing the holes shown in Figure 3.1 model in order to find the best void structure and develop optimal configurations for different loads for any environment.

Figure 3.1:

Exploring new structure.

Process

Figure 3.2 shows a flowchart that describes the process used to create a unit cell and test it for structural loading. The flowchart is a representation of the standard finite element based tool developed by Sigmund (2001). It starts with inputs such as the perimeter area to create the design space. Boundary conditions and loads are prescribed. Input of the design variable is loaded. An even distribution of the design variable is made across the design space. The MATLAB code representation of the model is run as depicted in Figure 3.2.

Figure 3.2:

Standard topology optimization flow chart.

As shown in the flowchart, also discussed by Liu and Tovar (2014), the iterative process continues with termination occurring when a maximum number of iterations is reached or where the change in material distribution from iteration to iteration is within a specified tolerance; for example, if available material is equal to or less than 0.01. Convergence takes place because the design variable is less than the convergence criteria.

The finite element code makes use of the sparse option in MATLAB. The top left and top right node element numbers are used to insert the element stiffness matrix at the right places in the global matrix. Both nodes and elements are numbered column-wise from left to right. Each node has a designation of two degrees of freedom which are translational for horizontal and vertical of an element node (Figure 3.5). Fixed and free degrees of freedom are easily entered in the code. The Youngs modulus and Poisson’s ratio are inputs to the code.

Step one in the process first requires specification of the applied loads, support conditions, volume domain of structure to be constructed, and geometry of the void and solid areas. In developing a design domain, a volume fraction (φ) is defined as a ratio of the object volume versus the available space volume in the material domain (Ω).

Finite elements

The design domain is discretized into elements. A single element is shown (Figure 3.3). Each element has four nodes each with two degrees of freedom (dof); totaling to eight dof. K_e is the element stiffness matrix with a dimension of 8 × 8.

Figure 3.3:

Element showing eight dof.

A 2D-quad element illustrated in Figure 3.2 has four nodes. Each node has the potential to translate in the x-direction (horizontal) and in the y direction (vertical), when a force is applied to the element. The force applied is proportional to element stiffness and resulting displacement. The relationship is captured in Eq. (3.1).

where the applied force f, is equal to the product of element stiffness K and its displacement u. The stiffness matrix K depends on the vector ρ of the element-wise constant material densities in each element, e = 1 … N. Therefore, K can be written as

where K_e is the global element stiffness and p is the power which penalizes intermediate densities. A p value of three or greater is preferred (Bendsoe and Sigmund 2004; Liu and Tovar 2014) because a high p level makes it uneconomical to have intermediate densities in the optimal design. Further, a p value of three or greater will result in obtaining a true “0–1” design for an active volume constraint (Bendsoe and Sigmund 2004).

Structural model validation

Bendsoe and Sigmund (2004) discuss how search for the optimal lay-out of structure with a specified region requires applied loads, support conditions, volume of structure to be constructed and location and size of void solid areas. Mativo and Hallinan (2017) use this approach to validate extensively the structural tool. The following three validation examples are provided.

Effects of constraints, and material distribution

A design domain uses support conditions such as fixed elements or nodes that may restrict displacement or transmission of loads. The design function in the design domain is density. The design problem for the fixed domain is formulated as a sizing problem by modifying the stiffness matrix so that it depends continuously on a function which is interpreted as a density of material. Figure 3.4 a and b shown below depicts material distribution that show size, shape, and topology of structures. The volume fraction is symbolized by φ. The difference of the topologies are results of point load applied at two different positions of design domain to support conditions imposed at the bottom of both left and right of the design domain.

Figure 3.4:

Effects of points of load application and support restrictions.

(a) φ = 0.40, F = top middle (b) φ = 0.40, F = two-thirds on top right.

Compliant mechanism

Sigmund (2007) explores the design of compliant mechanism using topology optimization. One of the most important objectives in compliant mechanism synthesis is to be able to control ratios between output and input displacements or output and input forces which are described by the geometrical and mechanical advantages respectively. Two examples shown in Figures 3.5 and 3.6, initially presented in Bendsoe and Sigmund (2004), were replicated by Mativo and Hallinan (2017), to validate the ability of the tool to create compliant mechanisms both symmetric and non-symmetric.

Figure 3.5:

(a) A basic compliant symmetric mechanism – example of and load model. (b) Optimized topology of the basic compliant mechanism of (a) (Mativo and Hallinan 2017).

Figure 3.6:

(a) A basic compliant non-symmetric mechanism – example of and load model. (b) Optimized topology of the basic compliant mechanism of (a) (Mativo and Hallinan 2017).

The tool validation for structural designs is affirmed. The tool is however, modified to address the current study.

Optimizing for structural model

Table 3.1 lists the design parameters considered for the structural model. These parameters include: design domain, geometry, fixed void region, fixed solid region, variable solid region, fixed nodes, free nodes, uniformly distributed compression loads, shear displacement, adjoint displacement, and mechanical material properties.

The imposed boundary conditions are represented as follows:

1. Fixed end: x₀ = 0; x_f = 0; where x₀ is initial position of x, and x_f is final position of x.

2. Free end: x_L′ = 0; x_L″ = δ; where x_L′ is the initial position at the top of the element, and x_L″ is the final position at the top of the element, and δ is the lateral displacement.

Also, the length of the legs is fixed, but the width of the legs can vary subject to the ability to sustain the imposed structural loads.

A matrix determined analytically using symbolic manipulation software by Bendsoe and Sigmund (2004) is modified and used in this numerical experiment. Earlier, in Figure 3.3 an illustration of the elemental degrees of freedom was given. It follows that the matrix terms k(1), k(2), …, k(8) represent node stiffness. The k_ij conditions apply for each node, such that a force in the “i” direction to have a deformation of 1 in the “j” direction only. Table 3.1 provides physical properties of the Young’s modulus and Poisson’s ratio used in the matrix computation. A Young’s modulus (E) of 8.1 and Poisson’s ratio (ν) of 0.23 are used in the matrices.

where k_e is the global element stiffness matrix.

Objective function

A representative static load is applied in the design tool to seek a design optimized for the structural loadings. No vibratory loads are applied in the design tool. The process to achieve this desired outcome is to reduce structural material carefully without compromise of structural fidelity. The systematic material reduction is conducted by manipulation of the volume fraction for the TEG leg. Therefore, an objective function for maximization of the TEG leg displacement has the form in Eq. (3.3).

subject to the following constraints: 0≤ρ<1, i.e., the material distribution (design variable) between 0 and 100%

where ρ is the element-wise constant material densities in element, e = 1, … , N. The u^* is the calculated minimum displacement, Eq. (3.4), that is able to supports the loads. The stiffness matrix K depends on ρ. r is the residual in obtaining the structural equilibrium. For topology optimization the equilibrium condition, r = 0, is found using an iterative procedure. In this equation, u and f are displacement and load vectors, respectively the shear and adjoint loads are unit valued displacements at the boundary. The displacements can be applied in lieu of a force on a node (Porter 1994).

The preferred displacement u* should be ≥0.029. This calculated value of 0.029 is considered the minimum displacement from a 100% volume fraction of a unit cell. A displacement greater than u* would allow the TEG leg to absorb vibrations it currently cannot. The u* value is obtained using the following equation:

where P is a point load, L is the length, E is the modulus of elasticity, and I is the inertia.

The goal is to find a maximum displacement of TEG legs that will ensure that the mechanical loads are sustained (i.e., that there is sufficient compliance to tolerate the loads) – ultimately while maintaining the ability to generate maximal power.

Exploring acceptable topologies

The interest in this study is to create compliant TEGs. Therefore, the topology optimization process controls the placement of material in the two legs of the unit cell. A successful optimized structural model requires that developed TEG should be able to accommodate translation induced by the shear load. Each leg is expected to support translation. A leg with one point connected to the TEGs top cover would allow a minimal lateral displacement before it results to undesired potential rotation. Two points connecting the TEG leg to the top cover would allow a higher lateral translation with less induced friction than three or more points would cause. According to Alread and Leaslie (2007), there are three basic ways that structural members can be connected. Pin connection allow members to rotate but not translate, moment connections allow neither translation nor rotation relative to one another, and rollers that allow only translation. Our design follows the latter.

Table 3.2 shows the optimal TEG structures obtained for variable porosity. In the table, the topologies shown in columns 2–4 are created by performing two changes, first by varying the TEG leg width size, an action done by varying width of the street, and secondly by changing the volume fraction or the amount of the TEG material in the legs.

Table 3.2:

Optimal TEG considering only structural loadings alone for varied leg size.

Volume fraction	2 × Baseline street (2B)	Baseline street (B)	½ Baseline street (½B)	Displacement
100				2B (u = 0.1259)
				B (u = 0.0504)
				½B (u = 0.0286)
90				2B (u = 0.1259)
				B (u = 0.0555)
				½B (u = 0.0376)
80				2B (u = 0.1260)
				B (u = 0.0636)
				½B (u = 0.0494)
70				2B (u = 0.1630)
				B (u = 0.0784)
				½B (u = 0.0688)
60				2B (u = 0.1902)
				B (u = 0.1031)
				½B (u = 0.0962)

For our case, a baseline model B is a topology whose street width size is equal to the width size of each leg as shown in column 3. Topologies with a double street width size 2B, halves the leg width size, are shown in column 2, while those with half the baseline street width size ½B are presented in column 4.

The volume fraction is varied from 100 to 60%. This effort is done to examine topologies whether topologies created with a lower volume fraction of TEG material can achieve greater compliance (u) while supporting the structural loads. The volume fraction is presented in column 1, while the displacement is shown in column 5.

Results of the effort presented in Table 3.2 show topologies that are suitable for a vibratory environment and those that are not. A suitable topology is a leg with two points that connect to the top plate. The two points indicate ability for the leg to translate without causing a moment. From the table, and enlarged topologies in Figure 3.7 a and b, it is observed that the baseline models presented in column 3 with a volume fraction around 80% have suitable topologies. Other leg topologies have one point touching the top cover which would restrict desired translation and possibly initiate rotation.

Figure 3.7 a:

Suitable leg with two distinct points [1 & 2] touching the top cover.

Figure 3.7 b:

Unsuitable leg with one point touching the top cover.

Figure 3.7 a and b are topologies illustrating legs with one, or two points touching the top cover.

Table 3.3 shows the optimal topology as a function of the porosity of the TEG for fixed leg width. Here, the porosity ranges from 0 to 0.8. Each of the shown cases represents the maximum transverse translation at a given porosity. The associated shear translation respectively ranges from 0.0504 to 0.6079. As the porosity increases, the shear translation increases. Thus, more porosity increases the compliance to shear loadings. The largest displacements occurred between volume fractions between 0 and 30%. This displacement is however not considered in this research due to minimal material to even converge and form TEG legs.

Table 3.3:

Percent material, displacements, and topology.

Volume fraction	Displacements (u)	Structural topology
100	0.0504
90	0.0555
80	0.0636
70	0.0784
60	0.1031
50	0.1290
40	0.1420
30	0.1778
20	0.2340
10	0.6079
0	0.6079

Of particular interest is the leg pattern at the top. Topologies with a volume fraction about 80% have two points on each leg. Two points are important because they indicate the ability to accommodate lateral movement at the top of the TEG unit cell which would allow use in a vibratory environment. While Table 3.2 depicts exploration of topologies and displacements of three different leg sizes with variable volume fraction, Table 3.3 presents further exploration of the optimized leg size that emerged from Table 3.2, in pursuit of an optimal displacement as it explores from maximum to minimum volume fractions for the TEG leg.

Figure 3.8 presents the transverse displacement as a function of volume fraction. It is clear that when the volume fraction exceeds 40%, the slope of the displacement reduces markedly. Thus, in a sense, the structural integrity stabilizes rapidly from this point on. In contrast, with a volume fraction less than 40%, the displacement slope increases rapidly. This may not be bad. If the environment calls for more compliance a higher void percentage is required.

Figure 3.8:

Plot of shear displacement versus volume fraction for structurally optimized TEG.

Table 3.4 shows the structurally optimized topologies for volume fractions ranging from 76 to 84%, with the corresponding topologies and displacements shown in columns 2 and 3. It can be observed that the topologies for volume fraction of 80% and above have legs with the same form. The leg topology at a volume fraction of 78% displays a slightly different form for each leg, which could be a concern for even translation. The optimal topology associated with a 76% volume fraction has legs connected only at one point at the top cove. Therefore, the topologies associated with volume fractions beneath 80% here are not feasible.

Table 3.4:

Structurally optimized topologies for volume fraction near 80%.

Volume fraction (%)	Topologies	Displacement (u)
84		0.0604
82		0.0620
80		0.0636
78		0.0662
76		0.0688

Thermal model

Optimal topologies for thermal loadings alone based upon the logic shown in the flowchart given by Figure 3.2 are likewise developed for variable volume fractions. A one-dimensional steady-state heat conduction model is assumed for the unit cell.

The heat transfer rate is governed by Fourier’s Law (Eq. (3.5))

where q_x is the x-component of the heat transfer rate, k is thermal conductivity of the material, A is the area normal to heat flow, and ∂T∂x is the temperature gradient. When the length, l, of the member or element is considered, Eq. (3.9) can be re-written as (Eq. (3.6))

Since we have a constant thermal conductivity, k, Laplace’s equation governs the change in temperature in the x-direction.

Integrating Eq. (3.7) twice results in Eq. (3.8).

To solve for c₁ and c_2, we assume two boundary equations as follows: at x = x₀, T = T₀, and at x = x_L, T = T_L. The solution to Eqs. (3.7) and (3.8) are as follows in Eqs. (3.9a) and (3.9b).

Similar to the structural model, the thermal model has a conductance matrix [K], a temperature matrix {T}, and the heat source matrix {q}. These matrices are created and assembled for each leg and computed to find thermal flow depicted as q = KT. Following the procedure of Bendsoe and Sigmund (2004), a 4 × 4 conductance matrix was developed to describe the thermal model.

Here k_e is thermal conductance.

The heat conductance in a global form can be written as K_thT = F_th. T is the temperature, and F_th is the thermal load vector which characterizes the heat flux at the boundary. K_th is the conductivity matrix and takes the form

where N is the number of elements, and ρ is the design variable. The elemental thermal conductivity K^e_th is given by:

Here, Q/t is the rate of thermal energy transferred through the material. The thermal conductivity is represented by k. The cross-sectional area is represented by A. The difference in temperature is represented by ΔT. The thickness or length through which heat transfers is represented by d.

The thermal design parameters shown in Table 3.5 include geometrical characteristics, material properties, and heat source and heat sink temperatures.

Ultimately, the maximum heat flow is dependent on the structure. Rowe (2006), states that the ideal geometry for thermo-elements used in TEGs should be long and thin (p. 1–11).

Objective function

In this model, the TEG legs are assumed to be made of homogenous material (Figure 3.9). The leg volume has a uniform cross-sectional area and perfect insulation is assumed for all sides.

Figure 3.9:

TEG unit cell subject to thermal loading.

The desired objective is to maximize the temperature gradient by maximizing heat flow through the TEG. The objective function for the quest is shown in Eq. (3.12).

subject to the following constraints: 0≤ρ<1, i.e., the material distribution (design variable) between 0 and 100%

where ρ are the element-wise constant material densities in element, e = 1, … , N.

Maximizing the heat flow through the TEG legs by varying volume fraction φ, subject to the equilibrium constant [K_thT – F_th = 0] and material volume constraint ∑e=1Nveρe≤V,0<ρmin≤ρe≤1, when e = 1, … , N are used in the analysis; and its conductivity captured in Eq. (3.10).

Cases examined

Optimal topologies were developed for varying volume fraction, as shown in Figure 3.10 below. Figure 3.7 a–c show the optimal topologies for respectively 100, 80, and 60% volume fractions.

Figure 3.10:

Optimal topologies for volume fractions of (a) 100%, (b) 80%, and (c) 60%.

The heat conductance patterns mirror Bejan’s Constructal Law (Bejan 2005, 2015; Bejan and Lorente 2008a; Bejan, Lorente, and Lee 2008b). These researchers argued that for a flow system to persist in time it must evolve in such a way that it provides easier access to its currents, as supported by Bejan’s Constructal Law of tree like flow patterns. Further, in his observation, Bejan suggested that nature is the law of configuration generation, or the law of design. Where others saw disorder, Bejan found a pattern. “One of the basic aims of design is to get the most amount of work done with the least amount of effort. What engineers try to accomplish through their plans, nature achieves on its own. Zane (2007) wrote about going with the flow and gave examples such as animate (people, trees) and inanimate (rivers, mud cracks) phenomenon, if given freedom, will organize themselves into patterns or shapes that help them get from here to there in an easier manner”.

TEG heat conversion to power

Three phenomena control the power generation within a TEG; namely the Seebeck effect, Joule effect, and Thomson effect.

The Seebeck effect describes how a temperature difference between two dissimilar conductors produces electromotive force between them. The p-type leg contains positive charge carrier (holes) while the n-type leg is a negative charge carrier (electrons). If the hot and cold junctions are maintained at different temperatures T_h and T_c and T_h > T_c, an open circuit electromotive force, V, is developed at the resistance load, R_L. This action is represented in Eq. (3.15):

where V is electromotive force or voltage, α is Seebeck coefficient, T_h is temperature at the heat source, and T_c is temperature at the sink.

The electromotive force is a catalyst for the thermally excited electrons to move. Although the temperature distribution along the legs is the same, the heat flux is higher in the n-type leg. Different charge carriers in the legs cause unique rate of vibrations resulting to current flow.

Joule heating occurs when current flows through a material offering resistance to the flow. The amount of heat produced is represented as shown in Eq. (3.16).

where Q_j is joule heating, I is current, and R is electrical resistance.

The Thomson effect relates to the rate of generation of reversible heat resulting from the flow of current along an individual conductor along which there is a temperature difference (Rowe 2006; Ruiz-Ortega and Olivares-Robles 2018; Saqr and Musa 2009). The Thomson effect is defined by Eq. (3.17).

where β is the Thomson coefficient, Q_t is the rate of reversible heat absorption, I is the current flow, and ΔT is the temperature difference between the ends of the conductor. This effect is normally minimal.

The performance of the thermoelectric generator is greatly influenced by the figure of merit (Z) depicted in Eq. (3.18) below.

Here α²σ is the electrical power factor and k is the thermal conductivity. When multiplied with temperature, Z becomes unitless and is represented as ZT, a dimensionless figure of merit. A higher ZT is associated with higher TEG power generation. Further, in order to obtain maximum figure of merit (Kanimba and Tian 2016), the geometry and material properties for the TEG should satisfy Eq. (3.19) (Chen, Meng, and Sun 2012; Kanimba and Tian 2016; Rowe 2006).

where A_p is cross sectional area positive leg, A_n is cross sectional area negative leg, L_n is length negative leg, L_p is length positive leg, k_n is thermal conductivity negative leg, k_p is thermal conductivity positive leg, ρ_p is electrical resistivity positive leg, and ρ_n is electrical resistivity negative leg.

The continuity equation for the steady-state current flow in a thermoelectric material can be represented as

where ∇⇀ is the differential operator with respect to length and j⇀ is the current density. The current density j⇀ and temperature gradient ∇⇀T affect the electric field E⇀. Differentiating Eq. (3.18) with respect to length results to an electric field shown in Eq. (3.21) (Lee 2013).

where ρ is electrical resistivity and α is Seebeck coefficient.

The heat flux is represented as (Eq. (3.22))

The heat diffusion for a steady state is given as (Eq. (3.23))

where q˙ is expressed as shown in Eq. (3.24)

Substituting Eqs. (3.22) and (3.24) into Eq. (3.23), results in Eq. (3.25)

where ∇⇀.(k∇⇀T) is thermal conduction, J2ρ is Joule heating, and TdαdT is the Thomson coefficient (Lee 2013).

Examining one dimension heat transfer

The governing equations given by Eqs.(3.20)–(3.25) simplify significantly with a one-D heat flow assumption, which is appropriate given the assumption of constant cross-sectional area and nearly constant heat flow in the y-direction. Thus, the energy balance given by Figure 3.11 can represent the heat and power flow through a thermoelectric leg with x = L being the hot side and x = 0 depicting the cold side. The schematic below represents the global energy balance on a leg.

Figure 3.11:

Illustration of equations used to evaluate heat conduction through a TEG leg.

A differential energy balance inclusive of the Thomson and Joule heating terms is thus given by Eq. (3.26) (Lee 2013):

Further, using Fourier’s Law and the definitions for Joule and Thomson heating given respectively by Eqs. (3.16) and (3.17), Eq. (3.26) can be expressed as

and re-written as

where T, k, A, μ, eJ, σ, and L are respectively temperature, thermal conductivity, cross section area, Thomson coefficient, electrical current density, electrical conductivity, and length of a TEG leg.

The boundary conditions given by Eqs. (3.29) and (3.30) are applied to Eq. (3.28),

Here the subscripts c, and h represent cold side (heat sink), and hot side (heat source), respectively; and then taking the limit as Δx goes to 0, the following differential Eqs. (3.31) and (3.32) result for the p- and n-type legs.

p-type leg

n-type leg

where I = eJA is the electrical current.

The heat flow through each of the legs is given respectively by Eqs. (3.31) and (3.32) for respectively the hot and cold sides of the TE.

where α is the Seebeck coefficient and represents α_ph – α_nh for x = L, and α_pc – α_nc for x = 0. These equations can be further simplified by using the following relationship for thermal and electrical conductance given by Eqs. (3.35) and (3.36).

Using Eqs. (3.29), (3.30) and (3.33)–(3.36), we can generate Eqs. (3.37) and (3.38)

The equation to determine power is given as

The equation used to calculate the thermal power efficiency of the TEG is

Power output in a thermoelectric generator is a function of electrical current. In this case, the interest is in flexibility of the TEG leg, and the power that could be generated in that state. In addition, although shown in Eqs. (3.37) and (3.40), the Thomson effect is neglected in this thermoelectric computation as it has been identified as having minimal influence (Rowe 2006).

The integrated model is a function of separately running the structural model matrices and thermal/power model matrices with their results normalized, weighted and combined to obtain optimized outcomes (Mativo and Hallinan 2017). The reason to run the two models separately is because of their uneven degrees of freedom (dof) that originate from differences in their discrete elements in the reference domain. The structural element node has the potential to translate (x, y, z) and the ability to rotate (x, y, z). This operation results to six dofs. The thermal element node can only translate (x, y, z) and hence has three dofs. All the boundary conditions stated for both structural and thermal models apply in the integrated model.

Integrated model

The combined model is used to examine what heat flow can be allowed subject to the constraint of meeting applied loads. Design parameters for both structure and thermal are presented in Table 3.8.

The objective function for the combined structural and thermal load is shown as maximizing displacement and maximizing heat flow as represented by:

subject to the following constraints:    0≤ρ<1, i.e., the material distribution (design variable) between 0 and 100%

where ρ is the element-wise constant material densities in element, e = 1, … , N. T_min is the temperature at the heat sink, while T_max is the temperature at the heat source. The u* is the minimum calculated displacement. The stiffness matrix K depends on the ρ. r is the residual in obtaining the structural equilibrium. For topology optimization the equilibrium r = 0 is found using an iterative procedure. u and f are displacement and load vectors, respectively.

The MATLAB tool used for the structural model and thermal models were significantly modified to integrate both models and calculate power. The flow chart illustrated in Figure 3.12 depicts the process. It shows a tool that simultaneously computes a TEG leg topology by distributing available material in the integrated design domain and creating corresponding heat path resulting to a reconfigured leg and power generation, respectively.

Figure 3.12:

Detailed flowchart of integrated TEG model with power extraction.

The MATLAB code for the combined model is divided into five sections. The first main section of the program defines the inputs and starts the distribution of the material evenly in the design domain. An initialization of design variables and physical variables is done. In this case, the design variable is the volume fraction limit of the material that is used to develop optimal topologies for special TE material volume fractions. The physical densities are initially assigned a constant uniform value, but are iteratively updated to form an optimal ability to support both structural and thermal loads. A power-law approach is used to penalize intermediate density values, driving them towards solid or void.

The second code section is the Optimality Criteria (OC) which updates the design variables depending on the used and remaining material. The OC is formulated on the conditions that if constraint 0 ≤ x ≤ 1 is active, then convergence is achieved with a positive move-limit of 0.2 and a damping factor of 0.3, for compliant designs (Bendsoe 1995; Liu and Tovar 2014; Vijayan and Thangavelu 2012). Material spread entirely on an element is considered to have a value of 1, while an element without any material has a value of 0. The iterative process continues with termination occurring when a maximum number of iterations is reached or where a tolerance of the material distribution is relatively small, for example if available material is equal to or less than 0.01 (Liu and Tovar 2014).

In the third section of the code, a mesh-independency filter whose function is to avoid numerical instabilities is applied. The filtered density defines a modified physical density that is now incorporated in the topology optimization formulation. During each iteration, the element sensitivity is modified. This process ensures that elements gain or lose material depending on the critical amount available and reaction loads, based on set filter settings.

The finite element code is the fourth section of the program. Here, two separate global stiffness matrices are simultaneously created by a loop over all elements and executed. A four node bi-linear element enables formation of an 8 × 8 stiffness matrix for structural loading while a 4 × 4 heat conduction matrix is formed for the mono-linear thermal element. The two matrices evaluation simultaneously is the uniqueness of the tool. Both matrices are weighted and normalized for optimized topologies before power extraction.

Finally, the fifth section of the code addresses calculation of the power generation. It should be noted that while the goal is to maximize power generation; the reality is that maximizing heat flux is the equivalent of this.

Results

Numerical experiments were conducted to study the integrated structural and thermal with power generation. Shear loading cases that were used in the structural only models are also applied to the integrated models and the resulting adjoint displacements calculated. The integrated models use optimized structural models and seek to maximize heat flux. In order to achieve an even structural and thermal density in the design domain, the baseline model at 100% volume fraction was used to determine the objective functions for both the structural and thermal models. For the structural model, both the mechanical loads and the constraints were applied resulting to an objective function of 57.43. For the thermal model, both the thermal loads for the heat source and heat sink; and constraints of insulated sides were applied resulting to an objective function of 79.91. Together their total was 137.34 which made the structural density to be at 41.8% and thermal to be 58.2%. To normalize the density to a 50–50 ratio, the density for the structure was multiplied by 1.2 while the thermal was multiplied by 0.86.

Table 3.9 shows the effects of optimization weighting on topologies. Weighting is done between 0 (thermal only) and 1 (structural only). It was observed that the higher the structural weighting, the less the displacement and the increase in power generated.

Table 3.9:

Effects of weighting (Ws) on the combined model.

Percent structural weighting (Ws)		Percent structural weighting (Ws)
0		0.2
0.4		0.6
0.8		1

Table 3.10 shows the optimal configurations for representative results for volume fractions ranging from 100 to 40%. Volume fractions between 30 and 0% have no material to create a leg to allow heat transport. As seen, unique topologies form when structural and thermal loads are simultaneously applied to varying volume fractions. Consistently, most of the heat transport occurs along the street where the most material is found. From the topologies, it is observed that the overall constraint is structural because in the end it provides a path to heat flow. The heat travels to the sink where it is readily distributed. A lower volume fraction yields more displacement at the cost of less power generated. In order to sustain shear load, less material is required, which means that the thermal resistance increases. Thus, for a given temperature difference, there is reduced heat transfer, and thus power generation potential.

Table 3.10:

Displacement, power, and shape of an integrated TEG for the vibratory environment.

Volume fraction	Displace-ment	Power generation at 110 °C	Power generation at 170 °C	Power generation at 230 °C	Shape
100	0.0504	0.1163	0.4467	1.0499
90	0.0555	0.1078	0.4319	0.9725
80	0.0636	0.0992	0.3972	0.8937
70	0.07841	0.0900	0.3604	0.8115
60	0.1031	0.0799	0.3231	0.7264
50	0.1290	0.0704	0.2824	0.6374
40	0.1420	0.0591	0.2362	0.5313

Power generation is favored highly with models that have more material. The highest power generation is shown to be from the model without a void while the least power generation has the least material. Power generated from models with volume fraction of 50% and less are not to be considered because they are not structurally capable of supporting all loads. Comparisons of maximum power generation for the baseline model, and the reconfigured integrated model are shown in Table 3.11.

The aim of this study was to search for TEG geometry that could be used in a vibratory environment. Based on the investigation, several models emerged with different potentials. Tradeoffs are seen in Table 3.11 and in Figure 3.13. Figure 3.13 presents a relationship between displacement and power generation. As can be seen, the less the displacement, the higher the power generation. This is expected because there is more material to carry heat.

Figure 3.13:

Chart display power generation versus displacement.

The less void in the leg, the higher the power and less shear displacement. For a balanced model that power generates and is structurally stable, a 80% volume fraction is recommended for it can allow a displacement of 0.0636 and produce about 80% baseline power.