Trung-Son Nguyen

Tung Le Duc

Thanh-Son Tran

Jean-Michel Guichon

Olivier Chadebec

Gérard Meunier

Adaptive Multipoint Model Order Reduction Scheme for Large-Scale Inductive PEEC Circuits

Keywords:

The model order reduction techniques based on the multipoint projection Krylov methods have become the methods of choice to generate large macromodels of multiport RLC circuits. A well-known difficulty with such methods lies in the need for clever point selection to attain model compactness and accuracy.

In this paper, we present an automatic methodology for optimizing sample point selection. This method is general, suitable for circuits obtained by the partial element equivalent circuit method coupled with an adaptive multilevel fast multipole method. Our algorithm has been validated on an industrial example to demonstrate the accuracy and robustness of the approach.

I. INTRODUCTION

P

explosion in memory by the number of terms in the full PEEC matrix. In order to save the memory and the computation time, we treat the interactions between far-away basis functions using the fast multipole method (FMM) [START_REF] Ardon | EMC modeling of an industrial variable speed drive with an adapted PEEC method[END_REF]. The corresponding matrix elements do not need to be explicitly stored, result in a significant reduction in required memory. By using FMM, the resulted model is available only in the frequency domain and the link with a conventional circuit representation is lost. To be included in time-domain simulation of circuit-solvers, this model needs to be replaced by a reduced-order model, which is much smaller but would produce a very good approximation of input-output behavior of the original device.

This problem can be done in two steps. The first step is the reduction order of frequency model returned by PEEC. We then obtained a simpler and lighter reduced-order model. In the second step, this reduced model is synthesized into a circuit model in netlist representation (RLC elements), which can be integrated in a circuit-solver [START_REF] Yang | RLC equivalent circuit synthesis method for structure-preserved reduced-order model of interconnect in VLSI[END_REF]. The pure RLC model is then simple enough to ensure the convergence of the resolution and accurate enough to return reliable results. The second step is presented in our recent work [START_REF] Nguyen | Circuit realization method for reduced order inductive PEEC modeling circuits[END_REF]. In this paper, we investigate the first step (reduction of frequency model): the model order reduction (MOR).

MOR methods for linear model reduction have greatly evolved and can be broadly characterized into two types: those based on subspace generation and projection (or momentmatching method), and those based on balancing techniques [START_REF] Antoulas | A survey of model reduction methods for large-scale systems[END_REF]. The balancing techniques can generate models with provable general error-bound [START_REF] Wolf | Gramian-based error bound in model reduction by Krylov subspace methods[END_REF], on the contrary, global error is unknown by projection methods [START_REF] Kamon | Generating nearly optimally compact models from Krylov-subspace based reduced-order models[END_REF]. However, the balancing techniques rely upon dense matrix computations, which seem unable to deal with large scale PEEC circuit compressed by FMM. In this case, due to a matrix compression algorithm, only projection method is suitable.

The projection method has led to construct the projection space from a rational or multipoint Krylov subspace [START_REF] Grimme | A rational Lanczos algorithm for model reduction II: Interpolation point selection[END_REF]. The multipoint approximants are usually more expensive to construct than single point interpolation as PRIMA [START_REF] Odabasioglu | PRIMA: Passive reducedorder interconnect macromodeling algorithm[END_REF], but they tend to have better quality models for given effort [START_REF] Wolf | Gramian-based error bound in model reduction by Krylov subspace methods[END_REF]. Multipoint-based approaches [START_REF] Feng | A fully adaptive scheme for model order reduction based on moment-matching[END_REF]- [START_REF] Rasekh | Multiorder arnoldi approach for model order reduction of PEEC models with retardation[END_REF] have recently gathered renewed attention due to their robustness, reliability.

Although very appealing, a major drawback of such methods in practical implementation is: how many interpolation points to choose, and how to place them. For large-scale systems, the computation cost of Krylov subspace at each expansion point ARTIAL ELEMENT EQUIVALENT CIRCUIT (PEEC) approach [START_REF] Ruehli | Equivalent circuit models for three dimensional multiconductor systems[END_REF] is known to be very suitable for the modeling of combined circuit and electromagnetic field problems. For low-frequency applications, capacitive effects can be neglected and the method can be restricted to its inductive formulation [START_REF] Ardon | EMC modeling of an industrial variable speed drive with an adapted PEEC method[END_REF], [START_REF] Le-Duc | Coupling between PEEC and magnetic moment method[END_REF]. With such approach, arbitrarily shaped 3-D conducting devices can be represented by equivalent circuits combining resistors, partial/mutual inductances, and current and/or voltage sources. The circuit-based model of electromagnetic devices allows its interface with SPICE-like circuit solvers. This approach is known to be very efficient for the electromagnetic modeling of power electronic modules for instance [START_REF] Ardon | EMC modeling of an industrial variable speed drive with an adapted PEEC method[END_REF], [START_REF] Le-Duc | Coupling between PEEC and magnetic moment method[END_REF] and finely model can be coupled with the external electrical circuit.

The increment of number of conductor elements (directly proportional with the complexity of the devices) causes an T.-S. Nguyen is with the Vietnam Fire and Rescue Police Department, Hanoi 10000, Vietnam (e-m ail: sonngtrung@ gm ail.com ).

T. Le Duc is with the Department of Electric Power Systems, Hanoi University of Science and Technology School of Electrical E ngineering, H anoi 10000, Vietnam (e-m ail: tung.leduc1@hust.edu.vn).

T.-S. Tran is with the Faculty of Electrical Engineering, Electric Power University, H anoi 100000, Vietnam (e-m ail: sontt@epu.edu.vn).

J.-M. Guichon, O. Chadebec, and G . Meunier are with the G2ELA B, G renoble 38000, France (e-mail: Jean-Michel.Guichon@g2elab.grenoble-inp.fr; Olivier. Chadebec@g2elab.grenoble-inp.fr; gerard.meunier@g2elab.grenoble-inp.fr). is expensive, which makes the choice of these points is very important.

Until now, this area of research has barely been touched upon and the problems above still remain open. An early method called Complex Frequency Hopping (CFH) is proposed in [START_REF] Chiprout | Analysis of interconnect networks using complex frequency hopping (CFH)[END_REF] and [START_REF] Achar | Simulation of high-speed interconnects[END_REF] to illustrate a principle of choosing multiple expansion points of the transfer function. By using a binary search algorithm, the expansion points are chosen with respect to the common poles contained in both circles of the neighboring expansion points. However, the choice of expansion points is not optimal, which can lead to a large number of expansion points. The ARMS technique in [START_REF] Villena | ARMS-automatic residue-minimization based sampling for multi-point modeling techniques[END_REF] proposed a fully automated procedure to choose expansion point, which is based on maximizing the subspace spawned by the sampling set in order to obtain a good approximation of the dominant subspace in the orthonormalization step (construction of Krylov subspace). However, this approach is not specifically directed to point placement but to choose the best future point in a set of initial point.

Obviously, strategies include random selection of points or uniform grinding of the domain. Both of these can be shown as uninformed algorithms about either the structure or the behavior of the system, and as such, they may not be optimal. Resampling plans have been proposed in [START_REF] Silveira | Resampling plans for sample point selection in multipoint model-order reduction[END_REF] to provide guidance for future sample point placement in linear MOR, based on variance measurements of the reduced transfer function. However, it requires the evaluation of several reduced models, which may be expensive, and, in addition, seem awkward to extend to multidimensional cases.

Intuitively, an efficient point selection algorithm should place new points, where the error or model uncertainty is large without evaluating the detailed model, which is expensive. There is no general answer to this question but we know some points could have large error: f min , f max , poles of the models (peak point in frequency response) [START_REF] Grimme | Krylov projection methods for model reduction[END_REF]. Unfortunately, these poles are not known a priori. The work in [START_REF] Gugercin | Rational Krylov methods for optimal H2 model reduction[END_REF] for SISO systems and a MIMO version [START_REF] Gugercin | H2 model reduction for large-scale linear dynamical systems[END_REF] shows how to compute these points by an iterative procedure and the algorithm is fast convergence. It aims at automatically finding the optimal ||H|| 2 -norm model for a fixed qth order, based on a refinement of the interpolation points. However, the algorithm needs to solve the original system at r-points at each refinement iteration, which makes the procedure expensive. The authors in [START_REF] Panzer | Model order reduction by Krylov subspace methods with global error bounds and automatic choice of parameters[END_REF], [START_REF] Wolf | H2 pseudo-optimal model order reduction[END_REF] proposed a cumulative reduction scheme (CURE framework), which offers more flexibility version of IRKA [START_REF] Chiprout | Analysis of interconnect networks using complex frequency hopping (CFH)[END_REF] without loss of precision. In this method, the authors reuse the previous constructed Krylov space to reduce computation cost. A recent work in [START_REF] Feng | A fully adaptive scheme for model order reduction based on moment-matching[END_REF] presented a fully adaptive scheme for Krylov-based MOR. Although, the proposed approach provides an interesting way of adaptive choosing moment (q), the computation time is rather high because at each expansion point we need to compute the original model and the reduced model. The author in [START_REF] Rewieński | Greedy multipoint modelorder reduction technique for fast computation of scattering parameters of electromagnetic systems[END_REF] also proposed a greedy multipoint MOR technique applied for second-order systems and a "goal-oriented" error estimator different than the system's impedance. Recently, there are many research related to MOR methods for finite element method [START_REF] Clenet | Reduction of a finite-element parametric model using adaptive POD methods-application to uncertainty quantification[END_REF], [START_REF] Sato | Generation of equivalent circuit from finiteelement model using model order reduction[END_REF] but it is difficult to apply for our research because of the fully dense matrix in PEEC method (FEM matrix is sparse).

The multipoint MOR method for delayed PEEC circuits including the capacitive effects [START_REF] Ferranti | Multipoint full-wave model order reduction for delayed PEEC models with large delays[END_REF] and parameterized MOR [START_REF] Ferranti | Interpolation-based parameterized model order reduction of delayed systems[END_REF] have been developed recently. In our research, the PEEC formulation is restricted on the inductive effects but it can be extended to general case as in [START_REF] Ferranti | Multipoint full-wave model order reduction for delayed PEEC models with large delays[END_REF] because the fully dense and sparse matrices are all treated in our study. An interesting approach for adaptive choice of expansion point is presented in [START_REF] Ferranti | Interpolation-based parameterized model order reduction of delayed systems[END_REF] but the validation process of the method needs the computation of original model. In our approach, the exact computation of original model is impossible because of FMM algorithm. We also avoid the high computation cost of full model and focused on memory and time saving for large-scale application.

In this paper, we will propose a novel methodology for automated selection of optimal expansion points for first-order systems. The procedure not only shares some ideas in [START_REF] Feng | A fully adaptive scheme for model order reduction based on moment-matching[END_REF] and [START_REF] Gugercin | Rational Krylov methods for optimal H2 model reduction[END_REF] but also exploits some properties of Krylov subspace and thick-restart techniques [START_REF] Freund | Recent advances in structure-preserving model order reduction[END_REF] to reduce the computation cost as much as possible. The originality of this study is its formulation coupled with the FMM method and a novel adaptive algorithm for selecting expansion points automatically in regard to minimize computation cost.

The paper is organized as follows: in Section II, the back ground state space equation of PEEC method is presented. In Section III, we discuss the implementation of our algorithm in presence of FMM and an adaptive point selection algorithm is introduced in Section IV. For the Section V, we show a computational experiment to illustrate the advantage of proposed approach and give several conclusions in Section VI.

II. PEEC RLMC CIRCUIT EQUATIONS

Considering a RLMC circuit (resistor-R, inductance-L, mutual inductance-M, capacitor-C) that is only excited by current sources. The incidence matrix of nodes/branches A, the branch currents i B , and nodal voltages v N can be written as follows:

A = [A R A C A L A I ]; i B = i R i C i L i I T v N = v R v C v L v I T (1)
where the subscripts R, L, M, C, and I are associated with the branches containing resistors, the partial and mutual inductances (calculated by PEEC integral method), the external capacitors, and the external current sources. The state equations of electrical circuits provide:

i I = -I t (t) ; i R = R EXT -1 • v R ;i C = C EXT d dt v C ;v L = L EXT d dt i L ( 2
)
where I t (t) is the vector of current sources, R EXT and C EXT are diagonal matrices, and L EXT is an inductance matrix from external circuit. For full matrix L PEEC (inductive coupling PEEC) and resistance diagonal matrix R PEEC , we have an equation restricted to the PEEC elements: In order to compute the impedance matrix, the final state space formulations are:

v lPEEC = R PEEC • i lPEEC + L PEEC d dt i lPEEC . (3)
E dx (t) dt = -G • x (t) + B • u (t) y (t) = B T • x (t) (4)
where

E = ⎡ ⎢ ⎢ ⎢ ⎣ A C • C EXT • A T C 0 0 L PEEC 0 0 L EXT ⎤ ⎥ ⎥ ⎥ ⎦ G = ⎡ ⎢ ⎢ ⎢ ⎣ A R • R EXT -1 • A T R A L -A T L R PEEC 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ (5)
B = A I 0 x = v N i L u = -i I . ( 6
)
By applying the Laplace transform to (4), we can calculate the matrix Z(s) or transfer function:

Z(s) = B T • (G + s • E) -1 • B. (7)

III. MULTIPOINT EXPANSION MOR COUPLED WITH FMM

With an expansion point s 0 , we can rewrite [START_REF] Wolf | Gramian-based error bound in model reduction by Krylov subspace methods[END_REF] as follows:

Z(s) = B T • (G + s 0 • E + (s -s 0 ) • E) -1 • B = B T • (I -(s -s 0 ) • D) -1 • F (8)
and

D = -(G + s 0 • E) -1 • E;F = (G + s 0 • E) -1 • B. (9)
In the conventional PRIMA method, only one expansion point (s = 0) is used to construct the Krylov subspace created by the method block Arnoldi [START_REF] Odabasioglu | PRIMA: Passive reducedorder interconnect macromodeling algorithm[END_REF]. This approach is equivalent to a Taylor series expansion of the last (8) at s 0 = 0 Fig. 1 K

m (D, F) ≡ span F, D • F, D 2 • F, ..., D q -1 • F = span {s 0 } (10)
with q order of Taylor series expansion (size of Krylov subspace).

Algorithm 1: Block Arnoldi for multipoint. 1.

K total = Ø 2.
For i = 1, 2 , . . . , r, Do 3.

Orthonormalisation F i with K total by modified Gram Schmidt process. If K total = Ø, orthonormalisation F i 4.

Construction K i by block Arnoldi method 5.

K total = ortho(K total , K i ) 6. EndDo

A. Multipoint Expansion Method

In this approach, several expansion points are selected for construction of Krylov subspace, the result is then more accurate Fig. 2.

Combination of several subspaces for each expansion point s j gives us the new Krylov subspace U m :

U m = span F, D • F, D 2 • F, ..., D q -1 • F points=0 , ...F j , D j • F j , D j 2 • F j , ..., D j q -1 • F j points =0
= span {s 0 , s 1 , ..., s j } [START_REF] Feng | A fully adaptive scheme for model order reduction based on moment-matching[END_REF] with

D j = -(G + s j • E) -1 • E;F j = (G + s j • E) -1 • B.
The Krylov subspace, therefore, consists of a first portion corresponding to the first point of expansion and a second portion corresponding to the second point and so on. Let us notice that the vectors constituting this subspace are orthogonal with each other so the blocks representing different points of expansion are orthonormal.

The following algorithm allows us to construct the Krylov subspace in [START_REF] Feng | A fully adaptive scheme for model order reduction based on moment-matching[END_REF] for the expansion point s0, . . . ,sr [START_REF] Ferranti | Multipoint full-wave model order reduction for delayed PEEC models with large delays[END_REF]. Recall that the Krylov subspace created by an expansion point s j by starting matrices D j and F j is K j . The set of vectors of the Krylov subspace is K total .

Using Algorithm 1, we have the Krylov subspace K total containing several expansion points. The linearly dependent vectors were deleted by deflation process [START_REF] Freund | SPRIM: Structure-preserving reduced-order interconnect macromodeling[END_REF]. Finally, we got the Krylov subspace U m . Reduced order matrices obtained with the congruent transformation are then: The transfer function Z R (s) is then reduced:

Ê = U T m • E • U m ; Ĝ = U T m • G • U m ; B = U T m • B. (12)
Z R (s) = BT • Ĝ + s • Ê -1 • B. ( 13
)
All the important properties of PRIMA method like the preservation of passivity and the preservation of moments are conserved. Consequently, the reduced order matrix Ĝ is then invertible [START_REF] Freund | SPRIM: Structure-preserving reduced-order interconnect macromodeling[END_REF].

For full PEEC matrix, the most consuming time of construction Krylov subspace is to solve the linear system of ( 9). Without a compression matrix algorithm, the computational cost (matrix-vector product) and the memory cost (for fully dense matrix storage) are at least O(N 2 ) (N is the size of the matrix) for a good implementation of iterative solver [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. An AMLFMM compression algorithm is presented in the next section in order to speed-up the computation time and to reduce memory consumption.

B. Fast Multipole Algorithm in Block Arnoldi Process

The main computational cost of this method is computing matrix-vector product (D x in the following) in the block Krylov method to obtain U m .

With the expansion point s = 0, D x = v is the product made by solving a linear system:

D • x = G -1 • E • x = v ⇔ G • v = E • x. (14)
The matrix G is a sparse matrix (which incorporates mainly the resistance); we use a sparse direct LU solver with AMD algorithm in order to conserve the sparse properties of LU decomposition and to enhance the performance [START_REF] Davis | Direct Methods for Sparse Linear Systems[END_REF]. In the case of circuits obtained with PEEC, the matrix E (especially L P E E C ) is compressed by an algorithm of type FMM [START_REF] Jazzar | Modélisation électromagnétique des grands systèmes: l'apport des méthodes intégrales dans l'étude du foudroiement des avions[END_REF] so the matrix-vector product E x is greatly accelerated from O(N 2 ) to O(NlogN) (N is the size of matrix D). An example of G matrix is shown in -Fig. 3(a). The white elements in low matrix block in Fig. 3(b) are compressed thanks to FMM. Without FMM, this part is fully dense matrix.

With the expansion point s ࣔ 0, we have a new linear system:

D • x = (G + s 0 • E) -1 • E • x = v ⇔ (G + s 0 • E) • v = E • x. (15)
Note that (G + s 0 E) is compressed by FMM, so that the system of ( 15) is no more sparse. To solve [START_REF] Rasekh | Multiorder arnoldi approach for model order reduction of PEEC models with retardation[END_REF], we need a solver compatible with FMM (only require product matrix-vector). In this case, iterative solver like GMRES is chosen by its efficiency. However, even if this technique is applied, we may not be able to obtain the solution because of the nonconvergence of GMRES algorithm. Therefore, we need to apply a preconditioning technique in order to ensure the convergence and/or to reduce the number of iterations. To do that, we modify the original system to the new one as [START_REF] Chiprout | Analysis of interconnect networks using complex frequency hopping (CFH)[END_REF], where P is called the preconditioner matrix:

P • (G + s 0 • E) • v = P • E • x. ( 16
)
The new system has a better condition number result in better chance to converge. In theory, the matrix P should be as similar as possible to (G + s 0 E) -1 to significantly improve the condition number and speed-up the convergence and the best matrix P is the matrix (G + s 0 E) -1 . However, to build a preconditioner matrix compatible with the FMM is somewhat different than an explicit matrix solver, because many of the matrix elements are no longer explicitly available. Therefore, the preconditioner has to be built from the near-matrix elements only. We decompose the matrix E into a near-field component and a far-field one thanks to FMM:

E = E near + E far . ( 17
)
The preconditioned matrix P is built from the total matrix G, s 0 , and E near as [START_REF] Grimme | Krylov projection methods for model reduction[END_REF], where only largest elements in each row of E near are contributed to the preconditioned matrix (selected thanks to a coefficient (coeff) range from 0 to 1):

P = [G + s 0 • E near (coeff)] -1 . ( 18
)
A value in the ith row of the matrix E near is retained if this value is bigger than the maximum value of this line (typically diagonal) multiplied by coeff. If coeff equals 0, the entire matrix E near is retained, if the coeff equals 1, only the diagonal of matrix E near is retained. The usage of coeff parameter allows us control the sparsity of matrix P. The sparsity of this matrix is in between matrix G and E near (see Fig. 3) so the direct sparse solver could be used effectively.

In the next section, we will present an algorithm for automated selection of expansion point for a fixed moment matching (qorder) at each point.

IV. ADAPTIVE MULTIPOINT SELECTION SCHEME

Up to now, IRKA algorithm [START_REF] Gugercin | Rational Krylov methods for optimal H2 model reduction[END_REF], [START_REF] Gugercin | H2 model reduction for large-scale linear dynamical systems[END_REF] or CURE framework [START_REF] Panzer | Model order reduction by Krylov subspace methods with global error bounds and automatic choice of parameters[END_REF], [START_REF] Wolf | H2 pseudo-optimal model order reduction[END_REF] give the best approximation of optimal ||H|| 2 -norm reduced model with a fixed q-order moment matching, among all reduced systems sharing the same set of poles [START_REF] Gugercin | H2 model reduction for large-scale linear dynamical systems[END_REF].

In this section, we present the basic idea of this algorithm and then we propose our new algorithm for adaptive selection of interpolation point based on IRKA algorithm.

In theory, for a dynamic system of classic form as [START_REF] Gugercin | Rational Krylov methods for optimal H2 model reduction[END_REF], the optimal reduced model might be expected to approximate the eigenvalues of the original system as [START_REF] Gugercin | H2 model reduction for large-scale linear dynamical systems[END_REF]:

dx (t) dt = -Ax (t) + Bu (t) y (t) = B T x (t) (19)
error = eig (A) -eig  < tol. (20)
As in IRKA algorithm for r expansion points, the author propose using r initial shifts randomly distributed within a region containing the mirror image of the numerical range of original system A as starting expansion point and use r eigenvalue of reduced system as a next expansion point in an iterative processes. The algorithm will stop if it reaches the condition as [START_REF] Panzer | Model order reduction by Krylov subspace methods with global error bounds and automatic choice of parameters[END_REF] for r values:

error = eig Âi+1 -eig Âi < tol. (21)
It was shown that this algorithm always converges after a few number of step and the resulted model is optimal in ||H|| 2 -norm sense [START_REF] Gugercin | Rational Krylov methods for optimal H2 model reduction[END_REF], [START_REF] Flagg | Convergence of the iterative rational Krylov algorithm[END_REF].

However, the value of r (number of expansion point) is not known. The r initial point impact on the iteration process and the choice of these values is important. It could be expensive if we use the eigenvalues of original system. For each refinement step, we need to compute Krylov subspace again and the subspace of previous steps is deleted. At each sample vector of Krylov subspace, we need to solve a large linear system of ( 16), which is expensive. From the computation point of view, the cost of refinement step is very high.

We should avoid computing the new vectors of Krylov subspace unless it is strictly necessary. For treating the problem of selecting initial points and limiting the number of these points, we propose a novel approach.

To choose a set that covers our space of interest, and because we do not know system frequency response, we select only two initial points f min and f max as expansion points. As the properties of Arnoldi process for q-order, it is easy to approximate the largest (dominant) eigenvalues of original system within a very small number of iteration. We expect that with these two points we have some information to compute (approximately) all the eigenvalues in the range [f min , f max ] thanks to the properties of multipoint expansion.

To do that, we generate the orthonormalized Krylov subspace at two points U m = spanf min , f max , then we compute the eigenvalues of reduced system:

λ R = eig Ĝ-1 • Ê . ( 22
)
In [START_REF] Kundur | Power System Stability and Control[END_REF], for classic dynamic system in [START_REF] Gugercin | Rational Krylov methods for optimal H2 model reduction[END_REF], the poles or peak points will be the imaginary value of eigenvalues:

poles = imag [eig (A)] . ( 23
)
For our dynamic system of the form as ( 4), the poles are computed as follows:

poles = imag eig E -1 • G = imag 1 eig (G -1 • E) . ( 24
)
Because the eigenvalue computation of original system is expensive, in the same fashion as IRKA algorithm, we compute the poles of reduced system:

poles R = imag ⎛ ⎝ 1 eig Ĝ-1 • Ê ⎞ ⎠ . ( 25
)
In [START_REF] Mark | The Arnoldi eigenvalue iteration with exact shifts can fail[END_REF], the authors indicated that we could find the largest eigenvalues of large scale matrix by Arnoldi method. The more difficulty is to find the smallest eigenvalues, so if we have a reduced system satisfy the condition in [START_REF] Davis | Direct Methods for Sparse Linear Systems[END_REF] for maximal and minimal eigenvalues, our reduced system is nearly optimal as ||H|| 2 -norm sense error = max min (poles)min (poles R ) , max (poles)max (poles R ) < tol.

(26) The eigenvalues of reduced system might be expected to approximate the eigenvalues of original system if the target error tolerance satisfies. In our practical applications, the maximum relative errors are found at the peak points, where the dominant resonance peaks appear, so the choice of expansion points near these peak points will reduce overall error [START_REF] Grimme | A rational Lanczos algorithm for model reduction II: Interpolation point selection[END_REF]. The only difficulty is how to compute these points precisely.

The IRKA algorithm is suitable for large-scale systems and it has always converged after a small number of steps while using eigenvalues as expansion points [START_REF] Gugercin | Rational Krylov methods for optimal H2 model reduction[END_REF]. Based on this idea, if we use the points min(poles R ) and max(poles R ) as expansions points, we will obtain a good approximation of original eigenvalues after few iteration steps. Since poles closest to the expansion point converge fastest, this idea is called thick-restart process in [START_REF] Freund | Recent advances in structure-preserving model order reduction[END_REF].

We propose an original adaptive three points selection scheme to construct reduced order model with automated selection of expansion points.

The main idea of the proposed algorithm is to keep building the projection basis by adding more expansion points in the range of interest until the reduced-order model becomes accurate enough in the neighborhood of the maximum, minimum, and the middle range eigenvalues (three values). The next expansion point is added based on the value of eigenvalue of reduced-order model. The accuracy of the model is assessed by computing error estimators at three points.

In this algorithm at each step, we add the previous subspace in the new subspace because the information from last step remains relevant from cycle to cycle. The reuse of Krylov subspace reduces the computation cost at each refinement step. The middle expansion point (e mid ) is used not only for testing the precision of our subspace at this point (which may be the maximum error), but also helps to enlarge subspace in order to make the computation of minimal and maximal eigenvalue is more accurate.

If three tolerances of low, high, and middle points are bigger than the target error tolerance (tol), the Algorithm 2 will start by adding more expansion points. First, it adds two expansion points at the step 2(a). In the case the tolerances of the new low Algorithm 2: Adaptive three points selection algorithm (A3PSA).

1. Initialization: given tol, q, s min = 2π × f min , s max = 2π × f max , tol h = 1, tol l = 1, tol mid = 1; e l = s min , e h = s max , e mid = 10 (log (el)+log (eh))/2 ; Krylov subspace Um = Ø 2.

While •

U m = U m ∪ spane mid → compute Ĝ, Ê, λ = eig( Ĝ-1 • Ê) • Computation as 2(b)
and high points are still bigger than the target error tolerance (tol), the new middle point will be added in the step 2(b).

The expansion point is automatically added at each step: 2 points for minimal and maximal eigenvalues and then 1 middle point if it is necessary. If the last three expansion points converge within the target error tolerance, the algorithm will stop and we have a nearly optimal reduced model in the sense of ||H|| 2 minimization. As in [START_REF] Freund | Recent advances in structure-preserving model order reduction[END_REF], this technique makes a smaller reduced model as PRIMA in the same accuracy.

As in [START_REF] Feng | A fully adaptive scheme for model order reduction based on moment-matching[END_REF], the two frequencies [f min , f max ] have been chosen for initial expansion points. Ideally, the optimal choice of new expansion point would be the point that causes the largest error but the authors in [START_REF] Feng | A fully adaptive scheme for model order reduction based on moment-matching[END_REF] point out that the choice of new expansion point as the midpoint of others two old expansion points in linear scale cannot warranty the largest errors at this point. In our experience in PEEC modeling, the next expansion point or midpoint should be chosen in logarithmic scale to minimize large error overall. The approach in [START_REF] Feng | A fully adaptive scheme for model order reduction based on moment-matching[END_REF] presents an adaptive scheme of choosing moments but many computation efforts are needed to compute the original model at new expansion point. In our case, the usage of compression algorithm (FMM) cause the exact value of original model even does not exist. In order to reduce the computation time and because of matrix compression property (FMM), in A3PSA algorithm, we use eigenvalues as the stopping criterion and error estimator.

Because of the properties of Krylov subspace and Taylor expansion, it is probable that all the eigenvalues in this range are converged within a target tolerance. This algorithm has the same property of IRKA such as: it converges fast after a few steps so the number of expansion points is small. The number of expansion point depends on the size of q-order. We will discuss about this value in the result section. We employ real expansion point s 0 for two reasons. First, the matrices in ( 15), E and G, to which the Krylov subspace method is applied, are real, thus avoiding the usage of complex arithmetic, which is four times as costly as real arithmetic. Second, when passive reduced-order models are constructed via the projection approach described in Section III-A, the projection matrix U m needs to be real in order to obtain real reducedorder matrices [START_REF] Villena | ARMS-automatic residue-minimization based sampling for multi-point modeling techniques[END_REF]. On the other hand, the usage of complex expansion points s 0 typically results in a significantly smaller state-space dimension of the reduced-order models, since s 0 can be placed closer to the frequency range of interest than any real s 0 [START_REF] Freund | SPRIM: Structure-preserving reduced-order interconnect macromodeling[END_REF].

In the next section, we will present the result for our test problem.

V. NUMERICAL RESULTS

We modeled an EMC filter in the software InCa3D [START_REF]InCa3D Software Available at: www[END_REF] (see Fig. 4) to determine the equivalent impedance.

The state space equation of PEEC discretization gives us an RLMC circuit consisting of 711 nodes, 4767 resistors, five capacitors, 4765 inductors, 4762 × 4762 inductive couplings. The frequency range of interest is from 1 Hz to 100 MHz.

The equivalent impedance is calculated by the InCa3D software and with the compression algorithm AMLFMM. Two results are compared in the Fig. 5, the maximal error relative is 1.3%.

In comparison in time and memory of computation with the InCa3D version without FMM, we have the Table I. Let us notice that we use a server with 16-GB memory and four processors Intel Xeon @3.0 GHz for InCa3D version (direct solver) while our reduction method with FMM was computed on Intel Core 2 Duo @2.66 Ghz with 2-GB memory. Matrix-vector products were computed by a parallel version of FMM developed in Java. Despite of using FMM, the computation time for this range is still too high (see Table I). The memories used is reduced by 5.4 times from 4900 MB (non-FMM) to 900 MB (with FMM), the runtime decreased from 22 601 s (direct solver LU) to 3500 s (iterative solver GMRES).

To eliminate the error caused by AMLFMM, we take the result by using compression algorithm as a new reference to compare with our algorithm A3PSA (see Algorithm 2). The accuracy of GMRES will be chosen at 10-12 for resolution of linear systems with point expansion s 0 = 0. We chose the coeff equals 0.05 to sparse the matrix E near .

We run our algorithm with different number of moment matching (q-order) at each expansion point, the tolerance 1% is chosen. The results later (see Table II) show the size of reduced model, the expansion points created for each Krylov subspace generated, the runtime for generating reduced model.

Let us notice that the bigger value of q is used, the smaller number expansion point is obtained. For example, q = 10, q = 8, we need three expansion points. But q = 5, we need five points and with q = 3, nine points are needed. In Table II, we can see the strategy of adding two points (in parenthesis) or one point. The peak points are well calculated. For all the result obtained, the maximal error is less than the fixed tolerance 1%. The highest computation time is 609 s at q = 7, it reduces with another size of expansion point, which is smaller than 3500 s as classical computation. Our approach with high accuracy can save not only the time of computation but also the memory requirement.

In the case of different tolerance (no tol, tol = 1E-2 and tol = 1E-4) and fixed q = 8, we obtained Table III and Fig. 6. If no tolerance is chosen, our algorithm use only two expansion points (f min and f max ), after 235 s of computation. If tol equals 1E-4, the A3PSA algorithm stops at six expansion points. It takes 841 s to obtain MOR size of 53.

At this value of q = 8, our algorithm is fast converged after a few steps and the resulted model is accuracy as fixed tolerance of choice. In general, our algorithm can generate rapidly a nearly optimal ROM for a given q-order. In comparison with FMM version, we reduce the runtime from 3500 to 411 s, while using the same amount of memory (900 MB).

VI. CONCLUSION

In this paper, we have proposed a multipoint rational approximation based on IRKA method for inductive PEEC large circuits compressed with the FMM in the frequency domain. An original algorithm for selecting expansion point automatically is also introduced. The size of reduced order model is nearly optimal for a fixed q-order. The reduced-order model is more accurate than one expansion point in PRIMA method for a wide spectrum of frequency. There are still some open problems like the adaptive choice of moment matching number or the choice of complex expansion points, which we are working on it. Regarding the PEEC formulation, one of perspective research is to introduce capacitive and delayed effects [START_REF] Ferranti | Multipoint full-wave model order reduction for delayed PEEC models with large delays[END_REF] in our approach to extend application domain.

To reuse themodel in the time-domain circuit simulation context, a method based on RLMCSYN approaches [START_REF] Nguyen | Circuit realization method for reduced order inductive PEEC modeling circuits[END_REF] is needed. The algorithm enables the preservation of PEEC RLM blocks structure and input and output incidence matrices. The final synthesized circuit can then be exported to a SPICE-like solver for the simulation in time domain [START_REF] Nguyen | Circuit realization method for reduced order inductive PEEC modeling circuits[END_REF].

M

  anuscript received Novem ber 9, 2016; revised Novem ber 29, 2016; accepted D ecember 7, 2016. This research is funded by Vietnam N ational Foundation for Science and Technology Development (NAFO STED) under grant num ber 102.01-2014.39. (Corresponding author: Dr. Tung Le Duc.)

Fig. 1 .

 1 Fig. 1. Taylor expansion around s = 0.

Fig. 2 .

 2 Fig. 2. Taylor expansion around several points s.

Fig. 3 .

 3 Fig. 3. Properties of system matrix (a) G matrix (b) E matrix with FMM compressed.

  (tol h > tol || tol l > tol || tol mid > tol) do (a) U m = U m ∪ spane l , e h → compute Ĝ, Ê, λ = eig( Ĝ-1 • Ê) (b) e lnew = min imag 1 λ ; e hnew = max imag 1 λ e m idnew = 10 lo g ( e l n e w ) + lo g ( e h n e w ) 2 tol l = | e lnewe l e lnew |; tol h = | e hnewe h e hnew |; tol m id = | e m idnewe m id e m idnew |; e m id = e m idnew ; e l = e lnew ; e h = e hnew (c) If (tol l < tol || tol h < tol)

Fig. 4 .

 4 Fig. 4. EMC filter structure in InCa3D software.

Fig. 5 .

 5 Fig. 5. Comparison of impedances calculated by InCa3D without FMM and FMM version (left) and relative error (right).

Fig. 6 .

 6 Fig.6. Relative error of different tolerance, q = 8.

  , f m a x = 10 8 ), f 1 = 4.309 × 10 5 , (f 2 = 7.02 × 10 4 , f 3 = 7.16 × 10 7 ), f 4 = 2.28 × 10 6 , (f 5 = 7.03 × 10 4 , f 6 = 8.08 × 10 7 ), f 7 = 2.35 × 10 6

				TABLE II
			Q-ORDER AND EXPANSION POINTS
	q	MOR size Time (s)	Expansion points (frequency)
	10	32	491	(f m in = 1, f m a x = 10 8 ), f 1 = 2.36 × 10 6
	8	26	411	(f m in = 1, f m a x = 10 8 ), f 1 = 2.37 × 10 6
	7	39	609	(f m in = 1, f m a x = 10 8 ), f 1 = 2.57 × 10 6 , (f 2 = 7.03 × 10 4 , f 3 = 7.96 × 10 7 )
	5	29	465	(f m in = 1, f m a x = 10 8 ), f 1 = 2.34 × 10 6 , (f 2 = 7.00 × 10 4 , f 3 = 7.91 × 10 7 )
	4	24	398	(f m in = 1, f m a x = 10 8 ), f 1 = 2.06 × 10 6 , (f 2 = 7.06 × 10 4 , f 3 = 7.94 × 10 7 )
	3	35	582	(f m in = 1

  in = 1, f m a x = 10 8 ), f 1 = 2.37 × 10 6 1E-4 53 841 (f m in = 1, f m a x = 10 8 ), f 1 = 2.57 × 10 6 , (f 2 = 7.03 × 10 4 , f 3 = 7.97 × 10 7 ), f 4 = 2.37 × 10 6

				TABLE III
			DIFFERENT TOLERANCE AND EXPANSION POINTS
	Tol	MOR size Time (s)	Expansion points (frequency)
	-	17	235	(f m in = 1, f m a x = 10 8 )
	1E-2	26	411	(f m

He is currently a CNRS Senior Researcher (directeur de recherche CNRS) and leads the "Models, Methods and Methodologies Applied to Electrical Engineering" research group of G2ELab, University of Grenoble Alpes, Grenoble. He has been a Visiting Scientist at the Technology Centre of the Federal University of Santa Catarina, Brasil, in 2012-2013. His research activities focus on computational electromagnetism (finite element and integral methods), equivalent magnetic sources identification by solving inverse problems and low magnetic fields measurements. He has coauthored more than 140 papers published in international journals and conference proceedings.

Gérard


Riccardo Manno
email: r.manno@bristol.ac.uk

Josiane Nguejio
email: josiane.nguejio@insa-rouen.fr

Fabrice Barbe
email: fabrice.barbe@insa-rouen.frd

Ivano Benedetti
email: ivano.benedetti@unipa.it

An investigation into the fracture behaviour of honeycombs with density gradients

In this study we perform an experimental and computational investigation about the fracture behaviour of polymer honeycombs presenting gradients in terms of lattice density. Such lattice relative density variations are introduced with the aim of mimicking the micro-morphology encountered in some natural materials, such as several kinds of woods, which seems related to the ability of the corresponding macro-material to delay the propagation of fracture under certain conditions. Starting from the conclusions of previous computational analyses, we perform a few experimental tensile tests on ABS model honeycombs obtained by additive manufacturing, with the aim of getting insights into their fracture behaviour and assessing the effect of the density gradients on the failure process, with respect to the behaviour observed in baseline homogeneous lattices. Following the performed tests, novel finite element analyses are performed, to help explain the observed failure processes and as preliminary calibration for further investigations addressed at maximising the lattice fracture toughness under tensile loading. With the emerging of more reliable and affordable additive manufacturing technologies, the present study contributes to exploring the possibility of modifying the effective properties of honeycombs and lattice materials through sensible small modifications of the micro-morphology.

INTRODUCTION

The development and employment of materials with suitable balances of sometimes contrasting properties, such as light weight and strength, stiffness and toughness, has always been inherently connected with the advancement of engineering and technology [START_REF] Fleck | Micro-architectured materials: past, present and future[END_REF]. In the last few decades however, a special sensitivity has arisen in some sectors, e.g. automotive and aerospace, towards the development of materials able to deliver their structural performance, possibly together with some other function, in a multi-functional perspective, at the lowest possible wight. In such fields, lower weight directly implies reduced fuel consumption at the same payload level, with obvious perceived benefits in terms of operational capability and environmental impact.

In the sketched background, composite materials provide a remarkable example, confirmed by an impressive rise in the number of academic and industrial studies and subsequent applications. For structural and multifunctional lightweight applications, increasing interest is also being attracted by the development or enhancement of other kinds of materials: metallic foams [START_REF] Banhart | Manufacture, characterisation and application of cellular metals and metal foams[END_REF] and micro-lattices [START_REF] Rashed | Metallic microlattice materials: A current state of the art on manufacturing, mechanical properties and applications[END_REF] offer an array of interesting properties such as low density, relatively high stiffness, high energy absorption, mechanical damping, acoustic and thermal absorption, tunable porosity. In 2015, the announcement of the creation of the "world's lightest metal", in fact a micro-lattice, by Boeing in collaboration with the HLR Laboratories made headlines [START_REF] Forbes | Boeing creates the world's lightest metal[END_REF], confirming the growing attention in the field. Despite their appealing features, a wider deployment of such materials in structural applications has been so far hindered by a certain scatter of their mechanical properties [START_REF] Ramamurty | Variability in mechanical properties of a metal foam[END_REF][START_REF] Smith | Steel foam for structures: A review of applications, manufacturing and material properties[END_REF], often due to a lack of control in the manufacturing processes, which poses a certain difficulty in the certification processes required for critical operations.

However, the recent fast developments in Additive Manufacturing (AM) technologies, with growing levels of quality control, have shown the possibility of manufacturing materials with multifunctional properties at a remarkable level of scale resolution [START_REF] Frenzel | Three-dimensional mechanical metamaterials with a twist[END_REF]. The use of such technologies, combined with structural optimisation tools [START_REF] Robbins | An efficient and scalable approach for generating topologically optimized cellular structures for additive manufacturing[END_REF], offers unprecedented opportunities to fill voids in the material design space, paving the way towards the developments of architectured materials [START_REF] Brechet | Architectured materials: Expanding materials space[END_REF] delivering sought performances as a consequence of their micro-architecture [START_REF] Frenzel | Three-dimensional mechanical metamaterials with a twist[END_REF].

While the use of such novel technologies may offer new opportunities for materials manufacturing, the materials design space to be explored, before delivering materials with the desired properties, might be too large and costly, even with the employment of optimisation tools. In this respect, Nature is often regarded as an invaluable source of inspiration, providing countless examples of materials able to deliver outstanding properties, thanks to an optimised hierarchical multi-scale arrangement of their elementary constituents [START_REF] Meyers | Biological materials science: biological materials, bioinspired materials, and biomaterials[END_REF], thus suggesting the possibility of mimicking such strategies in a bio-inspired materials design paradigm.

In this work, we present some preliminary results about the fracture performance of honeycombs exhibiting lattice density variations, similar to those observed in several kinds of woods and related to the ability of slowing crack propagation. The study extends the work performed in Refs. [START_REF] Manno | A computational study on crack propagation in bio-inspired lattices[END_REF][START_REF] Manno | Engineering the crack path in lattice cellular materials through bio-inspired micro-structural alterations[END_REF] with some experimental results about the tensile fracture behaviour of ABS specimens, obtained by fused deposition modelling, and their finite element interpretation.

TENSILE FAILURE EXPERIMENTAL TESTS

The investigation presented here moves from the computational analysis carried out in a previous study. In particular, an experimental campaign on ABS lattices is started with the aim of assessing the observations put forward in Ref. [START_REF] Manno | Engineering the crack path in lattice cellular materials through bio-inspired micro-structural alterations[END_REF] and, ultimately, evaluating the effects of lattice density gradients on the tensile failure behaviour of this kind of materials/structures.

Following preliminary computational indications, lattices identified by the ratio R = 0.6 are considered, as they should present the highest increment of fracture energy with respect to the reference homogeneous lattices. It is reminded that, in a lattice with alternating bands of different relative density, if l c is the width of the lower density bands and l f the width of higher density bands, the ratio R is defined by

R = l f l c + l f . (1)
From the geometrical point of view, besides selecting the value of R, the test specimens have been designed considering the constraints related with the printing technology, the loading device and the tensile mechanical conditions. In particular: a) the overall specimens dimensions, as well as the characteristic metrics of the lattice, i.e. the size of the cells and walls, the number of cells, the number of bands, have been set as a trade-off between the building plate dimensions and the printing resolution; b) in the printed specimens the lattice region is attached to two non-lattice full parts -see Fig. 1 -designed and printed so to provide suitable grips to the loading machine and to transmit quasi-uniform vertical displacements to the lattice upper and lower ends upon loading; c) the notch length has been determined, from a preliminary finite element analysis, as the minimal length ensuring that the stress in a vertical cell wall reaches a maximum value at the notch tip, thus avoiding crack initiation elsewhere in the lattice.

The geometry and loading of the specimens is further described in Fig. 2. Given the geometry and the position of the grips, the un-notched specimen deforms as shown in Fig. 2a, which results from the finite element simulation of the tensile test, assuming isotropic linear elastic properties for the material constituting the whole specimen. In this situation, due to the imperfect transmission of the displacement boundary conditions, higher values of axial stress are established in the central region of the specimen, as shown by the contour plot of σ yy in Fig. 2a. On the other hand, Fig. 2b shows the stress distribution in a notched specimen, illustrating how the notch length has been selected, so to ensure that fracture initiates from the notch tip, due to the higher value of σ yy in vertical walls located in its proximity.

From the manufacturing point of view, the specimens have been realised by ABS Filament Fusion Fabrication on a 3DGence Industry F340 3D printer. The infill density and the layer height have been set at 100% and 50 µm respectively.

To characterise the tensile properties of 3D-printed ABS, standard tensile test specimens have been printed in the same conditions as the lattice specimens and according to different orientations on the building plate, i.e. at 0°, 45°, 90°and 135°with respect to the front side of the plate. Following the procedure established in Ref. [START_REF] Dana | Polymer additive manufacturing of ABS Structure: influence of the deposition path on the tensile properties[END_REF] about the characterisation of 3D-printed-ABS properties, the tensile tests have been performed after 2-week storage in cool and dry atmosphere, to allow the relaxation of the material after processing, imposing a displacement rate of 1 mm/min. The tensile curves and the corresponding mean tensile properties are presented in Fig. 3. As classically observed concerning ABS, the stress response increases linearly with strain up to a maximum stress of about 30 MPa and then reaches a plateau eventually leading to fracture. As compared to the ultimate tensile stress, the elongation to fracture is highly dispersed, as also noted in Ref. [START_REF] Dana | Polymer additive manufacturing of ABS Structure: influence of the deposition path on the tensile properties[END_REF], and no significant effect of the building orientation could be observed, which suggests that the mechanical properties may reasonably be considered as isotropic at a macro-scale. Such isotropy can be explained by the fact that the filament deposition orientation is switched by an angle of 45°from layer to layer; this is programmed by the slicer software in order to favour isotropy -such as in a fibre-matrix laminatesand also to reduce the amount of porosity.

Summarising, lattice specimens with a 25 mm-long horizontal notch have been subjected to displacement controlled loading at the rate of 1 mm/min. Currently, 3 tests have been performed on notched specimens exhibiting density gradients and all the tests led to cracks propagating from the notch to the opposite side of the specimen, as shown in Fig. 4. The obtained experimental force-displacement curves are reported in Figs. [START_REF] Smith | Steel foam for structures: A review of applications, manufacturing and material properties[END_REF][START_REF] Frenzel | Three-dimensional mechanical metamaterials with a twist[END_REF] where they are also compared with the results of the finite element analysis discussed in the next section. It is observed that such curves show a good level of reproducibility, which is the necessary condition to exclude the possibility that printing defaults, necessarily present in the printed specimens, be the major structural characteristic governing the fracture propagation. In few of the manufactured specimens, undetectable printing defects initiated fracture events in regions far from the notch tip; such specimens have been thus excluded from the analyses.

A certain number of other tests, originally planned as part of this study, could not be performed as a consequence of the lockdown enforced throughout Europe to contrast the current outbreak of the CoViD-19 pandemic [START_REF] Sohrabi | World health organization declares global emergency: A review of the 2019 novel coronavirus (covid-19)[END_REF]. In particular, we could not perform tests on reference baseline homogenous lattices; for such a reason, in this study, only the response of honeycombs with density gradients has been experimentally investigated, while the response of homogeneous lattices has been estimated both analytically and numerically, without experimental validation. The presented results should then be intended as preliminary; however, the study, up to the forced interruption, has allowed the establishment a sound workflow for specimen manufacturing, testing and simulation, which will be useful in the planned continuation.

COMPUTATIONAL FAILURE ANALYSIS

The finite element analysis has been performed in ABAQUS/Explicit, identified as the suitable tool to investigate the highly non-linear processes of crack initiation, propagation and failure in polymeric lattices. The geometry of the investigated honeycombs, their mesh and the sets of appropriate boundary conditions have been generated using a dedicated Python pre-processor developed in-house. In particular, for reducing computational costs, only the geometry of the cellular portion of the specimens has been generated, while the full parts manufactured for easing the specimens clamping have been neglected. Hence, to mimic the test loading conditions, the boundary conditions have been applied enforcing a parabolic displacement profile at the edges of the virtual specimen. The geometry has been meshed using quad shell element with plane strain formulation (CPE4R) and element deletion has been enabled. The constitutive behaviour of the parent material has been initially calibrated using the experimental tests on the printed bulk ABS reported in Fig. 3. Specifically, the constitutive behaviour has been approximated in this stage as shown in Fig. 5, i.e. as an ideal material exhibiting a short flat plateau, after an initially elastic response, followed by a narrow softening segment.

The set of numerical results obtained from the modelling strategy described above are are reported in Figure 6, where they are compared with the experimental results available for the the tested honeycombs with density gradients and with the solution obtained from the analytical estimation of the stiffness of a honeycomb with homogeneous cellular density [START_REF] Fleck | Micro-architectured materials: past, present and future[END_REF]. It is observed that, while the slope of the initial linear elastic response is satisfactorily captured by FIGURE 6: Force VS displacement curves for honeycombs with R = 0 and R = 0.6. the numerical model, both the peak force and the failure strain are overestimated. The reasons behind such discrepancies might be related to a higher number of material defects in the printed honeycomb, with respect to the number of defects in the standard tensile test specimens used for inferring the material properties, and/or to the possible presence of a scale effect affecting the fracture behaviour of structural elements at different scales.

To get further insights into the influence of the parent material on the response of the honeycomb, a sensitivity analysis on the effect of the material constitutive properties is performed. Fig. 7 shows the load-displacement curves obtained when the parent material is modelled starting from the reduced constitutive curve plotted in Fig. 5. With such an assumption about the parent material, the slope of the elastic segment, the peak force and the strain at the onset of cracking are captured, which seems to confirm that the actual parent material is weaker that that inferred by the standard tensile tests. Furthermore, the crack path in the virtual specimen, reported in Fig. 8, is qualitatively compatible with what observed in the performed tests, see Fig. 4.

CONCLUSION

In this study, a sound workflow for the manufacturing, experimental tensile testing and computational analysis of honeycombs exhibiting density gradients has been established and some preliminary experimental and computational results about the tensile fracture behaviour of such cellular materials have been presented. Good experimental reproducibility of the fracture behaviour across the specimens has been observed: the curves are characterised by four main stages occurring at quite similar imposed displacements and leading to similar force evolution. Some discrepancy between the behaviour of the tested specimens and the simulated response has been identified, both in terms of peak load and failure strain, when the constitutive behaviour of the parent material is inferred starting from printed standard tensile specimens. On the basis of some sensitivity analysis on the influence of the parent material properties on the honeycomb failure, it has been hypothesised that such discrepancies may be related to the influence of inevitable manufacturing defects and/or the presence of a material scale effect. However, such aspects will be further investigated in the continuation of the present research.

FIGURE 1 :

 1 FIGURE 1: Example of specimen with density gradients for experimental analyses. The cell walls thickness is t = 0.54 mm. Left: side view. Right: front view.

FIGURE 2 :

 2 FIGURE 2: Spatial distribution of the axial (vertical) stress σ yy in the loaded non-homogeneous lattice as provided by FE simulation of the tensile test with assumed material isotropic linear elastic constitutive behaviour. Left: un-notched specimen. Right: notched specimen.

FIGURE 3 :

 3 FIGURE 3: Top: tensile stress-strain curves of different ABS samples printed with different orientations on the building plate. Bottom: material parameters extracted from the tensile tests.

FIGURE 4 :

 4 FIGURE 4: Post-failure configuration of ABS lattice specimens under progressive tensile loading; crack propagation initiated from the notch on the left hand side.

FIGURE 5 :

 5 FIGURE 5: Material constitutive model calibrated against the experimental tensile tests performed on printed bulk ABS.

FIGURE 7 :

 7 FIGURE 7: Force VS displacement curves for honeycombs with R = 0 and R = 0.6 using re-calibrated material properties.

FIGURE 8 :

 8 FIGURE 8: Crack path for a virtual specimen with R = 0.6.

ACKNOWLEDGMENTS

Ivano Benedetti acknowledges the support of INSA Rouen Normandie, where he has been invited as Foreign Visiting Professor between November and December 2019.