Porous structures have demonstrated their exceptional performances in various engineering fields Lakes_Nature1993 ; Sigmund_JMPS1997 ; Lu_Science1999 ; Liu_JAM2017 ; Krauss_Pqe1999 ; Liu_Pra2015 , and nowadays the manufacturability of porous structures is increasingly enhanced, thanks to the rapid development of additive manufacturing technology. This naturally gives rise to the demand for design tools for porous configurations, and the major challenge stems from the fact that a porous structure is generally associated with more than one length scale, normally a macroscopic scale characterised by its overall size and a microscale identified by the size of its constituting microstructural members. In consideration of such a multiscale feature, analysis of the behaviour of porous structures seems to be constrained down to the microscopic level, which, albeit a number of effective works Alexandersen_CMME2015 ; Liu_JAM2017 ; Liu_CMME2020 , often poses high demands on computation resources.
To resolve the issue of efficiency, homogenisation treatments for porous configurations are also under researchers’ attention. The aim is to decouple the original problem into two subproblems formulated on separate scales. The subproblem on the microscale, which is sometimes called the cell problem, is often localised of a macroscopic material point (Gauss point in numerical implementation) and designated to estimate the equivalent property. The macroscopic subproblem is coarse-grained, and solved with homogenised) properties fed from the solutions to the microscale subproblem. When the microscopic configuration of a porous material becomes spatially periodic, the solution based on the homogenised elasticity tensor sees its rigorousness via perturbation theory, i.e., the scale-separated solution asymptotically converges to that of the original multiscale problem as the length scale ratio tends to zero. This forms the backbone of the asymptotic-analysis-based homogenisation (AABH) methodBensoussan_Book1978 ; cioranescu_Book1999 ; Pavliotis_book2008 . Bendsøe and Kikuchi Bendsoe_CMME1988 pioneerly implemented the AABH method for the topology design of continuum structures. In their frameworks, any macroscopic point in the prescribed domain is assumed being composed of infinitesimal and periodically distributed microstructures. Structural optimisation is manifested through varying the microscopic parameters representing the geometry of the cell. By representing the microstructural configuration with a density field and assuming an explicit relation between the density and the assumed isotropic homogenised material property, the solid isotropic material with penalisation (SIMP) framework was introduced Bendsoe_so1989 ; Zhou_CMME1991 , then intensively developed and popularly used, because of its simplicity and effectiveness. Furthermore, Sigmund Sigmund_IJSS1994 proposed an inverse homogenisation approach to design the unit cell configurations with desirable elastic properties and carried out a series of microstructure design work, such as piezoelectric sensors Sigmund_JMR1998 , optimising thermal expansion coefficient Sigmund_JMPS1997 , and extremal elastic properties material Sigmund_2000JMPS (more details in the monograph Bendsoe_book2003 ). Then concurrent optimisation schemes Rodrigues_Smo2002 ; Coelho_SMO2008 ; LiuLing_CS2008 searching among both the macroscopic structural topologies and the microstructural configurations were proposed, and were then extended in several topics (e.g., Niu_SMO2009 ; Deng_SMO2013 ; Deng_SMO2017 , etc.). It is worth noting, however, that aforementioned works are devised based on the homogenisation results which are actually rigorous only for periodic configurations.
Compared with periodic porous structures, configurations infilled with spatially varying cells have greater design freedom. Early-stage attempts on the design and optimisation of graded porous structures have been made Zhou_JMS2008 ; Radman_JMS2012 ; Radman_CMS2014 ; Wang_CMME2017 ; Cheng_CMME2019 , but challenging issues persist. Firstly, only the interior layouts of the cells are permitted to vary, while the cell shape is kept fixed. Secondly, the microstructural variation must be in alignment with the directions that are orthogonal to one facet of the constituting cells. Thirdly, in most cases, spatial changes in cells are accomplished by (discretely) varying parameters controlling the cell structures. As a result, smooth connections between cells can not be ensured, as the parameter values are essentially piecewise constants.
Recently, further investigation into this subject has attracted wide attention, mainly based on conformal mapping Vogiatzis_CMME2018 ; Jiang_FME2019 ; Allaire_CMA2019 ; Xie_CAD2020 and the projection method Groen_IJNME2018 ; Groen_CMME2019 ; Groen_CMME2020 ; Perle_JCP2020 to generate configurations infilled with smoothly varying cells. The reason that conformal mapping is so attractive, as will be seen here, is that orthotropic porous configurations offer a natural choice for compliance optimisation Groen_IJNME2018 . Under these novel frameworks, the limiting issues outlined in the previous paragraph, such as limited cell orientation and smoothness, are properly addressed. Vogiatzis et al. Vogiatzis_CMME2018 considered using filling the unit cells with specific elastic properties to fill an arbitrary surface, but the issue of optimisation is not discussed. Jiang et al. Jiang_FME2019 discussed optimising porous structures with conformal mapping, but homogenised elasticity tensor is maintained to be isotropic. As a result, the rotational effect of the infilling cells is not considered. The projection-based method, which is also termed as a “de-homogenisation” method by some of its main contributors Groen_CMME2020 , sees its origin in the work by Pantz and Trabelsi Pantz_SIAM_JCO2008 . In this approach, the homogenisation-based topology optimisation is carried out on a coarse grid first. Then the optimised results are projected onto smooth graded and nearly orthotropic porous structures with a desired resolution, i.e. de-homogenisation. Here the projection method effectively offers a post-processing framework. The reason is briefed as follows. In order to fully resolve a porous configuration, de-homogenisation relies on the determination of a smooth mapping function, which indicates the spatial variation of the structure. But the design variables for compliance optimisation are the rotational angles of the matrix cells on pixels, which are related to the spatial derivatives of the mapping function. As a result, one has to integrate optimised results so as to finally describe the structure on the fine scale. However, with the conditions of total differential not met in general, the integrability of the results can not be ensured, and the projection method is then employed to “project” the derivatives onto a total differential form, so as to approximate the desired mapping function.
In parallel with the novel frameworks mentioned above, Zhu et al. Zhu_JMPS2019 proposed an AABH-plus-based optimisation scheme to design graded porous structures. The method consists of three modules: multiscale representation, AABH plus calculation, and compliance optimisation. For representation, a smooth and continuous function is introduced to map a graded porous structure to an artificial periodic structure. Then the moving morphable components/voids (MMC/MMV) approach Guo_JAM2014 ; Zhang_SMO2015 ; Zhang_CM2016 ; Zhang_JAM2017 ; Zhang_CMME2017 is adopted to describe the unit cell configuration. Thus the multiscale topology description function is formed through function composition Liu_JAM2017 ; Zhu_JMPS2019 . The AABH plus part is to (rigorously) derive the corresponding scale-separated form, by means of asymptotic analysis. For optimisation, the parameters controlling the (composite) topology description function (TDF) become the design variables. Upon modification, such treatments of representing graded microstructural configurations with composite TDFs actually provide a platform, on which existing descriptions of porous configurations get formulated in a unified manner Xue_CMME2020 . At its present stage, the AABH-based framework has yet realised its full potential due to the issues detailed as follows. Firstly, unlike periodic structures, cell configurations now differ point by point in the coarse-grained design domain. This means the microstructural cell problems should be computed as many times as the number of macroscopic finite elements, leading to a huge computational burden. This challenging issue is somehow mitigated by adopting a linearisation treatment Zhu_JMPS2019 or a zoning strategy Xue_Smo2020 . Secondly, the use of composite TDFs enhances the describability of graded porous configurations. Then how to manage such strong describability becomes a critical issue. For compliance optimisation, in particular, configurations infilled with orthotropic microstructures sometimes suffice in netting a nearly optimised solution for two-dimensional cases Groen_IJNME2018 . Therefore, the identification of effective porous configurations for optimisation is an issue worth further investigation.
To this end, a compliance optimisation framework among orthotropic porous configurations is proposed in this article. Here the macroscopic mapping function monitoring the gradual change of the matrix cell is restricted to holomorphic functions. Consequently, the porous configurations of interest are generated through gradually rescaling and rotating the matrix cell in space. As from the conformal mapping rule, the logarithm of the scaling factor denoted by here and the rotation angle denoted by here should just correspond to the real and imaginary parts of a holomorphic function. This greatly benefits the realisation of the resulting compliance optimisation, summarised as follows. Firstly, the porous configurations generated through holomorphic mapping functions maintain smooth connectivity naturally, and are orthotropic provided that the matrix cells are orthotropic. Secondly, the equivalent elasticity tensors should be computed from the cell problems only for the cases of the matrix cell, and the actual onsite (homogenised) elasticity tensors can be obtained by rotating the corresponding matrix cells implied by the angle . Hence the homogenised results here are expected (and are also shown) to be asymptotically consistent with that from fine-scale calculations. But the computational efficiency should be enhanced greatly, compared to early-stage AABH-plus-based schemes Zhu_JMPS2019 ; Xue_Smo2020 ; Xue_CMME2020 . It is shown that the optimisation for two-dimensional cases can be brought down roughly below 100 seconds on a desktop computer. Thirdly, the two design controllers and are interrelated through Cauchy-Riemamm equations. The number of the design controller, which is chosen to be here, is then reduced to 1. Fourthly, as must satisfy two-dimensional Laplace’s equation, is fully determined with its prescribed boundary values. Thus the actual design variables are evaluated at a collection of boundary points. Compared to the SIMP-based schemes, the number of design variables is much fewer. Fifthly, the design controller is a harmonic function, whose behaviour is guided by the maximum principle, i.e. the extremal values of must be attained on the domain boundaries. Thus the minimum size of microstructural members can be explicitly controlled through the design variables, the boundary values of . This is especially important when there are requirements on the minimal printable sizes in additive manufacturing. Finally, the rotational angle here is automatically continuous in space. This is in contrast with the projection-based approaches Groen_IJNME2018 ; Allaire_CMA2019 ; Xie_CAD2020 , where or serves as the actual design parameters, and a jump of angle between neighbouring pixels can not be fully avoided without imposing extra conditions.
To the present stage, features that distinguish the present work from others in the frontier, mainly the projection-based methods (e.g., Groen_IJNME2018, ; Groen_CMME2019, ; Groen_CMME2020, ; Allaire_CMA2019, ; Perle_JCP2020, ), are specified in brief, while further elaboration over this point is conducted in the main text. The projection-based methods seek the best cell orientation at each pixel in the design domain, but extra (and sometimes quite intensive) treatments to smoothly connect the neighbouring cells are needed. This is why the projection method Pantz_SIAM_JCO2008 is introduced. The reason, if viewed in a mathematical viewpoint, is a matter of total differential, i.e., knowing the partial derivatives at all pixel points does not mean one can naturally integrate them for a potential function, whose spatial gradients coincide with the prescribed partial derivatives. And the projection method is used to identify an approximate total differential form against the given partial derivatives. In contrast, the real and imaginary parts of a holomorphic function automatically satisfy the conditions for total differential, which are actually expressed by the Cauchy-Riemann equations. This is why no post-processing is needed for the proposed method based on conformal mapping.
The remainder of this article is organised as follows. In Section 2, the AABH-plus-based framework is briefly reviewed. This is followed by the introduction of using conformal mapping technique to derive the related formulation, and a novel optimisation framework towards graded orthotopic porous structures is established in Section 3. Then the issues concerning numerical implementation and corresponding sensitivity of the present approach are discussed in Section 4. In Section 5, numerical examples are presented to illustrate the performance of the proposed method, and conclusions and perspectives compose of Section 6 at the end of the article.
2 Review of AABH-plus-based optimisation framework
In this section, the AABH-plus-based optimisation framework is outlined first, and then key issues for its improvement are discussed.
2.1 AABH-plus-based optimisation framework
The AABH-plus-based framework Zhu_JMPS2019 mainly consists of three parts: the multiscale topology description function, AABH plus approach, and corresponding optimisation formulation.
2.1.1 Multiscale description of graded porous structures
The digital representation of a structure is usually carried out by a topology description function, say , defined by
where and denote the region occupied with solid materials and the whole design domain, respectively. A porous structure, as shown in Fig. 1, is normally associated with two length scales: a microscale length on which microstructures are resolved and a macro length on which the design domain is measured, and . Thus representing a porous structure pixel by pixel is costly. Here we can adopt a multiscale formulation to overcome this issue.
As shown in Fig. 1, a smooth macroscopic mapping function defined in the whole region is introduced, so as to map the graded porous structure of interest to a periodic configuration of period . Then the description of the resulting periodic structure can be confined within a single cell, say , which has been rescaled to be a unit cell, with being the dimensionality number. Suppose for is adopted to identify the material layout in the unit cell. Then the actual graded porous configuration can be described through function composition, i.e.,
Here the structured unit cell is termed as a “matrix cell”. Eq. (2) effectively represents a porous structure by gradually rescaling, rotating and (iso-volumetrically) distorting the matrix cell in space. The onsite operation on the unit cell is captured by the locally Jacobian matrix defined by
Eq. (2) makes use of the self-similar feature of graded microstructures. Hence fine mesh is only needed for resolving , while the mapping function is digitalised on the coarse-grained scale. Note that Eq. (2) can be generalised for cases with multiple matrix cells, and more details can be found in Xue et al. Xue_CMME2020 .
2.1.2 AABH plus formulation
The AABH plus formulation seeks for the scale-separated form of the multiscale formulation defined directly in the porous region, as the parameter
For a graded porous region, the equilibrium equation defined in it should read
for , where is the body force per volume; denotes the displacement field; referring to Eq. (2), is given by
representing the elasticity tensor field; is the elasticity tensor of the base material; a superscript is affiliated with a physical quantity to indicate it resolves the microstructural details in the design domain.
Given the structural periodicity when measured in coordinates, as illustrated in Fig. 1, the multiscale displacement field can be approximated by a series in terms of , which is formed by a set of scale-separated function with x now only measuring mean-field variations. Further asymptotic analysis leads to Zhu_JMPS2019
where is the homogenised displacement field; the third-order tensor is referred to as a first-order corrector satisfying Eq. (8), which is an artificially introduced intermediate quantity for computing the unit cell equivalent property; denotes the homogenised elasticity tensor; is given by
2.1.3 Optimsation framework
On taking the mapping function as the macroscopic design variables, the TDF of the matrix cell as the microscopic design variables, and the equivalent system compliance as the target function, a compliance structures optimisation framework for graded porous configurations can be formulated as follows
where and represent the boundary sections imposed with Dirichlet and Neumann conditions, respectively; denotes the surface force per area acting on the ;
is the outer normal vector of the design domain;is the volume fraction of the region part with solid materials; denotes the admissible upper bound of the volume fraction.
2.2 Challenging issues
Up to this stage, the potential of the AABH-plus-based framework has not been fully realised, mainly caused by the following reasons. Firstly, as pointed out earlier, since the Jacobian matrix depends on the macroscopic coordinate x, the first-order corrector should be solved for point-wisely. This is computationally expensive. To address this issue, actions were taken Zhu_JMPS2019 ; Xue_Smo2020 , but improvement is still in demand. Secondly, the mapping function is up to now restricted to polynomials, and it appears that the optimisation result is quite sensitive to the choice of basis functions. For compliance optimisation, high distortion of the matrix cell is not favoured, and the corresponding treatments for avoiding that are desired.
The present article is aimed to demonstrate that the use of conformal mapping appropriately addresses the challenging issues mentioned above. With the excellent properties of conformal mapping, the proposed framework manages to archive a delicate balance among computational efficiency, gradual varying freedom, and limiting distortion, for compliance optimisation of two-dimensional porous configurations.
3 AABH-plus-based optimisation framework combined with conformal mapping
In this section, we demonstrate how the concept of conformal mapping is integrated in the AABH-plus-based optimisation framework. Besides, the case of rectangular design domain where analytical solutions are available is discussed.
3.1 Properties of graded porous structures generated through conformal mapping
3.1.1 Problem set-up
In the original AABH-plus-based framework Zhu_JMPS2019 , the matrix cell can be rescaled, rotated and distorted in accordance with the onsite macroscopic mapping function . Now we restrict the porous configuration of our attention to be generated through conformal mapping, where the angles about any two crossing microstructural members stay unchanged after mapping. Thus the deformation paradigms of microstructures are limited to rescaling and rotation only. Here we use the symbol to emphasise that the corresponding mapping is conformal. Hence, the Jacobian matrix of can be specified to be
where R denotes a rotation matrix satisfying with I
being the identity matrix. For two-dimensional cases,can be further identified by
where and denote the scaling factor and the rotation angle of the matrix cell in a counterclockwise sense, respectively, as shown in Fig. 2. Here stands for the situation when one edge of the matrix cell is parallel to the -axis. Note that both and are functions of the macroscopic coordinate x, where holds for .
To ensure the existence of the mapping function , and must take exact differential forms, indicating that the components of should be interrelated by the conditions for total differential, i.e.,
which means that and satisfy the Cauchy-Riemann equations.
To this end, if we introduce a complex argument
where represents the imaginary unit, is conjugate to , i.e., they are just the real and the imaginary parts of a holomorphic function. As one conclusion from a holomorphic function, both and defined in are harmonic, that is,
3.1.2 Boundary conditions
Now we consider imposing proper boundary conditions for the two harmonic functions, and . The boundary of a domain, as shown in Fig. 3, may consist of several simple closed oriented curves, denoted by , where is set to be the outer boundary of the region; form the boundaries of its inner ”holes”; n and represent the outer normal vector and the tangent vector of the boundary, respectively. Here the direction of should be specified for a multiply connected region. Here is set counterclockwise on the outer boundary , and clockwise on all other parts, i.e. , as shown in Fig. 3.
Now we demonstrate that one only needs to impose boundary conditions for either or . Hence both and are almost fully determined (up to a constant) in the multiply domain definition. Here we assume the values of are given on . The directional derivative of along the tangential direction on the region boundary, which is also known, is expressed by
which indicates that the boundary information for is also known by means of Neumann boundary condition.
Therefore, a pair of definite problems for and are set to be
where denotes the values of on the boundary. Note that, if is one solution to problem (22b), so is , with being an arbitrary constant. Hence adding a supplementary condition
guarantees the uniqueness of the solution to problem (22b), where is a constant representing the mean value of the rotation angle.
Up to this stage, we manage to show that and within the design domain are fully determined with the boundary values as well as the mean notation angle . Hence we are enabled to just confine the design variables as on a collection of boundary points, as well as the averaged rotational angle .
where and are the Green’s functions corresponding to Dirichlet and Neumann conditions imposed for Laplace’s equation, satisfying
where is the Dirac-delta function.
It is worth noting that both and are independent of boundary conditions and just rely on the shape of . Hence we can get Green’s functions in an offline stage. When the boundary values are changed, one simply inserts them into Eq. (24), and and can be obtained without solving Laplace’s equation repeatedly.
3.1.3 Determination of the mapping function
To fully represent a porous configuration with Eq. (2), one still needs to express the mapping function . If is a simply connected region, it can be seen from Eq. (15) that the line integrals of and defined by Eq. (14) in are path independent. Therefore can be expressed by
where and are the values of and at , respectively, which cast no effects on the configuration of the resulting graded porous structure because of the scale separation.
3.2 Microstructural size control
As is harmonic in , the maximum principle reads
Therefore, the range for in the domain can be fully determined by simply controlling the values of on the boundaries, i.e. . Consequently, the upper and lower bounds of the size of microstructural members constituting a porous configuration can be explicitly imposed through the design variable, . This is of practical value when porous configurations are processed with additive manufacturing techniques, which often yield a minimal printable size. To actually impose the constraint, we first use to denote the (non-dimensioned) minimal size of the microstructures member in the matrix cell, with reference to Eq. (28b), the minimal size of the entire structure can be expressed by
where is recalled to be the characteristic length of the unit cells as indicated in Fig. 1. Therefore, as long as the minimum of the on satisfies
it can be ensured that the minimal size of the microstructure components falls be low the minimal printable size denoted by .
Compared with general cases, the porous structures considered here result only from the rescaling and rotation of the matrix cell. This also helps reduce the computational cost for homogenisation, if compared with the general case briefly reviewed in Sec. 2.1.2.
where and represent the first-order corrector and homogenised elasticity tensor of the matrix cell when and , respectively.
It is proved in Appendix Appendix B that when the base materials are elastically isotropic, the resulting first-order corrector and the homogenised elasticity tensor for arbitrary and can be linked with that of the unit matrix cell by
Eq. (34) is intuitive. Due to the conformality, each constituting cell of the porous structure is effectively obtained by rescaling and rotating the matrix cell. Furthermore, the equivalent properties of unit cells are independent of the scaling factor due to scale separation. Hence the equivalent elasticity tensor of each unit cell is determined by rotating the initial homogenised elasticity tensor counterclockwise with angle .
In the present case, therefore, one simply needs to calculate the homogenised elasticity tensor for a benchmarked case with and . This significantly reduces the computational cost for general cases where the cell problem has to be resolved point by point.
Based on the discussion above, an AABH-plus-based optimisation framework combined with conformal mapping is formulated to be
where represents the volume fraction of the solid parts constituting the porous configuration.
Note that we choose , rather than the rotational angle , as the microscopic design variables. This is because appears in the optimisation formulation by means of or . Whichever the case, ambiguities caused by a jump of in may appear. Note that this is a problem reported in other frameworks Groen_IJNME2018 . Such ambiguities can be circumvented by employing as the design variables.
3.5 Special cases of rectangular design domain
When the design domain takes a rectangular shape, where Laplace’s equation may be analytically solved, the computational burden for conducting optimisation here can be further levied. To see this, we consider a design domain as shown in Fig. 4, where ; and represent the values of correspondingly on . Note that, the compatibility conditions, , and should also hold.
Hence the governing equations of and outlined in the optimisation formulation (35) can be re-written by
With the use of the method of separation of variables, the solutions to problems (36) can be expressed in terms of series given by
where the detailed expressions for ,,, and , for as well as the derivations for Eqs. (37) can be found in Appendix Appendix C. Therefore the procedure of resolving Laplace’s equation can be skipped at each step during optimisation. Instead, one may just insert the design variables which are now and , into the expressions (37).
4 Issues on numerical implementation
In this section, several key numerical issues on implementing the present optimisation method, such as the digital representation of design variables, sensitivity analysis, are considered.
4.1 Digital representation of design variables
The design variables of the proposed optimisation framework are the logarithm of the scaling factor evaluated on the domain boundary, that is, and the TDF of the matrix cell. Now we need to determine appropriate digital representation for them.
4.1.1 Microscopic design variables
Following the preceding studies Zhu_JMPS2019 ; Xue_Smo2020 , the MMC framework Guo_JAM2014 ; Zhang_SMO2015 ; Zhang_CM2016 , is employed to represent the material layout within the matrix cell. In MMC framework, can be represented by the geometric parameters of components, and the microscopic design variables can be collected as
where , , and denote the central coordinates, the half-length, the half-width, and the oriented angle of th moving morphable component, respectively.
4.1.2 Macroscopic design variables
In this paper, interpolation is used to describe. Here we assume that the boundary of the design domain always consists of simple closed curves, each of which can be parameterised by the arc-length variable. Mathematically, this can be given by
for , where and are the arc length and the perimeter of , respectively. Furthermore, defined on the boundary can be expressed by
where for continuity. Now the arc length parameter is then linked with the boundary condition along the boundary tangent, i.e.,
Select points , , , equidistantly on , where and , and define as the values of at point , that is,
for , where denotes the arc-length difference between the neighbouring nodes on .
Through interpolation, is given by
for , where the are interpolation nodes on , and are corresponding basis functions.
With the treatments above, can be represented by the values taken on a finite set of points. Therefore, the macroscopic design variables are collected in
4.2 Sensitivity analysis
For any design variables, being either macroscopic or microscopic, say, , the adjoint method reads
Then we consider the derivatives of the homogenised elasticity tensor with respect to design variables . It is recalled from Eq. (34) that depends only on the rotation angle and the matrix cell . Therefore, the evaluation of
can be classified depending on the type of the design variables of interest, i.e.,