Abstract

This article presents a novel three-dimensional topology optimization framework developed for 3D spatial packaging of interconnected systems using a geometric projection method (GPM). The proposed gradient-based topology optimization method simultaneously optimizes the locations and orientations of system components (or devices) and lengths, diameters, and trajectories of interconnects to reduce the overall system volume within the prescribed 3D design domain. The optimization is subject to geometric and physics-based constraints dictated by various system specifications, suited for a wide range of transportation (aerospace or automotive), heating, ventilation, air-conditioning, and refrigeration, and other complex system applications. The system components and interconnects are represented using 3D parametric shapes such as cubes, cuboids, and cylinders. These objects are then projected onto a three-dimensional finite element mesh using the geometric projection method. Sensitivities are calculated for the objective function (bounding box volume) with various geometric and physics-based (thermal and hydraulic) constraints. Several case studies were performed with different component counts, interconnection topologies, and system boundary conditions and are presented to exhibit the capabilities of the proposed 3D multi-physics spatial packaging optimization framework.

1 Introduction

Effective design automation strategies are key to meeting the present and future needs for 3D physics-based systems involving spatial packaging problems. Systematic, flexible, and efficient design optimization methods are essential for achieving better system performance, compactness, and life-cycle cost, across a wide range of engineering industries. The optimal 3D spatial packaging of interconnected systems with physics interactions (or SPI2) consists of three individual sub-problems: component placement, interconnect routing, and physics-based performance evaluation [17]. Current practices treat these problems independently, relying heavily on packaging and routing design rules [8,9], design heuristics [1012], and expert intuition [13,14], and solution of these design problems often relies upon time-intensive direct human effort. Solving the SPI2 problem requires knowledge of the system requirements for manufacturing, assembly, testing, operation, maintenance, and repair. Often, engineers start with an existing system design and adapt it to meet new requirements. This process is not only time-consuming but also usually leads to designs that are sub-optimal in terms of density, system performance, and reliability. The associated design space involves many interrelated geometric constraints and physics-based considerations, which are cognitively difficult to navigate. Several inefficiencies, including designer bias and human error, lead to inferior designs; these inefficiencies may increase as product requirements grow more complex and required packaging volumes shrink. Therefore, there exists a critical need for design automation methods that support optimal solution of the combined component placement, routing, and physics evaluation problem for systems of increased complexity, and that reduce overall solution expense. Figure 1 shows a 3D representation of one of the many possible solutions to a simple SPI2 packaging problem.

Fig. 1
Three-dimensional design representation of a system with spatially interconnected components that shows one of the many possible solutions to a simple SPI2 packaging problem
Fig. 1
Three-dimensional design representation of a system with spatially interconnected components that shows one of the many possible solutions to a simple SPI2 packaging problem
Close modal

Here are some of the specific challenges that indicate the complexity of the 3D spatial packaging problem from various angles: (1) Both 3D component placement and interconnect routing are NP (nondeterministic polynomial time)-hard problems. Therefore, as the scale and complexity of the system increases, the number of possible solutions explodes combinatorially, increasing decision-making cost significantly. The 3D topological space is vast and challenging to navigate as there can be infinite design options depending on the tuning parameters. Therefore, it is essential to have sampling strategies that can cover the design space thoroughly and efficiently, (2) the 3D spatial packaging problem is a highly nonlinear optimization problem that simultaneously addresses component placement, routing, and physics performance evaluation. Therefore, there is a greater possibility of encountering local solutions with continuous spatial or parameter tuning when compared to design optimization of individual problem elements, (3) one key challenge in using gradient-based solution methods is that changes in interconnect spatial topology may impact the lumped-parameter system models (such as fluid loops) in ways that either prevent simulation of certain designs or at a minimum introduced non-smoothness, (4) creating design representations that can support topology, geometry, and physical aspects of the 3D spatial packaging problem in a unified way is one of the most challenging aspects. Conventional methods address at most a pair of them while solving multi-physics optimization problems. Previous work exists where all three aspects are included but they are specific to their applications and do not generalize, and (5) finally, human expertise solving these packaging problems is useful but can lead to bias or errors when the scale and complexity of the decision-making jas multiple attributes are involved. All these make the spatial packaging problem extremely challenging to tackle with.

Three-dimensional component placement (or packing) [1520] and 3D pipe routing or cable wiring [2140] problems are individually NP-hard problems; solving the packing and routing problem combined with multi-physics interactions and couplings between system elements is thus especially challenging. The 3D component placement, interconnect routing, and physics performance evaluation sub-problems, however, should be solved synchronously to achieve system-optimal designs. A sequential effort, such as place-then-route or vice versa, cannot fully exploit design coupling between all sets of decisions. The challenges grow in significance as system compactness and performance requirements intensify. For example, commercial aircraft engines a few decades ago had larger core diameters compared to current designs; consequently, modern aircraft engines have a much smaller core surface area than prior designs on which are package comparable externals (e.g., wires, pipes ducts, sensors, and actuators).

The authors have made recent advancements in simultaneous 2D spatial component placement and routing (PR) of fluid-thermal systems. A two-stage design framework [4144] was previously developed by the authors utilizing a 2D geometric projection-based topology optimization method [45] for solving SPI2 problems using a continuous design representation. The continuous representation enables the use of gradient-based methods to efficiently search the component packing and routing space. This methodology is centered on the use of simple geometric bars as basis elements to approximately model both the component geometry and routing paths. Two-dimensional bars have favorable geometric properties that can be exploited to represent both the packing and routing problem, and solve them simultaneously. This simple geometric representation was used in two stages to perform physics-based packing and routing for different performance objectives; the staged method ensures identification of a system-optimal solution. The success of this method has been demonstrated for simple 2D PR test cases in Refs. [41,46]. While the authors had success in the 2D design problem, the 3D problem is much more challenging due to the fundamentally different topological design space and many other features such as crossings between interconnects. This work focuses primarily on extending the prior 2D SPI2 design optimization capabilities to 3D systems by utilizing the 3D gradient-based geometric projection techniques [47].

In must be noted that the current scope of this proposed work only focuses on the demonstration of the 3D geometric projection method (GPM) applied through a gradient-based topology optimization (TO) approach. However, a more comprehensive design exploration process would involve efficient enumeration of distinct 3d layout topologies and optimizing each of them using a multi-start method with the proposed 3D GPM approach. This would help explore the design space more thoroughly and get a close system-optimal design solution through exploration of several local minima points. The authors have demonstrated such a comprehensive approach for 2d systems as discussed in Ref. [41], and the 3d topology enumeration algorithms were successfully demonstrated using spatial graphs representations in Ref. [48]. Future work would involve combining both the 3d topology enumeration and the 3D GPM topology optimization steps into a single-integrated design optimization problem.

1.1 Objectives and Contributions.

The main objective of this article is to demonstrate simultaneous multi-physics 3D component placement and interconnect routing optimization using gradient-based topology optimization methods. This work is a novel attempt to explore a geometric projection method for development of a 3D simultaneous packaging and routing algorithm with broad industrial applicability. The ultimate plan is to incorporate all kinds of geometric, topological, and physics-based constraints that support efficient optimization solutions. This initial work is focused on finding solutions that minimize the bounding box volume, including multi-physics capabilities such as thermal and hydraulic models, non-interference geometric constraints, physics boundary conditions, etc. This will pave way to develop a robust design optimization framework suitable for real world 3D SPI2 applications. Here is a list of the major contributions in this article:

  1. We present a novel technique that supports simultaneous multi-physics 3D optimization of component placement and interconnect paths, whereas existing methods treat component layout design and interconnect routing separately (e.g., optimal routing with fixed layout).

  2. Thermal and hydraulic physics-based performance objectives and constraints were incorporated into the optimization problem, subject to various boundary conditions and geometric constraints that prevent interference between components and interconnects. Both 1D lumped parameter and 2D finite element physics models are used within a single optimization problem to support physics-based evaluation.

  3. We use the 3D geometric projection method discussed in Ref. [49], which is an alternative to the well-established SIMP (solid isotropic material with penalization) method [50] for parameterizing and solving the design optimization problem [4].

  4. Design variables include component locations, orientations, and interconnect lengths, diameters, bend radius at joints, and piping trajectories. The physics-based constraints incorporated are maximum allowable device temperatures, maximum fluid flowrate, and pressure head loss within the fluid network.

  5. We demonstrated the effectiveness of the proposed method via the solution of several component placement-routing test cases that utilize multi-physics-based simulations.

The remainder of this article is organized as follows. Section 3 discusses the 3D GPM and the simultaneous component packing and routing (PR)-based optimization problem formulation. The sensitivity analysis performed to attain the first-order derivatives of the objective and constraint functions utilized by the optimization algorithm is derived in Sec. 4. Results from numerical case studies using the proposed optimization method are presented in Sec. 5. Finally, the conclusions and future work items are discussed in Sec. 6.

2 Geometric and Physics-Based Modeling

2.1 Geometric Projection.

Material distribution topology optimization (MDTO) methods optimally distribute material within a prescribed design domain to achieve optimal system performance, for example, structural compliance. Significant advancements have been made in the development of density-based topology optimization methods where the design space is parameterized by a discretized material density field. Here the local material properties are defined as differentiable functions of material density. The advantages of these MDTO methods are: (1) the ability to circumvent re-meshing due to changes in topology and shape of the structure during optimization and (2) support implementation of efficient gradient-based nonlinear programming algorithms by allowing computation of design sensitivities for optimization. Two examples of density-based TO methods are: (1) SIMP technique [51], which generates a solid-void design by penalization of intermediate density regions within the domain during the optimization process; and (2) level set TO methods [52,53] where the structural design is represented implicitly via the zero level set of a function. Level set TO methods do not require penalization of intermediate material densities as the level set representation inherently induces a solid-void design almost everywhere.

The above-mentioned density-based optimization approaches such as the level set TO and the SIMP-based techniques are quite popular as they generate complex, organic shapes that can be fabricated using additive manufacturing processes to generate structures that closely resemble the optimal structural topologies. However, in applications such as the proposed SPI2 design research where the designer chooses to utilize stock materials such as tubes, square or circular plates, bars, etc., it becomes very challenging to directly use the TO results from aforementioned methods due to limitations with the manufacturing processes. In such cases, the designer spends a significant amount of time in post-processing to achieve a design that simultaneously captures the intent of the result from the TO method and can be manufactured within the available processes. This in turn might lead to addition of weight on the structure and degradation of the overall structural performance. An alternative to this is to use high-level explicit design representations such as the GPM that renders designs based on solid geometric components. Thus, the GPM provides several advantages over conventional TO methods to represent complex geometries easily and optimize the locations and orientations of these objects with respect to multi-physics and thereby enabling easy post-processing for additive manufacturing such solutions. The GPM method is further described below.

In this GPM method based on geometric primitives, each component such as a bar (in 2D) or a plate (3D) is analytically represented by its location, dimension, orientation, and curvature radius. The GPM-based topology optimization method was developed further by Norato et al. [45] to generate designs that are manufacturable using conventional processes. This is done by projecting a discrete number of geometric shapes onto a continuous density field defined over a uniform finite element mesh. This enables the algorithms to determine the density distribution of the layout domain and model the physics (thermal) interaction between the different layout components (or devices). A routing segment is represented as a cylinder with curved edges and the density distribution of each elements is a function of the signed distance between the element and the routing bar. The signed distance function from an element centroid p to the boundary of the routing segment q is given as
(1)
where rb is the bar radius, and db is the distance from the element centroid to the medial axis of the bar. Using the notation in Fig. 2, the distance is given by
(2)
where xα/β=xαxβ, and
(3)
(4)
where ab is a unit vector along the bar’s medial axis:
(5)
(6)
Hence, the projected density distribution is given by
(7)
where H(ϕb/r) is the fraction of the mesh element that is contained within the bar:
(8)
where x is a substitute variable for the expression ϕb/r.
Fig. 2
Figure shows (a) the device and bar elements used to represent the components and interconnect routing, respectively, and (b) the geometric projection and parameters of the device and bar elements
Fig. 2
Figure shows (a) the device and bar elements used to represent the components and interconnect routing, respectively, and (b) the geometric projection and parameters of the device and bar elements
Close modal

2.2 Temperature Model.

The thermal interactions between the devices are modeled using a combination of finite element models and geometric projection. The system’s thermal model is assumed to be operating in steady-state conditions, and the thermal conductivity of all components and interconnects is assumed to be the same. Additionally, it is assumed that the temperature changes do not significantly affect the material’s properties of the components and the interconnects. The governing differential equation (GDE) for the heat conduction problem can therefore be expressed as
(9)
with boundary conditions:
(10)
(11)
where T(x) is the temperature distribution in the domain, κ is the thermal conductivity, Q is the heat generation per unit volume in the domain, q is the heat flux, n is the outward normal to the boundary, while ∂ΩT and ∂Ωq are the prescribed temperature boundary and heat flux boundary, respectively. This is considered the strong form of the GDE and we can obtain the weak form of the GDE using Galerkin’s weighted residual method by multiplying the weighting function and integrating it over the domain:
(12)
After integrating the above expression by parts, applying heat flux boundary condition, and discretizing the domain, the weak form of the GDE can be represented as a linear problem. The linear problem can be represented as
(13)
where T is the discretized temperature solution field and K and P are the global stiffness matrix and load vector, respectively. The element stiffness matrix [K] and the load vector {P} are given by
(14)
(15)
The second term in Eqs. (14) and (15) represents the contribution of heat convection along the boundaries of the domain to the stiffness matrix and the load vector, respectively. The first term in Eq. (14) represents the heat conduction term in the stiffness matrix, and the first term in Eq. (15) represents the heat generation term in the load vector. To model heat conduction only between the devices but not throughout the domain, we combine the finite element analysis (FEA) formulation with geometric projection by multiplying the elemental stiffness matrix and load vectors by a penalized element density. Therefore, the element stiffness matrix and load vector are represented as follows:
(16)
(17)
where ρe is the elemental density after projecting the entire layout, ρde is the element density after projecting the dth device to model its specific heat generation, and ρmin is used in the stiffness matrix to avoid numerical instability occurring from zero density regions.

2.3 Lumped-Parameter Pipe Flow Model.

A lumped-parameter model is developed to model the pressure throughout the flow in the routing pipes. This model is used as a computationally inexpensive alternative to a computational fluid dynamics model. The flow of fluids in a pipe can be described as
(18)
where Z1 and Z2 are the heights of the pipe inlet and outlet. The following assumptions are made in the lumped-parameter pipe flow model to simplify the problem:
  1. Flow is incompressible.

  2. All components are connected in series with no branches.

  3. There is no pump or turbine.

  4. The inlet and outlet cross-sectional area of each pipe are the same.

  5. Flowrate at the inlet is known.

  6. Specific weight is non-changing.

  7. Flow is turbulent everywhere.

Based on these assumptions, Eq. (18) equation can be simplified as
(19)
This approximate model indicates that head loss is proportional to pressure loss in a constant area flow pipe. Two principal sources of head loss are considered in this model: (1) losses due to surface friction in a straight pipe and (2) losses due to elbow fittings at pipe joints. The expression for induced turbulence head loss is given as
(20)
where K is the resistance or loss coefficient. Combining this above expression and the identity A1V1 = A2V2 with Eq. (19), we obtain the pressure head loss as
(21)
Pressure losses in an incompressible flow are additive; therefore, we can combine the individual pressure loss of each component, pipe, and elbow joint to obtain the overall pressure loss in the system, which is obtained as follows:
(22)
The loss coefficient for a straight segment of pipe with length li and diameter di is
(23)
and for an elbow with bend angle αi and bend radius ri, the loss coefficient is
(24)
Friction factor fi is a measure of resistance to flow in the pipe. Reference [54] defines the Blasius expression for friction factor as a function of the Reynolds number as
(25)
where the Reynolds number is
(26)
and ρm is mass density and μ is fluid viscosity.
(27)
Here θ is the angle at the joint between two connecting pipes. Let a and b represent the vectors representing the straight lines subtending an angle α as shown in Fig. 3 such that
(28)
(29)
The angle θ can be expressed as
(30)
where
(31)
A tolerance εc is used to avoid an indeterminate derivative at v = 1 [42]. The clipped length lc is given by
(32)
Fig. 3
Pipe elbow geometry
Fig. 3
Pipe elbow geometry
Close modal

3 Problem Formulation With the Three-Dimensional Geometric Projection Method

This study leverages the 3D GPM developed by Smith and Norato [47] to develop a 3D spatial layout and optimization algorithm. Sensitivity of the projection is provided to gradient-based optimization methods. The devices and connectors for industrial applications are represented by simple shapes such as cubes of fixed dimensions and 3D bars with rounded ends, respectively. Each bar is assigned three parameters: segment start and end points represented as x0 and xf and a thickness represented with parameter w. However, the current framework is based on the work by Jessee et al. [42] where the original formulations were modified to prevent bars from being removed. This way the flow network formed by the bars is preserved. The density of each element is evaluated by projecting the geometric shapes into a mesh. The devices are approximated by simple cubic shapes for this initial design study. Each device is defined by a reference point cd, rotation angles about the three coordinate axes from its center (θx,θy,θz), and a set of vectors bi pointing from the reference point to the vertices. The densities of the devices are calculated using the formulation as described in Ref. [42]. The center point and rotation angles are used as the device’s placement parameters in the optimization process. In this initial formulation, the component vertices are not considered as they can be expressed as a function of the device center points and rotation angles. It can be represented mathematically as
(33)
where Rθ is a rotation matrix defined by
(34)
The derivative of the rotation matrix with respect to the rotation angle about the x-axis is given by
(35)

The derivative with respect to the other axes can be obtained similarly. In our design representation for this work, the rotation matrix is used to rotate the device about its center by any angle between 180 deg (clockwise) and −180 deg (anticlockwise).

For the three-dimensional packaging and routing design framework, the objective is to minimize the bounding box volume of the optimal layout under certain geometric constraints. The system elements are predominantly defined as components and interconnects. The constraints as shown in Fig. 4 in the 3D system are enumerated below.

  1. Device-device constraint (gdd(x)): This constraint ensures that the devices defined within the system do not collide with one another.

  2. Device-interconnect constraint (gds(x)): This constraint ensures that the devices and the interconnects do not intersect one another.

  3. Interconnect-interconnect constraint (gss(x)): This constraint prevents the interconnects from intersecting one another.

  4. Minimum bar length constraint (gl(x)): This constraint ensures that the optimal interconnect segment lengths are not reduced to a point (or zero length). However, there could be cases where zero-interconnect lengths are common and very useful. For example, they represent direct connection of two devices, such that the outlet of the source device is directly connected to the inlet of the target device. For example, an engine component bolted directly to engine block.

Fig. 4
The 3D non-interference geometric constraints imposed between different system elements are shown in (a) between two devices (gdd), (b) between a device and a routing segment (gds), and (c) between two routing segments (gss). The mathematical constraint equations for each case are provided in Eqs. (38)–(40), respectively.
Fig. 4
The 3D non-interference geometric constraints imposed between different system elements are shown in (a) between two devices (gdd), (b) between a device and a routing segment (gds), and (c) between two routing segments (gss). The mathematical constraint equations for each case are provided in Eqs. (38)–(40), respectively.
Close modal

These constraints are provided in the two-dimensional design framework and are modified to be included in the three-dimensional framework. The system layout is shown in Fig. 5.

Fig. 5
Different components (devices and interconnects) in a 3D layout optimization representation
Fig. 5
Different components (devices and interconnects) in a 3D layout optimization representation
Close modal
Along with these constraints, a constraint on minimum interconnect bend angles is included. This helps to ensure that excessively sharp bends to not appear in the solution, even if physical constraints, such as pressure drops, are not included that could drive interconnects toward smooth bends. The objective function here, as mentioned previously, is the bounding box volume of the system design. Mathematically, the objective function and the associated constraints can be expressed as
(36)
The variable x represents the vector of design variables (placement and other continuous design decisions). The objective function that represents the bounding box volume (V) is evaluated using the following formula:
(37)
Here, the subscript Cmax represents the maximum value of the corresponding coordinate axis evaluated from the trial spatial layout of the system, and Cmin represents the corresponding minimum value. C represents the corresponding coordinate axis in 3-dimensions. Please note that there is one limitation with the way this objective function is evaluated. In the proposed approach, the bounding box volume calculations are based on alignment of the length, width, and height with the unit directions of the Cartesian coordinate system. Hence, if the components are diagonally oriented, then a small error occurs due to this calculation. This error is considered to be negligible when compared to the actual volume of the bounding box. In future work, a more accurate evaluation of the bounding box volume function will be implemented. The constraints gdd(x), gl(x), gl(x), and gl(x) are collectively referred to as the geometric constraints. The geometric constraints shown above are the three-dimensional versions of the planar formulations. The constraint between the ith and jth components is given by
(38)
where r is the maximum vertex distance from the device center location. The interference constraint between the ith device and jth routing segment is given as
(39)
where dij is the minimum distance between the device center point and the line segment representing the routing bar. The interference constraint between the ith and jth routing segment is given as
(40)
In this case, dij is the minimum distance between the two line segments representing the routing bars.
The angular constraint places a lower bound on the angle between two adjacent interconnections segments. This helps to smooth sharp acute angles that can otherwise be generated by the layout optimization algorithm, and helps prime the optimizer to generate feasible design solutions. The angular constraint can be written as
(41)
where
(42)
Here, θab represents the angle between two 3D vectors b and a. The parameter θ0 is a constant. The term θab is given by
(43)
where
(44)
(45)
As shown in Fig. 6, xa0 and xaf are the start and end coordinates of interconnection segment “a.” Similarly, xb0 and xbf are the start and end coordinates of interconnection segment “b.”
Fig. 6
Representation of the angle subtended between two interconnection bar segments
Fig. 6
Representation of the angle subtended between two interconnection bar segments
Close modal

The angular constraints were initially used as an alternative to physics constraints as a strategy to reduce interconnect bend sharpness. While physics-based constraints can help drive interconnects toward smoother designs (e.g., to reduce pressure losses due to sharp bends), the underlying combined finite element analysis and lumped-parameter models can be computationally expensive to evaluate. For some early-stage studies it may be beneficial to solve simplified problems where physics does not drive smoothness, but is still assured through direct angular constraints.

In this study, the temperature at the center of the devices is constrained not to exceed a specified critical temperature Tmax. Temperature constraints are represented as
(46)
(47)
where N is the shape functions of the element where the center lies and T is the nodal temperatures from the finite element formulation (KT = P). The pressure losses considered in our lumped-parameter model are losses due to friction between the fluid and inside walls of the routing bars and losses due to bends at the elbows. The head loss constraint is given as
(48)
(49)
where Ki is the loss coefficient. The loss coefficient due to surface friction is given as
(50)
The loss coefficients at elbow bends are expressed as
(51)
where the weight density ρw, friction factor fi, and weight flowrate w are all specified properties of the fluid or system being modeled.

4 Sensitivity Analysis

The goal of performing design sensitivity analysis (DSA) is that it discusses “how” and “how much” changes in the parameters of an optimization problem modify the optimal objective function value and the point where the optimum is attained and how each of the constraints are affected by any of the change in the respective parameters. In reality, the solutions of an optimization problem alone are not sufficient to study how effective the design method is, but performing DSA helps in in-depth study of how the different parameters behave with change in the inputs, boundary conditions, and their consequential effect on the output functions. The optimization algorithm “sequential quadratic programming method” provided in Matlab’s fmincon requires the first-order derivatives of the objective and constraint functions with respect to the design variables to determine the path of descent. These first derivatives are often referred to as sensitivities. The sensitivities of the physics constraints defined in Eq. (36) can be obtained through the formulations provided in Ref. [42]. The sensitivities of the geometric constraints with respect to the design variables are given as
(52)
(53)
(54)
The sensitivities of the angle constraint between two adjacent interconnection bars with respect to one of their expanded design variables (xa0x) are given by Eq. (55). The sensitivities with respect to the other design variables can be found similarly (xa0y,xa0z,xafx,xafy,xafz,xb0x,xb0y,xb0z,xbfx,xbfy,xbfz).
(55)
The derivative of the dot product with the x coordinate bar end point can be written as
(56)
The derivatives with respect to the other coordinate directions (for the y and z directions) and the remaining bar end points can be derived similarly. The derivative of the dot product with respect to the x coordinate of the point where the two bars meet can be computed as follows:
(57)
(58)
Here, the term a × b(2) represents the second component of the vector obtained by evaluating a × b. Similarly, the derivative with respect to the common node point of the two bars can be computed as
(59)
The sensitivity of the temperature constraints is derived using the adjoint method, where the constraints are expressed as
(60)
where the residual vector R is the zero vector such that:
(61)
and ΨT is the adjoint vector and is given by
(62)
This adjoint vector is obtained by setting the total derivative with respect to the state variables y to zero. The state variable vector y is the set of the unknown variables in the finite element problem and they are made up of the temperature at the free nodes, Tf, and the heat flux at nodes where the temperature boundary conditions are prescribed, Pp, such that y = [Tf, Pp]T. The sensitivity can therefore be obtained as
(63)
such that
(64)
(65)
(66)
The derivative of the residual with respect to the extended design variables is evaluated as
(67)
The stiffness matrix K and load vector P have been combined earlier with geometric projected densities; therefore, the total derivative of the temperature constraints is obtained as follows:
(68)
while p is the derivative of the element densities with respect to the device center location, rotation angles, interconnect bar ends, and widths. The sensitivity with respect to the bar ends and width is provided in Ref. [47], while the derivative with respect to the dth device center and rotation angle is given as
(69)
(70)
where
(71)
(72)
The terms ρ~e/xe0 and ρ~e/xef are the derivatives of the density with respect to the device edges, which are represented as bar ends in the formulation and their expression is provided in Ref. [47]. From Eq. (33), vertex coordinates xe0 and xef can be represented as
(73)
where bi is the vector from the device center to xe0 or xef. Their derivatives are therefore computed as
(74)
(75)
The derivative with respect to the other axes of rotation (y and z) can be found similarly. The pressure drop in the system is constrained to not exceed a specified maximum allowed system pressure drop, and the associated constraints are expressed as
(76)
The first derivative is taken with respect to the extended variables and is derived using the chain rule. The sensitivity is obtained as
(77)
The extended variable vector consists of device reference points, device rotation angles, bar end coordinates, widths, and elbow radii. The partial derivative with respect to each design variable is zero, except the bar widths, and it is given by
(78)
The derivative of the pressure drop with respect to the “friction factor” for each pipe elements is
(79)
The derivative of the friction factor in Eq. (77), d Ki/dx', is different for a straight pipe and an elbow section. For a straight section, the design variables are segment end points and segment diameter, xi0, xif, and di, and their sensitivities are given as
(80)
(81)
(82)
(83)
For an elbow section, the design variables are the four end points of connected segments, xa0, xaf, xb0, and xbf, diameter, di, and radius of the elbow, ri. The sensitivity of the elbow loss coefficient with respect to pipe diameter is
(84)
where the partial derivatives are given as
(85)
(86)
The sensitivity of the elbow loss coefficient with respect to radius is calculated as
(87)
The sensitivities obtained with respect to the four segment end points are as follows:
(88)
where the partial derivatives are evaluated as follows:
(89)
The derivatives of angle α are
(90)
where
(91)
The derivatives of the clipped length are
(92)
(93)
where
(94)
The derivatives, dv/dxaf, dv/dxa0, dv/dxbf, dv/dxb0, can then be calculated as
(95)
(96)

5 Results and Discussion

This section presents five numerical case studies that were performed to demonstrate the proposed 3D spatial layout optimization method described in Sec. 3. All reported numerical simulations were performed using a workstation with an Intel Xeon E5-2660 HDD CPU @ 2.00 GHz, 64 GB DDR4-2400 RAM, Windows 10 64-bit, and Matlab 2021a.

The design domain for all the case studies is chosen as a box of size 0.5 × 0.5 × 0.25 cubic units (which equals 6.25 × 10−2 cubic units). This 3D volume is the entire design space that is available for the components and the interconnect network to continuously move during the optimization process. For all case studies, the components are represented as cubes having sides of length 0.06 units. For every case study, the initial design layout has an initial bounding box volume, i.e., the tightest bounding box enclosing the system components and interconnect network evaluated using Eq. (37). For example, in case study 1, the initial design layout shown in Fig. 7(a) has an estimated bounding box volume of 2.379 × 10−2 cubic units, which corresponds to the total initial space occupied by the three components and their interconnect network. The bounding box volume minimization was considered as the optimization objective for all the case studies. The optimization design variables are component center location, orientation, interconnect thickness, length, interconnect elbow bend angle, and bend radius.

Fig. 7
Case study 1: three-component system with geometric constraints—(a) initial layout (volume: 2.379 × 10−2 cubic units) and (b) final optimal layout (volume: 1.044 × 10−4 cubic units). We observe a decrease in total system volume by 99.1%.
Fig. 7
Case study 1: three-component system with geometric constraints—(a) initial layout (volume: 2.379 × 10−2 cubic units) and (b) final optimal layout (volume: 1.044 × 10−4 cubic units). We observe a decrease in total system volume by 99.1%.
Close modal

5.1 Case Study 1: Simple Volume Minimization Using Geometric Constraints.

In this case study, we perform a simple bounding box volume minimization of a three-component (as shown in Fig. 7(a)) and an eight-component layout (as shown in Fig. 8(a)). Only the component rotation, and geometric constraints to ensure non-interference between the components and/or the interconnect network are considered in this study. The final optimal layouts for the three- the eight-component layouts are shown in Figs. 7(b) and 8(b), respectively. The initial and final bounding box volumes for both the layouts are depicted in Figs. 7 and 8. We observe a decrease in total volume by 95.6% and 93.4% for the three- and eight-component layouts, respectively. The corresponding 2D XY plane (or top) views of the initial and final optimal layouts in both the cases are shown in Fig. 9. The computational time taken for the three- and eight-component systems are 219 and 118 min, respectively. The eight-component system converged faster than the three-component system as more components and interconnects helped the system to attain a quicker solution before the constraints are violated as shown in Table 1.

Fig. 8
Case study 1: eight-component system with geometric constraints—(a) initial layout (volume: 2.534 × 10−2 cubic units) and (b) final optimal layout (volume: 1.668 × 10−3 cubic units). We observe a decrease in total bounding box volume by 93.4%.
Fig. 8
Case study 1: eight-component system with geometric constraints—(a) initial layout (volume: 2.534 × 10−2 cubic units) and (b) final optimal layout (volume: 1.668 × 10−3 cubic units). We observe a decrease in total bounding box volume by 93.4%.
Close modal
Fig. 9
Case study 1: (a) three-component system initial layout in the XY plane and (b) its final optimal layout, (c) eight-component system initial layout in the XY plane and (b) its final optimal layout
Fig. 9
Case study 1: (a) three-component system initial layout in the XY plane and (b) its final optimal layout, (c) eight-component system initial layout in the XY plane and (b) its final optimal layout
Close modal
Table 1

Summary of case studies: objectives and constraints

Case studyProblem descriptionConstraints under consideration (including with or without boundary (wall) conditions)
1Simple bounding box volume minimization with rotation and geometric non-interference constraintsComponent rotation, and geometric constraints with no boundary conditions (wall temperature or heat flux constraints)
2Finding optimal system volume with three components using angular constraints (the angle subtended between the interconnect segments is constrained within the given bounds)Angular and geometric non-interference constraints for component layout with no boundary conditions (wall temperature or heat flux constraints)
3Finding optimal system volume with same critical temperatures for components and heat dissipation ratesGeometric non-interference constraints, component rotation, physics-based constraints, same component temperature limits, and heat dissipation rates with boundary conditions (wall temperature or heat flux constraints)
4Finding optimal system volume with different critical temperatures and heat dissipation rates for componentsGeometric non-interference constraints, component rotation, physics-based constraints, varying temperature limits, varying heat dissipation rates with boundary conditions (wall temperature or heat flux constraints)
5Finding optimal volume layout with pressure head loss constraintsGeometric non-interference constraints, component rotation, physics-based constraints, temperature limits, heat dissipation rates, fluid pressure head loss limits with boundary conditions (wall temperature or heat flux constraints)
Case studyProblem descriptionConstraints under consideration (including with or without boundary (wall) conditions)
1Simple bounding box volume minimization with rotation and geometric non-interference constraintsComponent rotation, and geometric constraints with no boundary conditions (wall temperature or heat flux constraints)
2Finding optimal system volume with three components using angular constraints (the angle subtended between the interconnect segments is constrained within the given bounds)Angular and geometric non-interference constraints for component layout with no boundary conditions (wall temperature or heat flux constraints)
3Finding optimal system volume with same critical temperatures for components and heat dissipation ratesGeometric non-interference constraints, component rotation, physics-based constraints, same component temperature limits, and heat dissipation rates with boundary conditions (wall temperature or heat flux constraints)
4Finding optimal system volume with different critical temperatures and heat dissipation rates for componentsGeometric non-interference constraints, component rotation, physics-based constraints, varying temperature limits, varying heat dissipation rates with boundary conditions (wall temperature or heat flux constraints)
5Finding optimal volume layout with pressure head loss constraintsGeometric non-interference constraints, component rotation, physics-based constraints, temperature limits, heat dissipation rates, fluid pressure head loss limits with boundary conditions (wall temperature or heat flux constraints)

5.2 Case Study 2: Optimal Layout With Three Components Using Angular Constraints.

The second case study consists of a three-component layout. The initial design layout is shown in Fig. 10. The initial bounding box volume, i.e., the tightest bounding box enclosing the system components and interconnects, evaluated using Eq. (37) is approximately 2.379 × 10−2 cubic units.

Fig. 10
Case study 2: different views of the initial layout: (a) isometric view and (b) a view in the XY plane
Fig. 10
Case study 2: different views of the initial layout: (a) isometric view and (b) a view in the XY plane
Close modal

The optimization algorithm was run for about 350 iterations. The final layout obtained is shown in Fig. 11. Due to the computational expense associated with running a 3D geometric projection scheme, the optimization was terminated when the objective value approximately plateaued as shown in Fig. 12. The final optimal bounding box volume is approximately 0.018 cubic units. Thus for this simple three-component layout optimization case, an almost 60% reduction in the bounding box volume was achieved. There is not much reduction in the total volume here compared to case study 1 because of the angular constraints applied. This affects the overall bounding box volume as components cannot come closer than they otherwise would due to the restrictions caused by the angle subtended at the joints of their interconnects. Further, density improvement might be possible by increasing the number of segments representing an interconnect to help enrich the design space. It must be noted that the angular constraint might not work well for similar case studies but can be useful in complex case studies where system specifications need pipe bends at joints to be at fixed angles for design and manufacturing purposes. Furthermore, the increase in the bounding box volumes seen in this case study also motivated the inclusion of component rotational constraints and incorporating component orientation within the optimization framework presented here.

Fig. 11
Case study 2: different views of the final layout: (a) isometric view and (b) a view in the XY plane
Fig. 11
Case study 2: different views of the final layout: (a) isometric view and (b) a view in the XY plane
Close modal
Fig. 12
Bounding box volume objective function values over iterations during the optimization
Fig. 12
Bounding box volume objective function values over iterations during the optimization
Close modal

5.3 Case Study 3: Same Critical Temperatures and Heat Dissipation Rates for All Components.

Unlike earlier case studies, in this study we demonstrate the use of physics-based constraints, system boundary conditions, and geometric non-interference constraints. Figure 13 shows the 3D design domain containing a representative eight-component interconnected system with physics-based boundary conditions. We assume there is heat convection on the colder face (on the right at an ambient temperature, 0 °C) and a hotter face in the front of the layout (at a fixed temperature 100 °C). It is also assumed that all the other faces (boundaries) of the box are thermally insulated. The heat transfer coefficient across the colder face boundary is 35.80 W/m3.

Fig. 13
Three-dimensional design domain with physics boundary conditions used for case studies 3, 4, and 5. It is assumed that the hotter face (or boundary) is fixed at 100∘C, and the colder face has heat convection due to the ambient air at 0∘C. The remaining boundaries of the design volume are assumed to be thermally isolated. The heat transfer coefficient for thermal conduction is 35.80 W/m3.
Fig. 13
Three-dimensional design domain with physics boundary conditions used for case studies 3, 4, and 5. It is assumed that the hotter face (or boundary) is fixed at 100∘C, and the colder face has heat convection due to the ambient air at 0∘C. The remaining boundaries of the design volume are assumed to be thermally isolated. The heat transfer coefficient for thermal conduction is 35.80 W/m3.
Close modal

In case study 3, all the components have the same critical temperature of 30 °C and a heat dissipation rate of 3000 W/m3. Figures 14 and 15 show initial and final optimal layouts of the three- and eight-component systems, respectively. Figures 14(b) and 15(b) help visualize how all the components in both the three- and eight-component systems and their interconnect networks have moved toward the colder boundary and away from the hotter boundary during optimization in the final layout. This is a clear evidence that the physics-based optimization is successful as it tries to satisfy all the geometric and temperature constraints during bounding box volume minimization. Total reductions in bounding box volume of 98.8% and 84.99% were realized for the three- and eight-component systems, respectively.

Fig. 14
Case study 3: three-component system with same critical temperatures and heat dissipation rates for all components—(a) initial layout (volume: 2.379 × 10−2 cubic units), (b) final optimal layout (volume: 2.646 × 10−4 cubic units), and (c) closer view of the final layout. We observe an overall decrease in bounding box volume by 98.8%.
Fig. 14
Case study 3: three-component system with same critical temperatures and heat dissipation rates for all components—(a) initial layout (volume: 2.379 × 10−2 cubic units), (b) final optimal layout (volume: 2.646 × 10−4 cubic units), and (c) closer view of the final layout. We observe an overall decrease in bounding box volume by 98.8%.
Close modal
Fig. 15
Case study 3: eight-component system with same critical temperatures and heat dissipation rates for all components—(a) initial layout (volume: 2.534 × 10−2 cubic units), (b) final optimal layout (volume: 3.802 × 10−3 cubic units), and (c) closer view of the final layout. We observe an overall decrease in bounding box volume by 84.99%.
Fig. 15
Case study 3: eight-component system with same critical temperatures and heat dissipation rates for all components—(a) initial layout (volume: 2.534 × 10−2 cubic units), (b) final optimal layout (volume: 3.802 × 10−3 cubic units), and (c) closer view of the final layout. We observe an overall decrease in bounding box volume by 84.99%.
Close modal

5.4 Case Study 4: Different Critical Temperatures and Heat Dissipation Rates for All Components.

In case study 4, all the components have different critical temperatures and heat dissipation rates as shown in Tables 2 and 3 for the three- and eight-component systems, respectively. Figures 16 and 17 show initial and final optimal layouts of the three- and eight-component systems, respectively. Figures 16(b) and 17(b) help visualize how the components with lower critical temperatures and higher heat dissipation rates in both the three- and eight-component systems have moved toward the colder boundary and away from the hotter boundary in the final layout. The final component temperatures for both the layouts are shown in Tables 2 and 3. It can be observed that all of them have satisfied their respective critical temperature limits. Total reductions in bounding box volumes of 97.9% and 83.4% were obtained for the three- and eight-component systems, respectively. The objective function convergence plots for this case study are provided in Fig. 18.

Fig. 16
Case study 4: three-component system with different critical temperatures and heat dissipation rates—(a) initial layout (volume: 2.379 × 10−2 cubic units), (b) final optimal layout (volume: 4.897 × 10−4 cubic units), and (c) closer view of the final layout. Here we observe a reduction of the total bounding box volume by 97.9%.
Fig. 16
Case study 4: three-component system with different critical temperatures and heat dissipation rates—(a) initial layout (volume: 2.379 × 10−2 cubic units), (b) final optimal layout (volume: 4.897 × 10−4 cubic units), and (c) closer view of the final layout. Here we observe a reduction of the total bounding box volume by 97.9%.
Close modal
Fig. 17
Case study 4: eight-component system with different critical temperatures and heat dissipation rates—(a) initial layout (volume: 2.534 × 10−2 cubic units), (b) final optimal layout (volume: 4.217 × 10−3 cubic units), and (c) closer view of the final layout. Here we observe a reduction of the total bounding box volume by 83.4%.
Fig. 17
Case study 4: eight-component system with different critical temperatures and heat dissipation rates—(a) initial layout (volume: 2.534 × 10−2 cubic units), (b) final optimal layout (volume: 4.217 × 10−3 cubic units), and (c) closer view of the final layout. Here we observe a reduction of the total bounding box volume by 83.4%.
Close modal
Fig. 18
(a) Case study 4: three-component layout objective function convergence plot and (b) eight-component layout objective function convergence plot
Fig. 18
(a) Case study 4: three-component layout objective function convergence plot and (b) eight-component layout objective function convergence plot
Close modal
Table 2

Critical temperature constraints, heat dissipation rates, and final temperatures for components in case study 4A

Component (color)Heat dissipation rates (W/m3)Critical temperatures (°C)Final temperature (°C)
1 (blue)60007048.25
2 (cyan)30003023.14
3 (magenta)20002020
Component (color)Heat dissipation rates (W/m3)Critical temperatures (°C)Final temperature (°C)
1 (blue)60007048.25
2 (cyan)30003023.14
3 (magenta)20002020
Table 3

Critical temperature constraints, heat dissipation rates, and final temperatures for components in case study 4B

Component (color)Heat dissipation rates (W/m3)Critical temperatures (°C)Final temperature (°C)
1 (blue)50004013.48
2 (light blue)40005036.7
3 (green)30003015.38
4 (cyan)30003030
5 (magenta)20002017.979
6 (yellow)50005042.507
7 (black)30003028.617
8 (orange)10002017.034
Component (color)Heat dissipation rates (W/m3)Critical temperatures (°C)Final temperature (°C)
1 (blue)50004013.48
2 (light blue)40005036.7
3 (green)30003015.38
4 (cyan)30003030
5 (magenta)20002017.979
6 (yellow)50005042.507
7 (black)30003028.617
8 (orange)10002017.034

5.5 Case Study 5: Pressure Head Loss Constraints.

In this last case study, we have included the pressure head loss as an additional physics-based constraint during optimization. Once again, the three- and eight-component systems were considered here as shown in Figs. 19 and 20. The same critical component temperature constraints and heat dissipation rates are maintained as indicated in case study 4. The maximum head loss constraints are 20 m and 35 m for the three- and the eight-component systems, respectively. Sharper interconnect segment bends create greater pressure head loss in the pipes. Figures 19(b) and 20(b) show the final layouts after optimization. To satisfy the head loss constraint, the interconnects have smoother bends compared to the previous case studies. This suggests that the geometry at the interconnect joints is the dominant contributor to the head loss.

Fig. 19
Case study 5: (a) initial layout (volume: 2.379 × 10−2 cubic units) and (b) final optimal layout (volume: 1.558 × 10−3 cubic units) with a 2.5 m active head loss constraint. Here we observe a reduction of the total bounding box volume by 93.45%, and the pressure head loss is an active constraint at 20 m.
Fig. 19
Case study 5: (a) initial layout (volume: 2.379 × 10−2 cubic units) and (b) final optimal layout (volume: 1.558 × 10−3 cubic units) with a 2.5 m active head loss constraint. Here we observe a reduction of the total bounding box volume by 93.45%, and the pressure head loss is an active constraint at 20 m.
Close modal
Fig. 20
Case study 5: (a) initial layout (volume: 2.534 × 10−2 cubic units) and (b) final optimal layout (volume: 3.612 × 10−3 cubic units) with a 5 m active head loss constraint. Here we observe a reduction of the total bounding box volume by 85.74%, and the pressure head loss is an active constraint at 35 m.
Fig. 20
Case study 5: (a) initial layout (volume: 2.534 × 10−2 cubic units) and (b) final optimal layout (volume: 3.612 × 10−3 cubic units) with a 5 m active head loss constraint. Here we observe a reduction of the total bounding box volume by 85.74%, and the pressure head loss is an active constraint at 35 m.
Close modal

Total reductions in bounding box volumes of 93.95% and 85.74% were obtained for the three- and eight-component systems, respectively. To satisfy device temperature constraints, some of the routing interconnects touch the convection boundary. This conducts heat from the devices through the routing to the boundary where it can be dissipated. The optimization finds a balance between smooth bends and reducing pipe length to reduce head loss in a way that is best for system performance. However, the head loss constraints in this case study are active as there exists a trade-off between pressure head loss and volume minimization. For example, if the volume of the system has to decrease, the components should come closer which then requires that the interconnect elbows are sharper. Hence, the minimum bounding box volume is obtained when the head loss constraint remains active.

6 Conclusions and Future Work

A novel method for simultaneously optimizing the 3D placement of components (or devices) and their corresponding interconnect routing paths was presented in this work. The aim is to improve the spatial packaging density and physics-based system performance of a given 3D interconnected engineering system by minimizing the overall bounding box volume, satisfying spatial, physics-based (thermal, hydraulic, etc.), component orientation, and other geometric constraints that prevent interference between devices and/or the interconnect paths. A geometric projection of the various 3D system elements, as well as required sensitivity calculations, were performed that support optimization based on a 3D finite element mesh. A set of design variables was developed that links directly to the geometric functions of the system that represent the volume objective and the different types of constraints.

The proposed method was demonstrated using several case studies to perform simultaneous component packing and routing optimization with spatial, geometric, and physics-based constraints. In addition, component rotation was also demonstrated during the continuous gradient-based optimization process. This improved the spatial packaging density as components could reorient themselves during the optimization to attain more compact solutions. In case study 1, a simple bounding box volume minimization problem was solved with three-component and eight-component systems. Only the component rotation and the geometric constraints for non-interference between components and interconnects were considered for this case study. Case study 2 was performed to minimize the overall volume of a three-component system subject to geometric, rotational, and angular constraints as well. The angular constraints ensure that a specified minimum angle is subtended at interconnect pipe joints, and could be useful for spatially packaging complex 3D systems that must satisfy specific design and manufacturability constraints for pipe bends. Case studies 3, 4, and 5 include physics-based constraints and boundary conditions as part of the optimization problem. In case study 3, all the components have same critical temperatures and heat dissipation rates. In case study 4, each component has a unique critical temperature constraint and a heat dissipation rate. From this case study, it is evident that components having lower critical temperatures and higher heat dissipation rates move toward the colder boundary to satisfy the temperature constraints during volume minimization. The final case study included pressure head loss as an additional physics constraint along with temperature, and geometric constraints. Smoother pipe bends were obtained in the final optimal layouts from this study as sharper bends cause greater head losses. Hence, there exists a trade-off between system volume and pressure head loss.

Some major advantages of the proposed design optimization method are (1) the GPM method is an explicit design representation approach that helps capture the complex geometry of the components and the interconnects, (2) efficient continuous parameterization with a gradient-based optimization approach can be used with this representation when projected on a finite element mesh, (3) this approach allows integration of multi-physics model into the optimization problem formulation, and (4) the primitive geometric shapes of any 3D engineering system can be represented using this approach. Some of the limitations of the proposed approach are (1) it would be interesting how the proposed method would scale with more than 30 or 40 components and 100s of interconnects and (2) the proposed method is based on a gradient-based optimization method and hence only can perform local search. Hence, there is a need for a method that can help explore multiple discrete topologies and a multi-start approach would be more ideal for a comprehensive search.

In future work, we plan to incorporate additional physics models such as thermoelastic, structural, and electrical models for describing system performance objectives and constraints. The proposed design framework does not consider irregular devices shapes. New design representations and shape descriptors for modeling complex device geometries will be investigated as part of future work. Furthermore, we plan to demonstrate our design method on larger-scale industry-relevant 3D applications that include multi-port components with arbitrarily shaped component geometries, interconnects, and 3D spatial topologies. These features in addition to more comprehensive physics-based optimization approaches will be useful in designing a realistic 3D component packing and routing design optimization framework in the future. The method presented here is expected to lead to faster development of holistic SPI2 design automation methods which are robust, efficient, and provide better design solutions than purely manual or heuristic design methods.

Acknowledgment

This work was supported by the National Science Foundation Engineering Research Center for Power Optimization of Electro-thermal Systems (POETS) with cooperative agreement EEC-1449548. The authors are also thankful to industry partner liaisons from CU Aerospace Ltd., PC Krause and Associates, and Ford Motor Company for their invaluable technical feedback on this work.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Nomenclature

m =

meters

r =

maximum vertex distance from the device center location

w =

width (or thickness) of the interconnect segment

x =

design variable vector

z =

reduced design variable vector

N =

shape functions of the element where the center lies

T =

nodal temperatures from the finite element formulation (KT = P)

V =

bounding box volume

R =

rotation matrix

T =

temperature vector (state variables)

dij =

minimum distance between the device center location and line segment representing the routing bar

fi =

friction factor

li =

length of a component

nd =

total number of devices (or components)

ns =

total number of interconnect segments

xc =

X-coordinate of the component center

yc =

Y-coordinate of the component center

zc =

Z-coordinate of the component center

bi =

vector pointing from the reference point (or center) to the vertices of a component

cd =

component reference (or center) point

gheadloss =

head loss constraint

gtemp =

device temperature constraint

gθ =

angle constraint

pi =

port locations vector

x0 =

interconnect start point

xf =

interconnect end point

Ai =

area of cross section of the pipe

Ki =

loss coefficient due to surface friction

Tcd =

temperature at the center of the devices

Tmax =

specified critical temperature

z' =

expanded design variable vector

mm =

millimeters

f(x, T) =

objective function

g(x, T) =

constraint function

gphys(x, T) =

physics-based constraints—they depend both on design and on solutions to physics models

gdd(x) =

non-interference constraints between two components

gds(x) =

non-interference constraints between one routing segment and one component

gss(x) =

non-interference constraints between two routing segments

gl(x) =

minimum bar length constraint

EV =

electric vehicle

GPM =

geometric projection method

HVACR =

heating, ventilation, air-conditioning, and refrigeration

MMA =

method of moving asymptotes

PR =

packing and routing

SPI2 =

spatial packaging of interconnected systems with physical interactions

W/m3 =

Watt per cubic meters—used for component heat dissipation rate

K(x) =

global thermal stiffness matrix

P (x) =

global load vector

°C =

degree celsius—temperature units

α =

angle at the elbow bends between interconnect segments

θi =

angle of rotation along ith axis direction

ρw =

weight density

References

1.
Sakti
,
A.
,
Zeidner
,
L.
,
Hadzic
,
T.
,
Rock
,
B. S.
, and
Quartarone
,
G.
,
2016
, “
Constraint Programming Approach for Spatial Packaging Problem
,”
International Conference on Integration of Constraint Programming, Artificial Intelligence, and Operations Research (CPAIOR 2016)
,
Bannf, Canada
,
May 29–June 1
, Springer International Publishing, pp.
319
328
.
2.
Peddada
,
S. R. T.
,
Zeidner
,
L. E.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2021
, “
An Introduction to 3D SPI2 (Spatial Packaging of Interconnected Systems With Physics Interactions) Design Problems: A Review of Related Work, Existing Gaps, Challenges, and Opportunities
,”
ASME 2021 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference (IDETC-CIE)
,
Virtual, Online
,
Aug. 17–19
, p. V03BT03A034.
3.
Peddada
,
S.
,
Zeidner
,
L.
,
Ilies
,
H. T.
,
James
,
K.
, and
Allison
,
J. T.
,
2022
, “
Toward Holistic Design of Spatial Packaging of Interconnected Systems With Physical Interactions (SPI2)
,”
ASME J. Mech. Des.
,
144
(
12
), p.
120801
.
4.
Bhattacharyya
,
A.
,
Peddada
,
S. R. T.
,
Bello
,
W. B.
,
Zeidner
,
L. E.
,
Allison
,
J. T.
, and
James
,
K. A.
,
2022
, “
Simultaneous 3d Component Packing and Routing Optimization Using Geometric Projection
,”
AIAA 2022 SciTech Forum
,
San Diego, CA and online
,
Jan. 3–7
.
5.
Peddada
,
S. R. T.
,
Dunfield
,
N. M.
,
Zeidner
,
L. E.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2021
, “
Systematic Enumeration and Identification of Unique Spatial Topologies of 3D Systems Using Spatial Graph Representations
,”
ASME 2021 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference (IDETC-CIE)
,
Virtual, Online
,
Aug. 17–19
, p. V03AT03A042.
6.
Peddada
,
S. R. T.
,
2023
, “Automated Interference-Free Layout Generation Methods for 2D Interconnected Engineering Systems,” Tech. Rep., https://www.ideals.illinois.edu/items/126471
7.
Peddada
,
S. R. T.
,
2021
, “A Two-Stage Design Framework for Optimal Spatial Packaging of Fluid-Thermal Systems,” Ph.D. thesis, University of Illinois at Urbana-Champaign, Urbana, IL, May, http://hdl.handle.net/2142/110458
8.
Qiang
,
L.
, and
Chengen
,
W.
,
2011
, “
A Discrete Particle Swarm Optimization Algorithm for Rectilinear Branch Pipe Routing
,”
Assembly Autom.
,
31
(
4
), pp.
363
368
.
9.
Shao
,
X.-Y.
,
Chu
,
X.-Z.
,
Qiu
,
H.-B.
,
Gao
,
L.
, and
Yan
,
J.
,
2009
, “
An Expert System Using Rough Sets Theory for Aided Conceptual Design of Ship’s Engine Room Automation
,”
Expert Syst. Appl.
,
36
(
2, Part 2
), pp.
3223
3233
.
10.
López-Camacho
,
E.
,
Ochoa
,
G.
,
Terashima-Marín
,
H.
, and
Burke
,
E. K.
,
2013
, “
An Effective Heuristic for the Two-Dimensional Irregular Bin Packing Problem
,”
Ann. Oper. Res.
,
206
(
1
), pp.
241
264
.
11.
Tejani
,
G. G.
,
Savsani
,
V. J.
,
Patel
,
V. K.
, and
Savsani
,
P. V.
,
2018
, “
Size, Shape, and Topology Optimization of Planar and Space Trusses Using Mutation-Based Improved Metaheuristics
,”
J. Comput. Des. Eng.
,
5
(
2
), pp.
198
214
.
12.
Gulić
,
M.
, and
Jakobović
,
D.
,
2017
, “
Evolution of Vehicle Routing Problem Heuristics With Genetic Programming
,”
2013 36th International Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO)
,
Opatija, Croatia
,
May 20–24
, pp.
988
992
.
13.
Beckert
,
B. A.
,
2013
, “
When Engineering Intuition Is Not Enough
,”
Altair Eng. Res.
,
1
(
1
), pp.
1
4
.
14.
Bayrak
,
A. E.
,
Ren
,
Y.
, and
Papalambros
,
P. Y.
,
2016
, “
Topology Generation for Hybrid Electric Vehicle Architecture Design
,”
ASME J. Mech. Des.
,
138
(
8
), p.
081401
.
15.
Aladahalli
,
C.
,
Cagan
,
J.
, and
Shimada
,
K.
,
2006
, “
Objective Function Effect Based Pattern Search—Theoretical Framework Inspired by 3D Component Layout
,”
ASME J. Mech. Des.
,
129
(
3
), pp.
243
254
.
16.
Szykman
,
S.
, and
Cagan
,
J.
,
1995
, “
A Simulated Annealing-Based Approach to Three-Dimensional Component Packing
,”
ASME J. Mech. Des.
,
117
(
2A
), pp.
308
314
.
17.
Yin
,
S.
,
Cagan
,
J.
, and
Hodges
,
P.
,
2004
, “
Layout Optimization of Shapeable Components With Extended Pattern Search Applied to Transmission Design
,”
ASME J. Mech. Des.
,
126
(
1
), pp.
188
191
.
18.
Yin
,
S.
, and
Cagan
,
J.
,
2000
, “
An Extended Pattern Search Algorithm for Three-Dimensional Component Layout
,”
ASME J. Mech. Des.
,
122
(
1
), pp.
102
108
.
19.
Szykman
,
S.
, and
Cagan
,
J.
,
1997
, “
Constrained Three-Dimensional Component Layout Using Simulated Annealing
,”
ASME J. Mech. Des.
,
119
(
1
), pp.
28
35
.
20.
Yano
,
N.
,
Morinaga
,
T.
, and
Saito
,
T.
,
2008
, “
Packing Optimization for Cargo Containers
,”
2008 SICE Annual Conference
,
Tokyo, Japan
,
Aug. 20–22
, pp.
3479
3482
.
21.
Park
,
J.-H.
, and
Storch
,
R.
,
2002
, “
Pipe-Routing Algorithm Development: Case Study of a Ship Engine Room Design
,”
Expert Syst. Appl. (UK)
,
23
(
3
), pp.
299
309
.
22.
Betz
,
V.
, and
Rose
,
J.
,
1997
, “VPR: A New Packing, Placement and Routing Tool for FPGA Research,”
Field-Programmable Logic and Applications
,
Luk
,
W.
,
Cheung
,
P. Y. K.
, and
Glesner
,
M.
, eds.,
Springer
,
Berlin
, pp.
213
222
.
23.
Szykman
,
S.
, and
Cagan
,
J.
,
1996
, “
Synthesis of Optimal Nonorthogonal Routes
,”
ASME J. Mech. Des.
,
118
(
3
), pp.
419
424
.
24.
Guirardello
,
R.
, and
Swaney
,
R. E.
,
2005
, “
Optimization of Process Plant Layout With Pipe Routing
,”
Comput. Chem. Eng.
,
30
(
1
), pp.
99
114
.
25.
Lv
,
N.
,
Liu
,
J.
,
Xia
,
H.
,
Ma
,
J.
, and
Yang
,
X.
,
2020
, “
A Review of Techniques for Modeling Flexible Cables
,”
Comput. Aided Des.
,
122
, p.
102826
.
26.
Wehlin
,
C.
,
Persson
,
J. A.
, and
Ölvander
,
J.
,
2020
, “
Multi-objective Optimization of Hose Assembly Routing for Vehicles
,”
16th International DESIGN Conference (Proceedings of the Design Society)
,
Virtual, Online
,
Oct. 26–29
, Vol. 1, pp.
471
480
.
27.
Souissi
,
O.
,
Benatitallah
,
R.
,
Duvivier
,
D.
,
Artiba
,
A.
,
Belanger
,
N.
, and
Feyzeau
,
P.
,
2013
, “
Path Planning: A 2013 Survey
,”
2013 International Conference on Industrial Engineering and Systems Management (IESM)
,
Agdal, Morocco
,
Oct. 28–30
, pp.
1
8
.
28.
Hermansson
,
T.
,
Bohlin
,
R.
,
Carlson
,
J. S.
, and
Söderberg
,
R.
,
2016
, “
Automatic Routing of Flexible 1d Components With Functional and Manufacturing Constraints
,”
Comput. Aided Des.
,
79
(
10
), pp.
27
35
.
29.
Zhu
,
Z.
,
La Rocca
,
G.
, and
van Tooren
,
M. J. L.
,
2017
, “
A Methodology to Enable Automatic 3D Routing of Aircraft Electrical Wiring Interconnection System
,”
CEAS Aeronaut. J.
,
8
(
2
), pp.
287
302
.
30.
Belov
,
G.
,
Du
,
W.
,
Banda
,
M. G. D. L.
,
Harabor
,
D. D.
,
Koenig
,
S.
, and
Wei
,
X.
,
2020
, “
From Multi-Agent Pathfinding to 3D Pipe Routing
,”
Thirteenth Annual Symposium on Combinatorial Search (SOCS)
,
Online
,
May 26–28
.
31.
Rourke
,
P.
, and
Merritt
,
L.
,
2007
, “
Real-Time Semiautomatic 3D Pipe Routing
,”
J. Ship Prod.
,
23
(
03
), pp.
180
184
.
32.
Niu
,
Y.
,
Niu
,
W.
, and
Gao
,
W.
,
2017
, “Branch Pipe Routing Method Based on a 3D Network and Improved Genetic Algorithm,”
Civil, Architecture and Environmental Engineering
, 1st ed., Vol.
2
,
Taylor and Francis
,
Milton Park, Oxfordshire, UK
, pp.
1
6
.
33.
Zhou
,
Q.
, and
Lv
,
Y.
,
2020
, “
Research Based on Lee Algorithm and Genetic Algorithm of the Automatic External Pipe Routing of the Aircraft Engine
,”
Int. J. Mech. Eng. Appl.
,
8
(
1
), pp.
40
44
.
34.
Qu
,
Y.
,
Jiang
,
D.
, and
Yang
,
Q.
,
2018
, “
Branch Pipe Routing Based on 3D Connection Graph and Concurrent Ant Colony Optimization Algorithm
,”
J. Intell. Manuf.
,
29
, pp.
1647
1657
.
35.
Stanczak
,
M.
,
Pralet
,
C.
,
Vidal
,
V.
, and
Baudoui
,
V.
,
2020
, “
Optimal Pipe Routing Techniques in an Obstacle-Free 3D Space
,”
ICORES
,
La Valette, Malta
,
Feb. 22–24
.
36.
Suchorab
,
P.
, and
Kowalski
,
D.
,
2018
, “
Methods of Routing and Sizing of Water Supply Networks
,”
E3S Web Conf.
,
59
, p.
00024
.
37.
Agafonov
,
A. A.
, and
Myasnikov
,
V. V.
,
2018
, “
Vehicle Routing Algorithms Based on a Route Reservation Approach
,”
J. Phys.: Conf. Ser.
,
1096
(
1
), p.
012029
.
38.
Saleh
,
Y.
,
Tofigh
,
A.
, and
Zahra
,
A.
,
2014
, “
Transportation Routing in Urban Environments Using Updated Traffic Information Provided Through Vehicular Communications
,”
J. Transp. Syst. Eng. Inf. Technol.
,
14
(
5
), pp.
23
36
.
39.
Persson
,
L. G.
,
Santos
,
F. B.
,
Tavares
,
C. A. C.
,
de Andrade
,
A. E.
,
de Brito Alves
,
R. M.
, and
Biscaia
,
E. C.
,
2009
, “Methodology of Pipe and Equipment Layout for On-Shore Oil and Gas Industry,”
Comput. Aided Chem. Eng.
,
27
, pp.
1845
1850
.
40.
Peddada
,
S. R. T.
,
Herber
,
D. R.
,
Pangborn
,
H. C.
,
Alleyne
,
A. G.
, and
Allison
,
J. T.
,
2019
, “
Optimal Flow Control and Single Split Architecture Exploration for Fluid-Based Thermal Management
,”
ASME J. Mech. Des.
,
141
(
8
), p.
083401
.
41.
Peddada
,
S. R. T.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2020
, “
A Novel Two-Stage Design Framework for Two-Dimensional Spatial Packing of Interconnected Components
,”
ASME J. Mech. Des.
,
143
(
3
), p.
031706
.
42.
Jessee
,
A.
,
Peddada
,
S. R. T.
,
Lohan
,
D. J.
,
Allison
,
J. T.
, and
James
,
K. A.
,
2020
, “
Simultaneous Packing and Routing Optimization Using Geometric Projection
,”
ASME J. Mech. Des.
,
142
(
11
), p.
111702
.
43.
Bello
,
W. B.
,
Peddada
,
S. R. T.
,
Bhattacharyya
,
A.
,
Jennings
,
M.
,
Katragadda
,
S.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2022
, “
Underhood Spatial Packing and Routing of an Automotive Fuel Cell System (AFCS) Using 2D Geometric Projection
,”
AIAA SciTech 2022 Forum
,
San Diego, CA and Online
,
Jan. 3–7
.
44.
Peddada
,
S. R. T.
,
Rodriguez
,
S. B.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2020
, “
Automated Layout Generation Methods for 2D Spatial Packing
,”
ASME International Design Engineering Technical Conferences & Computers and Information in Engineering Conference (IDETC)
,
Virtual, Online
,
Aug. 17–20
.
45.
Norato
,
J.
,
Bell
,
B.
, and
Tortorelli
,
D.
,
2015
, “
A Geometry Projection Method for Continuum-Based Topology Optimization With Discrete Elements
,”
Comput. Methods Appl. Mech. Eng.
,
293
, pp.
306
327
.
46.
Peddada
,
S. R. T.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2020
, “
A Novel Two-Stage Design Framework for 2D Spatial Packing of Interconnected Components
,”
ASME 2020 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference (IDETC-CIE)
,
Online, Virtual
,
Aug. 17–19
, p. V11BT11A032.
47.
Smith
,
H.
, and
Norato
,
J. A.
,
2020
, “
A Matlab Code for Topology Optimization Using the Geometry Projection Method
,”
Struct. Multidiscipl. Optim.
,
62
, pp.
1579
1594
.
48.
Peddada
,
S. R. T.
,
Dunfield
,
N. M.
,
Zeidner
,
L. E.
,
Givans
,
Z. R.
,
James
,
K. A.
, and
Allison
,
J. T.
,
2023
, “
Enumeration and Identification of Unique 3d Spatial Topologies of Interconnected Engineering Systems Using Spatial Graphs
,”
ASME J. Mech. Des.
,
145
(
10
), p.
101708
.
49.
Smith
,
H.
, and
Norato
,
J. A.
,
2020
, “
A Matlab Code for Topology Optimization Using the Geometry Projection Method
,”
Struct. Multidiscipl. Optim.
,
62
(
3
), pp.
1579
1594
.
50.
Bendsøe
,
M. P.
,
1989
, “
Optimal Shape Design as a Material Distribution Problem
,”
Struct. Optim.
,
1
(
4
), pp.
193
202
.
51.
Rozvany
,
G.
,
2000
, “
The SIMP Method in Topology Optimization—Theoretical Background, Advantages and New Applications
,”
8th AIAA Symposium on Multidisciplinary Analysis and Optimization
,
Long Beach, CA
,
Sept. 6–8
.
52.
Yulin
,
M.
, and
Xiaoming
,
W.
,
2004
, “
A Level Set Method for Structural Topology Optimization and Its Applications
,”
Adv. Eng. Softw.
,
35
(
7
), pp.
415
441
.
53.
van Dijk
,
N. P.
,
Maute
,
K.
,
Langelaar
,
M.
, and
van Keulen
,
F.
,
2013
, “
Level-Set Methods for Structural Topology Optimization: A Review
,”
Struct. Multidiscipl. Optim.
,
48
(
3
), pp.
437
472
.
54.
Nakayama
,
Y.
,
2018
, “Chapter 7—Flow in Pipes,”
Introduction to Fluid Mechanics
, 2nd ed.,
Butterworth-Heinemann
,
Amsterdam, Netherlands
, pp.
135
161
.