Abstract

This work presents a method for the topology optimization of welded frame structures to minimize the manufacturing cost. The structures considered here consist of assemblies of geometric primitives such as bars and plates that are common in welded frame construction. A geometry projection technique is used to map the primitives onto a continuous density field that is subsequently used to interpolate material properties. As in density-based topology optimization techniques, the ensuing ersatz material is used to perform the structural analysis on a fixed mesh, thereby circumventing the need for re-meshing upon design changes. The distinct advantage of the representation by geometric primitives is the ease of computation of the manufacturing cost in terms of the design parameters, while the geometry projection facilitates the analysis within a continuous design region. The proposed method is demonstrated via the manufacturing-cost-minimization subject to a displacement constraint of 2D bar, 3D bar, and plate structures.

1 Introduction

Manufacturing cost is seldom considered in the topology optimization of structures, which is often based solely on structural criteria and weight. Consequently, the optimal design may exhibit superior mechanical performance but be costly to manufacture. While the amount of material in a structure may contribute significantly to its cost, there are also other costs associated with the fabrication process. In the case of the welded frames considered in this work, the manufacturing cost is also influenced by the cost of cutting, welding, and painting frame members.

The main obstacle to incorporating manufacturing cost in topology optimization lies in the difficulty to express it in terms of the design parameters. Conventional density-based and level-set topology optimization techniques employ a field representation of the structure, which endows the optimizer with substantial design freedom but makes it difficult to compute manufacturing cost.

Ground structure methods have been used to minimize the manufacturing cost of truss structures [1]. In these methods, which use 1D-elements to represent the structure and for analysis, it is easy to compute costs associated with, for example, the length of the truss elements and the number of struts. However, other cost components, such as the cost of welding, cannot be readily computed because they require a calculation of the welding length, which is difficult to compute from the 1D-representation. Moreover, the structure must be a topological subset of the initial design (the ground structure), thus ground structure methods cannot model arbitrary topologies. Also, a ground structure approach does not accommodate primitives like plates.

Some density-based topology optimization techniques for continua consider the manufacturing cost of additively manufactured single components [25] or assemblies of 2-dimensional stamped and spot-welded components [6,7]. For these manufacturing processes, it is possible to express the manufacturing cost in terms of quantities that can be computed from the field representation, such as surface area, bounding box, and the volume of the support material. Some topology optimization techniques (e.g., Ref. [8]) design multi-component structures in which the boundaries among components in the optimal design are intended to be joined via continuous welds. However, in these techniques, there is no computation of the weld length between components and no consideration for the welding cost.

It should be noted that some methods have incorporated various mechanisms to limit the number of structural members or holes in the optimization, which is an indirect way of controlling manufacturing cost. These works include the ground structure approach of Ref. [9] and the moving morphable components method of Ref. [10], which impose a constraint on the maximum number of members in the structure; and the work of Ref. [11], where a limit is imposed on the maximum number of holes in density-based topology optimization.

In this paper, we propose a method for the topology optimization of welded frame structures with regard to manufacturing cost. The frame is modeled as the union of bar primitives or plate primitives. The analysis is performed on a fixed finite element mesh as in density-based methods. To employ the primitive-based representation of the frame while performing the analysis on a fixed mesh, we employ the geometry projection (GP) method [12], whereby the geometric parameters of the primitives are smoothly mapped onto a density field that is subsequently used to interpolate material properties, just as in density-based methods. The GP mapping is differentiable, hence we can use efficient gradient-based methods for the optimization.

The proposed approach has several advantages. Each term of the manufacturing-cost function can be computed directly in terms of the geometric parameters of the primitives, the projected density field, or both. For example, the geometry representation allows us to determine the weld location among components and approximate the weld length. Re-meshing is circumvented as in conventional topology optimization techniques. Compared to ground structure methods that use 1D-elements for the analysis, the bars need not be connected during the optimization and the optimal design is not a subset of the ground structure, leading to increased design freedom and thus efficient structures that can be generated with only a few bars (the trade-off, however, is that the analysis in the GP method discretizes a continuum with 2D- and 3D-elements and thus is more computationally expensive). Moreover, the GP method captures intersections between 3D-primitives in a way that is not possible using 1D-elements. For instance, in methods that use truss elements to model cylindrical bars, two elements that come very close to each other may not intersect even if the corresponding cylindrical bars in 3D do intersect because of their dimensions. Finally, although stress constraints are out of the scope of the current work, we note that the continuum 2D- and 3D-representation of the structure can capture stress concentrations arising from the intersection of primitives, which cannot be captured by 1D-elements.

The remainder of the paper is organized as follows: in Sec. 2, we briefly introduce the GP method. Section 3 formulates the manufacturing-cost function based on the GP method for both bar and plate primitives. The optimization problem is stated in Sec. 4. A brief description of design sensitivities is given in Sec. 5. In Sec. 6, we demonstrate the proposed method in the design of 2D and 3D cantilever beams with bar primitives and of a Messerschmitt–Bölkow–Blohm (MBB) beam in 3D with plate primitives. Finally, we draw conclusions in Sec. 7.

2 Geometry Projection

The proposed technique employs the GP method to map the geometric primitives that make up the structure onto a fixed finite element mesh for analysis, thus avoiding re-meshing upon design changes [12,13]. The GP method consists of a differentiable map between the high-level parameters that describe the geometric primitives and a density field defined over the design region, which is subsequently discretized via a finite element mesh for analysis. The differentiability of the map ensures that efficient gradient-based techniques can be used for the optimization.

In this work, the projected density at point x corresponding to component c is computed as
ρc(x;R):=H~(ϕc(x,zc)R)
(1)
where
H~(x):={1ifx1(x+1)3(3x29x+8)16if|x|<10otherwise
(2)
is a C2 regularized Heaviside function, ϕc(x, zc) is the signed distance from x to the boundary of c, zc is the vector of design parameters defining component c, and R is the radius of the transition region of the Heaviside. In all the examples shown in Sec. 6, R is equal to the diagonal of the element, i.e., twice the size of the radius of the ball that circumscribes the element. We note that, as demonstrated in Refs. [12,13], R must be strictly smaller than the radius of the primitive to ensure well-defined sensitivities. Correspondingly, the element size must be at most half of the primitive radius.
The formulation of the signed distance ϕc depends on the specific representation of the geometric components. For the frame structures considered in this paper, we consider bar and plate primitives represented by offset solids. In the case of bars, the solid corresponds to all points within a distance rc of a line segment (the bar’s medial axis), and the corresponding vector of design parameters is zc = {x1c, x2c, rc, αc}, where x1c and x2c denote the endpoints of the medial axis and αc is a membership variable that will be introduced later in this section. In the case of plates, we use quaternions to represent the plate orientation to avoid issues like gimbal lock and 2π-periodicity that arise when using Euler angles [14]. The vector of quaternion components describing the orientation in space of the plate is denoted as qc={qcr,qcx,qcy,qcz}. The rotation matrix of plate c is given as
Rc=[ec1ec2ec3]ec1=[12qcy22qcz22qcxqcy2qcrqcz2qcxqcz+2qcrqcy]ec2=[2qcxqcy+2qcrqcz12qcx22qcz22qcyqcz2qcrqcx]ec3=[2qcxqcz2qcrqcy2qcyqcz+2qcrqcx12qcx22qcy2]
(3)
where ec1,ec2,ec3 are the basis vectors of plate c in the local coordinate system. Therefore, the vector of design parameters of the plate primitive is zc={xc0,lc1,lc2,qc,rc,αc}, where xc0 is the location of the center of plate, lc1 and lc2 are the dimensions of the rectangular medial surface, and rc is the semi-thickness of the plate. The design parameters are illustrated in Fig. 1.
Fig. 1
Geometric primitives used in this work with their respective geometric parameters: 2D bars (top), 3D bars (middle), and 3D plates (bottom). Dashed line in bars corresponds to the medial axis; rectangle in plates corresponds to medial surface.
Fig. 1
Geometric primitives used in this work with their respective geometric parameters: 2D bars (top), 3D bars (middle), and 3D plates (bottom). Dashed line in bars corresponds to the medial axis; rectangle in plates corresponds to medial surface.
Close modal

The convenience of offset solids is that the signed distance from any point x to the boundary of the primitive can be simply computed as ϕc(x) = rcdc(x), where dc(x) is the distance from x to the medial axis/surface. Moreover, dc can be computed in closed form as a function of zc. For details on the computation of the signed distance for offset bar and plate primitives and the corresponding design sensitivities, the reader is referred to Refs. [1416]. It should be noted, however, that the geometry projection is not limited to offset solids and can be used with any geometric representation, provided a signed distance to the boundary of the primitive can be computed.

In addition to the geometric parameters that determine the dimensions, position, and orientation of the component, we also ascribe a continuous membership variable αc to each primitive. When αc = 1, it signifies the component is part of the design, and an interior point of the component will have the elastic properties of the material that the component is made of. Conversely, when αc = 0, the component is removed from the structure, and the elastic properties at an interior point of the component will not be influenced by the component, regardless of the component’s dimensions. This membership variable enables the optimizer to remove geometric components by penalization. The penalized densityρ^c(zc,x) that incorporates with the membership variable αc is given by
ρ^c(zc,x):=μ(αcρc(zc,x))
(4)
where μ is a penalization function that renders intermediate values of αc structurally inefficient. This penalization ensures that the optimal design has bars with 0/1 membership variables and is essentially the same penalization used by density-based topology optimization technique. For instance, it may correspond to the power law of the solid isotropic material with penalization interpolation scheme (cf. Ref. [17]) given by μ(x) = xP, with P ≥ 3.
We combine multiple primitives via their Boolean union. Since the projected and penalized densities are ultimately implicit representations of the primitives, the Boolean union corresponds to their maximum (cf. Ref. [18]). The maximum function, however, is not differentiable, which precludes the use of efficient gradient-based optimizers. To circumvent this, we use a softmax differentiable approximation in which the interpolated material properties are given by [16]
C(x)=Cvoid+ρ~(Z,x)(CsolidCvoid)
(5)
where
ρ~(Z,x)):=c=1Ncwc(ρ^(Z,x))ρ^c(Z,x)
(6)
is a combined density that corresponds to the Boolean union of all components, and the weights wc are given by
wc(ρ^(Z,x)):=softargmaxc(ρ^(Z,x);p)=epρ^c(zc,x)j=1Ncepρ^j(zj,x)
(7)
In these expressions, Csolid denotes the elasticity tensor of solid material that each component is made of, ρ^ is the vector of penalized densities for all the components, and Cvoid is the elasticity tensor of a relatively weak material to prevent an ill-posed analysis. Nc is the number of components. Z={z1,,zNc} denotes the vector of design parameters for all components. As the parameter p → ∞, the weights wc in the softargmax function approach a one-hot vector, i.e., wc → 1 for the component with the largest penalized density and wc → 0 for all other components. In the finite element analysis, we assume for simplicity an element-uniform combined density ρ~e, which is computed at the element centroid; consequently the ersatz elasticity tensor of (5) is also element-uniform.

While the softmax material interpolation of (5) admits a different elasticity tensor for each component, which accommodates multi-material structures and anisotropic materials (cf. Ref. [16]), in this work we consider that all components are made of a single material with elasticity tensor Csolid.

3 Manufacturing-Cost Function

We formulate the manufacturing-cost function based on the fabrication cost presented in Ref. [19]. The manufacturing-cost function consists of two parts: material cost Cm and fabrication cost Cf. We take four stages of fabrication into consideration, namely preparation, cutting, welding, and painting. The relative weights of the terms making up this function correspond to monetary cost per unit (e.g., length, area) for each stage. Different supply chains would have different values of these weights and therefore may lead to different designs. The manufacturing cost of the structure is defined as
M(Z):=Cm(Z)+Cf(Z)=Cm(Z)+iCi(Z)
(8)
where Ci, i ∈ {1, 2, 3, 4}, is the cost for each manufacturing stage. It is important to note that M is an estimate, since an accurate assessment of the manufacturing cost of a structure can only be made once a detailed design of the structure and detailed knowledge of the cost structure for the material and the manufacturing process are available. Therefore, it is essential to emphasize that the place of the proposed technique in the structure’s design workflow is in the conceptual design stage, and that its aim is to produce designs with significantly improved manufacturing cost as compared to that of an optimal structure driven purely by mechanical criteria.

In order to tie the manufacturing cost to the geometric description of the primitives, each of the terms in (8) must be expressed as a function of the geometric parameters of the components. The remainder of this section details these terms. Some of them are computed directly in terms of the geometric parameters, while other terms are more easily computed in terms of the projected densities. In this sense, the manufacturing-cost function introduced here takes advantage of the dual representation (geometric parameters/densities). It is worth noting that even for terms that are computed in terms of projected densities, their calculation is possible because there is a projected density field ρc associated with each component c and ρc serves as a (fuzzy) point classification (i.e., an inside/outside test), which makes it possible to compute, for example, the weld length of Sec. 3.4.

3.1 Material Cost.

The material cost Cm is simply computed as the weight of the material cost wm times the mass of the whole structure ϱV:
Cm=wmϱV(Z)
(9)
where wm depends on the type of the material and the detailed supply chain, ϱ is the material density (which, in this work, we assume to be 1.0 for all the examples), and V is the volume of the structure. V can be computed as the sum of the volumes of all the elements in the mesh:
V:=e=1Neρeve
(10)
where Ne is the number of elements in the mesh, ve is the volume of a single element e, and ρe is the projected density of element e.

3.2 Preparation Cost.

Before components are fully welded, they usually need to be prepared, assembled, and pre-positioned by tack welding. The cost of preparation stage mainly depends on the number of the components and the weight of each component. According to Ref. [19], the cost for preparation, assembly, and tacking is defined as
C1:=w1(κϱV)1/2
(11)
where w1 is the corresponding weight of C1 and κ is the number of existing components in the structure. In the GP method, κ can be expressed as the sum of the membership variable αc of all the components:
κ:=c=1Ncαc
(12)

It should be noted that κNc. While Nc corresponds to all the components available to design the structure, κ represents the components actually used for the design.

3.3 Cutting Cost.

Cutting and edge grinding can be made with different technologies, such as flame cutting, plasma cutting, and laser cutting. The corresponding cost C2 is simply proportional to the cutting area of the component:
C2:=w2Acut
(13)
where the weight w2 of the cutting cost depends on the detailed manufacturing condition, and Acut is the total cross-sectional area to be cut for all the components. By using the GP method, the area of cutting for different components is computed as
Acut={αc(2rctc)for 2D bar primitivesαc(πrc2)for 3D bar primitivesαc(2rc(lc1+lc2))for 3D plate primitives
(14)
where tc is the out-of-plane thickness of bar c. These expressions assume each bar component only needs to be cut at one end, while plate components are cut at all edges. The presence of the membership variable αc in (14) ensures no cutting area is computed for components that have been removed from the design. Also, we note that this is an estimation, as an accurate calculation of cutting cost requires a detailed design.

3.4 Welding Cost.

The cost of welding C3 incorporates not only the welding process itself but also all the additional fabrication steps associated with welding. These additional steps include changing the electrode, deslagging, and chipping. C3 depends on the welding length lw of the whole structure:
C3:=w3lw
(15)
where the weight w3 of the welding cost depends on, e.g., the technology of welding. The length of welding can be computed as follows: the norm of the gradient of the projected density ρc(x) is nonzero for a point x on the boundary of the component and zero elsewhere. For two components i and j, the quantity
lij(x)={ρi(x)ρj(x)in 2Dρi(x)ρj(x)in 3D
(16)
is 1 if x belongs to the weld between the two components and 0 otherwise. The total weld length is subsequently given by
lw:=Ωi=1Ncj>iNcαiαjlijdΩ
(17)
where Ω denotes the design region. The condition j > i for the innermost sum of (17) ensures that the intersection between two components is counted only once, and it also prevents the situation in which i = j, for which Ωlii would render the perimeter in 2D and the surface area in 3D for component i. One could alternatively use the condition j < i to achieve the same purpose. Since we assume uniform projected densities and gradients within each element, the integral in (17) is computed as a sum over the elements (with the integrand multiplied by the element volume). Figure 2 shows an example of bar intersections plotting the integrand of (16).
Fig. 2
Plot of the quantity lij for j > i in (16) indicating the portion of the intersection between primitives counted towards the weld length lw: 2D bars (left), 3D bars (middle), and 3D plates (right). The weld regions shown at the intersections between primitives correspond to an iso-surface with lij ≥ 0.01.
Fig. 2
Plot of the quantity lij for j > i in (16) indicating the portion of the intersection between primitives counted towards the weld length lw: 2D bars (left), 3D bars (middle), and 3D plates (right). The weld regions shown at the intersections between primitives correspond to an iso-surface with lij ≥ 0.01.
Close modal
We note that from (1) and using the chain rule, we have that
ρc=1RH~(ϕc/R)ϕc
where H~ denotes the first derivative of H~. Since the distance function satisfies the eikonal equation ϕc=1, the norm of the gradient of the projected density for component c can simply be computed as
ρc=1RH~(ϕc/R)
(18)
which can be readily obtained from (2).

Note that even though the weld length is computed in terms of the projected densities and not directly in terms of the geometric parameters of the components, this calculation is only possible because we have distinct structural members and can compute ‖ρc‖ and αc for each component.

3.5 Painting Cost.

The fourth term of the fabrication cost in (8) includes the cost of painting and surface preparation. The surface preparation includes cleaning (e.g., sand-spray), grinding, and application of a top coat. This cost of painting C4 is proportional to the surface area of the whole structure
C4:=w4As
(19)
with the proportionality constant weight w4. For frames made of 2D bars, the surface area is computed as
As=c=1Ncαc((2πrc+2lc)tc+2(πrc2+2rclc))
(20)
For structures consisting of 3D bars, the surface area is given by
As=c=1Ncαc(2πrclc+4πrc2)
(21)
and for plates it is given by
As=c=1Ncαc(4πrc2+2πrc(lc1+lc2)+2lc1lc2)
(22)

4 Optimization Problem

We consider minimization of the manufacturing cost subject to a constraint that the displacement at a point p does not exceed a specified value:
minZM(Z)

(23)
Subject to:g(up(Z))0
(24)
K(Z)u(Z)=f
(25)
zi_zizi¯,i=1,,NdNc
(26)
The constraint of (25) represents the discretized finite element equations for the linear elasticity problem, with K, u, and f the global stiffness matrix, displacement vector, and force vector, respectively. We assume that the applied force is design-independent.
In the displacement constraint of (24), up = (uTCu)1/2 is the magnitude of the displacement at a specified point p, with C a square matrix such that Cii = 1 for the degrees-of-freedom i corresponding to point p, and Cij = 0 for every other component. For all the examples, the function g provides scaling to aid convergence. Since the displacement magnitude can attain extremely large values for disconnected design during the optimization, the following log-scale version of the constraint ensures that values in consecutive iterations do not change drastically:
g(up):=1log2logup+u¯2u¯
(27)
where u¯ is the maximum allowable displacement at p. For a positive displacement in double precision, we can assume that up ∈ [0, 1016]. Therefore, when up = 0, g(up) = −1; when up=u¯, g(up) = 0; when up = 1016, g(up)(53.151.443log2u¯). Thus, the displacement at a point constraint is bounded as g[1,(53.151.443log2u¯)], which is in the range of values recommended for the method of moving asymptotes (MMA). For the same reason, we linearly scale down the value of the manufacturing cost by 100 times. For all the examples in this paper, p is the point at which the load is located.
We impose lower and upper bounds, zi, zi¯, on the individual design variables zi in (26), with Nd the number of design variables per component (six for 2D bars, eight for 3D bars, and 11 for 3D plates). The endpoints xic of the medial axes of bars and the center points xc0 of plates are bound to lie inside the design region. The membership variables are bound as αc ∈ [0, 1], and the quaternion components as qci ∈ [−1, 1]. The dimensions of the components are given different bounds for the examples in Sec. 6. As in prior implementations of the GP method (cf. Ref. [12]), all design variables are scaled as
z^i:=ziz_iz¯iz_i[0,1]
(28)
and a move limit m is imposed at each iteration I on the scaled variables to improve convergence:
max(0,ziI1m)ziImin(1,ziI1+m)
(29)

5 Sensitivity Analysis

To employ efficient gradient-based optimizers, we compute design sensitivities of the manufacturing-cost function and the displacement constraint. For a given function J(z):=J~(ρ~(z),z), its sensitivities with respect to a design variable zi are computed as
Jzi=e=1NeJ~ρ~eρ~ezi+J~zi
(30)
where ρ~e is the combined density for element e. For functions that depend on the solution u to the equilibrium equation of (25), J~/ρ~e is computed using adjoint differentiation and is similar to the sensitivity computed in density-based topology optimization. For functions that do not depend on the analysis, J~/ρ~e and/or J~/zi can be computed directly, as is the case for all the terms of the manufacturing-cost function. The term ρ~e/zi is obtained from the geometry projection; its derivation is here omitted for brevity, and the reader is referred to Refs. [1416] for details. However, it is worth noting that the sensitivity of the norm of the projected density gradient is computed from (18) as
ρczi=1R2H~(ϕc/R)ϕczi
(31)
where H~ denotes the second derivative of H~ which is readily obtained from (2). The term ∂ϕc/∂zi can be found in the aforementioned references. The derivatives of all the terms in the manufacturing-cost function can be obtained from their respective expressions.

6 Numerical Examples

In this section, we present three examples to demonstrate the proposed method. The computer implementation of the method is done in matlab. The finite element analysis is performed using a regular mesh of bilinear quadrilateral elements. The linear system of equations arising from the finite element discretization is solved using the conjugate gradient method with a multigrid preconditioner. A modulus of Emin = 10−6 is used for the weak material with elasticity tensor Cvoid of (5). The solid material in all examples (i.e., Csolid in (5)) has a Young’s modulus of 1 and a Poisson’s ratio of 0.3

Each design update in the optimization is obtained as follows. An arbitrary initial design composed of a finite number of primitives is specified. For this initial design or any given design Z, we compute the projected densities of (1) for each component at the centroid of each element in the mesh. The combined density of (6) and the ersatz elasticity tensor of (5) are subsequently computed for each element. This is followed by the usual finite element assembly and solution, where each element stiffness matrix is computed using the corresponding ersatz material. The solution of the finite element analysis is used to compute the relevant objective and constraint functions in the optimization. Where applicable (e.g., when a displacement constraint is imposed), adjoint analysis is also performed to compute design sensitivities. The objective and constraints, together with their sensitivities, are passed to a gradient-based optimization algorithm, which produces a new design. The foregoing process is repeated until a stopping criterion is satisfied. For simplicity, we stop the optimization when a maximum number of iterations has been completed. The optimization problem is solved using MMA with the default parameters suggested in Ref. [20].

To demonstrate the effectiveness of the proposed approach, the examples considered in this work compare the manufacturing-cost-minimization problem of (23) to a compliance-minimization problem subject to a volume constraint. The latter problem can be stated as
minZC(Z):=u(Z)TfSubject to:v(Z)v¯K(Z)u(Z)=fzi_zizi¯,i=1,,NdNc
(32)
where C is the compliance of the structure, v is the volume fraction, and v¯ is the maximum allowable volume fraction.

6.1 Two-Dimensional Cantilever Beam With Bar Primitives.

In this first example, we design a cantilever beam with 2D bar primitives. The initial design, loading and boundary conditions are depicted in the top portion of Fig. 3. The dimensions of the design region are 60 × 10. The mesh size is 360 × 60. The left-hand edge of the design domain is fixed, and the load F = 0.1 is applied at the middle point of the right-hand edge. The initial design consists of 42 bars with radius and membership variables of 0.5. The bottom portion of Fig. 3 shows the weld locations (i.e., the integrand of αiαjlij in (17)) in the initial design. We bound the radius of each bar to [0.25, 2], and prescribe out-of-plane thickness of the 2D bar as tc = 1. The move limit is m = 0.025.

Fig. 3
2D cantilever beam example. Top: initial design, loading and boundary conditions. Bottom: weld locations shown.
Fig. 3
2D cantilever beam example. Top: initial design, loading and boundary conditions. Bottom: weld locations shown.
Close modal

To obtain the maximum allowable displacement u¯ for the displacement constraint in the manufacturing-cost-minimization problem of (23), we first perform the compliance minimization problem of (32) with a maximum material volume fraction v¯=0.3. The final design and its weld locations are shown in Fig. 4. We deem this design as a reference for comparison, and we choose the point where the load is applied as point p for the displacement constraint in the manufacturing-cost-minimization problems. Therefore, we subsequently assign u¯=uref, where uref denotes the displacement at point p in the reference minimum-compliance design.

Fig. 4
Minimum-compliance (reference) design for the 2D cantilever beam example. Top: bars in optimal design. Bottom: weld locations shown.
Fig. 4
Minimum-compliance (reference) design for the 2D cantilever beam example. Top: bars in optimal design. Bottom: weld locations shown.
Close modal

In practice, the weights for the terms in the manufacturing-cost function of (8) depend on many factors, including the size of the structure, the fabrication capabilities of the manufacturer, and the number of structures to be fabricated. It also depends on the type of structure; for instance, the cutting and preparation for a tubular frame for a racing car, which requires careful coping of the tubes, may be more expensive than the welding; on the other hand, the opposite is generally the case for a truss made of structural profiles joined to plate gussets. Here, we assign weights that are realistic for the manufacture of a steel frame truss whose dimensions are in inches. We assign wm = 1.3333 $/in3, w1 = 0.05 $/in3/2, w2 = 0.344 $/in2, w3 = 0.2143 $/in, and w4 = 0.0104 $/in2. These values are estimated based on some typical costs for preparation, cutting, welding, and painting in the United States. The manufacturing cost is minimized with these weights, and the resulting optimized design is shown in Fig. 5. We refer to this design as the baseline design.

The cost corresponding to each term in the manufacturing-cost function and other relevant measures for the reference and baseline designs are listed in Table 1. The comparison between the reference (minimum-compliance) and baseline (minimum-manufacturing-cost) designs clearly demonstrates that the manufacturing-cost function can effectively reduce every cost term compared to the reference design. The manufacturing cost of the baseline design is approximately 59% of that of the reference design, with a similar compliance C. We also observe that the manufacturing-cost-driven design employs significantly fewer bars (as indicated by κ in Table 1) and it substantially reduces the welding length lw and surface area As compared to the reference design.

Fig. 5
Minimal-manufacturing-cost (baseline) design for the 2D cantilever beam problem. Top: bars in optimal design. Bottom: weld locations shown.
Fig. 5
Minimal-manufacturing-cost (baseline) design for the 2D cantilever beam problem. Top: bars in optimal design. Bottom: weld locations shown.
Close modal
Table 1

Comparison between the reference, baseline, and weighted designs for the 2D cantilever problem

DesignMCmC1C2C3C4CgκAcutlwAs
Reference468.73239.993.749.51189.7225.7617.02170.2131.2327.65885.302476.97
Baseline278.28253.713.006.885.639.0417.10171.0618.9620.0026.29869.73
wm=5wm0289.09255.313.167.0713.0610.4617.10171.0020.9820.5660.981006.42
w1=5w10282.96253.463.286.739.569.9117.10171.0422.6619.5844.63953.62
w2=5w20274.30249.463.006.026.729.0717.10171.0119.3517.5131.39873.06
w3=5w30303.57275.773.078.176.689.8617.10171.0418.2723.7531.21948.19
w4=5w40304.29276.393.509.106.199.1017.10171.0023.6426.4628.91875.29
DesignMCmC1C2C3C4CgκAcutlwAs
Reference468.73239.993.749.51189.7225.7617.02170.2131.2327.65885.302476.97
Baseline278.28253.713.006.885.639.0417.10171.0618.9620.0026.29869.73
wm=5wm0289.09255.313.167.0713.0610.4617.10171.0020.9820.5660.981006.42
w1=5w10282.96253.463.286.739.569.9117.10171.0422.6619.5844.63953.62
w2=5w20274.30249.463.006.026.729.0717.10171.0119.3517.5131.39873.06
w3=5w30303.57275.773.078.176.689.8617.10171.0418.2723.7531.21948.19
w4=5w40304.29276.393.509.106.199.1017.10171.0023.6426.4628.91875.29

Note: The weights for the manufacturing-cost function in the baseline design are denoted by wi0.

We also note that there are co-linear bars in the optimal design of Fig. 5 that could be converted into a single bar (as long as they have the same radius), which would remove the weld between them and consequently decrease the welding cost. This could be achieved by a heuristic strategy that detects this situation and replaces the two bars with a single one, such as the one presented in Ref. [10]. However, generalizing this strategy to other geometric primitives, such as the plates shown in Sec. 6.3, is not straightforward. Therefore, this possibility is deferred to future work.

An interesting and important aspect of considering manufacturing cost as an optimization function is that we expect different supply chains (i.e., different values of the weights for the cost terms) to render different designs. In other words, the optimal design is a function of the supply chain. To demonstrate this, we repeat the manufacturing-cost-minimization with different weights for the terms of the cost function. For each optimization run, we assign the weight for one term to be five times the corresponding weight for the baseline design, keeping all other weights the same. The resulting designs are shown in the radar plot of Fig. 6 and the values of each cost terms and other measures are listed in Table 1.

Fig. 6
Radar plot of the modified manufacturing-cost-minimization results for the 2D cantilever beam problem, showing optimal designs (top) and the corresponding weld locations (bottom)
Fig. 6
Radar plot of the modified manufacturing-cost-minimization results for the 2D cantilever beam problem, showing optimal designs (top) and the corresponding weld locations (bottom)
Close modal

It is worth noting that all these designs exhibit the same structural performance, as the displacement constraint g is active in all these results. As expected, assigning different weights to the individual terms of the manufacturing-cost function renders different designs. Although the results corresponding to the increased weights for individual terms tend to render a lower value of the corresponding term, this is not necessarily the case because the entire manufacturing-cost function is minimized. In other words, if we view the manufacturing-cost-minimization problem as a multi-objective optimization, whereby each term of M is an objective, then we expect that each of the optimal designs with these modified weights corresponds to a non-dominated design in the Pareto front.

It should also be mentioned that the values of these weights cannot be arbitrary. If the weight for a particular term is too large, then the optimization can render designs that are poor in terms of their structural performance.

6.2 Three-Dimensional Cantilever Beam With Bar Primitives.

The second example is the design of a 3D cantilever beam with bar primitives. The initial design, depicted in Fig. 7, consists of 43 bars with a unit radius and membership variable of 0.5. The loading and boundary conditions are depicted in Fig. 8. To reduce the analysis time, we exploit the problem’s symmetry with respect to the xy plane and model only half of the domain. We note, however, that it is possible that an asymmetric design that is better than the design obtained with symmetry boundary conditions may be attained, as has been shown for frame structures in Refs. [21,22] and as is discussed in Ref. [23] for feature-mapping methods. The dimensions of the design domain are 60 × 10 × 5 and the mesh size is 180 × 30 × 15 elements. A load F = 0.1 is applied at the point (60, 0, 0). We bound the radius of each bar to [0.5, 1]. The move limit for this example is m = 0.05.

Fig. 7
Initial design for the 3D cantilever beam example consisting of 43 bars. The design is reflected to show the entire beam. Top: bars. Middle: iso-surface of combined density. Bottom: weld locations shown.
Fig. 7
Initial design for the 3D cantilever beam example consisting of 43 bars. The design is reflected to show the entire beam. Top: bars. Middle: iso-surface of combined density. Bottom: weld locations shown.
Close modal
Fig. 8
Rear view showing the loading and boundary conditions for 3D cantilever beam example. Only the half that is modeled is shown. The two square regions of size 2 × 2 indicate domains on which zero-displacement boundary conditions in all directions are imposed. A symmetry boundary condition is imposed on the x–y plane.
Fig. 8
Rear view showing the loading and boundary conditions for 3D cantilever beam example. Only the half that is modeled is shown. The two square regions of size 2 × 2 indicate domains on which zero-displacement boundary conditions in all directions are imposed. A symmetry boundary condition is imposed on the x–y plane.
Close modal

Similar to the previous example, we minimize the compliance subject to a maximum material volume fraction constraint limit of v¯=0.15. As before, we use the displacement at the point of application of the load for this minimum-compliance design to define the maximum allowable displacement u¯ for the displacement constraint of the manufacturing-cost-minimization problem. The final design and its corresponding welding locations are shown in Fig. 9. The values of different terms for the optimal design are listed in Table 2.

Fig. 9
Minimum-compliance (reference) design for the 3D cantilever beam example. The design is reflected to show entire beam. Top: bars. Middle: iso-surface of combined density. Bottom: weld locations shown.
Fig. 9
Minimum-compliance (reference) design for the 3D cantilever beam example. The design is reflected to show entire beam. Top: bars. Middle: iso-surface of combined density. Bottom: weld locations shown.
Close modal
Table 2

Comparison between the minimum-compliance and minimum-manufacturing-cost designs for the 3D cantilever beam example with 3D bar primitives

DesignMCmC1C2C3C4CgκAcutlwAs
Min C1214.19600.616.2822.83550.9933.468.4385.0835.1166.382571.123217.61
Min M766.18656.745.8521.2164.5217.838.4284.9927.8461.67301.111715.17
DesignMCmC1C2C3C4CgκAcutlwAs
Min C1214.19600.616.2822.83550.9933.468.4385.0835.1166.382571.123217.61
Min M766.18656.745.8521.2164.5217.838.4284.9927.8461.67301.111715.17

We perform the manufacturing-cost-minimization, using the same weights as the previous example. The final design is shown in Fig. 10. The cost corresponding to each term in the manufacturing cost and other relevant measures for the reference design and the minimum-manufacturing-cost design are listed in Table 2. The history of the manufacturing-cost objective and displacement constraint for this problem are shown in Fig. 11.

Fig. 10
Minimum-manufacturing-cost design for the 3D cantilever beam example. The design is reflected to show the entire beam. Top: bars. Middle: iso-surface of combined density. Bottom: weld locations shown.
Fig. 10
Minimum-manufacturing-cost design for the 3D cantilever beam example. The design is reflected to show the entire beam. Top: bars. Middle: iso-surface of combined density. Bottom: weld locations shown.
Close modal
Fig. 11
Objective function (scaled) and constraint value histories for 3D cantilever beam example
Fig. 11
Objective function (scaled) and constraint value histories for 3D cantilever beam example
Close modal

As shown in Table 2, the material cost is slightly increased in the minimum-manufacturing-cost design, but the fabrications costs are drastically decreased. Also, the number of existing components κ decreases from 35 to 27. For this example, the cost of welding dominates the optimization. As can be observed in Fig. 9, there is substantial overlap between bars in the minimum-compliance design, leading to a large weld length. In contrast, the weld length of the minimum-manufacturing-cost design of Fig. 10 is substantially shorter. The history plots of Fig. 11 show that the optimization exhibits good convergence towards the optimum. This example demonstrates that the proposed methodology can be readily applied to 3D bar structures.

6.3 Three-Dimensional Messerschmitt–Bölkow–Blohm Beam With Plate Primitives.

Our last example corresponds to an MBB beam design with the loading and supports depicted in Fig. 12. This example demonstrates the proposed method that can be readily applied with plate geometric primitives—something that cannot be easily done with ground structure approaches. Similar to the 3D cantilever beam example, we take advantage of symmetry and only model a quarter of the design region, with dimensions 60 × 10 × 5. The mesh size is 180 × 30 × 15 elements. The corners of the design domain are fixed only in the y-direction. The load of magnitude F = 0.1 is applied at the top center point along the negative z-direction. The initial design, consisting of 18 plates of dimension 5 × 5 with a fixed (non-designable) semi-thickness rc = 1.2R = 0.6928, is depicted in Fig. 13. The membership variable of all plates in the initial design is 0.5. The lengths of the plates are bounded to [0, 20]. The move limit for this example is m = 0.05. It is worth noting that a large plate semi-thickness rc may cause the optimizer to render designs in which plates collapse into bars (i.e., one of the dimensions of the rectangular medial surface approaches zero), as this minimizes the surface area of the structure.

Fig. 12
Rear view showing the loading and boundary conditions for the 3D MBB beam example with plate primitives. Only the quarter region that is modeled in the analysis is shown. Symmetry boundary conditions are imposed on the x–y plane and y–z plane.
Fig. 12
Rear view showing the loading and boundary conditions for the 3D MBB beam example with plate primitives. Only the quarter region that is modeled in the analysis is shown. Symmetry boundary conditions are imposed on the x–y plane and y–z plane.
Close modal
Fig. 13
Initial design of the 3D MBB beam example with 18 plate primitives (shown here reflected twice)
Fig. 13
Initial design of the 3D MBB beam example with 18 plate primitives (shown here reflected twice)
Close modal

As before, we first solve the reference minimum-compliance problem and set the corresponding displacement of the point where the load located as the displacement constraint limit for the minimum-manufacturing-cost problem. The weight of material cost for plate primitives is set to wm=1.45$/in3 based on a typical unit price of plate stock material in United States. The minimum-compliance and minimum-manufacturing-cost designs are shown in Figs. 14 and 15, and the corresponding cost of each term in M and other relevant measures are listed in Table 3.

Fig. 14
Minimum-compliance (reference) design for the 3D MBB beam example with plates. The design is reflected twice to show the entire beam. Top: plates. Middle: iso-surface of combined density. Bottom: weld locations shown.
Fig. 14
Minimum-compliance (reference) design for the 3D MBB beam example with plates. The design is reflected twice to show the entire beam. Top: plates. Middle: iso-surface of combined density. Bottom: weld locations shown.
Close modal
Fig. 15
Minimum-manufacturing-cost design for the 3D MBB beam example with plates. The design is reflected twice to show the entire beam. Top: plates. Middle: iso-surface of combined density. Bottom: weld locations shown.
Fig. 15
Minimum-manufacturing-cost design for the 3D MBB beam example with plates. The design is reflected twice to show the entire beam. Top: plates. Middle: iso-surface of combined density. Bottom: weld locations shown.
Close modal
Table 3

Comparison between the minimum-compliance and minimum-manufacturing-cost designs for the 3D MBB beam example with 3D plate primitives

DesignMCmC1C2C3C4CgκAcutlwAs
Min C5631.83599.444.494841.67180.146.0710.4799.0717.9714074.64840.60584.47
Min M3584.14849.414.952682.2343.913.6110.02100.2616.777797.20204.92347.73
DesignMCmC1C2C3C4CgκAcutlwAs
Min C5631.83599.444.494841.67180.146.0710.4799.0717.9714074.64840.60584.47
Min M3584.14849.414.952682.2343.913.6110.02100.2616.777797.20204.92347.73

As we can see from Table 3, the manufacturing cost of the minimum-manufacturing-cost design is only 63.54% of the cost of the reference design. The cost of cutting dominates the manufacturing cost of structure constructed with plates. From Fig. 15, we observe that the optimal design looks like a boxed beam, which helps reduce the cutting cost by reducing the perimeter of the plates.

To visualize designs with their welding locations for the two 3D problems, we employ the visualization software paraview [24,25]. The structure in the figures corresponds to a combined density value above 0.4, and the welding locations correspond to a value of αiαjlij above 0.05.

7 Conclusions

This work introduced a geometry projection technique for the topology optimization of welded frame structures made of bar or plate components with regard to manufacturing cost. The examples demonstrate the effectiveness of the proposed method for 2D and 3D problems. In particular, the manufacturing-cost-driven designs attain a manufacturing cost that is significantly lower than that of the minimum-compliance designs.

It is important to recall that the manufacturing-cost function employed in this work is a rough estimate to be used in the concept design stage, as an accurate estimation of cost requires a significant amount of additional information that is only available when a detailed design of the structure is available. Nevertheless, the proposed method is useful in rendering designs that are more manufacturing-cost effective and is a tool to incorporate manufacturing-cost considerations early in the design. In particular, this technique makes it possible to explore different topologies that may be better for different manufacturing supply chains with different cost structures, i.e., for which different components of the manufacturing cost have different costs.

The dual representation of high-level geometric parameters and projected densities for each primitive makes it possible to express the components of the manufacturing-cost function in terms of design variables, something which is not possible with density-based and level-set representations of the entire structure. The geometry projection method, in which each individual primitive is first mapped onto its own density field and all component densities are subsequently combined for the material interpolation, makes it possible to capture, for instance, the welding length between components. Moreover, the geometry projection enables the analysis with a fixed grid throughout the optimization, and the components need not be connected in a predefined ground structure. Since the geometry projection is differentiable, it is possible to use efficient gradient-based methods to perform the optimization.

We note that, as in any gradient-based technique, the proposed method will in general converge to a local minimum. While we did not observe poor local minima in our experiments, it is possible that they may occur, particularly since the design representation is more compact, as demonstrated in Ref. [23]. Finally, we note that stress constraints and other structural criteria are important in the design of frame structures and are deferred to future work.

Acknowledgment

The authors thank the United States National Science Foundation, Award CMMI-1751211, for support of this work. We are also grateful to Prof. Krister Svanberg for providing his MMA matlab optimizer to perform the optimization.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article may be obtained from the corresponding author upon reasonable request.

References

1.
Asadpoure
,
A.
,
Guest
,
J. K.
, and
Valdevit
,
L.
,
2015
, “
Incorporating Fabrication Cost Into Topology Optimization of Discrete Structures and Lattices
,”
Struct. Multidiscipl. Optim.
,
51
(
2
), pp.
385
396
.
2.
Liu
,
J.
,
Chen
,
Q.
,
Liang
,
X.
, and
To
,
A. C.
,
2019
, “
Manufacturing Cost Constrained Topology Optimization for Additive Manufacturing
,”
Front. Mech. Eng.
,
14
(
2
), pp.
213
221
.
3.
Ryan
,
L.
, and
Kim
,
I. Y.
,
2019
, “
A Multiobjective Topology Optimization Approach for Cost and Time Minimization in Additive Manufacturing
,”
Int. J. Numer. Methods Eng.
,
118
(
7
), pp.
371
394
.
4.
Sabiston
,
G.
, and
Kim
,
I. Y.
,
2020
, “
3D Topology Optimization for Cost and Time Minimization in Additive Manufacturing
,”
Struct. Multidiscipl. Optim.
,
61
(
2
), pp.
731
748
.
5.
Ulu
,
E.
,
Huang
,
R.
,
Kara
,
L. B.
, and
Whitefoot
,
K. S.
,
2019
, “
Concurrent Structure and Process Optimization for Minimum Cost Metal Additive Manufacturing
,”
ASME J. Mech. Des.
,
141
(
6
), p.
061701
.
6.
Zhou
,
Z.
,
Hamza
,
K.
, and
Saitou
,
K.
,
2014
, “
Decomposition Templates and Joint Morphing Operators for Genetic Algorithm Optimization of Multicomponent Structural Topology
,”
ASME J. Mech. Des.
,
136
(
2
), p.
021004
.
7.
Zhou
,
Y.
, and
Saitou
,
K.
,
2018
, “
Gradient-Based Multi-component Topology Optimization for Stamped Sheet Metal Assemblies (MTO-S)
,”
Struct. Multidiscipl. Optim.
,
58
(
1
), pp.
83
94
.
8.
Li
,
G.
,
Xu
,
F.
,
Huang
,
X.
, and
Sun
,
G.
,
2015
, “
Topology Optimization of an Automotive Tailor-Welded Blank Door
,”
ASME J. Mech. Des.
,
137
(
5
), p.
055001
.
9.
Torii
,
A. J.
,
Lopez
,
R. H.
, and
Miguel
,
L. F. F.
,
2016
, “
Design Complexity Control in Truss Optimization
,”
Struct. Multidiscipl. Optim.
,
54
(
2
), pp.
289
299
.
10.
Zhang
,
W.
,
Zhou
,
J.
,
Zhu
,
Y.
, and
Guo
,
X.
,
2017
, “
Structural Complexity Control in Topology Optimization Via Moving Morphable Component (MMC) Approach
,”
Struct. Multidiscipl. Optim.
,
56
(
3
), pp.
535
552
.
11.
Zhang
,
W.
,
Liu
,
Y.
,
Wei
,
P.
,
Zhu
,
Y.
, and
Guo
,
X.
,
2017
, “
Explicit Control of Structural Complexity in Topology Optimization
,”
Comput. Methods Appl. Mech. Eng.
,
324
, pp.
149
169
.
12.
Norato
,
J.
,
Bell
,
B.
, and
Tortorelli
,
D. A.
,
2015
, “
A Geometry Projection Method for Continuum-Based Topology Optimization With Discrete Elements
,”
Comput. Methods Appl. Mech. Eng.
,
293
, pp.
306
327
.
13.
Zhang
,
S.
,
Norato
,
J. A.
,
Gain
,
A. L.
, and
Lyu
,
N.
,
2016
, “
A Geometry Projection Method for the Topology Optimization of Plate Structures
,”
Struct. Multidiscipl. Optim.
,
54
(
5
), pp.
1173
1190
.
14.
Smith
,
H.
, and
Norato
,
J.
,
2022
, “
Topology Optimization of Structures Made of Fiber-Reinforced Plates
,”
Struct. Multidiscipl. Optim.
,
65
(
2
), pp.
1
18
.
15.
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
.
16.
Smith
,
H.
, and
Norato
,
J. A.
,
2021
, “
Topology Optimization With Discrete Geometric Components Made of Composite Materials
,”
Comput. Methods Appl. Mech. Eng.
,
376
, p.
113582
.
17.
Bendsøe
,
M. P.
, and
Sigmund
,
O.
,
1999
, “
Material Interpolation Schemes in Topology Optimization
,”
Arch. Appl. Mech.
,
69
(
9
), pp.
635
654
.
18.
Shapiro
,
V.
,
2002
, “Solid Modeling,”
Handbook of Computer Aided Geometric Design
,
G.
Farin
,
J.
Hoschek
, and
M. S.
Kim
, eds.,
Elsevier
,
Amsterdam
, pp.
473
518
.
19.
Jármai
,
K.
, and
Farkas
,
J.
,
1999
, “
Cost Calculation and Optimisation of Welded Steel Structures
,”
J. Construct. Steel Res.
,
50
(
2
), pp.
115
135
.
20.
Svanberg
,
K.
,
2007
,
MMA and GCMMA – Two Methods for Nonlinear Optimization
.
21.
Stolpe
,
M.
,
2010
, “
On Some Fundamental Properties of Structural Topology Optimization Problems
,”
Struct. Multidiscipl. Optim.
,
41
(
5
), pp.
661
670
.
22.
Rozvany
,
G. I.
,
2011
, “
On Symmetry and Non-uniqueness in Exact Topology Optimization
,”
Struct. Multidiscipl. Optim.
,
43
(
3
), pp.
297
317
.
23.
Wein
,
F.
,
Dunning
,
P. D.
, and
Norato
,
J. A.
,
2020
, “
A Review on Feature-Mapping Methods for Structural Optimization
,”
Struct. Multidiscipl. Optim.
,
62
(
4
), pp.
1597
1638
.
24.
Ahrens
,
James
,
Geveci
,
Berk
, and
Law
,
Charles
,
2005
, “ParaView: An End-User Tool for Large Data Visualization,”
The Visualization Handbook
, Vol.
717
,
C.
Hansen
,
C.
Johnson
, eds.,
Elsevier
,
Burlington-Oxford
, pp.
717
732
.
25.
Ayachit
,
U.
,
2015
,
The Paraview Guide: A Parallel Visualization Application
,
Kitware, Inc.
,
Clifton Park, NY
.