Full-Scale Isogeometric Topology Optimization of Cellular Structures Based on Kirchhoff–Love Shells

Mingzhe Huang,Mi Xiao,Liang Gao,Mian Zhou,Wei Sha and Jinhao Zhang

State Key Lab of Intelligent Manufacturing Equipment and Technology,Huazhong University of Science and Technology,Wuhan,430074,China

ABSTRACT Cellular thin-shell structures are widely applied in ultralightweight designs due to their high bearing capacity and strength-to-weight ratio.In this paper,a full-scale isogeometric topology optimization(ITO)method based on Kirchhoff–Love shells for designing cellular tshin-shell structures with excellent damage tolerance ability is proposed.This method utilizes high-order continuous nonuniform rational B-splines(NURBS)as basis functions for Kirchhoff–Love shell elements.The geometric and analysis models of thin shells are unified by isogeometric analysis (IGA) to avoid geometric approximation error and improve computational accuracy.The topological configurations of thin-shell structures are described by constructing the effective density field on the control mesh.Local volume constraints are imposed in the proximity of each control point to obtain bone-like cellular structures.To facilitate numerical implementation,the p-norm function is used to aggregate local volume constraints into an equivalent global constraint.Several numerical examples are provided to demonstrate the effectiveness of the proposed method.After simulation and comparative analysis,the results indicate that the cellular thin-shell structures optimized by the proposed method exhibit great load-carrying behavior and high damage robustness.

KEYWORDS Cellular thin-shell structures;isogeometric analysis;full-scale topology optimization;Kirchhoff–Love shells

Nomenclature

Acronyms

FEM finite element method

IGA isogeometric analysis

ITO isogeometric topology optimization

LSM level set method

MMA method of moving asymptotes

NURBS nonuniform rational B-splines

SIMP solid isotropic material with penalization

1 Introduction

Thin-shell structures are widely applied in the aerospace[1],automobile[2]and civil engineering[3]fields due to their excellent load-carrying capacity.The mechanical behavior of a thin-shell structure can be described by shell elements based on surface models because the thickness of the thin-shell structure is much smaller than its size in the other directions.Reissner–Mindlin [4] and Kirchhoff–Love [5] models are the two most widely used mathematical models for describing shells.Reissner–Mindlin theory considering the transverse shear force is primarily used to describe thick shells.It is very suitable for the finite element method(FEM)based on the C0-continuous Lagrangian basis function since the shell elements in this theory need to ensure only C0-continuity.Therefore,the Reissner–Mindlin shell model is widely used in most current industrial software.In contrast,Kirchhoff–Love theory describes thin shells based on the assumption that the shell is thin enough to neglect transverse shear deformation.The rotation of the shell is described by the first derivative of displacement;hence,at least the C1-continuity of the shell elements must be guaranteed in Kirchhoff–Love theory [6].Because the FEM has difficulty constructing C1-continuous elements,Kirchhoff–Love theory is rarely adopted in current mainstream industrial software.

Compared with FEM,isogeometric analysis (IGA) [7] can avoid remeshing via the same basis function to describe the geometric structure and construct the analysis grid,hence eliminating geometric approximation errors by unifying geometric and analysis models of structures and attracting widespread concern [8].NURBS basis functions are the most widely used functions in IGA due to their geometric flexibility and high-order continuity.Because nonuniform rational B-splines(NURBS)is convenient for constructing high-order continuous mesh elements,which aligns with Kirchhoff–Love theory’s requirement for the C1-continuity of shell elements,IGA is suitable for use with the Kirchhoff–Love shell theory.Kiendl et al.[9] developed a Kirchhoff–Love shell element based on IGA,which was applied in thin shells under large rotation by a geometrically nonlinear formulation.Guo et al.[10]proposed a variationally consistent weak coupling method for thin-walled shell patches based on IGA,which ensured the effective transfer of the displacement and bending moment between multiple NURBS patches.This method enabled a blended coupling of shells based on different mathematical models,such as Kirchhoff–Love and solid-like shell models.Li et al.[11] presented a novel method coupling the meshfree method and IGA for geometrically nonlinear analysis of thin-shell structures based on Kirchhoff–Love theory.Thanh et al.[12] further extended the above method to static and free-vibration analyses of cracks in thin-shell structures,and Zhang et al.[13]coupled the isogeometric–meshfree method with the peridynamic method for static and dynamic crack propagation.Hirschler et al.[14] proposed a dual domain decomposition algorithm to address the inconsistencies between multiple NURBS patches for accurately analyzing nonconforming multipatch Kirchhoff–Love shells.Zareh et al.[15]developed isogeometric Kirchhoff–Love shell elements based on Cr smooth rational triangular Bézier splines.These elements overcame limitations on complex geometries requiring multiple NURBS patches,offering a more efficient way to handle complex geometric models.Miao et al.[16] developed an isogeometric Bézier dual mortar method coupling multipatch Kirchhoff–Love shells,which strengthened the continuity of the solution at the boundaries of patches.Peng et al.[17] developed a novel continuum-based fast projection scheme for cloth simulation based on the isogeometric Kirchhoff–Love shell model.Current methods for Kirchhoff–Love shell analysis based on IGA have applications in various emerging research areas,such as metamaterials[18],hyperelastic materials[19],biomechanics[20],lattice shells[21],and optimization of stiffeners for thin shell structures[22].

Topology optimization[23],as an effective structural design method,can enable the performance potential of structures to be fully explored by iteratively finding the optimal distribution of materials in the design domain under given design constraints and objectives.Through topology optimization,unprecedented thin-shell structures that extend beyond the experience and intuition of designers can be generated to further improve the strength-to-weight ratio,hence achieving ultralightweight designs and obtaining broad application prospects [24].Recently,the combination of IGA and topology optimization,i.e.,isogeometric topology optimization(ITO),has attracted great attention.Hassani et al.[25] proposed an ITO method,in which a solid isotropic material with penalization(SIMP) method based on control points was employed.Wang et al.[26] proposed an ITO method based on the level set method (LSM),in which NURBS basis functions were used to construct physical fields and to parameterize level set functions.Gao et al.[27]developed an ITO method based on a smooth density distribution function,in which the Shepard function was used to construct a smooth density distribution function,which was employed to ensure that the structural boundary of the optimized result was clear and smooth.Zhu et al.[28] presented an explicit ITO based on graded truncated hierarchical B-splines.Meanwhile,some studies on the applications of ITO to shell structures have also been conducted.Kang et al.[29] proposed an ITO method for shell structures,in which trimmed surface analysis was used to handle complex topologies.Zhang et al.[30]presented an ITO method for a 3D shell under stress constraints,in which the moving morphable voids method was adopted to explicitly describe the topological configuration.

Additionally,multiscale topology optimization is widely employed for cellular structure design to achieve extraordinary performances,such as excellent dynamic behavior[31],high stiffness-to-weight ratio[32,33]and fast heat dissipation[34].The majority of multiscale topology optimization methods utilize the homogenization theory based on scale separation and periodic assumptions to establish the relationship between macro and micro scales [35].The geometry of the micro structures filling in the cellular structure is generally limited to regular shapes,such as squares and parallelograms.As a result,few investigations have been conducted on the topology optimization design of cellular structures with curved edges or surfaces by homogenization theory[36].Another strategy for designing cellular structures is the full-scale topology optimization method[37].Full-scale topology optimization generates many bone-like holes by imposing local volume constraints within the design domain without using homogenization theory.Therefore,this full-scale topology optimization method can be used to design irregular cellular structures with curved edges or surfaces.

Based on these considerations,a high-order continuous Kirchhoff–Love shell model based on IGA is proposed,and a cellular thin-shell structure with high damage robustness is designed via fullscale ITO.NURBS is used to describe the geometry of the shell and establish its analysis model to improve the analysis accuracy.A density field is built on the control mesh to describe the topological configuration of the shell.Local volume constraints,which are aggregated into an equivalent global constraint by thep-norm function,are introduced into ITO to obtain bone-like thin-shell structures.The simulation and comparison of several numerical examples show that the cellular thin-shell structures obtained by the proposed method have great bearing capacity and damage robustness.

The remainder of the paper is organized as follows:Section 2 describes the full-scale ITO model based on the Kirchhoff–Love shell in detail.Section 3 provides the numerical implementation and framework.Section 4 provides several numerical examples to verify the effectiveness of the proposed method and the excellent performances of the structures optimized by the proposed method.Section 5 summarizes and discusses the research work of this paper.

2 Full-Scale Isogeometric Topology Optimization Based on the Kirchhoff–Love Shell

2.1 Local Volume Constraint

To generate a bone-like cellular structure on a thin shell,the local volume constraint method is introduced[37].The average of the effective densities in the designated neighborhoodncof each control point is taken as the local volume fractionLcof the current control pointc.Lcis calculated as follows:

whereρirepresents the density of the control point in the designated neighborhoodnc,wherencrepresents the set of adjacent control points of the current control point.The distance between adjacent control points and the current control point must be less than the specified influence radiusr,as follows:

wherePcis the coordinate of the current control point andPiis the coordinate of the adjacent control point.The local volume constraint is represented as follows:

wheregcrepresents the local volume constraint of thecth control point,Lmaxdenotes the specified upper limit of the local volume fraction,andNcdenotes the total number of control points.Imposing a local volume constraint on each control point results in many constraints,and implementing these constraints is difficult.These constraints can be represented by an aggregated constraint max(Lc)≤Lmax.Since the aggregated constraint is nondifferentiable,thep-norm function is used to approximate the max function as follows:

Whenptends to infinity,||L||pis completely equivalent to max(Lc).Ifpis not infinite,the aggregated constraint can be written as:

The size ofpdetermines how strict the local volume constraint is.

2.2 Full-Scale ITO Model

In the full-scale design framework,an ITO model for the cellular thin-shell structure is established with the minimum compliance as the objective function and the global and local volume fractions as constraints,as follows:

whereρiis the effective density at theNccontrol point,namely,the design variable.Jdenotes the objective function,that is,the minimum compliance.ueis the element displacement,andΩis the design domain of the thin shell.Kis the global stiffness matrix,uis the displacement vector,andfis the force load vector.The equilibrium equation is constructed by the Kirchhoff–Love model based on IGA.G(ρi)represents the global volume constraint,Vmaxrepresents the specified upper limit of the global volume fraction,andνrepresents the element volume.g(ρi)represents the aggregated local volume constraint.ρminrepresents the lower limit of the design variable.represents the interpolated element stiffness matrix[23]and is calculated as follows:

wherex(ρi)is the effective density field constructed by the effective densityρion the control point and the corresponding NURBS basis functionRi.x(ρi)is expressed as follows:

whereτdenotes the penalty factor,andnrepresents the number of control points,which is the number of effective densities.φi(ρi)represents the Heaviside step function[38],written as:

whereσrepresents the projection parameter to control the sharpness of the Heaviside step function.The Heaviside step function is adopted to promote the generation of the 0–1 solution.The larger theσvalue is,the more the effective density field obtained by the Heaviside step function tends toward being a 0–1 distribution.However,directly applying a largeσcan result in highly nonlinear equations.To ensure the stability of the iteration process,σis initialized asσ=1,and value is doubled after every 50 iterations untilσ=16.This is a commonly used technique for promoting convergence in topology optimization [38].The effective density of any point on the shell can be obtained through the effective density fieldx(ρi),which reflects whether each point on the shell corresponds to solid material or void.The topological configuration of the shell is optimized by optimizing the reasonable distribution of effective density.

2.3 NURBS Basis Functions

To construct NURBS basis functions,first,the B-spline basis function with a monotonically increasing open knot vectoris expressed by the recursive formula [8],as follows:

wherenandprepresent the number and degree of B-spline basis functions,respectively.The knot vector represents the coordinates of the element nodes in the parametric domain,and the knot spans represent the elements.Then,introducing the weightωiyields the NURBS basis functions as follows:

whereωiis the weight.Furthermore,a NURBS curve can be written as:

wherePiis the coordinate of the control point corresponding to the NURBS basis function.Similarly,a NURBS surface can be described by the tensor product of the NURBS basis functions in directionsξandη,written as:

wherenandmare the number of control points in two parametric directions,andpandqare the corresponding polynomial degrees.Specifically,is written as follows:

whereNi,p(ξ)andNj,q(η)represent B-spline basis functions in two parametric directions.

2.4 Kirchhoff–Love Shell Model

In this paper,the IGA-based Kirchhoff–Love model is used to analyze the mechanical behavior of the thin shell.First,the three-dimensional shell is described with the curvilinear coordinate(ξ1,ξ2,ξ3),whereξ1andξ2represent the in-plane directions,andξ3represents the shell thickness direction.As shown in Fig.1,the position vectorsxandXof any point in the shell in the current and reference configurations,respectively,are represented as follows[39]:

whereaαandAαrepresent the covariant basis vectors of the mid-surface in the current and reference configurations,respectively.Here,the Greek indexα=1,2.The corresponding inverse basis vectorsaαandAαare defined as follows:

Figure 1:Geometric description of the Kirchhoff–Love shell.The mid-surface is highlighted in blue,red points represent control points,and red dashed lines represent the control mesh.The knot vectors are{0,0,0,0,1/3,2/3,1,1,1,1}in two parametric directions

The displacement vectoruof the mid-surface between the current and reference configurations is calculated as follows:

The covariant basis vectors at any material point in the current and reference configurations are defined as follows:

Similarly,the corresponding inverse basis vectorsqαandQαare defined as follows:

The covariant metric coefficients and the inverse metric coefficients at any material point in the current and reference configurations are expressed as follows:

Based on Kirchhoff–Love’s hypothesis,the Green-Lagrange strain tensor[40]can be written as:

Bringing Eqs.(22)into(24)yields

42.Promised:Cinderella promises to leave the ball before midnight but ultimately breaks this promise with her late departure. The breaking of the promise gives Cinderella a slight hint of imperfection and humanity. It also shows how much she is enamored with the prince.Return to place in story.

whereεαβandκαβrepresent the membrane strain and bending strain,respectively,as follows:

Next,the strain energy density per unit area of the Kirchhoff–Love shell can be expressed as:

wheretis the shell thickness.Cis the material tensor,written as:

whereEis Young’s modulus andνis Poisson’s ratio.

2.5 Governing Equations and Discretization in IGA

The potential energy of the Kirchhoff–Love shell can be calculated by the following formula:

whereδεandδκrepresent the first variation of membrane and bending strains.nandmrepresent the membrane force resultant tensor and bending moment tensor,respectively,and are expressed as:

Eq.(30) is solved iteratively.Assumingui+1=ui+Δu,the variational Eq.(30) is linearized as follows:Find Δu∈such that

whereKmandKbare the membrane stiffness and bending stiffness,respectively,and are expressed as:

Eq.(32)can be converted into a discrete equation by numerical integration:

whereKis the positive definite stiffness matrix,uis the displacement vector,andfis the force load vector.Eq.(34)is also the equilibrium equation in Eq.(6).Furthermore,in IGA,uis expressed by the linear combination of the displacementui,jat the control point and the corresponding NURBS basis function,as follows:

Similarly,the coordinatesr(ξ1,ξ2)of any point on the mid-surface of the shell can be expressed as:

The stiffness matrixKis composed of the element stiffness matrixKe.Keis composed of the element membrane stiffness matrixand element bending stiffness matrix,as follows:

whereCeis the element material tensor,Bmis the membrane strain-displacement matrix,andBbis the bending strain-displacement matrix.3 × 3 Gauss-quadrature points are selected in each element to calculate the element stiffness matrix.ωijis the weight of the Gauss-quadrature point in each element.J1is the Jacobian matrix mapped from the parametric domain to the physical domain,andJ2is the Jacobian matrix mapped from the parent domain to the parametric domain.Here,Gauss-Legendre quadrature rules are employed for numerical integration in IGA.

3 Numerical Implementation and Framework

3.1 Sensitivity Analysis

Since the above full-scale ITO model involves multiple constraint functions,method of moving asymptotes (MMA) algorithm is adopted to update the design variables [41].The MMA method requires calculating the sensitivity of the objective function and constraints.The sensitivity of the objective functionJ(ρi)with respect to the design variableρiis calculated as follows:

Using the chain rule,the sensitivity of local volume constraintg(ρi)with respect to design variableρiis calculated as follows:

Similarly,the sensitivity of the global volume constraintG(ρi)with respect to the design variableρiis as follows:

3.2 Optimization Design Framework

The optimization design framework of the cellular thin-shell structure is shown in Fig.2.The specific steps of the framework are as follows:

Step 1:Construct the design domain of the thin shell based on the NURBS basis function and set boundary conditions;

Step 2:Construct the effective density field and local volume fraction field;

Step 3:Compute the displacement field of the thin-shell structure by IGA;

Step 4:Calculate the objective function and constraints and conduct sensitivity analysis;

Step 5:Update effective densities by the MMA method;

Step 6:If the maximum change in effective density between two iterations is less than 1%or the number of iteration steps reaches 800,stop the iteration;otherwise,repeat steps 3-6;

Step 7:Output the optimized result to obtain the final cellular thin-shell structure.

4 Numerical Examples

In this section,several numerical examples are provided to demonstrate the effectiveness of the proposed method.Unless otherwise specified,stainless steel material with Young’s modulusE=200 GPa and Poisson’s ratiov=0.3 is used in all examples.The upper limits of the global and local volume constraints are set to 0.5 and 0.6,respectively.Both the penalty factor and the filter radius are set to 3.

4.1 Scordelis–Lo Roof

First,a benchmark example is used to verify the accuracy of the proposed method.Fig.3 displays the design domain of the Scordelis–Lo Roof,a very popular and extremely severe example.Fig.4 shows the IGA mesh of the Scordelis–Lo Roof.The thickness of the shell is 0.25 mm.Two curved edges are fixed,and a uniform vertical load F=90 N is imposed on the entire shell.To align with the standard description of this benchmark problem,the material properties are set to Young’s modulusE=432 MPa and Poisson’s ratio v=0 here.The detailed data of the NURBS on the Scordelis–Lo Roof are displayed in Table 1.

Table 1:Detailed data of NURBS on the Scordelis–Lo Roof

Figure 3:Design domain of the Scordelis–Lo Roof

Figure 4:IGA mesh of the Scordelis–Lo Roof

The vertical displacement field obtained by the proposed method is shown in Fig.5.The maximum vertical displacement of the middle side calculated by the proposed method is 0.3077,which is close to the reference solution of 0.3086 in[42]and the solution of 0.308 in[15].Hence,the accuracy of the proposed method is demonstrated through this benchmark example of the Scordelis–Lo Roof.

Figure 5:Vertical displacement field of the Scordelis–Lo Roof obtained by the proposed method

4.2 Blunt Cone Shell

A design domain of a blunt cone shell is shown in Fig.6.The four corners of the blunt cone shell are fixed,and a vertical downward force load F=100 N is applied at its center point.The thickness of the thin shell is 1 mm.The design domain of the blunt cone shell is divided into 200 × 200 IGA elements for topology optimization,as shown in Fig.7.The parameters of the corresponding IGA mesh are provided in Table 2.In the optimization model,the influence radius of the local volume fraction is set to 10.

Table 2:Detailed data of NURBS on the blunt cone shell

To verify the analysis results of IGA,the analysis result of the blunt cone shell obtained by the commercial FEM software COMSOL is compared with that obtained by IGA.By importing the STL file of the blunt cone shell into COMSOL,the same geometric modeling can be obtained.Fig.8 depicts the displacement fields obtained by IGA and COMSOL.Here,the number of elements of IGA is set to 50×50 to demonstrate the analysis accuracy of IGA,while the number of elements in COMSOL is 7942.Comparison reveals that the contour lines of the two displacement fields by IGA and COMSOL are almost identical.Meanwhile,the maximum displacements of the analysis results of the two methods are very similar.Therefore,this comparison illustrates that IGA requires fewer elements to achieve the same level of accuracy as the FEM.

Figure 6:Design domain of the blunt cone shell

Figure 7:IGA mesh of the blunt cone shell,where the mesh resolution is set to one-tenth of the actual resolution for better visualization

Figure 8:Displacement fields of the blunt cone shell:(a)by IGA;(b)by COMSOL

The iteration curves of the blunt cone shell are shown in Fig.9.The iteration curves fluctuate significantly in the first 200 steps due toσin the Heaviside function doubling every 50 steps,as described in Section 2.2.After 200 iterations,the iteration curves are stable,and the compliance reaches the minimum while meeting the volume constraint.Fig.10 displays the final optimized result,namely,the optimized effective density distribution.Many small voids are produced in the optimized results due to local volume constraints.The material is primarily distributed along the main force transmission path from the load application point at the center of the blunt cone shell to the four fixed corners.

Figure 9:Iteration curves of the blunt cone shell

Figure 10:Optimized result of the blunt cone shell:(a)optimized cellular thin-shell structure;(b)top view of the cellular thin-shell structure

The influences of different parameters on the optimized result are also discussed through this numerical example.Fig.11 shows four optimized results under different upper limitsLmaxof the local volume constraint and influence radiusr.The upper limit of the local volume constraint controls the local porosity,while the influence radius affects the size of the holes.The larger the influence radius is,the larger the size of the holes.Therefore,the size of the holes in the optimized result can be controlled by adjusting the upper limit of the local volume constraint and the influence radius.

Figure 11:(Continued)

4.3 Free-Form Surface

Fig.12 shows the design domain of a free-form surface,where F1=50 N and F2=10 N.The design domain is discretized into 200×200 IGA elements,as displayed in Fig.13.The thickness of the free-form surface is 1 mm.The parameters of the corresponding NURBS patch are provided in Table 3.The influence radius of the local volume fraction is set to 16.

Table 3:Detailed data of NURBS on the free-form surface

Figure 12:Design domain of the free-form surface

Figure 13:IGA mesh of the free-form surface,where the mesh resolution is set to one-tenth of the actual resolution for better visualization

Fig.14 depicts the iteration curves of the free-form surface.The iteration curves become stable after 200 iterations and the convergence condition is satisfied at the 653rd iteration step.Fig.15 shows the optimized results obtained by the proposed method.A free-form surface with many holes is generated,and the material is mainly distributed along the main force transmission path.Fig.16 shows the 3D-printed sample obtained by additive manufacturing for the optimized result of the free-form surface.

Figure 14:Iteration curves of the free-form surface

Figure 15:Optimized results of the free-form surface: (a) optimized cellular thin-shell structure;(b)top view of the cellular thin-shell structure

Figure 16:The 3D-printed sample of the optimized result of the free-form surface

In Fig.17,the optimized structure obtained by the proposed method and the thin-shell structure with uniformly distributed holes are imported into COMSOL for simulation comparison.The same material properties and loading and boundary conditions are set in COMSOL.The optimized structure obtained by the proposed method has a 31.3% lower compliance and a 48.0% lower maximum displacement than the thin-shell structure with uniformly distributed holes.These results demonstrate that the proposed method can effectively improve the mechanical properties of thin-shell structures.

Figure 17:The free-form surfaces and corresponding displacement fields for(a)and(b):the proposed method,and(c)and(d):uniformly distributed holes

4.4 Cylindrical Shell

Fig.18 displays a design domain of the thin cylindrical shell,where the tangential force F=100 N.The thickness of the cylindrical shell is 1 mm.The design domain is meshed with 300×120 IGA elements,as shown in Fig.19.The parameters of the corresponding IGA mesh are listed in Table 4.In the initial setting,the influence radius of the local volume fraction is set to 12.

Table 4:Detailed data of NURBS on the cylindrical shell

Figure 18:Design domain of the cylindrical shell

Figure 19:IGA mesh of the cylindrical shell,where the mesh resolution is set to one-tenth of the actual resolution for better visualization

The iteration curves of the cylindrical shell are shown in Fig.20.Due to changes inσ,the iteration curves still exhibit slight oscillations in the first 200 iterations.Then,the iteration curves become stable and reach the minimum compliance at the 800th iteration step.Fig.21 displays the optimized result of the cylindrical shell by the proposed method.The material is mainly distributed along the main force transmission path from the load application points at the top of the cylindrical shell to the fixed bottom.Since the applied forces are all tangential forces with the same magnitude and clockwise direction,the optimized structure is highly symmetrical.Meanwhile,many bone-like cellular structures are present in the cylindrical shell due to the influence of local volume constraints.Fig.22 shows the 3D-printed sample obtained by additive manufacturing for the optimized result of the cylindrical shell.

Figure 20:Iteration curves of the cylindrical shell

A major advantage of the cellular thin-shell structure optimized by the proposed method is its excellent damage tolerance ability.Even if the structure is damaged to a certain extent,the performance of the structure will not be greatly impacted.Figs.23a and 23b show the optimized results by the proposed method and classical topology optimization,respectively.Unlike the proposed method,the classical topology optimization method does not set local volume constraints.The material in the blue area is removed to simulate local failure.Comparing the compliances of the two structures before and after local failure reveals that before the local failure,the compliance of the optimized structure obtained by the proposed method is 20.3%greater than that by the classical topology optimization;however,after the local failure,the compliance of the optimized structure obtained by the proposed method is 39.4%less than that by the classical topology optimization.In addition,the compliance of the optimized structure obtained by the proposed method is 21.2%greater after local failure,while the compliance of that by the classical topology optimization is 140.6%greater.These results demonstrate that local failure has less influence on the mechanical properties of the structure optimized by the proposed method than on the structure optimized by classical topology optimization.

Figure 21:Optimized results of the cylindrical shell: (a) optimized cellular thin-shell structure;(b)elevation view of the cellular thin-shell structure

Figure 22:The 3D-printed sample of the optimized result of the cylindrical shell

The two structures after local failure are then imported into COMSOL for simulation.The same material properties and loading and boundary conditions are applied in COMSOL.The displacement fields of the two structures are shown in Figs.23c and 23d.The maximum displacement of the optimized structure obtained by the proposed method is 80.1% smaller than the maximum displacement of that by the classical topology optimization.Therefore,the cellular thin-shell structure optimized by the proposed method has excellent robustness to local failure.

5 Conclusions

This paper presents a full-scale ITO method to realize the design of cellular thin-shell structures with excellent bearing capacity and damage tolerance.High-order continuous NURBS is adopted to build the geometric and Kirchhoff–Love models of thin shells.The displacement field of a thinshell structure is computed by IGA with high computational accuracy.A full-scale ITO framework is established by introducing local volume constraints,which can generate bone-like structures on thin shells.The size of the holes in the cellular thin-shell structures can be adjusted by changing the local volume constraints and their influence on radii.Local volume constraints are aggregated by thep-norm function to facilitate an effective optimization process.The effectiveness of the proposed method is verified by several numerical examples.These numerical examples show that the proposed full-scale ITO method can steadily generate cellular structures in various irregular design domains of thin shells.The comparative discussion demonstrates that the optimized structure obtained by the proposed method has high bearing capacity and great robustness to local damage,enabling its use in practical engineering applications with high requirements for safety factors.

However,importantly,the optimized results obtained by the proposed method may have buckling or stress constraint problems due to the presence of overly thin substructures.To address these issues,future work will consider incorporating buckling and stress constraints into the optimization process to optimize the mechanical properties of cellular thin-shell structures more comprehensively.

Acknowledgement:Thanks to the help of three anonymous reviewers and journal editors,the logical organization and content quality of this paper have been improved.

Funding Statement:This research was supported by the National Key R&D Program of China(Grant Number 2020YFB1708300),China National Postdoctoral Program for Innovative Talents (Grant Number BX20220124),China Postdoctoral Science Foundation(Grant Number 2022M710055),the New Cornerstone Science Foundation through the XPLORER PRIZE,the Knowledge Innovation Program of Wuhan-Shuguang,the Young Top-Notch Talent Cultivation Program of Hubei Province and the Taihu Lake Innovation Fund for Future Technology(Grant Number HUST:2023-B-7).

Author Contributions:The authors confirm contribution to the paper as follows:study conception and design:Mingzhe Huang,Mi Xiao,Liang Gao;data collection:Mingzhe Huang,Mian Zhou;analysis and interpretation of results:Wei Sha,Jinhao Zhang;draft manuscript preparation:Mingzhe Huang,Mi Xiao.All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials:The data and models generated and optimized in this paper are available from the corresponding author on reasonable request.

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