Development of Compliant Thermoelectric Generators (TEGs) in Aerospace Applications Using Topology Optimization

TEGs Application and Failure

Thermoelectric generators (TEGs) are made to be mounted on flat surfaces free of vibrations. They are ill-suited for non-planar areas that lend themselves to bending stresses. According to Custom Thermoelectric (2011), TEGs are typically mounted under compression not to exceed 200 psi with thermal grease or flexible graphite. Vibrations and non-uniform thermal and/or mechanical loading can also induce unintended stresses (Laird Technologies 2014; Global Thermoelectric 1992). Research by McCarty (2008) indicates that “the rigid nature of TEGs makes them particularly sensitive to principal stresses. Brittle materials under plane-stress conditions will fail if any point within the material experiences principal stresses exceeding the ultimate normal strength of the material as explained by the maximum normal stress theory. Failure starts with a brittle fracture whereby a rapid run of cracks through stressed material occurs. Choi, Seo, and Choi (2011) observe the fact that dislocations mostly occur on the edge of TEGs confirms that the failure is often because of shear stresses induced by thermal expansion. No plastic deformation is visible, and in many cases no special pattern on their fractured surface is observed.

Bismuth Telluride Properties

Bismuth Telluride (Bi₂Te₃) is the most common material for TEG legs, and is used widely in thermoelectric generators at temperatures of up to about 250 °C (Rowe 2006; Tong et al. 2010; Custom Thermoelectric n.d; Thermal Electric Corporation n.d). Its attractive properties include electrical conductivity of 1.1×10⁵ S.m/m², a very low lattice thermal conductivity of 1.20 W/mK, and an electrical resistivity of 10µΩ/m. Rowe (2006) observes that these properties make Bismuth Telluride suitable for attaining a high figure-of- merit (Z), which is the ratio of the electrical power factor (σα ²) and the thermal conductivity (k), where σ is electrical conductivity, and α is Seebeck coefficient. Materials which possess a Z > 0.5 are usually regarded as thermoelectric materials (Rowe 2006). Further, the crystalline structure of Bi₂Te₃ lends to its anisotropic nature, with a density of 7.8587 g/cm³, a modulus of elasticity that ranges from 45 GPa to 62.8 GPa, and Ultimate Tensile Strength (UTS) of 7.4 GPa (Tong et al. 2010). Williams (2010) and Adams (2010) state that as Bi₂Te₃ tensile strength is lower than its compressive strength, it will show brittle behavior. Therefore, legs made of Bismuth Telluride exhibit brittle behavior and result in forming rigid TEGs (McCarty 2008; Rowe 2006). For this reason, TEGs applications have not been feasible for use on non-planar surfaces and those that experience tensile or shear stress.

Search of Flexible TEGs

In previous and current research on TEGs optimization (Hannan et al. 2014; Quan et al. 2013), work has been done on ways to develop flexibility in energy harvesting. Further, attempts to develop flexible TEGs have been made in the recent past. Four types of flexible TEG designs have been developed. These have employed flexible foil structures, wavy-slit technology, carbon nanotubes (CNT), and graphene nanoribbons (GNRs). In their research, Qu, Plötner, and Fischer (2001) determined that flexible foil substrate technology relies on embedding thermo-elements in epoxy. This design is constrained by epoxy thickness. Foil substrates are typically made of flexible epoxy film categorized as thin or thick, with an average 50 µm for thin and about 190 µm for thick. Additionally, thermocouple strips capable of generating voltage have been embedded in the epoxy film. Glatz, Muntwyler, and Hierold (2006) argued that because of their limited thickness, thin film deposited materials have to be laid out laterally rather than vertically, inducing thermal losses through the supporting material and limiting the integration density. They further observed that placing a thermocouple onto a thin membrane reduces thermal losses but does not allow for effective thermal contacting the cold and hot side via the top and bottom surface of the thin device. He and his team therefore suggested and developed a thermoelectric wafer in a 190 µm thick flexible polymer mold formed by photolithographic patterning. Their preliminary efforts, together with Saqr and Musa (2009) led to a proposal of a model of vertical micro thermoelectric generators. Since the epoxy thickness required would be necessarily small, a TEG of this design tends to have low power capacity.

Shiozaki et al. (2004) proposed a flexible thermopile generator with slits (FTGS) to permit application of TEGs to non-planar surfaces. Devices of this nature have thermocouples placed on a polyimide sheet. Each thermocouple is placed at 45^o angle vertically, effectively separating p from n thermocouples. The cold junctions are formed by bending the thermopile sheet to a wavy form. The design by Shiozaki et al. (2004) forms the wavy and slit flexible thermopile generator. Using an approach similar to Shiozaki’s, Lon E. Bell registered patent #6,700,052 B2, in March 2, 2004, in which he claimed “a flexible thermoelectric comprising: a plurality of thermoelectric elements; and first and second substrates sandwiching the plurality of thermoelectric elements and having electrical conductors that interconnect ones of the plurality of thermoelectric elements, wherein at least one of the first and second substrates is constructed of a substantially rigid material, said substrates configured to flex in at least one direction.” One of the challenges facing the wavy-slit design is the reconfiguration of thermal and electrical continuity from floating elements. This design adds weight because of rerouting continuity components.

Carbon Nanotubes (CNTs) pose a means for developing composite thermoelectric devices. CNTs are mechanically strong and light weight. However, their high thermal and electrical conductivity pose a challenge to their integration into thermoelectric devices. Koplow et al. (2008) observed that for highly efficient devices, efficient generators should consist of materials with high Seebeck coefficients to provide significant voltages, low electrical resistivities to minimize internal losses, and low thermal conductivity to minimize heat losses. The solution to CNT use in TEGs lies in doping them in order to effect desired properties. Dragoman, Dragoman, and Plana (2007) observed that the Seebeck coefficient, α, is strongly dependent on the CNT conductance, G, e. g., the transmission coefficient carrier through the CNT. The mobility of CNT is μ=(G)l_fp/Ne, where l_fp is the mean free path of the carriers and Ne is the charge density. Therefore, the Seebeck coefficient and the mobility are related through conductance G. Dragoman, Dragoman, and Plana (2007) established that in CNTs, the mobility decreases with temperature, while α increases rapidly at low temperatures and increases slowly in the 200–300K range.

In Graphene Nanoribbons (GNR) technology, thin strips are increasingly being explored for use in TEGs. Their high electrical and thermal conductances place them in a close category with CNTs. However, a unique difference exists such that thermal conductivity is significantly decreased under tensile strain, but is insensitive to compressive and torsional strains (Wei et al. 2011).

TEGs in Aerospace

Thermal management has increasingly been identified as a significant challenge for both military and civilian aircraft, according to Wissler (2009). Recently, the Defense Advanced Research Projects Agency (DARPA) identified thermal management as a main obstacle to further enhancement of Department of Defense capacity (Cohen 2012; Keller 2013). Martinez (2014) indicates commercial aircraft have also experienced increased heat generation through advancement in avionics, passenger capacity, and passenger service expectations. Although TEGs could be considered as part of thermal management problems, the aircraft environment is unfortunately highly problematic for their application. The brittle nature of TEGs makes them highly susceptible to mechanical failure in the environment within aircraft where they might be employed, such as in the engine, aircraft skin, and landing gear. Further, the rigid structures of these devices restrict their use to flat surfaces. Were non-planar thermoelectric devices capable of being developed, their suitability for aircraft applications would be enhanced considerably. To this end, this study investigates a method of developing more compliant thermoelectric configurations tailored to the unique structural and thermal loadings they might be subjected to in an aircraft environment and be applied in locations where a non-planar configuration is necessary.

Thermoelectric generators could be used to reduce aircraft weight associated with thermal management by reducing the amount of heat that needs to be managed as a result of conversion of thermal to electrical energy; improve system efficiency when paired with solar cells to produce more power; reduce costs such as fuel by engine waste heat harvesting; and reduce carbon emissions as suggested by Aljazeera and agencies (2015), Huang (2009), and Callier (2010). Further, aircraft maintenance, the third highest aircraft expense after labor and fuel, can be significantly reduced if TEGs are used to power aircraft health monitoring sensors (Samson et al. 2010). In addition, Huang (2009) and Fleming, Ng, and Ghamaty (2004) suggest that the micro air-vehicle (MAV) challenge of developing a light weight propulsion system could be enhanced by converting waste heat into useful electricity through the use of a TEG.

Table 1 shows the extreme thermal environments present in various aircraft systems. Heat generation from aircraft avionics, the more electrical aircraft (MEA), and the landing gear, provides opportunities for energy harvesting using TEGs, as the operating temperature of these systems falls within the acceptable range for TEGs (Hufford 2014). In avionics applications, the Institute for Interconnecting and Packaging Electronic Circuits (IPC) suggests some limiting cases for different thermal parameters, which should be maintained in commercial aircraft. For example, the recommended worst case thermal conditions for commercial aircraft are a minimum of −55 °C and a maximum of 95°C. However, the actual thermal profile experienced by avionics systems tends to go beyond these limits. The temperature extremes are much lower than the range specified, with a high temperature of about 55 °C. According to Das (1999) the temperature difference in cycles and the number of cycles are much higher than the limits set by IPC. Choi, Seo, and Choi (2011) and FerroTec Thermoelectric Technical Reference (n.d.) reveal four additional factors that relate to failure rate in thermal cycling include (1) the total number of cycles, (2) the total temperature excursion over the cycle, (3) the upper temperature limit of the cycle, and (4) the rate of temperature change.

Increased electrical power demands in MEA, has established a need for more battery power. A response to such power demand has led to the creation of more powerful batteries such as the lithium-ion cell LVP65-8-402 battery used in the Boeing 787 fleet. According to the National Transportation Safety Board (National Transport Safety Board 2013), a recent testing showed that the heat generated inside the battery during the heaviest current loading condition of a full auxiliary power unit (APU) start could expose a cell to temperatures exceeding the maximum approved operating temperature of the battery (158ºF [70 °C]) without detection by the battery’s monitoring system.

Brake temperatures on aircraft landing gears can reach much higher temperatures. For example, a Boeing 767 fitted with carbon brakes realizes a maximum temperature of 427 °C while the MD-11 is rated at 550 °C as reported in Boeing (1990). High temperatures of this nature used with TEGs can offer significant auxiliary power for aircraft inspections and as supplement power on the ground or in flight. Although the TEGs based upon Bismuth Telluride are limited to a maximum temperature of 250 °C, other TEG materials could work at these temperatures. Thermal Electric Corporation (n.d) provides an example of, Hybrid BiTe – PbTe can operate for temperatures up to 360 °C, SnSe – PbSnTe up to 600 °C, Calcium Manganese (CMO) up to 800 °C, and CMO cascade with Bi₂Te₃ stacked up to 600 °C.

Although the aerospace industry has made inroads using TEGs as its “space battery” for deep space applications (Das 1999), the aircraft sector has lagged behind. The main challenges the aircraft sector faces in the use of TEGs are vibrations and non-planar surfaces. According to Adams (2010), aircraft vibrations at large amplitudes for an extended period of time cause fatigue stress cycles that could in turn induce shear stresses. Such vibrations are often experienced in wings and landing gears, which are both prime candidates for use of TEG technology. Rigid thermoelectric devices placed on such a wing or landing gear will fail due to their inability to withstand such shear stresses (Laird Technologies 2013). In addition, the landing gear incorporates many non-planar components which pose a challenge to clamping rigid TEGs for energy harvesting.

For purposes of reporting type of vibrations on aircraft, Carbaugh et al. (n. d.) suggest two categories of which one is high frequency tactile vibration typically, more than 25 Hz, and the other is low frequency, typically less than 20 Hz. High frequency vibrations could be associated with sound that related to a small-mass component acting on the frame, examples being a loose door, access panel or fairing. On the other hand, low frequency vibrations relate to large-mass components acting on the airframe, examples being the rudder, horizontal stabilizer, or elevator. Table 2 describes the types of vibratory loads present in various aircraft and flight conditions.

Goals

The authors have not identified any similar studies of TEG design using topology optimization techniques. There are countless examples in the literature on compliant mechanisms (Sigmund 2007; Bruns and Tortorelli 2001; and Venkiteswaran, Turkkan, and Su 2017); structural designs (Ming et al. 2017; Wu and Hu 2017; Priya and Inman 2009); thermal designs (Kumar et al. 2013; Koga et al. 2013; LaGrandeur et al. 2006); and even combined structural/thermals designs (Li et al. 1999; Yang and Li 2014; Hanser and Bumgarner 2011), but the application to TEGs appears to be unique. The goals for this research are to: (1) develop a general framework for optimizing the topology of a TEG subject to structural and thermal loads (any combination); and (2) test that framework on several case studies representing conditions applicable to aircraft in an effort to validate the framework. To achieve these goals a tool that will be developed to allow design and development of a compliant minimal weight Bismuth Telluride TEG capable of generating maximum power tailored to the structural and thermal environment it is present in. In order to address the second goal, it is necessary to first establish a design space for thermal and structural loads in aircraft.

Design Domain

The SIMP approach considers a design domain that is discretized into finite elements (Figure 1). Each element is assigned a density amount ω, which is treated as the design variable. The algorithm, driven by a sensitivity filter, determines changes in material distribution element to element. This material distribution is dependent on load and boundary conditions applied to the design domain. An ideal structural solution results with full dense material, ω=1, or voids, ω=0. In thermal design, ω=1 represents conductive material while ω=0.001 represents a non-conductive material.

Figure 1:

Discretized form elements (e).

A SIMP approach precludes that a constant load applied to a material with a high stiffness (K) will result in a minimum displacement (min(u)) which can be written as min(fu). The applied load and displacement is spread across the material and is accurately represented as min(∫(fu)dΩ) where dΩ is material domain. In a basic study, such as illustrated in Figure 1, shows the discretized design domain that experiences small displacements moving to the right due to the imposed shear force on the left. The design domain has a unit thickness depth.

In a structural design, a finite element mesh given by Figure 1 can be used for two structural fields of interest, i. e., stiffness (E) and deflection (u), and if E can be described as a constant in each element, we can write the discrete form of structural optimization as

max fTu

where u and f are displacement and load vectors, respectively. The stiffness K depends on the stiffness E_e in element e, numbered as e = 1, …, N. K can be written as

where K_e is the global element stiffness.

An application to a thermal design introduces modifications to accommodate thermal heating conditions. The material phases have isotropic conductivities of fully conductive material and non-conductive. A heat source and heat sink are added to the design domain.

Constraints

The geometrical domain is assumed to be fixed. Initial material distribution within the fixed geometric domain is assumed to be spread evenly to each element. Imposed initial conditions dictate material distribution within the geometry. Applied mechanical and thermal loads to the geometry with imposed boundary conditions cause new material distribution arrangement. However, the new material arrangement is dependent on the amount of available material and extent to which it should be distributed. The ratio of material amount versus the geometry is known as volume fraction (φ), which ranges from completely filled geometry at 100 % to completely void 0 %. Structural loads can include pin loads or distributed loads. Together with boundary conditions, applied loads determine the size, shapes, and topology of objects. Thermal load sources and sinks determine heat path.

Mathematically, these constraints are expressed as:

1. v ≤ V , resulting volume (v) is less or equal to the original design domain volume (V) within the design domain;

2. Applied structural loading accounting for shear and compressive loads at x = 0, (f) f=Ku, is equivalent to the stiffness (K) and displacement (u) of a material; or in a reaction it could be formulated as R=Ku−F ; and

Tool and Tool Validation

A standard topology optimization tool that includes an optimizer and Finite Element subroutine was created in MATLAB to enable building models that would solve linear compliance minimization problems (Bendsoe and Sigmund 2004). The flowchart showing the program is shown as Figure 2.

Figure 2:

Basic flow chart.

The authors modified the tool, as seen in Figure 3, to study more complex cases seeking compliant and integrated solutions. This basic tool is divided into four parts namely, a main program; an Optimality Criteria based optimizer; a mesh-independency filter; and a finite element subroutine. Where the nomenclature is:

Figure 3:

Detailed flowchart of integrated TEG model.

lk – global stiffness matrix for thermal

lkk – global stiffness matrix for structural

FEtn – Finite Element Analysis for thermal

FEs – Finite element analysis for structural

Utn – global temperature vector

Us – global displacement vector

u – displacement

c – optimization of the objective function

dc – derivative of c

oc – optimality criteria subroutines looks for passive densities and sets their density to minimum 0.001

check – sensitivity analysis by a call to mesh-independency filter

change – main loop is terminated if change in design variable is less than 1 %

According to Bendsoe and Sigmund (2004), to search for the optimal lay-out of a structure with a specified region only requires applied loads, support conditions, volume of structure to be constructed and location and size of 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 (Ω). The available space is considered to be the 100 % volume where the designed object is created within and is treated as a continuous variable for a given φ of the design domain (Figure 4). The volume is reduced by removal of material hence creating a minimized shape that will withstand the highest applied load. The result of this process creates the lightest possible structural.

Figure 4:

Different volume fractions correspond to design domain occupied space.

In general the passive elements are set to a density equal to 0.001. To ensure a complete distinction between the object and unused volume, a penalty is imposed to distinguish them. A penalty of 1 makes intermediate densities unfavorable compared to the cost (volume). A penalization power (penal) of ≥ 3 adequately distinguishes the two (Figure 5). The basic density –penalization relationship is ρ^p where ρ is the “density (thickness of material in Ω) and p is the penalty. The penalty has to be equal or greater than 0. For ρ of 0 and p of 10 the result is zero, meaning void. For ρ of 0.5 and p of 10 the result is 0.00098 meaning very little material. For ρ of 1 and p of 10 the result is 1 meaning full thickness. Recommended penalty is given to be 3.

Figure 5:

Penalization powers.

Finally, a filter size (rmin) divided by the element size (more like a mesh size) is administered, typically between 1 and 1.6, where 1 is coarse filter and 1.6 is high filtering, as shown in Figure 6. High filtering takes considerable computing time but provides better outcome.

Figure 6:

Filter sizes.

The basic MATLAB tool (Figure 2) starts in its main program by distributing the material in the design domain evenly. After initialization, it starts to call to finite element subroutine which returns the displacement vector. The element stiffness matrix is called only once as the stiffness matrix for solid material is the same for all elements. The objective function follows and loads are applied. The sensitivity analysis follows the global displacement vector after which mesh-independency filter and the optimizer follow, respectively. The material volume is a monotonically decreasing function of the Lagrange multiplier. The material volume can be found by bi-sectioning algorithm by guessing a lower and upper bound for the Lagrangian multiplier. This is done until its size is less than the convergence criteria. Only those square elements with a side length of two times round (rmin) around are considered elements and a searched by the filter. A filter set less than one in the call routine is equal to the original sensitivities and makes it inactive.

The finite element code makes the 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. Fixed and free degrees of freedom are easily entered in the code. Young modulus and Poisson’s ratio can be entered in the code as well.

Effects of Constraints

Limitations and restrictions are constraints imposed on the design domain. Support conditions such as fixed elements or nodes in the domain may restrict displacement or transmission of loads. Figure 7 shows a point load applied at two different positions of design domain and support conditions (fixed nodes) imposed on the bottom of both left and right of the design domain.

Figure 7:

Effects of points of load application and support restrictions.

As can be seen in Figure 7, all variables are held constant except load application locations and the results show altered topology for load transmission.

Unit Cell Design

A unit cell design is established that is comprised of two elements (legs), a street, heat spreader, and a heat sink (Figure 8). In application one or many more unit cells would be present. The unit cell shown in this figure considers mechanical boundary conditions that include the applied compressive loads from the top, and shear loads from the side on the upper end of the left leg, and fixed (zero deflection) or anchored surface, at the bottom surface. Additionally, the geometrical domain is assumed fixed. The thermal loading applied on each TEG leg is subjected to a distributed heating on the top surface. The sides are assumed insulated. The temperature of the bottom surface is set to a constant value, controlled by the heat sink to which heat is being rejected.

Figure 8:

A general unit cell for a TEG subject to mechanical and thermal loadings.

The configuration shown in Figure 8 could represent a situation like an aircraft skin, with heat input from the fuel (used for heat exchange to the external environment) or a number of other combined thermal and structurally loaded environments. In this unit cell, S₁ and S₂ are the unknown leg shape, size, with topology of thermoelectric p-type and n-type legs (e. g., contiguous p- and n-type materials) separated by an insulated space between them known as the street. The S₁ and S₂ legs carry both structural and thermal loads. In this figure they are shown as rectangular elements, however, here they are intended only to represent the geometrical domain where a composite TEG system can be developed. The street serves as insulation space between the legs and is configured as fixed void space in the unit cell design model. To counter the applied shear load indicated in Figure 8, an adjoint load is introduced to the right of the model. The unit cell model design is complex and requires a compliant mechanism approach.

Compliant Mechanisms

Compliant mechanisms attain their mobility from flexibility of their constituents rather than the use of hinges, bearings, and sliders as their rigid body counterparts (Bendsoe and Sigmund 2004). Nishiwaki et al. (1988), Bruns and Tortorelli (2001), and Lopes and Novotny (2016); discuss several approaches to designing compliant mechanism with constraints. Further, 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 (Figures 9 and 10) are provided to validate the ability of the tool to create compliant mechanisms both symmetric and non-symmetric.

Figure 9:

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

Figure 10:

(a & b) Design domain with boundary conditions (a) and output (b).

The input actuator is assumed to be a linear strain based thus modelled by a spring with a stiffness K_in and a force f_in. Strain actuated principles are found in several systems such as, electro-thermal heating, and shape memory, which are characterized by their blocking force f_in and their displacement u_in. Figure 9 only ensures that the component of the output displacement along the desired direction is maximized.

The example in Figure 10 portrays the ability to control output displacement which can be achieved through a mix of boundary conditions and forces.

The constraints in Figure 10 manage to control the output displacement as desired. This is possible by a maximum cross-sensitivity between two perpendicular output directions.

Structural

The goal of developing a compliant mechanism, such as depicted in Figure 11, is to maximize displacement (u_out) performed on a work-piece modeled by a spring with stiffness k_out. A low k_out will allow a large deflection and vice versa. Of particular interest in this study is to allow a horizontal displacement.

Figure 11:

A depiction of a TEG complaint mechanism for a large deflection.

Assuming the material domain is isotropic and assuming the input actuator is a linear strain, it can be modeled by a spring with a stiffness k_in and a force f_in.

The optimization problem for this system can be written as

max

where r is the residual in obtaining the structural equilibrium. For topology optimization the equilibrium r=0 is found using an iterative procedure. For this case, r is the finite element residual for the analysis problem with the applied load f_in.

Thermal

A two-dimensional steady-state heat conduction model is assumed for the unit cell is employed.

The purpose of integrating thermal and structural is to generate power. Power is a function of TEG legs, figure of merit (ZT), thermal conductivity, and electrical resistivity. In traditional thermoelectric designs, TEG legs have constant cross-sectional which lend them to be applied into thermoelectric equations easily. The search for optimized TEG legs produced uneven cross-sectioned TEG legs, hence the traditional thermoelectric equations are not useful in this case. Therefore, this paper will not be considering electrical energy extraction in the TE. Electrical energy extraction will be considered in a future paper.

Structural

A Marlow TG12-6 TEG was used as a baseline dimension for the design and study of the unit cell. Therefore, a unit cell model of dimensionally proportional design domain of 60 units in the x- direction (nelx) and 20 units in the y-direction (nely) was created. A fixed solid region and a fixed void region were created within the design domain. The fixed solid region was created to depict the top and bottom framework regions that would allow both mechanical and thermal loads to be applied, and also a place that served as an anchor to the unit cell. The fixed solid region occupied a total of 0.2 space in of the vertical dimension equally divided for the top and bottom ends. The size of the fixed solid region was optimized to carry a shear load. The fixed void region was set in the middle of the design domain and acted as the street or insulator.

Baseline: Design Domain Capacity

Material available for distribution is expected to cover the entire design domain including the fixed void regions if present. However, in cases with fixed void regions, total available material for distribution becomes excess with only an adequate fraction of it disbursed to the solid region the domain. The fixed solid regions get filled first, and then the rest of the solid regions get the remaining material.

As a baseline for the case considered, 100 % material for distribution was available for the 60×20 design domain. A fixed void region of 20×16 (26.7 %) precluded that a maximum of 73.3 % of available material was to be used to fill the entire remaining domain. Further, it is worth noting that the framework in the design domain formed a fixed solid region of 2×2×60 (27 % of solid space) required a φ of 0.20 material.

Numerical experiments for structural (mechanical) loading only were conducted to examine potential topologies that yielded large displacements. A compression load was distributed across the top of the model while a shear load was applied at the top left side of the model (Figure 6). Results indicated large displacements were obtained within φ range from 0.588 to 0.590. Figure 12 depicts selected topologies from these ranges of φ and displacement (u).

Figure 12:

Selected topologies from baseline model.

Topologies between of φ between 0.580 and 0.586, and 0.592 and 0.594, appeared to have partially fragmented right legs that seemed to affect their displacements. The effect is shown on the line graph in Figure 13 as a sharp change to compliance reduction.

Figure 13:

Highest displacement area of the baseline model.

An analysis to establish peak displacement indicated its occurrence right before a φ of 0.6, and the results are shown in Figure 13. Initial calculation indicated that a φ of less than 0.58 did not provide structural stability and therefore the study pegged its starting analysis from this point to a φ of 0.64 as the top end which resulted in minimal displacement.

Displacement, Space ratios, and Volume Fractions

A relationship between volume fractions, displacement, and space ratios exist. Volume fraction (φ) is the amount of material distribution availed for the design domain; displacement refers to a shift caused by both the compression and shear loads applied; and space ratios are areas in the design domain occupied by fixed void regions and solid regions. Figure 13 shows all three variables in the study.

Displacements shown in the line graph in Figure 14 indicate the effects of volume fractions and space ratios indicate that a higher displacement is experienced by models with less solid material than those with more material, for example, basing from a baseline (bl) of Marlow TEG 12-6 configuration, a model with 25 % less solid (ls) has a displacement greater than 0.45 while a model with a 25 % more solid (ms) has a displacement less than 0.1.

Figure 14:

Volume fraction; displacement; and ratio of fixed void vs. solid.

A baseline for space ratios was based on the Marlow TG12-6 which saw equal thirds of the design domain framework designating two solid regions and a fixed void region. Also, the graph shows a general high displacement between a φ of 0.5 and 0.6 for all models except for the one with a 25 % more solid. A minimum φ of 0.2 was and a maximum of 0.9 were used in the study exploration. The minimum φ would fill the top and bottom fixed solid regions. A maximum φ of 0.733 would completely fill the solid spaces of the model. As observed from Figure 14 the line graphs after a φ of 0.733 are basically flat indication no movement at all.

Three main observations from the study were:

1. The space ratio affected the displacement greatly. A higher volumetric fraction, meaning higher density, led to a higher ratio of solid to void region and resulted to a less displacement. This observation translated that more material caused a higher degree of stiffness.

2. A lower ratio of solid region required a minimum φ of 2 to fill solid fixed regions. The most φ used in the numerical experiments with this configuration was 0.638 which resulted to minimum displacement.

3. The baseline model showed a better balance of space ratio, u and φ than the model with 25 % more solid.

Thermal

A thermal study to validate the tool was investigated using compliance minimization to heat conduction and the thermal flow was studied. First, the study was set to understand thermal flow effects and patterns. This information would aid in learning of about integrated structures that had both mechanical and thermal loads. Second, the study would help establish relationships between volume fraction (φ), thermal conductivity (k), and displacement (u). The following aspects of material and geometry were made:

1. Material phases include 1 for fully conductive and 0.001 for insulation.

2. The plate in compliance minimization study was evenly heated (constant heat source in all nodes)

3. TEG legs were heated at the top with a heat sink placed at the bottom end of each leg

Heat Flow Validation

Numerical experiments on topology optimization of heat conduction with two different sink locations are shown in Figures 15a and b. The first topology’s sink is located at the center left while the second topology had its sink across the bottom. Heat flow show distinct patterns corresponding to where heat sinks are placed.

Figure 15:

(a) Heat conduction with sink at center right, and (b) With sink at the bottom edge.

Close observations of the heat flow optimized topologies for heat conduction suggest tree-like patterns. It should be noted that the seemingly void spaces are not empty, they indicate locations where heat flow occurs first before heading to the trunk and finally passing through the base (sink). This heat conductance flow pattern would be continuous in a steady-state condition. While the tool works for thermal compliance minimization due to heat generation on the domain modeled as a heat plate (Figure 15), it required modification to be adapted to heat flow in thermoelectric leg transmitting heat from source at the top of the leg to its sink at the bottom. Illustration of this modification is shown in Figure 16.

Figure 16:

Topology of thermal flow from top heat to sink at bottom.

Boundary Conditions for Figure 16 include a heat source at the top, heat sink at the bottom, and insulation on the left and right sides. These boundary conditions are similar to the TEG Model.

TEG Thermal Baseline

A review of thermal flow in a Marlow T12_6 TEG geometry with 100 % volume fraction (φ) is shown in Figure 17a. Illustrations of TEG geometries with φ’s of 80 % and 60 % are shown in b and c. Boundary conditions for the model include distributed heat source at the top and sink at the bottom edge, insulation at both outside sides and both inside sides of the TEG unit cell.

Figure 17:

(a, b, & c) Heat conduction through a TEG unit cell.

The heat conductance patterns favor Bejan’s Constructal law (Bejan 2015; Bejan, Lorente, and Lee 2008; Bejan and Lorente 2008; and Bejan 2005). He 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. Tree like flow patterns in Figures 15 and 16 support Constructal law. 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”. The tree flow system is viewed as a transport system that can shed some light on the thermal flow in a TEG device.

Exploring Heat Conductance Topologies

Heat conductance numerical experiments were conducted to determine relationship between thermal conductivity and heat flux when holding volume fraction constant and varying thermal legs. First, a baseline unit cell with equal size of legs and street was used, followed by cutting the street size in half and concluding by doubling the street size. A volume fraction (φ) of 0.5 was used in the numerical experiments. The outcomes of the numerical experiments shown in the cluster of Figure 18 indicate flow patterns generated. Heat flux is shown to increase with a minimized element width. This phenomenon is acceptable because all variables being equal, the rate of heat flow will increase with a reduced element size. In addition, it is evident that heat conductance flow patterns show no change regardless of conductance levels. However, a high thermal conductivity (k) correlates with a higher heat flux (q) as shown.

Figure 18:

a) Heat conductance patterns for a baseline full unit cell model with geometry [20-20-20] where first 20 is for left leg, next 20 is for street, and last 20 is for right leg. (b) Heat conductance patterns for half Street width size [25-10-25]. (c) Heat conductance pattern for 50 % increased Street width size [15 30 15].

Referring to results in Figure 18 a, b, and c, we can see a clear pattern emerging:

1. With all variables held constant, a high k yields a high q

2. With φ held constant, thinner legs produce a higher q

Topologies

The size, shape and optimized TEG legs are very different from the current designs. These new legs are created to provide flexibility while carrying compression, shear, and thermal loads. Traditional legs are optimized for compression and thermal loads only. Though the new legs provide a certain amount of flexibility, they pose a new challenge of narrow links shown in the topologies. Numerical analysis did not indicate this as a problem however; a possibility of accumulation of heat in the narrow links could be a concern. Authors suggest identifying polymers with similar thermal properties to Bismuth Telluride and adding them to the spaces around the optimized legs. This addition would do two things: (1) it will increase stiffness to the leg and reduce the likelihood of breaking, and (2) it will allow even flow of heat from heat source to the sink.