Microstructural Topology Optimization for Periodic Beam-Like Structures Using Homogenization Method

1  Introduction

Beam structures are widely used in many engineering fields, including aviation engineering, civil engineering like bridge and building construction. The microstructures determine the macroscopic performances of beams, and how to find proper microstructures to meet the design requirements becomes a key problem. Compared to size and shape optimizations, topology optimization offers greater design freedoms and serves as an effective tool for achieving the desired structural performance. With nearly 40 years of development, various topology optimization methods have been proposed, such as the homogenization method [1], the Solid Isotropic Material with Penalization (SIMP) method [2], the Bi-directional Evolutionary Structural Optimization (BESO) method [3], and the Level-set method (LSM) [4,5]. Among them, SIMP is widely used for its simple concept.

Sigmund conducted pioneering work on the optimization of extreme material properties using the asymptotic homogenization method (AHM) to evaluate equivalent properties of periodic unit cells [6–8]. This methodology has been extended to design high-performance functional materials, such as maximum bulk modulus [9], improved thermal conductivity [10], and optimized complex modulus (viscoelastic materials) [11]. Zhai et al. [12] promoted a heat conduction-driven method for controlled microstructural boundary expansion to ensure property continuity in functionally graded materials. Isotropic porous materials demonstrate significant engineering potential due to their direction-independent properties. However, naturally occurring isotropic porous materials are rare. Xiang et al. [13] constructed a topological adjustment strategy to derive isotropic porous structures from anisotropic configurations. The AHM-driven topology optimization method can obtain metamaterials with properties difficult to attain naturally, including zero or negative thermal expansion coefficients [14,15], chiral characteristics [16], and negative Poisson’s ratio [17,18]. Long et al. [19] found that incorporating negative Poisson’s ratio material into positive Poisson’s ratio material can enhance the microstructure’s stiffness. Wang et al. [20] advanced the design of mechanical cloaking, demonstrating that strategic assembling of various microstructures enables precise manipulation of elastic responses around shielded objects, thereby rendering them undetectable within their surroundings.

Over the past decade, two-scale structural optimization has remained a prominent research focus, offering enhanced design flexibility compared to single-scale microstructural design paradigms [21]. Within the two-scale optimization framework built by Liu et al. [22], the AHM established the linkage between macroscopic structural performance and micro unit cell configuration. This combining method reduces the computational burden a lot since it transforms the original optimization analysis of a large heterogeneous structure into an analysis of an equivalent homogeneous structure. Building upon this optimization framework, subsequent investigations were extended to address thermoelastic problem [23] and dynamic problem [24], demonstrating its versatility across multi-physics domains. To fully exploit two-scale topology optimization capabilities, researchers conducted systematic investigations into the design of functionally graded microstructures tailored for heterogeneous macro structures [25–27]. The accelerated maturation of additive manufacturing technologies demands high manufacturability of topology optimization solutions. To integrate with computer-aided design, scholars evolved two-scale optimization frameworks incorporating with moving morphable component/bar methods [28,29]. Stiffened shells have extensive applications in aerospace engineering. Zhou et al. [30] applied a two-scale technique to maximize the critical buckling load of grid-stiffened cylindrical shells. By implementing a p-norm constraint, Wu et al. [31] succeed in bone-inspired porous structures. With growing interest in nonlinear optimization [32], Jia et al. [33] established a crashworthiness optimization framework for porous architectures that eliminates the need for sensitivity analysis. By using de-homogenization techniques, Groen et al. [34] realized high-fidelity multiscale structural optimization on desktop computing platforms.

The Giavotto beam theory, which accounts for both warping and shear deformations, has been widely used in beam section optimization [35–37]. In consideration of beam-column connections, Grubits et al. [38] optimized steel I-beams by factoring in geometric and material nonlinear influences. To reduce computational cost in non-uniform beam optimization, Liu et al. [39] developed a mapping strategy that extrapolates sectional properties from a single reference profile to all geometrically similar cross-sections. In two-scale optimization combined with AHM, macroscopic structure analyses often rely on three-dimensional (3D) elastic theory. However, using solid elements for structures with numerous beams greatly increases computational cost without significantly improving accuracy. Moreover, the AHM, constrained to full-directional periodicity [40], is unsuitable for beam structures with limited scale separation along width and thickness dimensions. Yi et al. [41,42] established a AHM for periodic beams and achieved the microstructural design with the specified bending stiffness and extreme torsion stiffness. However, this method ignored the shear effects. Xu and Qian [43] presented a novel method for beam microstructure topology optimization based on the relaxed Saint-Venant solution, addressing the material separation problem (shear effects) with a constraint related to the ratio between two kinds of strain energies. However, considering the influence of shear deformation through the Saint-Venant solution, rather than via the shear stiffness, makes the optimization process less intuitive and complicates the analysis.

We propose a topology optimization framework integrating with fully coupled homogenization theory [44] for periodic microstructure design in beam structures. Unlike conventional AHM-based method [22], our method employs Timoshenko beam theory for macroscopic analysis. This treatment resolves the inaccuracies in equivalent constitutive modeling for 3D homogenized analyses of beam structures while improving computational efficiency through macroscale dimensionality reduction. Compared to previous beam-theory based studies, our method considers deformations not only tension, bending, and torsion, but also for shear. The optimization objectives focus on the extreme tension (compression)–torsion property, the maximization of minimum critical buckling load, and the minimization of structural compliance. Both numerical and experimental results validate the effectiveness of the proposed optimization framework. Notably, the experimental results reveal the effects of beam length on twisting behavior and the energy absorption mechanism of the optimized chiral structure. It is of significance in accuracy and efficiency, and has a great prospect for the optimization design of practical engineering structures with numerous beam components.

The remaining parts of this paper are structured as follows: First, the homogenization or effective stiffness prediction method for periodic beam-like structures is reviewed in Section 2; In Section 3, the topology optimization formulations are presented; In Section 4, several examples are given to show the topology optimization results by this proposed approach; Section 5 gives the conclusions of this work.

2  Homogenization Method for Periodic Timoshenko Beams

In this section, the homogenization method [44] for the effective stiffness prediction of periodic beam-like structures will be briefly reviewed. The generalized constitutive relationship of the equivalent Timoshenko beam as:

[N1M3M2T1Q2Q3]=[D11D12D13D14D15D16D21D22D23D24D25D26D31D32D33D34D35D36D41D42D43D44D45D46D51D52D53D54D55D56D61D62D63D64D65D66][ε1κ3κ2κ1γ12γ13](1)

where [N1,M3,M2,T1,Q2,Q3] and [ε1,κ3,κ2,κ1,γ12,γ13] denote the generalized forces and strains of beams, including the tension, the bending in two directions, the torsion, and the shearing in two directions; Dαβ(α,β=1,2,…,6) is the component of the effective stiffness matrix D.

To achieve the relationship in Eq. (1), the homogenization method considers the following unit cell problems as:

{∂σijε[α]∂xj=0(α=1,2,3,4,5,6)inVσijε[α]nj=0onSnp(2)

where i, j and later appeared m, n are equal to 1, 2, 3 for 3D problems; V is the unit cell domain; Snp and the latter appeared Sp respectively represent the non-periodic faces (or free boundary faces) and the periodic faces, as shown in Fig. 1; a variable with a superscript (ε) indicates that this variable is an actual field with microscopic oscillations; σijε[α] is the actual stress as:

{σijε[α]=Cijmnεmnε[α]εmnε[α]=12(∂umε[α]∂xn+∂unε[α]∂xm)(3)

Figure 1: Sketch of boundary surfaces

In Eq. (3), Cijmn is the fourth-order elastic tensor; εmnε[α] is the actual strain in terms of the actual displacement umε[α], which is separated into two parts as:

umε[α]=um0[α]+um1[α](4)

where um0[α] and um1[α] respectively denote the homogenized displacement and the warping (or the perturbed) displacement. The related strains with um0[α] and um1[α] can be calculated as:

εmn0[α]=12(∂um0[α]∂xn+∂un0[α]∂xm),εmn1[α]=12(∂um1[α]∂xn+∂un1[α]∂xm)(5)

With Eqs. (3)–(5), Eq. (2) can be re-expressed as:

{∂∂xj(Cijmn12(∂um1[α]∂xn+∂un1[α]∂xm))=−∂∂xj(Cijmnεmn0[α])inVσijε[α]nj=0onSnp(6)

where εmn0[α] (or ε3D0[α]) is the homogenized elastic strain corresponding to the α-th generalized unit strain state, and they are listed as:

ε3D0=[ε110ε220ε3302ε2302ε1302ε120],ε3D0[1]=[100000],ε3D0[2]=[x200000],ε3D0[3]=[x300000],ε3D0[4]=[0000−x2x3],ε3D0[5]=[000001],ε3D0[6]=[000010](7)

With the given εmn0[α], um1[α] can be solved by Eq. (6) with the following constraint:

{u11[α]|Sp−=u11[α]|Sp+=0u21[α]|Sp−=u21[α]|Sp+u31[α]|Sp−=u31[α]|Sp+⟨um1[α]⟩=0⟨u2,31[α]−u3,21[α]⟩=0(8)

where ⟨•⟩=1|V|∫V•dV is the average operator defined over the unit cell domain V.

After obtaining um1[α], the actual stress σijε[α] or strain εmnε[α] can be achieved. Further, the effective stiffness Dαβ (α,β=1,2,…,6) can be predicted through

Dαβ=1l∫Vεijε[α]Cijmnεmnε[β]dV=1l∫V(εε[α])TCεε[β]dV(9)

where l denotes the length of a unit cell.

3  Topology Optimization Formulations

3.1 Material Interpolation Scheme

Following the modified SIMP interpolation rule [45], the Young’s modulus of element e, i.e., Ee, can be written as:

Ee=Emin+ρ~ep(E0−Emin)(10)

where e=1,2,…,N, and N denotes the element number in a discretized unit cell; ρ~e∈[0,1] is the physical density of element e, and ρ~e=0 or 1 represents the element is void or solid; Emin and E0 respectively denote the Young’s moduli of void and solid materials; p represents the penalization power with the value of p=3. Emin is set to 10−9 for avoiding the singularity of the stiffness matrix. Further, we denote Ce as the elastic constitutive matrix of element e, and it can be expressed in terms of Ee as:

Ce=EeC0(11)

where C0 is the elastic constitutive matrix with the Young’s modulus E=1.

The physical density ρ~e has the following relationships with the design variable ρe. First, to avoid the checkerboard phenomenon, the density filtering [46,47] is applied to ρe as:

ρ¯e=1∑l∈NeHel∑l∈NeHelρl(12)

where ρ¯e is the filtered density; Ne denotes the set of elements whose center–to–center distance Δ(e,l) to element e is less than the filter radius rmin; Hel is a weighting factor as:

Hel=max(0,rmin−Δ(e,l))(13)

The filtering, however, can lead to the increase of gray elements between solids and voids. To get a clear white–and–black topology configuration, the Heaviside projection [48] is applied to ρ¯e as:

ρ~e=tanh⁡(ξη)+tanh⁡(ξ(ρ¯e−η))tanh⁡(ξη)+tanh⁡(ξ(1−η))(14)

where η is the threshold value and ξ controls the sharpness of the projection. In this work, η is set to be 0.5, and ξ starts at 1 for the first 60 steps, then doubles every 40 subsequent steps until reaching 64 by the maximum iteration of 300.

The above treatments of design variables, including the density filtering and Heaviside projection, can well solve the checkerboard problem and obtain a clear 0–1 topology.

3.2 Optimization Problem and Sensitivity Analysis

Taking the volume fraction as the constraint, the microstructural optimization of a beam-like structure can be formulated as:

{find ρ=(ρ1,ρ2,…,ρN)Tmaxorminf(ρ)s.t.vf=∑e=1Nρ~eve/∑e=1Nve≤vf∗0≤ρe≤1,e=1,…,N(15)

where f(ρ) is the optimization objective. In micromechanical beam optimization, two principal methods exist for microstructural design: one directly employs equivalent microstructural properties as objectives [42], while the second evaluates structural performance of beams containing heterogeneous microstructures as design targets [43]. Our work systematically addresses both paradigms. For microstructural extreme stiffness optimization, we demonstrate the method’s effectiveness through tension (compression)–torsion metamaterial design, as other extreme stiffness optimization problems follow similar principles. For microstructure design targeting macroscopic performance, our framework integrates two critical objectives: stiffness maximization, which remains the predominant objective in structural topology optimization, and buckling resistance enhancement, a key requirement for slender beam applications. ve is the volume of element e; vf and vf∗ respectively denote the present and the upper limit of the volume fraction of a unit cell. Furthermore, through Eqs. (12) and (14), the sensitivities of objective function f(ρ) can be arrived through the chain rule

∂f(ρ)∂ρe=∂f(ρ)∂ρ~e∂ρ~e∂ρ¯e∂ρ¯e∂ρe(16)

Case 1: Optimization of extreme stiffness of microstructure

As in Eq. (1), the diagonal stiffnesses Dαα (α=1,2,…,6) respectively represent the tension stiffness (D11), the two bending stiffnesses (D22andD33), the torsion stiffness (D44), and the two shearing stiffnesses (D55andD66), while Dαβ (α,β=1,2,…,6) for α≠β denote their coupling stiffnesses. When optimization problems targeting the extreme stiffness of microstructures, the objective function f(ρ) is formulated as Dαβ(ρ) or their combinations guided by specific design requirements. For example, in optimizing maximum bending stiffness, the objective function can be expressed as f(ρ)=aD22+bD33. To achieve balanced stiffness in both directions, coefficients a and b should both be set to 0.5. When intentionally improving stiffness in a designated orientation, the corresponding coefficient should be assigned a higher value compared to its orthogonal counterpart.

Using Eq. (9), the sensitivity of Dαβ to physical variable ρ~e is derived as:

∂Dαβ∂ρ~e=1l∫V(εε[α])T∂C∂ρ~eεε[β]dV+1l∫V∂(εε[α])T∂ρ~eCεε[β]dV+1l∫V(εε[α])TC∂εε[β]∂ρedV(17)

Since the last two terms are equal to zero [42], Eq. (17) becomes:

∂Dαβ∂ρ~e=1l∫V(εε[α])T∂C∂ρ~eεε[β]dV(18)

Further, with Eqs. (10) and (11), Eq. (18) transforms to:

∂Dαβ∂ρ~e=pρ~ep−1(E0−Emin)l∫Ve(εε[α])TC0εε[β]dVe(19)

where Ve is the region of element e. For objective function in a combinational form of Dαβ, such as f(ρ)=aD22+bD33, its sensitivity can be dealt with via the linear combination of the sensitivities of D22 and D33. As for the extreme tension (compression)–torsion stiffness, the objective function is selected as f(ρ)=D14.

Case 2: optimization of global performance of beam

(1) Maximization of minimum critical buckling load

The beam, especially for a slender beam, tends to buckle when the axial force reaches or exceeds a certain limit. This instability problem may lead to system damage. Therefore, it is of great importance to maximize the minimum critical buckling load to prevent the structure from entering unstable states. Note that only global buckling for the macroscale structure is considered in the present work.

For linear buckling problems with distinct eigenvalues, the optimization objective in Eq. (15) is rewritten as the minimum or the first-order critical buckling load λ1(ρ), which satisfies the relationship as:

(KE+λ1KG)φ1=0(20)

where φ1=[v(1),w(1),θ3(1),θ2(1),…,v(k),w(k),θ3(k),θ3(k)]T is mode vector of the macroscale structure, of which k denotes the node number in the equivalent beam model; KE and KG respectively denote the assembled elastic and geometric stiffness matrix related to the bending and shearing in x1–x2 and x1–x3 planes, and can be expressed as follows [49].

KE=∫Ω(B¯sub)TDsubB¯subdΩ(21)

KG=∫Ω[∂Nv∂x1∂Nw∂x1]T[∂Nv∂x1∂Nw∂x1]dΩ(22)

where Ω denotes the domain of a beam structure; the submatrix Dsub, extracted from the equivalent stiffness matrix D, exclusively considers the stiffness in two principle planes; B¯sub is the generalized strain–displacement matrix; Nv and Nw are the shape function matrices for interpolations of v and w. As shown in Eq. (22), KG is uncorrelated with the design variables in 1D (beam) problem but becomes correlated in both 2D and 3D problems.

Here, for the convenience of the sensitivity derivations, the λ1 can be expressed in another form as:

λ1=−φ1TKEφ1φ1TKGφ1(23)

Then with Eq. (23), the sensitivity of λ1 to ρ~e is derived as:

∂λ1∂ρ~e=−φ1T∂KE∂ρ~eφ1+2φ1TKE∂φ1∂ρ~e+λ1(φ1T∂KG∂ρ~eφ1+2φ1TKG∂φ1∂ρ~e)φ1TKGφ1=−φ1T∂KE∂ρ~eφ1+λ1φ1T∂KG∂ρ~eφ1+2φ1T(KE+λ1KG)∂φ1∂ρ~eφ1TKGφ1=−φ1T∂KE∂ρ~eφ1φ1TKGφ1=−φ1T(∑e=1∫Ωe(B¯sub)T∂Dsub∂ρ~eB¯subdΩe)φ1φ1TKGφ1(24)

where Ωe is the region of the element e in the domain Ω.

For linear buckling problems involving repeated eigenvalues, sensitivity analysis cannot be directly applied via Eq. (24) due to the eigenvectors are non-unique. To handle this problem, one can solve the following eigensystem as:

{ΦT(∂KE∂ρ~e+λ1∂KG∂ρ~e)Φ+∂λ1∂ρ~eI}Ψ=0(25)

where Φ consists of the original orthogonal eigenvectors that correspond to the repeated minimum eigenvalue λ1, and I is the identity matrix with a dimension corresponding to the multiplicity of the repeated minimum eigenvalues. The eigenvalue of Eq. (25) is the desired sensitivity of repeated eigenvalues and Ψ is the corresponding eigenvector [50,51]. The correctness of the sensitivity derivation, including both distinct and repeated eigenvalues, has been validated using the differential sensitivity method.

(2) Minimization of structural compliance

For an optimization problem concerning the minimum structural compliance, it aims to find a microstructural configuration that maximizes the global stiffness of the periodic beam. The objective function adopts the strain energy C(ρ)=UTKU/2, where KU=F governs the equilibrium condition. The U, K and F are the displacement, assembled stiffness matrix, and load vectors of the beam structure, where

U=[u(1),v(1),w(1),θ1(1),θ3(1),θ2(1),⋯,u(k),v(k),w(k),θ1(k),θ3(k),θ2(k)]T(26)

K=∫Ω(B¯)TDB¯dΩ(27)

and B¯ is the generalized strain–displacement matrix in consideration of all deformation modes.

Then, the sensitivity of strain energy C to ρ~e is derived as:

∂C∂ρ~e=12UT∂K∂ρ~eU+UTK∂U∂ρ~e=−12UT∂K∂ρ~eU=−12UT(∑e=1∫ΩeB¯T∂D∂ρ~eB¯dΩe)U(28)

Using the computed sensitivities, we employ the Method of Moving Asymptotes (MMA) [52] to update design variables. This gradient-based optimization algorithm can accommodate both single- and multi-constraint problems. The iterative process terminates upon reaching the prescribed iteration limit. Fig. 2 presents the key workflow for microstructural topology optimization in periodic beam structures.

Figure 2: Flowchart of microstructure topology optimization in periodic beam structures

4  Optimization Results and Validations

In this section, the Young’s modulus and the Poisson’s ratio of the solid material are set as E0=1 and v=0.3. The dimensions of the unit cell are 24 × 24 × 24 for the optimization problems of tension (compression)–torsion stiffness and structural compliance. For the optimization of critical buckling loads, we consider two unit cells with dimensions of 24 × 24 × 20 (without repeated minimum eigenvalue) and 24 × 24 × 24 (with repeated minimum eigenvalues). All unit cells are discretized by eight-node hexahedron elements with a uniform size of 1 × 1 × 1. For those optimizations involving macro structural analyses, the beam is discretized with elements of size 24. The filter radius rmin is set to 2.5 in this work.

4.1 Optimization of Extreme Tension (Compression)–Torsion Stiffness

To obtain a microstructure with extreme tension (compression)–torsion stiffness, the objective function is selected as f(ρ)=D14, under the volume constraint of vf∗=25%. As illustrated in Fig. 3, the objective value increases rapidly at the beginning of the iterations and then maintains the convergence during each ξ except for the oscillations that occurred with a sudden change of ξ. Compared to the objective iterations, the volume fraction variations are relatively stable and can satisfy the target of 25% finally.

Figure 3: Iteration histories and optimized microstructure for extreme tension (compression)–torsion stiffness

The optimized microstructure, comprising of inclined rods and square loops, demonstrates axial tension (compression)–torsion characteristic, with the resultant stiffness reaching D14=373.6. The calculated stiffness matrix D exhibits symmetry, with other stiffness values as follows: D11=77, D12=D13=D15=D16=D23=D24=D25=D34=D36=D45=D46=D56=0, D22=D33=5834.6, D26=D35=186.1, D44=2898.6 and D55=D66=11.9.

In the following, the compression–torsion and energy absorption performances of the acquired chiral microstructure are investigated through the quasi-static compression test. To analyze the effect of length, three groups of chiral structures with 3, 5, and 8 unit cells along the axial direction are fabricated with nylon P9000. The tests were performed on Instron 2345 with a velocity of 1 mm/min. In Fig. 4a, the torsion degree of the specimen ends is released with the aid of the thrust bearing. Moreover, the circular scale is fixed on the platform, and the torsion angle can be obtained by recording the pointer rotation angle with the video camera [53]. The experimental test would be stopped when the torsion angle of the chiral structure does not change anymore.

Figure 4: Compression-torsion experimental setup and test specimens: (a) experiment setup; (b) test specimens with different lengths

Fig. 5a illustrates deformation progression across four temporal phases for 3 unit cells, while Fig. 5b,c shows the corresponding force–displacement and torsional angle–displacement curves, respectively. In Fig. 5b, with increasing of unit cells from 3 to 8, the structural stiffness and peak force monotonically decrease due to the longer beams being more prone to buckling, of which the forces for 5 and 8 unit cells drop rapidly after reaching the peak force due to global buckling. The whole deformation process with 3 unit cells is complicated and can be roughly divided into two stages, i.e., the twist-domination stage and the compression-domination stage. In Stage 1, corresponding to the four deformation states in Fig. 5a, we know that the twisted deformation is successive, i.e., the bottom cell first begins to twist, followed by the middle cell, and finally the upper cell. After that, the deformation mode comes into Stage 2, in which the twist effect weakens and deformation is gradually dominated by the compression until the structure is to be dense. As shown in Fig. 5c, the maximum torsion angle decreases with the increasing length of chiral structures. And the trends of torsion angle vs. displacement can also be divided into two stages. In Stage 1, the torsion angle changes quickly, while the slope becomes smaller in Stage 2. Moreover, it is noted that the slope transformation point varies with the length of beams. As to the chiral structure with 3 unit cells, the slope becoming smaller in Stage 2 is due to the potential of structural torsion being fully exploited, vs. the same phenomenon that happens in chiral structures with 5 and 8 unit cells is due to the global buckling.

Figure 5: Experimental results of the chiral structures under compressive loading: (a) deformation states at different times; (b) force–displacement curves; (c) degree–displacement curves. Note: color lines in (b) keep the same meaning as in (c)

For comparative assessment of the chiral structure’s performance, a square-tube specimen with equivalent mass and length (matching the chiral structure with 3 unit cells) was also tested. The energy absorption (EA) is calculated via the force–displacement curves before the densification [54]. The comparison results in Fig. 6 show that the peak force of the chiral structure is much lower than that of the square tube with the comparable EA capacity, which indicates that the present chiral structure with tension (compression)–torsion property has a strong buffering capacity and can be used in structural protection.

Figure 6: Comparison results between the chiral and square-tube structures: (a) force–displacement curves; (b) peak force and EA capacity

4.2 Optimization of Global Beam Performances

It is known that macroscopic performance mainly depends on the microstructures. In this part, the global buckling and overall stiffness of cantilever beam problems are respectively discussed for the guidance of the microstructure design.

4.2.1 Maximization of Minimum Critical Buckling Load

This example aims to design microstructures for maximizing the minimum critical buckling loads of beams with different lengths (composed of 5 and 20 unit cells) and different sections (rectangular 24 × 20 and square 24 × 24). The initial design domain of unit cell with dimensions of 24 × 24 × 24 is shown in Fig. 7, in which the outer frames with a thickness of 2 are the undesignable areas during the optimization. Similarly, the unit cell with rectangular sections (24 × 20) also has undesignable frames with a thickness of 2. The volume constraint is vf∗=40%.

Figure 7: Design domain of a unit cell for maximizing the minimum critical buckling load

From Table 1, one can observe that the optimized minimum critical buckling load decreases with the increase of the structural length, illustrating that long beams are more prone to buckling. Moreover, the materials are observed to tend to move from the surface centers to the eight corners with the increase of the structure length. This is caused by the decreasing shear effects in long beams. Besides, the optimized microstructures are also different for unit cells with different sections. For the unit cell with rectangular sections of 24 × 20, the optimized topologies for the upper–lower and front–back surfaces are different. This is because of that the stiffnesses in the two principal planes are different. Thus, to make the beam have similar resistance to buckling in different principal planes, the optimized topologies are inconsistent. For the unit cell with square sections of 24 × 24, the same topologies hold for the four (upper–lower and front–back) surfaces due to the consistency of the two principal planes.

Taking the square beam with 20 unit cells as an example, the predicted minimum bending stiffness by the homogenization method [44] is D22=1.75×104. The analytical minimum critical buckling load of the cantilever beam, calculated using the formula π2EI/(2l)2, is determined as 0.19 with EI=1.75×104, consistent with the value in Table 1. This agreement between numerical and analytical results validates the accuracy of the computational method in predicting buckling behavior.

4.2.2 Minimization of Structural Compliance

In this case, we consider a cantilever beam structure consisting of 20 unit cells under a shear force F=0.2 acting on the right end. The material volume constraint is set as vf∗=40%.

As illustrated in Fig. 8, when only with volume fraction constraint, we can see that the materials are distributed on the upper and lower surfaces with D66=10.1, and there are no materials in the middle. The optimization result is reasonable for this slender beam since its deformations are mainly caused by bending, and bending stresses are maximum at the top and bottom. The material separation, however, is unreasonable in terms of engineering application. Hence to avoid material separations, besides the volume constraint, the shear stiffness constraint is introduced as D66≥D66∗. With the additional constraint D66∗=20 and D66∗=35, the materials firstly connect the middle parts of the lateral faces, then extend to the whole lateral faces, reflecting the shear stiffness contributions within a unit cell. Furthermore, with the shear stiffness increasing from 10.1 to 35, the compliance monotonically increases with the specific values of C=35.3, 36.3, and 39.2. Though the optimized structures with additional shear constraints lost some stiffness performance, the material connections are ensured, further demonstrating the necessity and effectiveness of shear constraints.

Figure 8: Optimized microstructures for minimum structural compliances with or without shear constraint

5  Conclusion

This study proposes a topology optimization framework for periodic beam microstructure design that simultaneously addresses extreme microstructural stiffness and enhanced macrostructural performance while incorporating shear stiffness constraints. Through the proposed method, the optimization problems of tension (compression)–torsion stiffness, critical buckling load, and structural compliance are systematically investigated. Experimental results elucidate the effects of beam length on both twisting behaviors and EA characteristics in tension (compression)–torsion beams. Notably, the optimized chiral structure exhibits promising potential for impact-resistant applications. Through the discussions on the critical buckling load problems, the effects of beam lengths and sections on the optimization results are revealed. The shear stiffness predominantly governs material distributions on lateral faces, effectively eliminate the material separation of slender beams in structural compliance optimization.

Acknowledgement: Not applicable.

Funding Statement: This work was supported by the National Natural Science Foundation of China (grant number 11902015), the Open Fund of Deceleration and Landing Laboratory of the Fifth Academy of Aerospace Science and Technology Group (grant number EDL19092138), and the Ministry of Education Chunhui Plan (HZKY20220014).

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Jiao Jia; methodology, Jiao Jia; software, Jiao Jia and Xin He; validation, Xin He and Zhenchen Liu; data curation, Xin He; writing—original draft preparation, Jiao Jia; writing—review and editing, Jiao Jia; visualization, Xin He; supervision, Jiao Jia; project administration, Jiao Jia and Shiqing Wu; funding acquisition, Jiao Jia, Zhenchen Liu and Shiqing Wu. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author, Jiao Jia, upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest to report regarding the present study.

References

1. Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Comput Meth Appl Mech Eng. 1988;71(2):197–224. doi:10.1016/0045-7825(88)90086-2. [Google Scholar] [CrossRef]

2. Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Arch Appl Mech. 1999;69(9):635–54. doi:10.1007/s004190050248. [Google Scholar] [CrossRef]

3. Querin OM, Young V, Steven GP, Xie YM. Computational efficiency and validation of bi-directional evolutionary structural optimisation. Comput Meth Appl Mech Eng. 2000;189(2):559–73. doi:10.1016/S0045-7825(99)00309-6. [Google Scholar] [CrossRef]

4. Allaire G, Jouve F, Toader AM. A level-set method for shape optimization. Comptes Rendus Mathématique. 2002;334(12):1125–30. doi:10.1016/s1631-073x(02)02412-3. [Google Scholar] [CrossRef]

5. Wang MY, Wang X, Guo D. A level set method for structural topology optimization. Comput Meth Appl Mech Eng. 2003;192(1–2):227–46. doi:10.1016/S0045-7825(02)00559-5. [Google Scholar] [CrossRef]

6. Sigmund O. Materials with prescribed constitutive parameters: an inverse homogenization problem. Int J Solids Struct. 1994;31(17):2313–29. doi:10.1016/0020-7683(94)90154-6. [Google Scholar] [CrossRef]

7. Sigmund O. Tailoring materials with prescribed elastic properties. Mech Mater. 1995;20(4):351–68. doi:10.1016/0167-6636(94)00069-7. [Google Scholar] [CrossRef]

8. Sigmund O, Torquato S. Design of materials with extreme thermal expansion using a three-phase topology optimization method. J Mech Phys Solids. 1997;45(6):1037–67. doi:10.1016/S0022-5096(96)00114-7. [Google Scholar] [CrossRef]

9. Sigmund O. A new class of extremal composites. J Mech Phys Solids. 2000;48(2):397–428. doi:10.1016/S0022-5096(99)00034-4. [Google Scholar] [CrossRef]

10. Radman A, Huang X, Xie YM. Topological design of microstructures of multi-phase materials for maximum stiffness or thermal conductivity. Comput Mater Sci. 2014;91(17):266–73. doi:10.1016/j.commatsci.2014.04.064. [Google Scholar] [CrossRef]

11. Huang X, Zhou S, Sun G, Li G, Xie YM. Topology optimization for microstructures of viscoelastic composite materials. Comput Meth Appl Mech Eng. 2015;283(3):503–16. doi:10.1016/j.cma.2014.10.007. [Google Scholar] [CrossRef]

12. Zhai X, Wang W, Chen F, Wu J. Topology optimization of differentiable microstructures. Comput Meth Appl Mech Eng. 2024;418(2):116530. doi:10.1016/j.cma.2023.116530. [Google Scholar] [CrossRef]

13. Xiang J, Chen J, Zheng Y, Li P, Huang J, Chen Z. Topological design for isotropic metamaterials using anisotropic material microstructures. Eng Anal Bound Elem. 2024;162(2):28–44. doi:10.1016/j.enganabound.2024.01.025. [Google Scholar] [CrossRef]

14. Guo Y, Wang H, Wang W, Chen C, Wang Y. Topology optimization for metamaterials with negative thermal expansion coefficients using energy-based homogenization. Adv Eng Softw. 2024;198(10):103794. doi:10.1016/j.advengsoft.2024.103794. [Google Scholar] [CrossRef]

15. Li Z, Gao W, Wang MY, Luo Z. Design of multi-material isotropic auxetic microlattices with zero thermal expansion. Mater Des. 2022;222:111051. doi:10.1016/j.matdes.2022.111051. [Google Scholar] [CrossRef]

16. Chen W, Huang X. Topological design of 3D chiral metamaterials based on couple-stress homogenization. J Mech Phys Solids. 2019;131:372–86. doi:10.1016/j.jmps.2019.07.014. [Google Scholar] [CrossRef]

17. Han Z, Zhou Y, Xu Z, Wei K, Zhao J, He Z, et al. High stable auxetic metamaterials developed through feature-control topology optimization and additive manufacturing. Thin Walled Struct. 2025;213:113305. doi:10.1016/j.tws.2025.113305. [Google Scholar] [CrossRef]

18. Zhang W, Qin L, Wang J, Xu W. A pixel design method for mechanical metamaterials based on topology optimization. Mech Adv Mater Struct. 2024;31(8):1777–85. doi:10.1080/15376494.2022.2142711. [Google Scholar] [CrossRef]

19. Long K, Du X, Xu S, Xie YM. Maximizing the effective Young’s modulus of a composite material by exploiting the Poisson effect. Compos Struct. 2016;153(2):593–600. doi:10.1016/j.compstruct.2016.06.061. [Google Scholar] [CrossRef]

20. Wang L, Boddapati J, Liu K, Zhu P, Daraio C, Chen W. Mechanical cloak via data-driven aperiodic metamaterial design. Proc Natl Acad Sci U S A. 2022;119(13):e2122185119. doi:10.1073/pnas.2122185119. [Google Scholar] [PubMed] [CrossRef]

21. Wu J, Sigmund O, Groen JP. Topology optimization of multi-scale structures: a review. Struct Multidiscip Optim. 2021;63(3):1455–80. doi:10.1007/s00158-021-02881-8. [Google Scholar] [CrossRef]

22. Liu L, Yan J, Cheng G. Optimum structure with homogeneous optimum truss-like material. Comput Struct. 2008;86(13–14):1417–25. doi:10.1016/j.compstruc.2007.04.030. [Google Scholar] [CrossRef]

23. Yan J, Cheng GD, Liu L. A uniform optimum material based model for concurrent optimization of thermoelastic structures and materials. Int J Simul Multidisci Des Optim. 2008;2(4):259–66. doi:10.1051/ijsmdo/2008035. [Google Scholar] [CrossRef]

24. Liu Q, Chan R, Huang X. Concurrent topology optimization of macrostructures and material microstructures for natural frequency. Mater Des. 2016;106:380–90. doi:10.1016/j.matdes.2016.05.115. [Google Scholar] [CrossRef]

25. Jia J, Da D, Loh CL, Zhao H, Yin S, Xu J. Multiscale topology optimization for non-uniform microstructures with hybrid cellular automata. Struct Multidiscip Optim. 2020;62(2):757–70. doi:10.1007/s00158-020-02533-3. [Google Scholar] [CrossRef]

26. Gao J, Luo Z, Li H, Gao L. Topology optimization for multiscale design of porous composites with multi-doma in microstructures. Comput Meth Appl Mech Eng. 2019;344:451–76. doi:10.1016/j.cma.2018.10.017. [Google Scholar] [CrossRef]

27. Sivapuram R, Dunning PD, Kim HA. Simultaneous material and structural optimization by multiscale topology optimization. Struct Multidiscip Optim. 2016;54(5):1267–81. doi:10.1007/s00158-016-1519-x. [Google Scholar] [CrossRef]

28. Liu C, Du Z, Zhang W, Zhu Y, Guo X. Additive manufacturing-oriented design of graded lattice structures through explicit topology optimization. J Appl Mech. 2017;84(8):081008. doi:10.1115/1.4036941. [Google Scholar] [CrossRef]

29. Nguyen MN, Hoang VN, Lee D. Multiscale topology optimization with stress, buckling and dynamic constraints using adaptive geometric components. Thin Walled Struct. 2023;183(3–4):110405. doi:10.1016/j.tws.2022.110405. [Google Scholar] [CrossRef]

30. Zhou Y, Tian K, Xu S, Wang B. Two-scale buckling topology optimization for grid-stiffened cylindrical shells. Thin Walled Struct. 2020;151:106725. doi:10.1016/j.tws.2020.106725. [Google Scholar] [CrossRef]

31. Wu J, Aage N, Westermann R, Sigmund O. Infill optimization for additive manufacturing-approaching bone-like porous structures. IEEE Trans Vis Comput Graph. 2018;24(2):1127–40. doi:10.1109/TVCG.2017.2655523. [Google Scholar] [PubMed] [CrossRef]

32. Habashneh M, Cucuzza R, Aela P, Movahedi Rad M. Reliability-based topology optimization of imperfect structures considering uncertainty of load position. Structures. 2024;69:107533. doi:10.1016/j.istruc.2024.107533. [Google Scholar] [CrossRef]

33. Jia J, Zhang JC, Long K, Wu SQ, Zhao ZL. Crashworthiness topology optimization of bone-like porous structures. Eng Struct. 2025;336(3–4):120411. doi:10.1016/j.engstruct.2025.120411. [Google Scholar] [CrossRef]

34. Groen JP, Stutz FC, Aage N, Bærentzen JA, Sigmund O. De-homogenization of optimal multi-scale 3D topologies. Comput Meth Appl Mech Eng. 2020;364(6):112979. doi:10.1016/j.cma.2020.112979. [Google Scholar] [CrossRef]

35. Blasques JP, Stolpe M. Multi-material topology optimization of laminated composite beam cross sections. Compos Struct. 2012;94(11):3278–89. doi:10.1016/j.compstruct.2012.05.002. [Google Scholar] [CrossRef]

36. Blasques JP. Multi-material topology optimization of laminated composite beams with eigenfrequency constraints. Compos Struct. 2014;111(7):45–55. doi:10.1016/j.compstruct.2013.12.021. [Google Scholar] [CrossRef]

37. Liu S, An X, Jia H. Topology optimization of beam cross-section considering warping deformation. Struct Multidiscip Optim. 2008;35(5):403–11. doi:10.1007/s00158-007-0138-y. [Google Scholar] [CrossRef]

38. Grubits P, Cucuzza R, Habashneh M, Domaneschi M, Aela P, Rad MM. Structural topology optimization for plastic-limit behavior of I-beams, considering various beam-column connections. Mech Based Des Struct Mach. 2025;53(4):2719–43. doi:10.1080/15397734.2024.2412757. [Google Scholar] [CrossRef]

39. Liu J, Li Q, Liu S, Tong L. Concurrent optimization design of axial shape and cross-sectional topology for beam structures. Struct Multidiscip Optim. 2019;59(6):2287–302. doi:10.1007/s00158-019-02195-w. [Google Scholar] [CrossRef]

40. Terada K, Hirayama N, Yamamoto K, Muramatsu M, Matsubara S, Nishi SN. Numerical plate testing for linear two-scale analyses of composite plates with in-plane periodicity. Int J Numer Meth Eng. 2016;105(2):111–37. doi:10.1002/nme.4970. [Google Scholar] [CrossRef]

41. Yi S, Xu L, Cheng G, Cai Y. FEM formulation of homogenization method for effective properties of periodic heterogeneous beam and size effect of basic cell in thickness direction. Comput Struct. 2015;156:1–11. doi:10.1016/j.compstruc.2015.04.010. [Google Scholar] [CrossRef]

42. Yi S, Cheng G, Xu L. Stiffness design of heterogeneous periodic beam by topology optimization with integration of commercial software. Comput Struct. 2016;172:71–80. doi:10.1016/j.compstruc.2016.05.012. [Google Scholar] [CrossRef]

43. Xu L, Qian Z. Microstructural topology optimization of periodic beam structures based on the relaxed Saint-Venant solution. Struct Multidiscip Optim. 2021;63(4):1813–37. doi:10.1007/s00158-020-02778-y. [Google Scholar] [CrossRef]

44. Huang Z, Xing Y, Gao Y. A new method of stiffness prediction for periodic beam-like structures. Compos Struct. 2021;267(4):113892. doi:10.1016/j.compstruct.2021.113892. [Google Scholar] [CrossRef]

45. Sigmund O. Morphology-based black and white filters for topology optimization. Struct Multidiscip Optim. 2007;33(4):401–24. doi:10.1007/s00158-006-0087-x. [Google Scholar] [CrossRef]

46. Bourdin B. Filters in topology optimization. Int J Numer Meth Eng. 2001;50(9):2143–58. doi:10.1002/nme.116. [Google Scholar] [CrossRef]

47. Bruns TE, Tortorelli DA. Topology optimization of non-linear elastic structures and compliant mechanisms. Comput Meth Appl Mech Eng. 2001;190(26–27):3443–59. doi:10.1016/S0045-7825(00)00278-4. [Google Scholar] [CrossRef]

48. Wang F, Lazarov BS, Sigmund O. On projection methods, convergence and robust formulations in topology optimization. Struct Multidiscip Optim. 2011;43(6):767–84. doi:10.1007/s00158-010-0602-y. [Google Scholar] [CrossRef]

49. Cook RD. Concepts and applications of finite element analysis. New York, NY, USA: John Wiley & Sons, Inc.; 2007. [Google Scholar]

50. Chen TY. Design sensitivity analysis for repeated eigenvalues in structural design. AIAA J. 1993;31(12):2347–50. doi:10.2514/3.11934. [Google Scholar] [CrossRef]

51. Rong JH, Xie YM, Yang XY. An improved method for evolutionary structural optimisation against buckling. Comput Struct. 2001;79(3):253–63. doi:10.1016/S0045-7949(00)00145-0. [Google Scholar] [CrossRef]

52. Svanberg K. The method of moving asymptotes—a new method for structural optimization. Int J Numer Meth Eng. 1987;24(2):359–73. doi:10.1002/nme.1620240207. [Google Scholar] [CrossRef]

53. Ma C, Lei H, Hua J, Bai Y, Liang J, Fang D. Experimental and simulation investigation of the reversible bi-directional twisting response of tetra-chiral cylindrical shells. Compos Struct. 2018;203(4):142–52. doi:10.1016/j.compstruct.2018.07.013. [Google Scholar] [CrossRef]

54. Zhang XC, An CC, Shen ZF, Wu HX, Yang WG, Bai JP. Dynamic crushing responses of bio-inspired re-entrant auxetic honeycombs under in-plane impact loading. Mater Today Commun. 2020;23:100918. doi:10.1016/j.mtcomm.2020.100918. [Google Scholar] [CrossRef]