Designing modular 3D printed reinforcement of wound composite hollow beams with semidefinite programming

TL;DR

This paper presents a convex semidefinite programming approach to design modular 3D-printed reinforcement structures within composite hollow beams, significantly improving their vibrational stability and manufacturability.

Contribution

It introduces a novel optimization formulation for internal lattice structures that enhance stability and manufacturability of composite beams using 3D printing.

Findings

Fundamental free-vibration eigenfrequency nearly doubled.

Internal reinforcement improves structural stability under shear and buckling.

Design process is fully automated and validated through prototype testing.

Abstract

Fueled by their excellent stiffness-to-weight ratio and the availability of mature manufacturing technologies, filament wound carbon fiber reinforced polymers represent ideal materials for thin-walled laminate structures. However, their strong anisotropy reduces structural resistance to wall instabilities under shear and buckling. Increasing laminate thickness degrades weight and structural efficiencies and the application of a dense internal core is often uneconomical and labor-intensive. In this contribution, we introduce a convex linear semidefinite programming formulation for truss topology optimization to design an efficient non-uniform lattice-like internal structure. The internal structure not only reduces the effect of wall instabilities, mirrored in the increase of the fundamental free-vibration eigenfrequency, but also keeps weight low, secures manufacturability using conventional three-dimensional printers, and withstands the loads induced during the production process. We showcase a fully-automatic procedure in detail for the design, prototype manufacturing, and verification of a simply-supported composite machine tool component, including validation with roving hammer tests. The results confirm that the 3D-printed optimized internal structure almost doubles the fundamental free-vibration eigenfrequency, allowing to increase working frequency of the machine tool, even though the ratio between elastic properties of the carbon composite and the ABS polymer used for 3D printing exceeds two orders of magnitude.

Figures

Tables

Table 1. Table 1 : Material properties of the wound composite beam laminae. E 1 subscript 𝐸 1 E_{1} and E 2 subscript 𝐸 2 E_{2} stand for the Young moduli in the fiber and transverse directions, respectively; G 12 subscript 𝐺 12 G_{12} denotes the shear modulus, ν 12 subscript 𝜈 12 \nu_{12} and ν 23 subscript 𝜈 23 \nu_{23} are Poisson’s ratios. θ 𝜃 \theta constitutes the angle between the 1 1 1 -direction and x 𝑥 x , rotating around the beam surface normals. Finally, ρ 𝜌 \rho and t 𝑡 t denote the density and thickness of the plies.

Layer	$E_1$ [GPa]	$E_2$ [GPa]	$G_{12}$ [GPa]	${\nu }_{12}$ [-]	${\nu }_{23}$ [-]	$\theta$ [$\mathrm{deg}$]	$\rho$ [kg/m3]	$t$ [mm]
$1$	$128.2$	$5.0$	$3.4$	$0.34$	$0.35$	$89.3$	$1,428$	$0.25$
$2$	$421.9$	$3.7$	$3.2$	$0.37$	$0.35$	$0.0$	$1,680$	$1.25$
$3$	$130.9$	$5.0$	$3.4$	$0.34$	$0.35$	$26.9$	$1,458$	$0.18$
$4$	$130.9$	$5.0$	$3.4$	$0.34$	$0.35$	$-26.9$	$1,458$	$0.36$
$5$	$130.9$	$5.0$	$3.4$	$0.34$	$0.35$	$26.9$	$1,458$	$0.18$
$6$ (casing)	$2.0$	$2.0$	$0.7$	$0.37$	$0.37$	$0.0$	$1,040$	$0.80$

Table 2. Table 2 : Comparison of model prediction of eigenfrequencies f FEM subscript 𝑓 FEM f_{\mathrm{FEM}} and measured natural frequencies f EXP subscript 𝑓 EXP f_{\mathrm{EXP}} using the roving hammer test. Accuracy of individual measurements is denoted by A 𝐴 A and the deviation of the model from the experiment by D 𝐷 D .

Eigenmode	$f_{\mathrm{FEM}}$ [Hz]	$f_{\mathrm{EXP}}$ [Hz]	$A$ [Hz]	$D$ [$\%$]
First shear	$600.1$	$658$	$2$	$-8.8$
First bending $y$	$747.0$	$714$	$2$	$+4.6$
First bending $z$	$748.2$	$724$	$2$	$+3.3$
Third buckling	$682.3$	$846$	$4$	$-19.3$
First torsion	$833.7$	$864$	$4$	$-3.5$
Fourth buckling	$708.9$	$896$	$4$	$-20.9$
Fifth buckling	$741.4$	$942$	$6$	$-21.3$
Sixth buckling	$790.6$	$1004$	$6$	$-21.3$
Second bending $z$	$1020.0$	$1102$	$6$	$-7.4$

Equations

$$\underset{\mathbf{x}}{\min }$$

$$s.t.$$

$$\mathbf{X}⪰0,$$

$$\underset{\mathbf{a},{\mathbf{u}}_{\mathrm{cm}},\mathbf{u}}{\min }$$

$$s.t.$$

$$\left({\mathbf{K}}_{\mathrm{cm}}^{\mathrm{IS}}(\mathbf{a})+{\mathbf{K}}_{\mathrm{cm}}^{\mathrm{C}}\right){\mathbf{u}}_{\mathrm{cm}}$$

$$-\overline{u}1\le {\mathbf{u}}_{\mathrm{cm},\mathrm{disp}}$$

$$0\le \mathbf{a}$$

$$\overline{\lambda }=4{\pi }^2{\overline{f}}^2.$$

$${\mathbf{K}}_{\mathrm{fv}}^{\mathrm{IS}}(\mathbf{a})=\sum _{e=1}^{n_{\mathrm{b}}}{\hat{\mathbf{K}}}_{\mathrm{fv},e}^{\mathrm{IS}}{a}_e,{\mathbf{M}}_{\mathrm{fv}}^{\mathrm{IS}}(\mathbf{a})=\sum _{e=1}^{n_{\mathrm{b}}}{\hat{\mathbf{M}}}_{\mathrm{fv},e}^{\mathrm{IS}}{a}_e,$$

$${\mathbf{K}}_{\mathrm{fv}}^{\mathrm{IS}}(\mathbf{a})+{\mathbf{K}}_{\mathrm{fv}}^{\mathrm{C}}-4{\pi }^2{\overline{f}}^2\left({\mathbf{M}}_{\mathrm{fv}}^{\mathrm{IS}}(\mathbf{a})+{\mathbf{M}}_{\mathrm{fv}}^{\mathrm{C}}\right)⪰0,$$

$$\left(\begin{array}{cc}{c}_{\mathrm{cm}}& -{\mathbf{f}}_{\mathrm{cm}}^{\mathrm{T}}\\ -{\mathbf{f}}_{\mathrm{cm}}& {\mathbf{K}}_{\mathrm{cm}}^{\mathrm{IS}}(\mathbf{a})+{\mathbf{K}}_{\mathrm{cm}}^{\mathrm{C}}\end{array}\right)⪰0,$$

$$c_{\mathrm{cm}}=c_{\mathrm{cm},0}\frac{\overline{u}}{\max \left\{|{\mathbf{u}}_{\mathrm{cm},\mathrm{disp}}|\right\}},$$

$$c_{\mathrm{cm},0}={\tilde{\mathbf{f}}}_{\mathrm{cm}}^{\mathrm{T}}{\left({\tilde{\mathbf{K}}}_{\mathrm{cm}}^{\mathrm{C}}\right)}^{-1}{\tilde{\mathbf{f}}}_{\mathrm{cm}}.$$

$$\underset{\mathbf{a}}{\min }$$

$$s.t.$$

$$\left(\begin{array}{cc}{c}_{\mathrm{cm}}& -{\mathbf{f}}_{\mathrm{cm}}^{\mathrm{T}}\\ -{\mathbf{f}}_{\mathrm{cm}}& {\mathbf{K}}_{\mathrm{cm}}^{\mathrm{IS}}(\mathbf{a})+{\mathbf{K}}_{\mathrm{cm}}^{\mathrm{C}}\end{array}\right)$$

$$1\overline{a}\ge \mathbf{a}$$

$$[box=]\mathbf{K}(\mathbf{a})\mathbf{u}=\mathbf{f}$$

$$[box=]\left(\begin{array}{cc}{\mathbf{K}}_a(\mathbf{a})& {\mathbf{K}}_b\\ {\mathbf{K}}_b^{\mathrm{T}}& {\mathbf{K}}_c\end{array}\right)\left(\begin{array}{c}{\mathbf{u}}_a\\ {\mathbf{u}}_b\end{array}\right)=\left(\begin{array}{c}{\mathbf{f}}_a\\ {\mathbf{f}}_b\end{array}\right),$$

$$[box=]{\mathbf{u}}_b={\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{f}}_b-{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{K}}_b^{\mathrm{T}}{\mathbf{u}}_a$$

$$[box=]\left[{\mathbf{K}}_a(\mathbf{a})-{\mathbf{K}}_b{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{K}}_b^{\mathrm{T}}\right]{\mathbf{u}}_a={\mathbf{f}}_a-{\mathbf{K}}_b{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{f}}_b.$$

$$[box=]c={\mathbf{u}}_a^{\mathrm{T}}{\mathbf{f}}_a+{\mathbf{u}}_b^{\mathrm{T}}{\mathbf{f}}_b.$$

$$[box=]c={\mathbf{u}}_a^{\mathrm{T}}\left[{\mathbf{f}}_a-{\mathbf{K}}_b{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{f}}_b\right]+{\mathbf{f}}_b^{\mathrm{T}}{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{f}}_b,$$

$$[box=]\left(\begin{array}{cc}c& -{\mathbf{f}}^{\mathrm{T}}\\ -\mathbf{f}& \mathbf{K}(\mathbf{a})\end{array}\right)⪰0$$

$$[box=]\left(\begin{array}{cc}c-{\mathbf{f}}_b^{\mathrm{T}}{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{f}}_b& -{\mathbf{f}}_a^{\mathrm{T}}+{\mathbf{f}}_b^{\mathrm{T}}{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{K}}_b^{\mathrm{T}}\\ -{\mathbf{f}}_a+{\mathbf{K}}_b{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{f}}_b& {\mathbf{K}}_a(\mathbf{a})-{\mathbf{K}}_b{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{K}}_b^{\mathrm{T}}\end{array}\right)⪰0.$$

$$[box=]\begin{array}{r}{\mathbf{K}}_a(\mathbf{a})-{\mathbf{K}}_b{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{K}}_b^{\mathrm{T}}-\left(-{\mathbf{f}}_a^{\mathrm{T}}+{\mathbf{f}}_b^{\mathrm{T}}{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{K}}_b^{\mathrm{T}}\right)\\ {\left(c-{\mathbf{f}}_b^{\mathrm{T}}{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{f}}_b\right)}^{-1}\left(-{\mathbf{f}}_a+{\mathbf{K}}_b{\left({\mathbf{K}}_c\right)}^{-1}{\mathbf{f}}_b\right)⪰0.\end{array}$$

$$[box=]\left(\begin{array}{cc}{\mathbf{K}}_a(\mathbf{a})-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_a(\mathbf{a})& {\mathbf{K}}_b-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_b\\ {\mathbf{K}}_b^{\mathrm{T}}-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_b^{\mathrm{T}}& {\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\end{array}\right)⪰0.$$

$$[box=]{\mathbf{K}}_c{\mathbf{u}}_b-\lambda {\mathbf{M}}_c{\mathbf{u}}_b=0,$$

$$[box=]\begin{array}{r}{\mathbf{K}}_a(\mathbf{a})-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_a(\mathbf{a})-\left({\mathbf{K}}_b-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_b\right)\\ {\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right)}^{-1}\left({\mathbf{K}}_b^{\mathrm{T}}-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_b^{\mathrm{T}}\right)⪰0.\end{array}$$

$$[box=]\begin{array}{r}{\mathbf{K}}_a(\mathbf{a})-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_a(\mathbf{a})-\left({\mathbf{K}}_b-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_b\right)\\ {\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right)}^{†}\left({\mathbf{K}}_b^{\mathrm{T}}-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_b^{\mathrm{T}}\right)⪰0,\end{array}$$

$$\begin{array}{r}\left[\mathbf{I}-\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right){\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right)}^{†}\right]\\ \left({\mathbf{K}}_b^{\mathrm{T}}-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_b^{\mathrm{T}}\right)=0,\end{array}$$

$$[box=]\begin{array}{r}[\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right)-\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right)\\ {\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right)}^{†}\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right)]\mathbf{C}=0.\end{array}$$

$$[box=]\text{span}{\left\{{\mathbf{u}}_b:\left({\mathbf{K}}_c-4{\pi }^2{\overline{f}}^2{\mathbf{M}}_c\right){\mathbf{u}}_b=0\right\}}^{\perp }.$$

Full text

††journal: arXiv.org

1 Introduction

Offering excellent stiffness-to-weight ratios, high damping, and a low sensitivity to fatigue and corrosion, carbon fiber reinforced polymers (CFRPs) are employed in high-tech applications, including bodies of racing cars in automotive [1], propulsors and turbines in naval [2, 3], wing boxes and ailerons in aerospace [4], and rocket bodies in the space industry [5]. Considerable attention has been therefore paid to methods for optimizing structural performance of these laminates, particularly the laminate layup [6, 7]. To concurrently maximize bending stiffness and keep weight low, the outer dimensions of these structures tend to be maximized, while the wall thickness is minimized. The “thin-walledness” of the resulting structures, combined with their anisotropy, renders them highly sensitive to shear and wall buckling instabilities manifested in low fundamental free-vibration eigenfrequencies.

Below, we review common approaches to topology optimization that reduce wall instabilities by designing an internal structure, Section 1.1. Section 1.2 provides a brief introduction to semidefinite programming and highlights several applications in structural optimization. Finally, Section 1.3 reveals the merits of designing internal structures using semidefinite programming.

1.1 Topology optimization

Topology optimization techniques [8] provide the means for reducing wall instabilities when designing sufficiently stiff yet lightweight structures. In the simplest setting—beam cross section optimization—we search an optimal two-dimensional cross-sectional shape or a stiffening of structures whose outer shape is predefined [9]. Blasques [10], for example, maximized the fundamental free-vibration eigenfrequency while accounting for the mass and shear center position constraints; Nguyen et al. [11] optimized cross-sections of prismatic beams to maximize their buckling loads.

The design of optimal core sandwich structures, whose skins are stiffened by a thick core is a related challenge. For honeycomb, solid, truss, and foam rectangular panels under in-plane compression or shear loads, optimal periodic topologies can be found analytically by considering the optimality criterium of all failure modes occurring simultaneously [12]. For complex boundary conditions, parametric shape-optimization studies are usually performed. Wang and McDowell [13] studied the geometry of a metal honeycomb sandwich beam core under torsion and bending and Xu and Qiu [14] optimized the lattice core of a composite sandwich panel to increase the fundamental eigenfrequency while accounting for uncertainties in the model. They concluded that bending eigenfrequencies increase with increasing strut thicknesses, with an increase in the elastic and shear modulus of the composite, and with a decrease in density. Although Daynes et al. [15] optimized spatially-graded lattice structures within a single sandwich panel domain, surprisingly, almost no prior research seems to have stepped beyond parametric intuition-based designs [16, 17], the rare exception being the multi-scale topology optimization approach investigated by Coelho and Rodrigues [18].

Questioning whether, where, and how to stiffen already engineered designs in order to further improve their structural performance constitutes the central question of the reinforcement problem [19, 20], superseding the former dimensional reduction and periodicity assumptions. Initial studies in this area have considered maximization of the fundamental eigenfrequency [20] and improving the structural frequency response of plane elastic structures [21] using the homogenization and optimality criteria methods, respectively.

Using the ground structure approach for topology optimization of truss structures, Bendsøe et al. [22] fixed cross-sectional areas of a set of bars and searched for their stiffest truss reinforcement, a (non-smooth) convex quadratic programming formulation. Alternatively, the effect of a fixed boundary structure has been approximated by an appropriate application of nodal forces to the ground structure [23, 24], but this choice influences, however, the optimized design.

In the setting of continuous topology optimization, Luo and Gea [25] developed a systematic optimization approach for the topology and orientation design of composite stiffeners of plates and shells in both static and dynamic settings, and Wang et al. [26] optimized the overall structural rigidity of an automobile body through a maximization of the fundamental eigenfrequency. In aerospace applications, Maute and Allen [27] optimized a wing’s internal structure, subjected to fluid-surface interactions; Aage et al. [28] performed an extremely large-scale optimization of the internal structure of a Boeing 777 wing, while avoiding the traditional rib and spar designs [29]. In military applications, topology optimization was the basis for the design of additively-manufactured lattice-reinforced penetrative warheads [30] and for optimizing the layout weight of stiffeners in composite submarines subjected to nonsymmetric wave slap loads [31].

Other methods relevant to internal structure design have arisen in conjunction with recently introduced coating and infill optimization problems. Clausen et al. [32] developed a formulation for the optimization of (uniformly) coated structures, wherein a base material, infill, was surrounded by another material at the interfaces, finding a porous, complex infill significantly improves both structural buckling resistance and robustness to local perturbations when compared to optimized solid structures of equal weight and similar stiffnesses [33, 34]. In three dimensions, optimized designs further exploit the merits of closed shell surfaces through the sandwich effect [34].

Inspired by natural, bone-like microstructures, Wu et al. [35] optimized a spatially non-uniform porous infill, Wang et al. [36] developed a sequential approach for generating graded lattice mesostructures, and Zhu et al. [37] introduced a novel asymptotic-analysis-based homogenization approach. All these methods automatically design stiff yet porous infills for additive manufacturing products while superseding the traditional pattern-based designs [38]. Finally, Wu et al. [39] extended their approach to the ultimate setting of a concurrent optimization of coated structures and porous infills, and Groen et al. [40] have developed a homogenization-based method to accelerate solutions.

1.2 Semidefinite programming

It has been shown in recent decades that several structural optimization problems can be modeled as semidefinite programs. Linear semidefinite programming (SDP) is a subset of convex optimization of the form

[TABLE]

and involves minimization of a linear function (1a) over a spectrahedron, which is an intersection of an affine space (1b) with the cone of symmetric positive semidefinite matrices (1c). In (1c), the notation “$\succeq\mathbf{0}$” enforces positive semidefiniteness of the left hand side. Due to the linear dependence of $\mathbf{X}$ on $\mathbf{x}$ (1b), (1c) is commonly referred to as a linear matrix inequality (LMI).

Applications of semidefinite programming to structural design were pioneered by Ben-Tal and Nemirovski [41], de Klerk et al. [42], and Vandenberghe and Boyd [43] who developed formulations for minimum-compliance and weight truss topology optimizations. The main added value of SDP lies in its ability to effectively avoid the non-differentiability of multiple eigenvalues for free-vibrations [44, 45] and buckling [46, 47], robust optimization [48], and bounds improvement for optimization problems in a discrete setting [49]. Semidefinite programming has also found applications in optimal materials design, the Free Material Optimization approach [50], or in the limit analyses [51].

1.3 Aims and novelty

In this contribution, we consider an industrial problem of designing the least-weight internal structure of a thin-walled filament-wound composite machine tool component prone to shear and buckling wall instabilities. The beam laminate was designed for bearing dynamic loads, allowing us to describe the wall instabilities naturally in terms of free-vibrations eigenfrequencies.

In current production process, the wall instabilities are reduced by inserting a uniform foam core structure into the beam interior, an uneconomical and labor-intensive process. Conversely, we have aimed to automatically design a structurally-efficient internal structure which can easily be manufactured using conventional low-cost $3$D printers.

To this goal, we extended the convex (linear) semidefinite programming formulation introduced by Ohsaki et al. [44] and Ben-Tal and Nemirovski [48] to design globally-optimal least-weight lattice-like internal structures and apply it to increasing the fundamental eigenfrequency and decreasing the compression-molding compliance of a thin-walled composite beam prototype. Note that Achtziger and Kočvara [45] avoided prescribed structural elements but allowed for a non-structural mass and Ohsaki et al. [44] did not consider prescribed mass or stiffness.

After introducing the case study of a simply-supported CFRP beam design in Section 2, we develop its finite element representation in Section 3.1. For this representation, a semidefinite programming formulation for truss topology optimization of internal structures is developed in Section 3.2. Having designed the optimal internal structure, we post-process the optimization outputs and export, in a fully-automated way, the internal structure for additive manufacturing in Section 3.3. During manufacturing, the internal structure serves as the support for carbon fibers in the filament-winding production phase, and a prototype is created. Section 4 describes verification and experimental validation of the prototype and concludes that its response agreed well with the model prediction.

2 Case study

As the basic structure, we consider a prismatic, laminated composite beam $1,000$ mm long, with a $80\times 80$ mm thin-walled cross-section $2.2$ mm thick, Fig. 1(b). According to current manufacturing technology, beam production consists of several steps, in which a supporting structure made of manually processed high-density foam is wound biaxially with a combination of ultra high modulus (UHM) and high modulus (HM) carbon fibers saturated with epoxy resin. The supporting structure prevents cross-section distortions induced by compression-molding loads as shown in Fig. 1(c). Subsequently, the beam is cured, the supporting structure is pulled out, and the beam outer surface is finalized.

The final product is exposed primarily to loads that induce bending. For this purpose, most of the carbon fibers are aligned with the beam’s longitudinal axis (layer $2$ in Table 1), denoted by $x$ in Fig. 1(a), whereas the remaining layers reduce the susceptibility to delamination. See Table 1, where all layers are listed by their orientations relative to the beam’s longitudinal axis, $\theta$. This layered composition reliably transmits the design forces to the supports, and is thus fully sufficient in this sense.

Attributed to transversely isotropic material properties, the beam’s walls are, however, prone to elastic wall instabilities under shear and buckling, which also manifests in free-vibration modes and frequencies of the non-reinforced beam. Figure 2 confirms that the first fundamental eigenmode with a frequency of $128.5$ Hz corresponds to shear wall instabilities, whereas the second eigenmode combines bending with buckling; all higher eigenmodes (not shown) exhibit similar wall instabilities. Because the fundamental eigenfrequency limits the maximum working frequency of the machine part, its increase is of considerable interest.

Although the effect of these instabilities can be reduced by additional laminate layers or by also keeping the uniform foam structure for operational loads, the added weight, decrease in the bending eigenfrequencies, and labor-intensive production process render these approaches both time-inefficient and uneconomical.

3 Optimal design of internal structure

The aim of this section is to cast the optimal internal structure design problem in the form of a linear semidefinite program (1). The internal structure has to withstand compression molding loads with a maximum deflection bound, while the internal structure is temporarily supported by a steel mandrel passing through the beam interior, Fig. 1(c). Most importantly, the internal structure is supposed to increase the beam fundamental eigenfrequency via reduction of wall instabilities.

In this section, we first describe the finite element model of the composite beam. This finite element model serves then as the basis for establishing the optimization problem formulation, yielding an optimal internal structure design. The section is concluded by discussing post-processing steps necessary to maintain manufacturability of the design.

3.1 Finite element model

The outer composite beam surface is discretized with shell elements which are supplied with the material properties from Table 1. The beam internal structure is modeled by bar (truss) elements, with the isotropic Acrylonitrile Butadiene Styrene (ABS) material properties [52]: elastic modulus $E_{\mathrm{ABS}}=2$ GPa, Poisson ratio $\nu_{\mathrm{ABS}}=0.37$, and density $\rho_{\mathrm{ABS}}=1,040$ kg/m3.

Special care needs to be paid to establishing a rigid connection between the internal structure and the carbon composite. The so-called casing, see Fig. 3, which is a $0.8$ mm thin layer of printed beam walls, further prevents leaking of the epoxy resin into the beam’s interior. Casing is modeled as the bottom layer of the laminate composition, recall Table 1.

The finite element model for the optimization part was developed in Matlab. In this model, the outer laminate was modeled with four-node Mitc4 elements [53]. The composite beam interior was discretized into the ground structure [54], a set of admissible truss111Based on comparative simulations (not shown), modeling internal structure with trusses or beams leads to an insignificant difference in the structural response which enabled us to employ truss topology optimization approaches in Section 3.2. elements, whose cross-sections we search in the optimization part. The ground structure was constructed from $47\times 4\times 4$ modular building blocks shown in Fig. 4, to guarantee manufacturability of the entire internal structure with $3$D printing. Note that the bars placed within the location of the steel mandrel were removed from the ground structure and that the shell element nodes coincided with the ground structure nodes, resulting in a rather coarse discretization of the outer layer.

3.2 Formulation of the optimization problem

3.2.1 Non-convex formulation

Adopting the previously described discretization, our goal is to find the cross-sectional areas $\mathbf{a}$ of $n_{\mathrm{b}}$ bars in the minimum-weight (or volume) ground structure, such that the fundamental eigenfrequency exceeds the user-defined lower threshold $\overline{f}$, taken as $300$ Hz in what follows, while exhibiting limit displacements $\overline{u}$ of the reinforced structure during the compression molding load case. This leads to the following optimization problem

[TABLE]

with

[TABLE]

In this formulation, the vector $\bm{\ell}$ appearing in the objective function (2a) collects the bar lengths in the truss ground structure. The Rayleigh quotient in (2b) involves stiffness, $\mathbf{K}^{\mathrm{C}}_{\mathrm{fv}}$, and mass, $\mathbf{M}^{\mathrm{C}}_{\mathrm{fv}}$, matrices of the outer shell structure for free-vibration analysis together with stiffness, $\mathbf{K}^{\mathrm{IS}}_{\mathrm{fv}}(\mathbf{a})$, and mass, $\mathbf{M}_{\mathrm{fv}}^{\mathrm{IS}}(\mathbf{a})$, matrices of the internal structure. The design-dependent contributions of the internal structure are obtained as

[TABLE]

where $\hat{\mathbf{K}}_{\mathrm{fv},e}^{\mathrm{IS}}$ and $\hat{\mathbf{M}}_{\mathrm{fv},e}^{\mathrm{IS}}$ stand for the stiffness and mass matrix of individual bars in the free-vibration (fv) setting, respectively; $a_{e}$ is the $e$-th component of $\mathbf{a}$ and $\overline{\lambda}$ the limit fundamental free-vibrations eigenvalue.

The constraints (2c) and (2d) address the compression-molding (cm) load case, recall Fig. 1(c). Specifically, (2c) introduces the generalized nodal displacements $\mathbf{u}_{\mathrm{cm}}$ in response to the generalized load vector $\mathbf{f}_{\mathrm{cm}}$ corresponding to the compressive load, and $\mathbf{u}_{\mathrm{cm,disp}}$ denotes the displacement components of $\mathbf{u}_{\mathrm{cm}}$. The stiffness matrix corresponding to this load case consists again of the design-independent, $\mathbf{K}_{\mathrm{cm}}^{\mathrm{C}}$, and design-dependent, $\mathbf{K}_{\mathrm{cm}}^{\mathrm{IS}}(\mathbf{a})$, parts; the latter is obtained as in (4). The symbol $\mathbf{1}$ denotes a column vector of all ones. Notice that the stiffness matrices in (2b) and (2c) differ because of different boundary conditions in the operational, Fig. 1(a), and manufacturing, Fig. 1(c), load cases. The constraint (2d) requires the displacement components of $\mathbf{u}_{\mathrm{cm}}$ to remain smaller than the user-defined limit value $\overline{u}$, considered to be $0.5$ mm in this study. Finally, (2e) requires the cross-sectional areas of the bars to be non-negative and smaller than $\overline{a}=200$ mm2, a value set by the additive manufacturing constraints.

A closer comparison of the optimization problem of Eq. (2) and that of Eq. (1) reveals that the problem of Eq. (2) lacks the structure of a semidefinite program. Namely, the objective function (2a) and the matrices in the constraints depend affinely on the design variables, $\mathbf{a}$. However, the constraints (2b) and (2c) are non-convex as the stiffness and mass matrices may become singular when the zero lower-bound for cross-sectional areas is attained in (2e). Moreover, (2b) might become non-differentiable when an eigenvalue with multiplicity higher than one is encountered. Altogether, this renders the problem (2) extremely difficult to solve in its original form. In the following section, we show how to re-cast the problem of Eq. (2) as a linear semidefinite programming problem.

3.2.2 Convex semidefinite program

Similar eigenvalue constraints such as (2b) have already been studied in detail by Ohsaki et al. [44] and Achtziger and Kočvara [45]. Their results allow us to rewrite (2b) equivalently as a convex LMI

[TABLE]

where the left hand side expression is a linear function of $\mathbf{a}$. This constraint also avoids the non-differtiability of multiple eigenvalues, see, e.g., [45], and effectively eliminates the kinematic variables $\mathbf{u}$ from the problem formulation.

To attain convexity of the final formulation, the compression molding constraints (2c)–(2d) must be enforced only approximately in the form of the LMI [42, 43, 48]:

[TABLE]

in which $c_{\mathrm{cm}}$ denotes a prescribed upper bound on compliance (work done by external forces) of the compression molding load case. As found from parametric studies (not shown), an appropriate value of the bound is provided as

[TABLE]

where $c_{\mathrm{cm,0}}$ stands for the compliance of the non-reinforced structure:

[TABLE]

Here, $\tilde{\mathbf{K}}_{\mathrm{cm}}^{\mathrm{C}}$ and $\tilde{\mathbf{f}}_{\mathrm{cm}}$ are constructed from $\mathbf{K}_{\mathrm{cm}}^{\mathrm{C}}$ and $\mathbf{f}_{\mathrm{cm}}$, respectively, by application of appropriate boundary conditions. For this particular problem, this compliance bound resulted in a maximum deflection of $0.4$ mm.

The final linear semidefinite programming formulation eventually reads as

[TABLE]

This formulation now possesses the structure of the linear semidefinite program introduced in Section 1.2, and thus can be solved efficiently via modern interior-point methods.

For numerical solution, we adopted the state-of-the-art industrial optimizer Mosek [55]. After discretization, the problem in Eq. (9) has $10,216$ admissible bars in total, with the corresponding sizes of the linear matrix inequalities $5,154\times 5,154$ (free-vibration, Eq. (9b)) and $4,608\times 4,608$ (compliance, Eq. (9c)). After tweaking the optimization problem with the steps outlined in the following subsection, the optimization process itself required $13$ GB of memory, and terminated after $5.75$ core hours running on Intel${}^{\text{\textregistered}}$ Xeon${}^{\text{\textregistered}}$ Gold 6130 processors at the MetaCentrum222https://metavo.metacentrum.cz/ virtual organization cluster. The resulting distribution of the optimal internal structure is shown in Fig. 5. Note that the internal structure increased the original weight of the beam, $1,094$ g, by an additional $488$ g ($280$ g of casing and $208$ g of reinforcing bars).

Improving solver performance

To reduce the number of iterations and time per iteration to solve problem (9), we rescale the cross-sectional areas to obtain the optimal values of the order of $1.0$ mm. Second, to improve both the numerical stability and convergence of the algorithm considerably, we rescale Eqs. (9b) and (9c) with the square root of the Frobenius norm estimates of $\mathbf{K}_{\mathrm{fv}}^{\mathrm{C}}$ (Eq. (9b)), and $\mathbf{K}_{\mathrm{cm}}^{\mathrm{C}}$ (Eq. (9c)). Finally, using the static condensation, Appendix A, and Schur complement, Appendix B, decomposition techniques, the sizes of LMIs reduce to $3,426\times 3,426$ (free-vibration, Eq. (9b)) and $2,880\times 2,880$ (compliance, Eq. (9c)). Consequently, memory usage was decreased from $21$ GB to $13$ GB, and the solution process was accelerated by $71\%$ (from $19.5$ to $5.75$ core hours).

3.3 Post-processing

Manufacturing of the optimal design is preceded by three preprocessing steps addressing individual bars, segmentation into modules, and conversion to a solid model. Note that we checked that none of the steps led to the constraint violation and have a rather negligible impact on the objective function, i.e., after all post-processing steps, the internal structure volume increased from $168.3$ cm3 to $175$ cm3.

Bars post-processing

In the initial step of module post-processing, we assign square cross-sections to each bar with the square side length according to the optimal area, $d_{e}=\sqrt{a_{e}}$. Next, we check potential intersection of bars and place a node at each intersection, which subdivides them into two and defines the new element lengths. Third, for each bar, we set the cross-sectional size $d_{e}$ to at least $l_{e}/40$, because more slender bars are difficult to manufacture with the Prusa $3$D printers used in this study. In addition to the optimized bars, the internal structure is extended with short L-shaped beams that ensure mechanical interaction between the internal structure and the steel mandrel, thus defining an empty $20.05\times 20.05$ mm space along the beam longitudinal axis $x$ for its insertion, see Figs. 9 and 8a.

Segmentation

To enable parallel manufacturing with conventional 3D printers, we split the optimized internal structure into $48$ segments of approximately 20 mm in length, see Fig. 6 where ten selected segments are shown, to be assembled later on the steel mandrel. Such segmentation requires re-alignment of bars within each beam cross-section and along the beam longitudinal axis to ensure the correct external beam dimension and a clearly defined interface among adjacent modules, see Fig. 7 for an illustration. Note that segment production does not require any supporting material when printed along the beam longitudinal axis $x$, which would be impossible when printing the internal structure as a single-piece product.

Solid conversion

Axial model conversion is performed independently and in parallel for each node of the ground

structure. We determine first all bars attached to the considered node, elongate them by one half of their cross-sectional side lengths at both of their ends, and cut the more distant half of each of these bars off. These half-bars are then modeled by a mesh-based representation. Geometries of individual nodes then result from the mesh-boolean operations performed with the Cork333https://github.com/gilbo/cork library. Finally, the overall segment geometry consists of the union of all nodal geometries, see Fig. 9 for typical topologies of post-processed segments, and can be readily exported to patch-based Stl file format, for example.

4 Results

4.1 Manufacturing

After the automated export of the optimized internal structure into Stl format, the part was additively manufactured using the Fused Deposition Modeling method with Prusa i3 MK3 printers. Printed segments were inserted on a $20\times 20\times 1,200$ mm steel mandrel of $1.5$ mm wall thickness, with its surface lubricated with Vaseline to simplify the pull-out process, and connected with acetone etching and a thin layer of epoxy glue.

The prototype beam was produced by CompoTech Plus company using the filament winding technology with axial fiber placement. This technology relies on the positioning of the tows of carbon fibers impregnated by the epoxy resin on the casing, placed in specified directions and specified quantity to reach expected dimensions and mechanical properties of the final product. The casing defines the internal shape of the beam and acts as an internal mold. After the fiber placement operation, the product (with still liquid resin) is placed into the press, the outer shape is formed, and the composite is consolidated. In the press, the product hardens at the room temperature. Finally, the prototype, Fig. 10, is postcured at the elevated temperature of $90^{\circ}$C.

The successful manufacturing process was followed by inspection of the prototype using an endoscope camera. Video and photograph sequences, see Fig. 8, revealed that the internal structure successfully withstood the compression molding pressure without any significant visual defects. The minor deviations from the assumed model reside in a small amount of Vaseline residue and slight leakage of the epoxy resin through casing interfaces. Another difference appeared in the increased outer dimensions of the beam, $0.32$ mm on average, caused by an insufficiently closed press cover. The total prototype weight of $1.768$ kg was therefore $186$ g higher than the model predictions due to the additional epoxy resin.

4.2 Verification

Recall that in the optimization we required increasing the fundamental free-vibration eigenfrequency above $300$ Hz. To check this value, we employed an independent model in Ansys. Compared to the model used for optimization, this model employs the dimensions measured in-situ, more refined discretization of the outer shells (element type Shell181), and models the internal structure with beam elements (Beam188) instead of trusses. Besides, the composite shells are supplemented with an additional layer of epoxy resin to account for the increased epoxy content. As a result, the model predicts that the beam fundamental eigenfrequency was increased by $92$% from $128.5$ Hz to $246.7$ Hz, compare Figs. 2 and 11. The effect of wall instabilities was reduced jointly in all the remaining eigenmodes (not shown).

Even though the fundamental eigenfrequency did not exceed the limit value, we find these results satisfactory because of two reasons: First, we attribute this discrepancy mainly to the manufacturing imperfections, which can be attributed to the prototype character of the manufacturing process and can be easily resolved in serial production. Second, the constraint violation is comparable to the difference between numerics and experiments as shown in the next section.

4.3 Validation

Dynamic response was validated with the roving hammer test in the free-free-vibration setting because it eliminates the need to reproduce the simply supported kinematic boundary condition in the experiment. To this goal, the beam was suspended at one of its ends, three piezoelectric acceleration transducers Type 4507B005 Brüel&Kjaer were placed on the beam’s outer surface, two of which were located in the middle of adjacent sides of the beam’s cross-section at one-eighth of the beam’s length, and the third one was placed at the corner of the beam, Fig. 12. Two adjacent sides of the beam surface were marked with a regularly spaced grid of $54$ points, $27$ on each side. These points then served as the excitation points for the impact hammer Type 8206 Brüel&Kjaer equipped with a force transducer.

Measurement was realized using data acquisition front-end hardware Type 3560B Brüel&Kjaer. The frequency response functions (FRFs) were evaluated from the recorded response (acceleration) and excitation (force) using the Fast Fourier Transform for all 54 points. The natural frequencies and mode shapes were evaluated from the FRFs with MEscope software developed by the Vibrant Technology company.

Experimentally determined natural modes and the values of natural eigenfrequencies, Fig. 13 top, were compared with the results of numerical simulations, Fig 13 bottom. Direct comparison in Table 2 reveals sufficient agreement of up to $9\%$ for eigenfrequencies of shear, bending, and torsional eigenmodes. In the case of buckling, we failed to measure the first and second buckling natural modes, and for the higher eigenmodes, the model predictions underestimate the natural frequencies by more than $20\%$. We attribute these deviations to the overall difficulty of measuring the buckling natural modes and to the manufacturing defects discussed in the previous section.

5 Summary and outlook

This contribution introduces and investigates a unique, fully-automatized procedure from an idea to prototyping, with applications to the manufacturing of thin-walled structural composite hollow beams. In particular, the considered prototype product is stiffened with a low weight internal structure designed by an efficient convex linear semidefinite programming formulation. This formulation increased the fundamental free-vibration eigenfrequency above a specified threshold value while avoiding the traditional issue of non-differentiability of multiple eigenvalues [45], and limited structural compliance of a compression-molding load case. The optimization output of the non-uniformly distributed lattice-like internal structure was further automatically post-processed and converted into a solid model ready for support-less additive manufacturing.

Our methodology was verified by designing and producing the simply-supported CFRP beam prototype. Optimization yielded an internal structure of $488$ g which increased the fundamental eigenfrequency by $92\%$ and limited the effect of wall instabilities. Moreover, the deflections within the compression-molding load case were limited to $\pm 0.5$ mm.

After a successful prototype production, the structural response was validated using the roving hammer test, which showed that bending, torsional, and shear eigenmodes exhibited good agreement with model predictions. For the wall buckling eigenmodes, however, the finite element model underestimated the natural frequencies by almost $22\%$. We attribute this to difficulties in measuring these natural modes and to manufacturing defects associated with compression-molding deformations of the casing.

Improving the structural response with a material more than two orders of magnitude more compliant when compared to CFRP suggests concentrating on substituting ABS with high-stiffness continuous carbon fiber in future studies. Another essential future enhancement resides in accelerating the optimization algorithm by exploiting the range-space sparsity [56] associated with the segment-based internal-structure decomposition.

Acknowledgments

We thank Edita Dvořáková for providing us with her implementation of the Mitc4 shell elements [57], and Ondřej Rokoš and Stephanie Krueger for a critical review of the initial versions of this manuscript.

The work of Jan Novák and Robin Poul was supported by the Technology Agency of the Czech Republic, through the project TAČR TH02020420. Marek Tyburec, Jan Zeman, and Matěj Lepš acknowledge the support of the Czech Science Foundation project No. 19-26143X.

Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum provided under the program ”Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042), is greatly appreciated.

Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time due to legal or ethical reasons.

Appendix A. Static condensation of static LMI

Consider the equilibrium equation

[TABLE]

split into two sets of equations

[TABLE]

such that only the principal submatrix $\mathbf{K}_{a}(\mathbf{a})$ depends affinely on $\mathbf{a}$444In the context of this article, the matrix $\mathbf{K}_{a}(\mathbf{a})$ comprises the degrees of freedom of the truss ground structure, $\mathbf{K}_{c}$ contains the remaining (rotational) degrees of freedom, and $\mathbf{K}_{b}$ is the coupling term.. Assuming that the system (10) is solvable uniquely for some $\mathbf{a}$, i.e., it holds that $\exists\mathbf{a}\geq\mathbf{0}:\mathbf{K}(\mathbf{a})\succ\mathbf{0}$, where “$\succ\mathbf{0}$” denotes positive definiteness of the left hand side. Note that $\mathbf{a}=\mathbf{1}$ is sufficient for verification that no rigid movement within the structure can occur. Because $\mathbf{K}_{c}$ is therefore invertible, the degrees of freedom $\mathbf{u}_{b}$ can be expressed from the second row in terms of $\mathbf{u}_{a}$

[TABLE]

and inserted back into the first row,

[TABLE]

Structural compliance (work done by external forces) is expressed as

[TABLE]

After inserting (12) and acknowledging that $\mathbf{K}_{c}^{-1}$ is Hermitian, we obtain

[TABLE]

i.e., compliance of the condensed problem (13) and a constant term. Because the compliance of the condensed problem is positive by definition, the constant term represents a non-negative lower bound on compliances achievable by the internal structure design.

Finally, the LMI

[TABLE]

is equivalent to a smaller LMI

[TABLE]

Further, if $c>\mathbf{f}_{b}^{\mathrm{T}}\left(\mathbf{K}_{c}\right)^{-1}\mathbf{f}_{b}$ is a prescribed constant (i.e., not a variable), then (17) is further reducible, using the Schur complement lemma, e.g., [58, Proposition 16.1], to a yet smaller LMI

[TABLE]

Appendix B. Reducing size of free-vibration LMI

In the case of the free-vibration constraint, we need to directly apply the (generalized) Schur complement lemma. Beginning with reordering of rows and columns, we split the symmetric LMI (9b) such that only the $\mathbf{K}_{a}(\mathbf{a})$ and $\mathbf{M}_{a}(\mathbf{a})$ matrices are functions of $\mathbf{a}$, and the other blocks are constant,

[TABLE]

For the (standard) Schur complement trick we require $\mathbf{K}_{c}-4\pi^{2}\overline{f}^{2}\mathbf{M}_{c}\succ\mathbf{0}$ [58, Proposition 16.1]. Since $\mathbf{K}_{c}\succ\mathbf{0}$ (boundary conditions exclude rigid motions), and $\mathbf{M}_{c}\succ\mathbf{0}$ by definition, we only need to secure that the fundamental eigenfrequency $f_{0}$ of the generalized eigenvalue problem

[TABLE]

with $\lambda=4\pi^{2}f^{2}$, is strictly greater than $\overline{f}$.

Let us therefore first assume that $0\leq\overline{f}<f_{0}$. Then, the inverse of $\mathbf{K}_{c}-4\pi^{2}\overline{f}^{2}\mathbf{M}_{c}$ exists and (19) can be rewritten equivalently using the Schur complement lemma into a smaller-sized LMI

[TABLE]

Second, consider that $f_{0}<\overline{f}$. Because the matrix $\mathbf{K}_{c}-4\pi^{2}(f_{0}+\varepsilon)^{2}\mathbf{M}_{c}$ is indefinite for any $\varepsilon>0$, which renders the original LMI (19) infeasible, the eigenfrequency $f_{0}$ constitutes an upper bound for achievable fundamental eigenfrequencies of the reinforced structure. From the mechanical point of view, the eigenmodes $\mathbf{u}_{b}$ associated with $f_{0}$ excite degrees of freedom not reinforced by the internal structure, and therefore the associated eigenfrequencies can not be increased by any admissible internal structure design (given the specific discretization).

In the case $\overline{f}=f_{0}$, reduction of (19) relies on the generalized Schur complement lemma [58, Theorem 16.1], so that (19) is equivalent to

[TABLE]

where $\left(\bullet\right)^{\dagger}$ denotes the Moore-Penrose pseudo-inverse of $\bullet$, and $\mathbf{I}$ is the identity matrix. The second condition (22b) holds iff the columns of $\mathbf{K}_{b}^{\mathrm{T}}-4\pi^{2}\overline{f}^{2}\mathbf{M}_{b}^{\mathrm{T}}$ are in the image of $\mathbf{K}_{c}-4\pi^{2}\overline{f}^{2}\mathbf{M}_{c}$. Indeed, (22b) can then be rewritten to

[TABLE]

with the columns of $\mathbf{C}$ being the coefficients of linear combinations of the columns of $\mathbf{K}_{c}-4\pi^{2}\overline{f}^{2}\mathbf{M}_{c}$, making the term in the square brackets vanish [58, Lemma 14.1].

Because $\mathrm{Im}(\mathbf{K}_{c}-4\pi^{2}\overline{f}^{2}\mathbf{M}_{c})=\mathrm{Ker}(\mathbf{K}_{c}-4\pi^{2}\overline{f}^{2}\mathbf{M}_{c})^{\perp}$ by [58, Lemma 13.1], it is spanned by

[TABLE]

Clearly, $\overline{f}=f_{0}$ might be achieved iff the columns of the coupling term $\mathbf{K}_{b}^{\mathrm{T}}-4\pi^{2}\overline{f}^{2}\mathbf{M}_{b}^{\mathrm{T}}$ are orthogonal to the eigenmodes occurring in (20) at $f_{0}$. From the mechanical point of view, induction of these eigenmodes would result in a decrease of the associated eigenfrequencies. Note that in practice, equation (22b) can be verified numerically, but it does not guarantee a feasible solution to (22a), because other (higher) eigenfrequencies associated with eigenmodes of (20) may decrease below $f_{0}$ due to the coupling term.