Pengfei Li

Yi Wu

Julien Yvonnet
email: julien.yvonnet@univ-paris-est.fr

A SIMP-Phase field topology optimization framework to maximize quasi-brittle fracture resistance of 2D and 3D composites

Keywords: Fracture resistance, Topology optimization, SIMP, Phase field method, two-material structures, 3D applications

A SIMP-Phase field topology optimization framework to maximize quasi-brittle fracture resistance of 2D and 3D composites

Introduction

Recently, design of composites with enhanced resistance to damage and fracture has emerged as a new challenging and exciting topic in computational mechanics. With the quick development of 3D printing and additive manufacturing techniques [START_REF] Tuan D Ngo | Additive manufacturing (3d printing): A review of materials, methods, applications and challenges[END_REF], new technological solutions which were not possible a few years ago can now be considered, like on-demand geometries of multi-materials structures and microstructures. More specifically, 3D printed bi-materials [START_REF] Anish Ravindra Amin | Dynamic response of 3d-printed bimaterial structure using drop weight impact test[END_REF][START_REF] Kao | Low-velocity impact response of 3d-printed lattice structure with foam reinforcement[END_REF][START_REF] Kao | Loading-unloading cycles of 3d-printing built bi-material structures with ceramic and elastomer[END_REF][START_REF] Kao | Bending behaviors of 3d-printed bi-material structure: Experimental study and finite element analysis[END_REF][START_REF] Wang | 3d printing of polymer matrix composites: A review and prospective[END_REF] offer new exciting possibilities such as designing composites with non-trivial periodic microstructures and ad-hoc functionalities. Among them, particle-matrix or skeleton-filled matrix composites able to increase the fracture resistance as compared to existing composites is of industrial and technological critical importance, for applications in aircraft, automotive or biomechanics, among many others.

One central ingredient for this task is the use of Topology Optimization (TO). TO has been widely studied and employed in both academic and engineering applications since it was firstly proposed by Bendsøe and Kikuchi [START_REF] Philip | Generating optimal topologies in structural design using a homogenization method[END_REF]. It is generally accepted that topology optimization methods can be mainly categorized into three kinds of families, namely the Solid Isotropic Material with Penalization (SIMP) method [START_REF] Martin P Bendsøe | Optimal shape design as a material distribution problem[END_REF][START_REF] Martin | Material interpolation schemes in topology optimization[END_REF][START_REF] George | Generalized shape optimization without homogenization[END_REF], the level set method [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Yu | A level set method for structural topology optimization[END_REF], and the Evolutionary Structural Optimization (ESO) method [START_REF] Yi | A simple evolutionary procedure for structural optimization[END_REF]. Review of these methods can be found in state-of-the-art papers [START_REF] Joshua | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF][START_REF] George | A critical review of established methods of structural topology optimization[END_REF][START_REF] Nico P Van Dijk | Level-set methods for structural topology optimization: a review[END_REF][START_REF] Zhu | Topology optimization in aircraft and aerospace structures design[END_REF]). An interesting comparison review on these techniques, with advantages and drawbacks, can be found in [START_REF] Sigmund | Topology optimization approaches[END_REF].

Initially mainly applied to maximize the stiffness of structures under volume constraints, TO has since then been extended and applied to a broad variety of problems like mechanical and thermal loads of structures, fluid flow, dynamics, acoustics and biomechanics, among many others (see a survey paper in [START_REF] Joshua | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF]).

More recently, TO has been applied to enhance the mechanical resistance of structures and materials subjected to damage and cracking.

One pioneering work combining TO and fracture mechanics can be traced back to Challis et al. [START_REF] Vivien | Fracture resistance via topology optimization[END_REF], where the level-set TO method has been used to maximize fracture resistance of structures, defined as the elastic energy released by the crack that are in tension and was calculated using the virtual crack extension. However, in this work, no real crack tips with stress singularity in the context of linear fracture mechanics were modeled. Another related technique can be found in [START_REF] Kang | Topology optimization considering fracture mechanics behaviors at specified locations[END_REF], where pre-defined cracks were inserted and TO used to minimize the J-integral around crack tip singularity as a fracture criterion.

An important progress was made by optimizing the topology while taking into account the incremental damage response of the structure during a full load, from initiation to damaged/cracked structures. A first series of works have been proposed where damage mechanics was considered during the TO problem. In [START_REF] Amir | A topology optimization procedure for reinforced concrete structures[END_REF][START_REF] Amir | Reinforcement layout design for concrete structures based on continuum damage and truss topology optimization[END_REF], Amir and Sigmund used a gradient enhanced model to define the truss topology and optimal cross sections of reinforcement bars, and where the external work was maximized as an evaluation of the fracture energy. In [START_REF] Kato | Multiphase layout optimization for fiber reinforced composites considering a damage model[END_REF], Kato and Ramm investigated optimal placement and shape of reinforcement in composites with respect to damage criterion to optimize the structural ductility during a full incremental procedure. In [START_REF] Frédéric | Optimization of multiphase structures considering damage[END_REF], Hilchenbach and Ramm optimized the position and shapes of stiff inclusions during the nonlinear loading of a structure. Even though not implying TO, the work shares similarities with the procedures developed in TO in this context. In [START_REF] Kai | Failure mitigation in optimal topology design using a coupled nonlinear continuum damage model[END_REF], James and Waisman developed TO within a nonlocal damage framework using SIMP where TO was performed with respect to a maximum damage criterion during the whole nonlinear load of a structure while minimizing the quantity of material within the structure as an objective. In [START_REF] Liu | Multi-material topology optimization considering interface behavior via xfem and level set method[END_REF], Liu et al. investigated TO with respect to damage induced by debonding at the interfaces between two materials in a structure during an incremental loading with Level-set TO and XFEM description of interfaces location. In [START_REF] Li | Design of fracture resistant energy absorbing structures using elastoplastic topology optimization[END_REF][START_REF] Li | Topology optimization of energy absorbing structures with maximum damage constraint[END_REF][START_REF] Li | Failure resistant topology optimization of structures using nonlocal elastoplastic-damage model[END_REF], Li et al. proposed SIMP TO using an elastoplastic-damage model where TO was performed to find the optimal structural topologies with high energy absorption capacity while constraining the damage indicator.

An important recent step was to include brittle fracture propagation [START_REF] Seleš | Crack propagation prediction in heterogeneous microstructure using an efficient phase-field algorithm[END_REF][START_REF] Wu | Comprehensive implementations of phase-field damage models in abaqus[END_REF] within TO. In [69] Zhang et al. used TO with XFEM and took into account full crack propagation to optimize several indicators such as tensile stress during loading of a structure. One issue with XFEM as a crack propagation simulation tool is its complexity in 3D and its inability to initiate the cracks, as well as handling multiple, complex crack networks which may connect and merge during the process. The development of the variational approach to fracture [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Bourdin | The Variational Approach to Fracture[END_REF][START_REF] Gilles | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Tung Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF][START_REF] Tung Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF][START_REF] Zhou | Phase field modelling of crack propagation, branching and coalescence in rocks[END_REF], also called phase field method in the literature, offers several advantages in the context of crack propagation simulation, like: the possibility to initiate cracks from undamaged configurations; the possibility to handle arbitrary crack networks (including branching, merging, in both 2D and 3D) without specific treatment and use of classical finite elements; a variational framework allowing to include many models or mechanisms, and a mesh-independence due to an appropriate regularization process. This point is of special interest in TO approach where the use of a fixed mesh is required. In [START_REF] San | Optimization of carbon black polymer composite microstructure for rupture resistance[END_REF], San and Waisman combined phase field and genetic algorithms to find the optimal location of particles in order to maximize indicators such as the peak force, maximum deformation at failure point and maximum fracture energy during an incremental procedure. The first works to our knowledge combining phase field and TO was introduced in Xia et al. and Da et al. in [START_REF] Da | Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage[END_REF][START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF], where the BESO TO [START_REF] Huang | Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[END_REF] was used to optimize the fracture resistance of two-phase structures with respect to inclusions shapes, including cracks in both bulk and interfaces, and applied to periodic composites and multiple loads in [START_REF] Da | Topology optimization for maximizing the fracture resistance of periodic quasi-brittle composites structures[END_REF]. In [START_REF] Russ | A novel topology optimization formulation for enhancing fracture resistance with a single quasi-brittle material[END_REF][START_REF] Russ | Topology optimization for brittle fracture resistance[END_REF], Russ and Waisman developed a SIMP TO combined with phase field to optimize the fracture energy in one-phase material structures, and Wu et al. [START_REF] Wu | Level-set topology optimization for maximizing fracture resistance of brittle materials using phase field fracture model[END_REF] developed a Level-Set TO-phase field approach to optimize the fracture resistance of composites.

In this work, we extend our previous BESO framework [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF] to a combined SIMP and phase field for maximizing the fracture energy of two-phase composites. It is worth noting that the present framework shares many similarities with the recent framework of [START_REF] Russ | A novel topology optimization formulation for enhancing fracture resistance with a single quasi-brittle material[END_REF]: phase field and SIMP are combined for maximizing fracture resistance and the objective function used in both of these works include the external work computed incrementally during a full crack propagation simulation. The main contribution here is the application of the SIMP-phase field to two-material structures, where the objective here is not to minimize the total weight as in [START_REF] Russ | A novel topology optimization formulation for enhancing fracture resistance with a single quasi-brittle material[END_REF] but to modify the topology of a second phase material under the constraint of a volume fracture. Another contribution is a comparison of the present SIMP results with BESO formulations with respect to convergence and performance (value of objective function). Finally, the last contribution is applications of this framework to 3D 2-material structures which show the full potential of the approach.

The present paper is organized as follows. In Section 2, we review the phase field method for brittle fracture based on the variational framework and provide the details of the finite element discretization. In Section 3, we propose the SIMP topology optimization approach for the design of brittle composites to maximize the fracture resistance in two-phase composites. Section 4 provides the overall algorithms to practically implement the present framework. Finally, a series of 2D and 3D numerical examples are presented in Section 5 to show the efficiency and potential of the approach.

Phase field fracture formulation

In this work, we briefly review the phase field method for fracture which will serve as one main ingredient in the present topology optimization framework.

Let Ω ⊂ R D be an open domain with D = 2, 3, describing a cracked solid as depicted in Fig. 1. The external boundary of Ω is denoted by ∂Ω ∈ R D-1 . Cracks which may propagate within the solid are collectively denoted by Γ. In this work, we adopt the framework proposed in [START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via t-convergence[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF] for a regularized representation of discontinuities. In this regularized framework, the crack is approximately represented by a scalar phase field 0 ≤ d (x, t) ≤ 1 (see Fig. 1(b)), such that when d = 0 the material is undamaged, and when d = 1 the material is fully broken.

Phase field approximation of crack

The scalar phase field d (x, t) can be determined through solving the following boundary value problem subjected to Dirichlet boundary conditions d = 1 on the crack (see [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] for more details):

           d (x, t) -2 ∆d (x, t) = 0 in Ω, d (x, t) = 1 on Γ, ∇d (x, t) • n = 0 on ∂Ω, (1)
where ∆ (•) and ∇ (•) are the Laplacian and gradient operator respectively, is a length scale parameter that governs the width of the regularization zone and gives for → 0 the exact sharp crack in Fig. 1(a), and n the outward normal on ∂Ω. It can be shown that ( 1) is the Euler-Lagrange equation associated with the variational problem:

d = Arg inf d∈S d Γ d (d) , Γ d (d) = Ω γ d (d) dΩ, S d = {d | d (x) = 1, ∀x ∈Γ} , (2)
where Γ d (d) represents the total length of the crack in 2D and the total crack surface area in 3D, and γ d (d) is the crack surface density function per unit volume defined by:

γ d (d) = d 2 2 + 2 ∇d • ∇d, (3)
where the second term in γ d (d) penalizes high values of ∇d(x) and where d varies between 0 and 1. Note that in the absence of the second right-hand term in (3), a local damage model is found, with well-known related non-convergence issues with respect to the mesh discretization.

It must be noted that, does not represent physically the exact crack width, but a parameter which is used to regularize the discontinuities. It has been shown that this parameter can be treated as a material parameter related to the Young's modulus, the tensile strength, and the critical energy release rate of the material in [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Tung Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Zhang | Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale[END_REF]. In our previous work [START_REF] Tung Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF], an inverse approach was developed to identify this parameter, by combining simulations and experiments. More recently in [START_REF] Molnár | Toughness or strength? regularization in phase-field fracture explained by the coupled criterion[END_REF], G. Molnár et al. assessed its mechanical significance with the help of the coupled criterion.

Variational framework

We follow the framework presented in [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF] to construct the variational principle, which involves: irreversibility condition, stability condition and energy balance. In this framework, stability condition provides mechanical balance equation and damage criteria. The energy balance provides damage consistency. The total energy for a cracked body is defined as

W (u, d) = Ω ψ e (ε (u) , d) dΩ + Ω ψ d (d) dΩ - ∂Ω F F • udS - Ω f • udΩ, (4)
in which u is displacement field, ε (u) = 1 2 ∇u + ∇u T , f and F are body forces and prescribed traction over the boundary ∂Ω F , respectively. Above, ψ d is the damage dissipative potential defined as

ψ d = g c γ d (d) , (5)
where γ d (d) is defined in (3), and g c is the Griffith-type critical energy release rate. ψ e is the stored elastic energy density function defined as [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] ψ e (ε,

d) = (1 -d) 2 + k ψ + e (ε) + ψ - e (ε) (6)
in which k is a small numerical parameter to prevent loss of definite posedness of the elastic tensor in case of full damage, ψ + e and ψ - e are the tensile and compressive energies,

ψ ± e = 1 2 λ tr [ε] 2 ± + µtr ε ± 2 , (7)
with λ and µ the standard lame parameters. Note that only tensile damage degradation is taken into account in the elastic energy density (6) through a decomposition of the elastic strain ε into tensile/positive and compressive/negative parts:

ε = ε + + ε -with ε ± = 3 i=1 ε i ± n i ⊗ n i , (8)
where x ± = 1 2 (x ± |x|), ε i and n i are the eigenvalues and eigenvectors of ε.

Irreversibility condition

The irreversibility condition is imposed on the damage variable to avoid material regeneration. It can be expressed as

ḋ ≥ 0, 0 ≤ d ≤ 1 (9)
which is ensured by using an appropriate history function [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF].

First-order stability condition

The first order stability condition (see [START_REF] Mielke | Evolution of rate-independent systems[END_REF][START_REF] Mielke | A mathematical framework for generalized standard materials in the rate-independent case[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF]) is expressed by:

D δu W (u, d) + D δd W (u, d) ≥ 0, (10)
where

D v f (u) = d dh f (u + hv) h=0 . ( 11
)
is the directional derivative. Applied to (4), it yields:

Ω σ : ε (δu) dΩ + Ω ∂ψ e ∂d + ∂ψ d ∂d δddΩ - ∂Ω F F • δudS - Ω f • δudΩ ≥ 0 ( 12
)
where

σ = ∂ψ e ∂ε = (1 -d) 2 + k λ tr [ε] + 1 + 2µε + + λ tr [ε] -1 + 2µε - (13)
in which 1 is the second-order identity tensor and σ + /σ -are the undamaged tensile/compressive stress tensors. From ( 12), the following results stem out:

For δd = 0, find u ∈ S u , S u = u|u(x) = ū on ∂Ω u , u ∈ H 1 (Ω) such that: Ω σ : ε (δu) dΩ - ∂Ω F F • δudS - Ω f • δudΩ = 0, (14)
which is the weak form of the mechanical equilibrium equation, with δu ∈ S 0 u , S 0 u = δu|δu(x) = 0 on ∂Ω u , u ∈ H 1 (Ω) .

For δu = 0 we obtain:

Ω ∂ψ e ∂d + ∂ψ d ∂d δddΩ ≥ 0 (15)
which is the weak form of the damage criterion. In local form, it can be expressed as:

f d (d) = - ∂ψ e ∂d + ∂ψ d ∂d ≤ 0 in Ω. ( 16
)

Energy balance

The energy balance represents the need for the total energy to remain constant as the state variables evolve. Following a procedure analogous to the treatment of the stability condition, this condition leads to

Ω -σ : ε ( u) - ∂ψ e ∂d + ∂ψ d ∂d ḋ dΩ + ∂Ω F F • udS + Ω f • udΩ = 0. ( 17
)
For u = 0, and using ( 16), the damage consistency condition is obtained:

f d (d) ḋ = 0. ( 18
)

Alternate minimization

In this section, a staggered alternate minimization algorithm is applied, which naturally stems out from the energetic principles. This procedure takes advantage of the fact that although the global energy is non-convex, it is convex with respect to u and d individually [START_REF] Alessi | Variational approach to fracture mechanics with plasticity[END_REF][START_REF] Li | An extension of the phase field method to model interactions between interfacial damage and brittle fracture in elastoplastic composites[END_REF]. With the total energy (4) at hand, the alternate minimization follows.

Minimization with respect to the displacement field:

D δu W (u, d) = 0 ( 19
)
which leads to

R 1 = Ω σ : ε (δu) dΩ - ∂Ω F F • δudS - Ω f • δudΩ = 0 ( 20
)
which corresponds to the weak form of the mechanical problem to be solved for u, given d.

Minimization with respect to the damage field:

D δd W (u, d) = Ω ∂ψ e ∂d + ∂ψ d ∂d δddΩ = 0. ( 21
)
which corresponds to the global problem to be solved to find the field d(x) (phase field problem), given u.

Governing equations

The associated Euler-Lagrange equations to [START_REF] Da | Topology optimization for maximizing the fracture resistance of periodic quasi-brittle composites structures[END_REF] are given by:

           ∇ • σ + f = 0 in Ω, u = ū on ∂Ω u , σn = F on ∂Ω F . (22)
Using ( 21) and the property:

(∆d) δd = ∇ • (∇dδd) -∇d • ∇(δd) (23)
as well as the divergence theorem and ∇d • n = 0, we obtain the weak form of the phase field problem as:

Ω -2(1 -d)ψ + e + g c d δd + g c ∇d • ∇(δd) dΩ = 0. ( 24
)
To prescribe irreversibility, we employ the technique introduced in [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids[END_REF] which consists in substituting the above weak form with:

Ω {-2(1 -d)H(ε) + 2ψ c d} δd + 2ψ c 2 ∇d • ∇(δd) dΩ = 0, (25)
Table 1: Governing equations of the phase field model.

Irreversibility ḋ ≥ 0, 0 ≤ d ≤ 1 Mechanical balance ∇ • σ + f = 0 in Ω u = ū on ∂Ω u , σn = F on ∂Ω F Constitutive law σ = (1 -d) 2 + k λ tr [ε] + 1 + 2µε + + λ tr [ε] -1 + 2µε - Damage criterion ψ c d -2 ∆d -(1 -d) H(ε) ≥ 0 H (ε) = max s∈[0,t] ψ + e ( ε, s) -ψ c + Damage consistency ψ c d -2 ∆d -(1 -d) H(ε) ḋ = 0
in which

H (ε) = max s∈[0,t] ψ + e (ε, s) -ψ c + (26)
and ψ c is a specific fracture energy density, which can be further related to a critical fracture stress σ c by:

ψ c = 1 2E σ 2 c , (27)
where E is the Young's modulus.

The corresponding Euler-Lagrange equations to [START_REF] T Guillén-Hernández | In-situ strength effects in long fibre reinforced composites: A micro-mechanical analysis using the phase field approach of fracture[END_REF] are given by:

           ψ c d -2 ∆d = (1 -d) H (ε) ∇d • n = 0 on ∂Ω, d = 1 on Γ. (28)
The different equations of the model are summarized in Table 1. The FEM discretization details are provided in Appendix.

SIMP Topology optimization formulation

In this section, we present a SIMP topology optimization framework for maximizing the fracture resistance of a two-phase (composite) structure.

Optimization problem statement

The topology optimization problem is conducted with respect to a density variable ρ(x) which is associated with the inclusion phase. In other words, ρ(x) = 1 corresponds to inclusion phase and ρ(x) = 0 corresponds to the matrix phase.

For stability considerations, here displacement-controlled loading is adopted. For a prescribed displacement load, the fracture resistance maximization is equivalent to the maximization of the mechanical work. Recalling that the fracture problem is quasi-static, we introduce a pseudo time t associated with the external load evolution, with t ∈ [0, t max ], where t max denotes the maximum loading time corresponding to the maximum prescribed displacement u max at the failure step.

The optimization problem is then defined as follows:

Maximize : J(ρ, u, d) (29)
subject : R 1 (ρ, u(t),

d(t)) = 0, ∀t ∈ [0, t max ] (30) R 2 (ρ, u(t), d(t)) = 0, ∀t ∈ [0, t max ] ( 31
)
f inc = V (Ω inc ) V (Ω) = Ω ρ(x)dΩ V (Ω) (32) 0 ≤ ρ(x) ≤ 1 (33) u(t) ∈ S u (34)
d(t) ∈ S d , (35)
where V (Ω inc ) is the inclusion volume, V (Ω) is the total volume of the domain, and

J = t max 0 F ext (t) • u(t)dt, (36)
where R 1 and R 2 are given by ( 20) and ( 59), F ext is the external force response at the load point and f inc is the prescribed volume fraction of the inclusion phase. Following [START_REF] Mp Bendsøe | Theory, methods and applications[END_REF], the material interpolations for the two-phase material are defined as

   E(x) = (ρ(x)) p E inc + (1 -(ρ(x)) p ) E mat , ψ c (x) = (ρ(x)) p ψ c,inc + (1 -(ρ(x)) p ) ψ c,mat , (37)
where E and ψ c are the Young's modulus and the fracture energy density. (•) inc and (•) mat are the parameters corresponding to the inclusion and the matrix phase, respectively. The Poisson's ratios of the two material phases are assumed identical. Above, p is the penalty coefficient to enforce solutions close to ρ = 0 or 1. Following [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF], we choose p = 3. Fig. 2 shows the difference of Young's modulus interpolation for SIMP and BESO [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF]. The continuous material interpolations for SIMP ensures that the optimization problem is smooth and that the objective function is differentiable.

Discrete topology optimization problem

The design domain Ω is discretized into N e finite elements and each element e is assigned with a topology design variable, or element density ρ e , which is allowed to vary continuously in the range [0, 1]. Here, a density of ρ e = 1 corresponds to an element completely filled with the inclusion phase, whereas ρ e = 0 corresponds to an element completely filled with the matrix phase. We define the vector {ρ} = ρ 1 , ρ 2 , ..., ρ Ne containing the discrete values of densities in elements of the mesh.

The discrete form of ( 29)-( 35) is then defined as: Maximize :

0 0.2 0.4 0.6 0.8 1 E mat E inc SIMP BESO data3
J ∆u (ρ, u, d) (38)
subject :

K n u u n -F n u = 0, n = 1, 2, ..., n load (39)
u n ∈ S u , n = 1, 2, ..., n load (42)

d n ∈ S d , n = 1, 2, ..., n load , (43)
and J ∆u in ( 38) is approximated by:

J ∆u = n load n=1 ∆J n ≈ 1 2 n load n=1 F n ext + F n-1 ext T ∆u n , (44)
where ∆u n denotes the prescribed load increment at load n. Above, v e is the volume (area in 2D) of the e-th element and f inc is the target inclusion volume which is prescribed during the optimization process. The stiffness matrix K n u at the n-th load increment is constructed following [START_REF] Yi | A simple evolutionary procedure for structural optimization[END_REF].

It should be noted that in this work the continuous topology design variable ρ e ∈ [0, 1] allows using well-proven gradient-based optimization update approaches (e.g., Optimality Criteria methods (OC) [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF][START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF], the Method of Moving Asymptotes (MMA) [START_REF] Svanberg | The method of moving asymptotes a new method for structural optimization[END_REF] and so on), in turn ensuring algorithmic convergence within a reasonable number of topological iterations (from 10 to 1000 iterations) [START_REF] Sigmund | Topology optimization approaches[END_REF].

Another important remark is that in the present paper, we do not explicitly use the second constraint [START_REF] Kao | Loading-unloading cycles of 3d-printing built bi-material structures with ceramic and elastomer[END_REF] in the discrete problem ( 38)- [START_REF] Mielke | Evolution of rate-independent systems[END_REF] and in the following sensitivity analysis. In fact, this constraint is taken into account implicitly as the regularized fracture problem is solved in a staggered solving procedure. When the mechanical problem [START_REF] Da | Topology optimization for maximizing the fracture resistance of periodic quasi-brittle composites structures[END_REF] is solved, the damage variable d(x) is assumed to be known from previous iteration and then the phase field problem equation R 2 in ( 59) is supposed to be verified. This assumption considerably simplifies the sensitivity analysis presented next.

Sensitivity analysis

In order to solve the optimization problem ( 38)-( 43), the sensitivity of the objective function J corresponding to change in the design variable must be determined. The derivation of the sensitivity requires using the adjoint method (e.g., [START_REF] Buhl | Stiffness design of geometrically nonlinear structures using topology optimization[END_REF][START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF]). Assuming that the mechanical problem [START_REF] Wu | Comprehensive implementations of phase-field damage models in abaqus[END_REF] has been solved, we introduce the Lagrangian:

J ∆u ≈ J = 1 2 n load n=1 F n u + F n-1 u T ∆u n + (λ n 1 ) T R n + (λ n 2 ) T R n-1 (45)
in which R n and R n-1 are the residuals of (65) at n-th and (n-1)-th load increments, respectively. λ n 1 and λ n 2 are Lagrange multipliers which have the same dimension with the displacement vector u.

Notice that for displacement-controlled loading problem, displacement components at the boundary nodes and force components at the free nodes are fixed, hence, they are independent of the current value of ρ. Here we introduce a division of all degrees of freedom into essential (index E; associated with Dirichlet boundary conditions) and free (index F) nodal values. For a vector v and a matrix M we have

v ∼   v E v F   and M ∼   M EE M EF M FE M FF   . ( 46
)
We then have these unknowns at the n-th load increment

∂u n ∂ρ e =   0 ∂u n F ∂ρ e   , ∂∆u n ∂ρ e =   0 ∂∆u n F ∂ρ e   , F n u =   F n u,E 0   , ∂F n u ∂ρ e =   ∂F n u,E ∂ρ e 0.   (47)
Here for arbitrary load increment indices n = 1, ..., n load , m = 1, ..., n load , we have

∂ ∂ρ e (F m u ) T ∆u n = ∂F m u ∂ρ e T ∆u n + (F m u ) T ∂∆u n ∂ρ e = ∂F m u ∂ρ e T ∆u n . ( 48
)
With the above property at hand, the derivative of the objective function J in ( 45) is given by

∂J ∂ρ e = 1 2
n load n=1 ∂ F n u + F n-1 u T ∂ρ e ∆u n + (λ n 1 ) T ∂R n ∂ρ e + (λ n 2 ) T ∂R n-1 ∂ρ e . ( 49
)
The derivatives of R m at the equilibrium of the m-th load increment with respect to ρ e can be expanded as

∂R m ∂ρ e = ∂F m u ∂ρ e - ∂K m u ∂ρ e u m -K m u ∂u m ∂ρ e . (50)
Using ( 47) and ( 50), ( 49) can be reformulated as

∂J ∂ρ e = 1 2 n load n=1 ∂F n u,E ∂ρ e T ∆u n E + λ n 1,E + ∂F n-1 u,E ∂ρ e T ∆u n E + λ n 2,E -(λ n 1 ) T ∂K n u ∂ρ e u n + K n u ∂u n ∂ρ e -(λ n 2 ) T ∂K n-1 u ∂ρ e u n-1 + K n-1 u ∂u n-1 ∂ρ e . (51)
In order to eliminate the unknowns

∂F n u,E ∂ρ e and ∂F n-1 u,E
∂ρ e in (51), we choose

λ n 1,E = -∆u n E and λ n 2,E = -∆u n E . (52)
Then we can re-write [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF] as

∂J ∂ρ e = - 1 2
n load n=1 (λ n 1 ) T ∂K n u ∂ρ e u n + K n u,FE λ n 1,E + K n u,FF λ n 1,F T ∂u n F ∂ρ e + (λ n 2 ) T ∂K n-1 u ∂ρ e u n-1 + K n-1 u,FE λ n 2,E + K n-1 u,FF λ n 2,F T ∂u n-1 F ∂ρ e . (53)
To eliminate the unknowns 53), we choose

∂u n F ∂ρ e and ∂u n-1 F ∂ρ e in (
λ n 1,F = K n u,FF -1 K n u,FE ∆u n E and λ n 2,F = K n-1 u,FF -1 K n-1 u,FE ∆u n E . (54)
Using ( 52) and ( 54), we can obtain the final objective derivative

∂J ∂ρ e = - 1 2
n load n=1 (λ n 1 ) T ∂K n u ∂ρ e u n + (λ n 2 ) T ∂K n-1 u ∂ρ e u n-1 . (55)
For each element e, (55) can be re-written as

α e = n load n=1 (∆α e ) n = - 1
in which (•) e is the element component of (•) and k n u,e is the element stiffness matrix for element e at the n-th load step. For the sake of simplicity, we use α e and α to represent element sensitivity and vector of all element sensitivities. It should be noted that in this work we need to evaluate the sensitivities of all elements, unlike our previous works [START_REF] Da | Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage[END_REF][START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF] which only compute the sensitivities of the inclusion elements and set directly the sensitivities of the matrix elements to zeros. It is worth noting that the present framework can easily include interfacial damage by using the extended phase field framework proposed in [START_REF] Tung Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]. Such an extension in the context of BESO can be found in [START_REF] Da | Topology optimization for maximizing the fracture resistance of periodic quasi-brittle composites structures[END_REF][START_REF] Da | Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage[END_REF]. Such formulation allows using a fixed regular mesh, which is advantageous for TO algorithms. In contrast, other formulation have been proposed to include interfacial damage within Phase Field by using cohesive elements e.g. in [START_REF] T Guillén-Hernández | A micromechanical analysis of inter-fiber failure in long reinforced composites based on the phase field approach of fracture combined with the cohesive zone model[END_REF][START_REF] T Guillén-Hernández | In-situ strength effects in long fibre reinforced composites: A micro-mechanical analysis using the phase field approach of fracture[END_REF], but require meshing explicitly the interfaces. This might not be well suited to TO problems, where 3D complex inclusion geometries can be involved, as shown in the next examples.

Overall algorithm

This section present the overall algorithms of the proposed method. We first introduce the algorithm used to simulate one crack simulation from the initial nucleation to complete failure, which will be used within one step of the topology optimization. During this simulation, the sensitivities are computed at the same time than the fracture evolution. Then, the overall topology optimization is presented.

In the fracture evolution problem, a staggered scheme is employed following [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], where at each load increment the crack phase field problem is solved for fixed displacement field which is known from the previous time step. The displacement problem is then solved for the obtained crack phase field. Finally the sensitivity analysis is implemented after the staggered scheme.

The flowchart for fracture evolution problem and sensitivity calculations is provided in Algorithm 1.

Algorithm 1: Fracture evolution problem and sensitivity calculations.

which is determined by the prescribed filter radius r min and the element center-to-center distance ∆ (e, j) between element e and j.

The optimization problem ( 29)-( 35) could be solved using several different approaches such as OC method [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF][START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF] and the MMA [START_REF] Svanberg | The method of moving asymptotes a new method for structural optimization[END_REF] as illustrated in section 3.1. In this work, the OC method is employed to update the design variables. The overall algorithm is illustrated in Algorithm 2.

Err = | k m=k-4 J m -k-5 n=k-9 J n | k-5 n=k-9 J n else Err = 1 end k = k + 1 end

Numerical examples

In this section we present several examples in both 2D and 3D to show the capability and convergence of the method, and compare the results with our previous work done by BESO method [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF]. In all 2D examples uniform meshes of quadrilateral bilinear elements with the plane strain assumption have been employed. Uniform meshes of eight-node cubic elements have been used for the 3D examples. Both damage and displacement fields are discretized with the same finite element meshes. The characteristic length scale parameter for the phase field problem in (3) and the filter radius r min in [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF] are both set to be twice the typical finite element size = r min = 2h e . For the sake of clear visualization, only the crack phase field with values higher than 0.4 in 2D examples and values higher than 0.95 in 3D examples are plotted. The material properties are shown in Table 2. The inclusion volume fraction f inc is defined by [START_REF] Kao | Bending behaviors of 3d-printed bi-material structure: Experimental study and finite element analysis[END_REF]. In the different simulations, the following features are noted.

For SIMP topology optimization simulations, there is no initialization of the inclusion geometry. The initial density ρ 1 e is set to be uniform and equal to f inc .

For BESO topology optimization simulations, an initial geometry of inclusion is required to match f inc at the first iteration. If a homogeneous design is used, the volume fraction has to be reduced at each iteration by setting ρ e = 0 in more elements, as ρ e can only be zero or one in each element. We refer to this initial geometry to "initial BESO design".

To evaluate the improvement of the fracture resistance in the present SIMP context which does not require an initial geometry, a "guess" design is defined in some examples to evaluate the improvement of the fracture resistance between optimized and guess designs.

The volume fraction f inc is fixed during the whole optimization process.

To avoid interfering with the topology optimization process, a region around the initial crack is defined where the design variable are enforced to ρ e = 0 (remain matrix material). This region is defined as embedding all nodes at a distance 2 from the initial crack surface.

The incremental loading process goes on until the reaction force is below a prescribed value, indicating that the structure is fully broken.

2D reinforced plate with one pre-existing crack notch

In the following example, a 2D plate with one pre-existing crack notch, as described in Fig. 3, is considered. The structure is composed of a matrix material, and we seek the shape of an inclusion which provides the maximal fracture resistance for the whole composite structure. The dimensions of the plate are 100 × 50 mm 2 , and the domain is uniformly discretized into 120 × 60 square bilinear elements. The boundary conditions are as follows: on the upper and lower ends, the y-displacement are fixed, while the x-displacement are free. On the left and right ends, the horizontal displacements are prescribed with incremental displacement loads with ∆U = 0.004 mm for the first five load increments and ∆U = 0.001 mm for the following load increments. The preexisting crack is simulated by prescribing Dirichlet conditions on the crack phase field with d = 1 along the crack. Additionally, the optimality convergence tolerance is set at 10 -5 (tol = 10 -5 ). Fig. 4 shows the evolution histories of inclusion topologies and their final crack patterns. Here, the inclusion volume fraction is set to f inc = 5%. As can be observed from Fig. 4, the fracture resistance of the composite structure increases with accumulation of the inclusion material around the crack paths obtained from previous design iteration and reduction of intermediate densities (sometimes called "grey zones" in the topology optimization studies), and then converges to an almost constant value of 18.8 mJ. Detailed propagation of the crack phase field corresponding to its load-displacement curve for the optimized design is shown in Fig. 5. The crack propagates vertically into the inclusion material and two other cracks initiate around the left and right corners of the inclusion pattern, and then continue to propagate until the structure is fully broken.

In order to illustrate the robustness of the method for other inclusion volume fractions and investigate the influence of f inc on the numerical results, two additional simulations are conducted: f inc = 2% and f inc = 10%. Fig. 6 and Fig. 7 show the evolution histories of inclusion topologies and their final crack patterns for cases f inc = 2% and f inc = 10%, respectively. Fig. 8 shows comparison of the optimized designs and corresponding final crack patterns for cases f inc = 2%, f inc = 5% and f inc = 10%. As can be observed, with the increase of the volume fraction, significant changes and more complex shapes of inclusions are obtained, in tandem with an increase of the fracture resistance. To investigate the effect of finite element mesh size on the numerical results, two additional simulations are performed for f inc = 5% with finer meshes, involving 180 × 90 square bilinear elements, and 240 × 120 square bilinear elements, respectively. The filter radius r min is fixed and is equal to 1.667 mm, and the length scale = 2h e is changed with mesh refinement. Fig. 9 shows the comparison of optimized inclusion designs and corresponding final crack patterns for the three different mesh models. The corresponding load-displacement curves are shown in Fig. 10. From these two figures, it can be noticed that the optimized design changes with finer meshes, but the overall response of the structure as well as the optimized fracture energy does not change significantly. This suggests that simulations performed with reasonably fine meshes provide usable designs for practical applications. For this reason, in the following, we use the 120×60 mesh model. In the following, we compare the proposed SIMP approach with BESO topology optimization solution [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF]. Three additional simulations using the BESO method with the same volume fraction f inc = 5% but different initial designs (see Fig. 11) are performed. To fully display the convergence histories of BESO method, we do not prescribe a convergence criterion for these three BESO simulations and let the algorithm run until a maximum number of iterations equal to 500 is reached. Fig. 12 shows the comparison of convergence histories for BESO and SIMP solutions. We can note that in the case of BESO, the initial designs have a strong influence on the optimized objective function and on the convergence rate: (i) using BESO with design B, C or SIMP, the same convergence value 18.8 mJ is reached, while for initial design A with BESO the convergence value is about 19.6 mJ; (ii) solutions by BESO with initial designs B, C and SIMP solution require a similar number of convergence steps of about 200, while for initial design A with BESO the convergence number is about 400. From Fig. 12, we can also observe spurious oscillations for initial design A and C with BESO after a stable convergence value, while this does not seem to occur with SIMP. Fig. 11 shows the influence of initial design on the optimized design with BESO method. Fig. 13 shows the comparison of load-displacement curves for the 7 optimization processes. From Fig. 12 and Fig. 13, we can conclude that even though SIMP has slightly better convergence properties with less oscillations, the convergence of both methods is comparable when an initial design is set such as the target volume fraction f inc is met at the first iteration. The final results are also very close. The only difference is that starting from the target volume fraction f inc is simpler with SIMP, as a value can be set in each element between 0 and 1, while in BESO an initial guess design matching the target volume fraction is required, which may be not trivial to define for more complex geometries. However, an initial homogeneous design is also possible using BESO, but the actual volume fraction will be changed at each iteration to reach f inc as ρ e in each element can only be 0 or 1. We investigate this case in Figs 12 and 13 (solution referred to as "homogeneous design, BESO": we can see that the convergence curve (Fig. 12) shows much more oscillations for BESO than SIMP for a homogeneous design, even though we can note in Fig. 13 that the optimized fracture energy is higher for BESO.

3D reinforced sample with one pre-existing crack notch surface

In this example, the aim is to illustrate the applicability and convergence of the present method for 3D problems. For this purpose, a 3D sample with one pre-existing crack notch surface, as described in Fig. 14 As in the 2D case of section 5.1, the pre-existing crack is modeled by prescribing Dirichlet conditions on the crack phase field with d = 1 along the crack surface. On the upper and lower ends of the sample, vertical incremental displacement loads with ∆U = 0.005 mm are prescribed for the first four load increments and ∆U = 0.002 mm for the following load increments. Additionally, the optimality convergence tolerance is set at 10 -5 (tol = 10 -5 ).

In Fig. 14(b) we provide a guess design for the inclusion, with f inc = 5%, corresponding to a simple parallelepipedic domain which will serve as a comparison with the SIMP optimized topology solution. Fig. 15 shows the evolution histories of inclusion topologies and their final crack patterns using SIMP. For 3D visualization purpose, only values of ρ e ≥ ρ thr are plotted in constant blue color (ρ thr is a threshold value, defined to make sure elements with volume fraction around f inc appear). Here the SIMP solution converges in roughly 10 iterations and strictly in 77 iterations. As a comparison, BESO with initial homogeneous design takes 150 iterations to reach f inc = 0.05, showing the lower computational cost and better convergence of SIMP in this case. The final BESO topologies are not shown here to avoid too many figures. Different views of the optimized design using SIMP are shown in Fig. 16. Detailed propagation of the crack phase field corresponding to its load-displacement curve for the optimized design is shown in Fig. 17. It can be observed that the crack first initiates from the surface of the inclusion phase and then interacts with the preexisting crack. Next, it propagates along the surface of the inclusion phase until crossing the whole domain and leading to the failure of the structure. Fig. 18 shows comparison of load-displacement curves and final crack patterns for guess and optimized design. Here, the fracture resistance of the optimized design is 20% higher as compared with the guess design.

3D reinforced sample with two pre-existing crack surfaces and a parallelepipedic cavity

This final example demonstrates the potential of the method in complex 3D structural problems. The geometry and boundary conditions for this example are shown in Fig. 19. The dimensions of the 3D sample are 50 × 50 × 60 mm 3 . The domain is discretized into 100 × 100 × 120 eight-node cubic elements. As can be observed from Fig. 19, there are two pre-existing crack notch surfaces and a parallelepipedic cavity. The first pre-existing crack notch whose geometry is 26 × 10 mm 2 is embedded in the left surface of the whole sample, the second pre-existing crack notch whose geometry is 16 × 8 mm 2 is embedded in the left surface of the parallelepipedic cavity, and the parallelepipedic cavity whose geometry is 20 × 20 × 7.5 mm 3 is created by removing the elements at the position of the parallelepipedic cavity. Similar to the 3D case in section 5.2, the pre-existing cracks are simulated by prescribing Dirichlet conditions on the crack phase field with d = 1 along the crack surfaces. On the upper and lower ends of the sample, vertical incremental displacement loads with ∆U = 0.005 mm are prescribed with for the first four load increments and ∆U = 0.002 mm for the following load increments. Additionally, the optimality convergence tolerance is set at 10 -3 (tol = 10 -3 ). Fig. 20 shows the evolution histories of inclusion topologies and their final crack patterns. Here f inc is set to 5%. As can be observed from Fig. 20, the inclusion material tends to accumulate around the pre-existing crack notch surfaces and the parallelepipedic cavity to prevent crack propagation and interaction in the sample. From Fig. 20, we can note that the SIMP solution is converged in roughly 12 iterations and strictly in 23 iterations. As a comparison, BESO with initial homogeneous design takes about 150 iterations to reach f inc = 0.05, showing again much less computational costs and better convergence of SIMP in this case. Here again, BESO final designs are not shown for avoiding too many figures. For better visualization, different views of the optimized design are shown in Fig. 21. Detailed propagation of the crack phase field corresponding to its load-displacement curve for the optimal design is shown in Fig. 22. It can be observed that the cracks nucleate firstly at the two pre-existing crack notch surfaces and the corners of the cuboidal cavity, and then propagate and interact following the surface of the inclusion phase until reaching the fully broken state of the sample.

Computational times

Finally, a summary of the computational times for the different examples is reported in Table 3. In this work, a workstation with 4 cores, 16 Go Ram and 3.00 GHz processor was used for all 2D cases. For all 3D cases, a workstation with 24 cores, 768 Go Ram and 2.70 GHz processor was used. The present code has been implemented in Matlab .

Conclusion

In this work, we proposed a framework employing SIMP topology optimization and phase field method to fracture to maximize the fracture resistance of composites (two-phase materials) structures. The method allows taking into account the whole fracture process, from initiation to complete failure of the specimen. The continuous density representation of density obtained by the SIMP method allows a good convergence of the scheme and to improve the fracture resistance of a structure embedding a reinforcement phase (inclusion) for a fixed volume fraction.

We have observed that even though SIMP has slightly better convergence properties with less oscillations than BESO, the convergence of both SIMP and BESO is comparable when the initial volume fraction is set by an appropriate initial design in BESO. The final results are also very close. However, it is not always easy to define such initial design, especially in 3D. When using a homogeneous design with BESO, the convergence can be much slower (up to 6 times in some examples) than SIMP, as the actual volume fraction will be changed at each iteration by removing/adding material in each element. Then, starting with an initial homogeneous design is straightforward and leads to faster convergence with the present SIMP framework than with BESO.

This has been illustrated by conducting 3D applications in complex configurations, where defining an initial design with a given volume fraction may be non-trivial. Then, the present SIMP-phase field framework is a good candidate for reducing the computational times in designing materials and structures with enhanced fracture energy. An extension to interfacial damage interacting with bulk fracture has been proposed in [START_REF] Da | Topology optimization for maximizing the fracture resistance of periodic quasi-brittle composites structures[END_REF][START_REF] Da | Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage[END_REF] and could be straightforwardly applied to the present SIMP framework. Another exciting perspective for this work would be to produce such bi-materials with optimized geometries and validate experimentally the numerical predictions.

Figure 1 :

 1 Figure 1: Phase field approximation of a sharp crack discontinuity. (a) A sharp crack surface Γ embedded into the solid Ω. (b) The regularized representation of the crack by the phase field d (x).

Figure 2 :

 2 Figure 2: Young's modulus interpolations for SIMP and BESO.

ρ

  e v e /( Ne e=1 v e ) = f inc (40) 0 ≤ ρ e ≤ 1, e = 1, 2, ..., N e

  1, . . . , N e ,

Algorithm

Figure 3 :

 3 Figure 3: Plate with one pre-existing crack notch subjected to incremental traction load: geometry and boundary conditions.

Figure 4 :

 4 Figure 4: Evolution of inclusion topologies and associated final crack patterns (f inc = 5%).

Figure 5 :

 5 Figure 5: Load-displacement curve and crack propagation for the optimized design: (a) U = 0 mm; (b) U = 0.023 mm; (c) U = 0.034 mm; (d) U = 0.038 mm; (e) U = 0.041 mm.

Figure 6 :Figure 7 :

 67 Figure 6: Evolution of inclusion topologies and associated final crack patterns (f inc = 2%).

Figure 8 :

 8 Figure 8: Comparison of optimized inclusion designs and corresponding final crack pattern at the failure load: (a) f inc = 2%, (b) f inc = 5% and (c) f inc = 10%.

Figure 9 :

 9 Figure 9: Comparison of optimized inclusion designs and corresponding final crack pattern for three different refined meshes: (a) 120 × 60, (b) 180 × 90 and (c) 240 × 120 (f inc = 5%).

Figure 10 :

 10 Figure 10: Load-displacement curves for different refined meshes (f inc = 5%).

Figure 11 :

 11 Figure11: Influence of initial design on the optimized inclusion topologies using BESO method (f inc = 5%). Crack and inclusions are denoted by red and blue color, respectively.

Figure 13 :

 13 Figure 13: Load-displacement curves for initial and optimized designs (f inc = 5%).

Figure 14 :

 14 Figure 14: 3D sample with one pre-existing crack subjected to uniaxial tension: (a) geometry and boundary conditions; (b) guess design (f inc = 5%) using SIMP. Crack and inclusions are denoted by red and blue color, respectively.

Figure 15 :

 15 Figure 15: 3D sample with one pre-existing crack: evolution of inclusion topologies and associated final crack patterns (f inc = 5%). Crack and inclusions are denoted by red and blue color, respectively.

  (a), is considered. The dimensions of the 3D sample are 20 × 20 × 60 mm 3 . The domain is discretized into 40 × 40 × 120 eight-node cubic elements. The dimensions of the pre-existing crack notch surface are 10 × 10 mm 2 .

Figure 16 :

 16 Figure 16: Different views of the optimized design of inclusion for the 3D sample with one pre-existing crack obtained by the SIMP formulation.

Figure 17 :

 17 Figure 17: Load-displacement curve and crack evolution for the optimized design in the 3D sample with one preexisting crack obtained by the SIMP formulation: (a) U = 0.03 mm; (b) U = 0.032 mm; (c) U = 0.034 mm; (d) U = 0.038 mm. Crack and inclusions are denoted by red and blue color, respectively.

Figure 18 :Figure 19 :

 1819 Figure 18: Load-displacement curves and final crack pattern for guess and optimized design in the 3D structure with one pre-existing crack obtained by the SIMP formulation. On the right figure, crack and inclusions are denoted by red and blue color, respectively.

Figure 20 :Figure 21 :

 2021 Figure 20: 3D sample with two pre-existing cracks: convergence of the topology optimization process, evolution of inclusion topologies and associated final crack patterns (f inc = 5%) obtained by the SIMP formulation. Crack and inclusions are denoted by red and blue color, respectively.

Figure 22 :

 22 Figure 22: Load-displacement curves and crack pattern evolution for optimized design in the 3D structure with two pre-existing cracks obtained by the SIMP formulation: (a) U = 0.03 mm; (b) U = 0.032 mm; (c) U = 0.034 mm; (d) U = 0.036 mm. Crack and inclusions are denoted by red and blue color, respectively.

2 :

 2 Overall topology optimization algorithm Initialize inclusion densities ρ 1 . %% Topology optimization iteration %% Set k = 1, Err = 1 and iterative tolerance tol while Err > tol do %% Regularized fracture problem and sensitivity calculations %% Compute regularized fracture problem, obtain sensitivity vector α k and total mechanical work J k with ρ

k (x) from Algorithm 1 %% Filtering scheme %% Obtain the smoothed α k with (57) %% OC method %% Update the design variable ρ k+1 with smoothed α k %% Compute convergence %% if k ≥ 10 then

Table 2 :

 2 Material parameters used in the numerical simulations[START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF]

	Name	Notations Inclusion	Matrix
	Young's modulus	E inc , E mat	52 GPa 10.4 GPa
	Poisson's ratio	ν	0.3	0.3
	Critical fracture stress σ inc c , σ mat c	0.03 GPa 0.01 GPa

  Figure 12: Convergence of the topology optimization scheme for BESO and SIMP schemes (f inc = 5%).

		34				1
							initial design A, BESO initial design B, BESO	0.9
		30				initial design C, BESO SIMP	0.8
	External work (mJ)	18 22 26				homogeneous design, BESO volume fraction with homogeneous design, BESO	0.3 0.4 0.5 0.6 0.7	Volume fraction
		14			correspoding optimized designs	0.2
							0.1
		10	0	50	100 150 200 250 300 350 400 450 500 0
						Design iteration number
					initial design A, BESO
			0.4	optimized design A, BESO initial design B, BESO
					optimized design B, BESO
					initial design C, BESO
					optimized design C, BESO
			0.3	optimized design, SIMP
					optimized design with
	Load [kN]	0.2	homogeneous design, BESO
			0.1
				0	0	0.01	0.02	0.03	0.04
						Displacement [mm]

Table 3 :

 3 Computational times for the different examples

	Problem	No. elements No. design	Average CPU	Total sim-
			iterations	time (s)	ulation
					time (h)
	2D reinforced plate (f inc = 5%) 7200	244	33.6	2.28
	3D reinforced sample in 5.2	192, 000	77	924	19.77
	3D reinforced sample in 5.3	1, 176, 000	23	5055	32.3

Acknowledgements

This work has benefited from the project "Multiscale Modeling and Experimental Investigation of damage in composite components obtained by additive manufacturing process" (MMELED), ANR-16-CE08-0044-04. Yi Wu gratefully acknowledged the support from China Scholarship Council (CSC No. 201906130024).

Input: Densities ρ k Output: the total mechanical work J k and sensitivity vector α k Initialize u 0 = 0, H 0 = 0, J 0 = 0 and α 0 = 0. Loop over load increments n load for n = 1, . . . , n load do %% Crack phase field problem %% Given u n-1 and H

In order to remove instabilities such as checkerboard patterns and to avoid mesh-dependency in topology optimization process, element sensitivities are smoothed by means of a filtering scheme [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF]

in which w ej is a linear weight factor

Appendix: Finite Element discretization

The weak form of the mechanical problem can be found in [START_REF] Da | Topology optimization for maximizing the fracture resistance of periodic quasi-brittle composites structures[END_REF]. Using [START_REF] Kai | Failure mitigation in optimal topology design using a coupled nonlinear continuum damage model[END_REF], we can re-write the associated weak form for the damage problem [START_REF] Da | Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage[END_REF] as: find

and δd(x) ∈ S 0 d , S 0 d = δd|δd(x) = 0 on ∂Ω, d ∈ H 1 (Ω) . In this work, we adopt the same finite element discretization for the approximation of the crack phase field d and the displacement field u. We can express the discretization of the phase field problem as:

where N d and B d are matrices of damage shape function and of damage shape function derivatives, respectively, and d e denote nodal damage in one element. The discretization of damage problem (59) results into the following discrete system of equations:

in which

and

where H is given in [START_REF] Frédéric | Optimization of multiphase structures considering damage[END_REF]. The displacement problem can be discretized as:

where u e , N u and B u are nodal displacement components in one element, a matrix of displacement shape function and a matrix of displacement shape function derivatives, respectively. Using the weak form [START_REF] Da | Topology optimization for maximizing the fracture resistance of periodic quasi-brittle composites structures[END_REF], we obtain the following discrete system of equations:

with the force vector

and the stiffness matrix

where [σ] and [ε] are the vector forms corresponding to the second order tensors of stress σ and strain ε. R ± and P ± are two operators for the decomposition of strain into the tensile and compressive parts (see [START_REF] Tung Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF] for more details) and P ± are the matrix forms corresponding to the fourth order projection tensor P ± = ∂ε ± ∂ε , which can be found in [START_REF] Ambati | Phase-field Modeling and Computations of Brittle and Ductile Fracture for Solids and Shells[END_REF][START_REF] Miehe | Algorithms for computation of stresses and elasticity moduli in terms of seth-hill's family of generalized strain tensors[END_REF].

Conflict of interest statement

On behalf of all authors, the corresponding author states that there is no conflict of interest.


Faten Sahel

Pascal Guilbault

Farouk Vallette

Sylvain Feruglio

A crosstalk modelling method between a power supply and a nearby signal in high-density interconnection PCBs

Keywords: signal and power integrity analysis, simulation flow, power supply induced crosstalk

In this paper, we suggest a simulation method that allows the early detection of crosstalk noise, caused by power supplies on nearby sensitive signals, in high density mixedsignals boards. This paper focuses on the PCB and power supply modelling and discusses the accuracy of ground planes modelling in S-parameters modules, as current return paths in ground planes can be a major means of coupling. Our approach is based on, and validated with a comparative study between simulation results and measurements on a real case study board developed by Atos for HPC applications.

I. INTRODUCTION

To meet the expectations of the ever-evolving electronic technology, today's electronic systems are expected to be faster, smaller, and more efficient than ever. Consequently, Printed Circuit Boards (PCB) are increasingly multi-layered, denser, and cluttered with miniaturized components. They contain mixed-signals circuits, some of which operate at high currents, and links of various data rates. This encumbrance of the PCB results in a closeness between some disruptive elements and sensitive ones, which undergo some disturbances, such as crosstalk. The crosstalk phenomenon occurs when the signal propagating in a channel, called 'aggressor', unintentionally causes disturbances on another neighbouring channel, called 'victim'. The disturbances, commonly known as noise, can lead to errors in the transmitted data, and in the worst case, may result in a malfunctioning circuit.

Among the disruptive modules on the board, electromagnetic interferences from power supplies are the most critical; they not only hinge on the quality of the power supplies themselves, depending on their characteristics, their associated filtering and decoupling circuits, but also on their operating mode and the power they provide [START_REF]Switch-mode power supply[END_REF].

Power Distribution Network (PDN) provides the circuit with stable delivery and return paths of the supply voltage [START_REF] Shringarpure | Printed circuit board power distribution network modeling, analysis and design, and, statistical crosstalk analysis for high speed digital links[END_REF]. Robust PDNs and strict routing rules normally minimize the possible perturbations caused by the switching of power supplies on sensitive signals. But with High-Density Interconnect (HDI) technology, it becomes harder to maintain the proper spacing between sensitive signals and disruptive ones [START_REF] Holden | Eight key HDI design principles[END_REF]. And as strict as the routing rules can be [START_REF] Zhang | High Speed Digital Design[END_REF], it is unlikely for PCB designers to predict and prevent signal coupling every time, what with the parasitic effects induced by the switching transistors and those caused by the current return paths in ground planes [START_REF] Kocubincki | Mid-frequency noise coupling between DC-DC converters and high-speed signals[END_REF][START_REF] Lin | The Effect of Return Path Current on Analog Sensor Readings in Mixed-Signal Embedded Applications[END_REF].

Therefore, virtual prototyping is commonly used to verify the signals' behaviour in transmission lines, detect the parasitic effects on sensitive signals due to crosstalk, and thus allowing necessary adjustments to be made to reduce the noise, prior to first prototypes built.

There are several EDA tools that achieve virtual prototyping, using different solvers to model PCBs and the behaviour of the signals, though the accuracy of these tools depends on the modelling mode of their solvers. Finite Element Method (FEM) 3D solvers are the most appropriate for a thorough modelling, as they analyse every aspect and variation of the system to model (lines, power and ground planes, vias, packages and materials). They also include the current in the vias since the electromagnetic waves propagations are considered on the X, Y and Z axes, and thus provide accurate results. However, FEM 3D solvers are memory resourceshungry and require very long simulation times. They are ideal for analysing small systems but would be too consequential if used for analysing crosstalk in ultra-dense PCBs. 2.5D solvers would be more fitting though, surface and density-wise, but their modelling mode is simpler and less precise than FEM 3D solvers, making the ensuing results less accurate [START_REF] Camposano | Moore meets maxwell[END_REF].

Consequently, modelling complex and ultra-dense boards proves to be a challenging matter, as they are not easy to accommodate to present solver methods. The boards that Atos develop for the BullSequana ™ Supercomputers range are destined for High Performance Computing (HPC) applications : they can reach 22 layers and have usually 2 processors functioning at very high currents (>200A) and low voltages, powered by fast switching power supplies, which creates coupling noise on nearby sensitive signals. We therefore aim to spot the coupling effects stirred by the power supplies on the neighbouring signals before the board manufacturing, rather than during the prototype testing. As there isn't currently a reasonable way to do so, our objective is to implement a simulation method, using existing EDA tools, that would allow us to visualize the disturbances that might occur on the sensitive signals near the power supplies.

In the remainder of this paper, we suggest a method consisting of a simulation flow using a 2.5D solver tool; ADS SI/PI pro ™ form Keysight Technologies Inc. ® , applied to a board developed by Atos, that enables the detection of a coupling phenomenon between a power supply and a nearby sensitive signal. Firstly, we introduce our case study, then we detail the procedure that led to our established simulation flow, and finally we compare the obtained results to the measurements of the case study.

II. CASE STUDY

The printed circuit board we choose to develop and validate our methodology is a small mezzanine board called CZBS (Compute Mezzanine BXI Single NIC), used for interconnections between HPC boards. the Bull eXascale Interconnect (BXI) network is a performance booster for communications in HPC applications; it manages the communications independently from the CPUs, leaving them to be fully dedicated to computational tasks [START_REF] Derradji | The BXI Interconnect Architecture[END_REF]. The BXI Network Interface Connector (NIC) is mounted on the mezzanine board, and interconnects with the mother board through a PCI Express interface.

The CZBS board (Fig. 1) is a 118mm x 130mm x 2mm 12-layer PCB (Fig. 2), made of 2 external layers, 4 internal routing layers, and 6 power and ground planes. It contains an Application-Specific Integrated Circuit (ASIC), mainly powered with a voltage of 0.9V, delivered by a buck converter from a 12V input voltage. The power supply, placed on the top layer, causes disturbances on some neighbouring signals, especially on the enable signal of a clock buffer on the first routing layer (L1), right under the power supply. As we can see in Fig. 3, this signal spreads under the transistors, the inductor, and near the switching node of the buck converter. For clarity purposes, only the elements involved in our case study appear in Fig. 3, the other components, vias and layers are not shown.

The choice of this board as a case study is motivated by the presence of crosstalk between the power supply and the clock enable signal, while being small enough to allow a reasonably fast modelling in an EDA tool. The aim of our study is to be able to simulate the coupling effect and observe the disturbances on the victim signal, through a refined modelling of the studied system, and correlate to measurements on the actual board.

III. SIMULATION PROCESS

To achieve our goal, we aspire to reproduce the interaction between the power supply and the clock enable signal through simulation, using as precise a modelling of the case study as possible.

A. Simulation overview

To ensure the accuracy of our depiction, we need to consider all the elements involved in the coupling phenomenon, meaning:

• The PCB layout, which contains the aggressor and victim signals, and whose properties greatly affect the quality of the signal transmission. The PCB is modelled using the ADS Signal Integrity (SI) simulation tool to extract the layout of the board, and therefore capture signals and power planes' behaviours in an S-parameters model.

• The power supply, which generates the perturbations on the victim signal, is emulated by its SPICE model.

Both models (PCB extraction and power supply models) are to be merged in a transient analysis.

B. Basic simulation method

Regarding the PCB modelling part, we run an SI analysis consisting of components' ports (the MOSFETs, inductor, decoupling capacitors, and load of the considered buck regulator, as illustrated in Fig. 4), and nets involved in our case study (the required signals for the power supply functioning, and the clock enable signal). In ADS SI/PI pro ™ , the nets are sorted into categories (ground, power, or signal), depending on their functions. We use a logarithmic analysis, with a 100Hz -10GHz frequency range, and 100 points/decade. The extractor takes about a day to simulate our design in a 64GB machine and supplies the results in an S-parameters format, containing 48 ports.

Concerning the buck converter modelling, we use the SPICE electrical models of the driver and MOSFETs provided by their manufacturer. For more accurate results, we include a lumped model of the PDN, obtained from the above-mentioned SI simulation, in the SPICE analysis of the power supply. The PDN parasitic effects being included, their impact on the converter performance are thus considered in generating realistic gate command signals.

The transient simulation we launch includes the results of both the layout SI modelling and the power supply SPICE analysis. It contains the S-parameters model of the layout, the components involved in the case study, and the gates' drive signals of the MOSFETs. Fig. 5 displays the simulation flow in a detailed diagram.

C. Limits of the standard approach

The obtained results using the SI tool default setup relate a very optimistic configuration, as the disturbances are barely apparent. We infer that the ground planes might be considered too ideal, and that the current return paths might not be modelled accurately enough. They indeed play a major role in our case study, as the power supply is isolated from the victim signal by a ground plane, which if deemed ideal, over-shields the two elements from the coupling between them. We therefore need to integrate the real behaviour of the ground planes to the PCB modelling, to get more accurate results.

The S-parameters block representing the PCB contains the ports of the components used in the modelling, and a reference node. Each port has a positive pin and a negative one, representing the positive and negative terminals of the components. However, as indicated in Fig. 6-a, the negative pins are not accessible to the user [START_REF] Brim | S-parameter Modeling and Simulation for Signal Intergrity Analysis[END_REF]. The measurable voltages at the positive pins relate to (Vpn+-Vref). The reference node is a mathematical node chosen amidst the ground planes, which form the reference planes in the analysis, and is connected to the circuit ground (node 0) [START_REF] Chitwood | Issues with Interfacing '2N' and 'N+ref' Behavioral Models[END_REF]. We assume that the lack of crosstalk noise on our victim signal is due to the ground planes being assigned as reference, and hence indirectly connected to a perfect ground (node 0) through the reference node.

We therefore seek to verify the signals' behaviours at each components' terminal. To be able to observe them, we need an S-parameters structure that has a port for each terminal, as represented in Fig. 6-b.

D. Suggested solution

To do so, we proceed by creating an additional virtual layer, above the bottom layer, on which we place a ground plane. We consider this layer virtual because it doesn't exist on the board, it is added solely for analysis purposes, since it is mandatory for the analysis to use a ground net as reference.

We choose to place it above the bottom layer, so it doesn't disturb the studied system in any way, for it is as far away from the case study as possible. Fig. 7 shows the cross section of the modified layout, with the added virtual layer, between the sixth power plane and the bottom layer.

We describe this plane as the only GND plane for the analysis, thus the only one being connected to the ground pin reference of the S-parameters block. Consequently, since it is not directly connected to the components, they aren't affected by the supposed ideality of the reference. The real ground planes of our board are simply set as power planes, and thereby conserve their characteristics.

Through this procedure, the actual GND planes are more accurately modelled (since they are not connected to the reference ground), and the S-parameters block has now a port for each positive and negative pins of each component, allowing us to visualize the actual signals circulating in the circuit, by measuring the difference between the ports matching the positive and negative terminals of the components (Vpn+-Vpn-). The new Sparameters structure contains 88 ports, since 40 ports out of the 48 ones from the first structure have now ports representing their positive and negative pins (namely the decoupling capacitors, the load, the power sources and the low MOSFET). The 8 other ports correspond to components' pins that are not connected to the ground in our board.

IV. RESULTS

Fig. 8 displays the evolution of the aggressor and victim signals of our case study. The aggressor (red), is the converter's switch signal, which presents a typical behaviour of a switch node, with a 7.5% duty cycle, and a 440kHz switching frequency. The victim (blue) is the clock enable signal, which should ideally be a 3.3V static signal. We can clearly see that the noise on the victim signal and the switching of the aggressor are synchronous.

The reference results in Fig. 8-a are obtained by measurement with a 1GHz bandwidth oscilloscope, which ensures the crosstalk noise isn't filtered. Fig. 8-b andc-are the simulation results from a transient analysis and represent the standard and revised simulation methods respectively.

As shown in Fig. 8-b, the standard simulation does not cover the crosstalk modelling between the power supply and the victim signal, considering that we can hardly distinguish any noise on it at all. On the other hand, Fig. 8-c shows that with our approach, the crosstalk disturbances are clearly visible, and match the measurements. A first order comparison indicates that the obtained disturbances indeed resemble the measured ones in shape, and in amplitude to a lesser extent.

Out of the possible reasons for this difference in amplitude, the aggressor's behaviour stands out. The switch node seems indeed to be stronger in simulation than in measurements, which could explain why the noise amplitude on the victim signal in the simulation is greater than expected. However, our simulation method remains efficient, since we are now able to spot the crosstalk noise on the victim signal, in a dense PCB, using a 2.5D solver.

The crosstalk noise level on our case study seems a bit excessive, we therefore modify the board by adding more decoupling capacitors and a gate resistor to the power supply system, in order to reduce the coupling. We apply our simulation process to the altered design layout and compare the results to the measurements on the revised board, as shown in Fig. 9. The aim of this step is not to correct the board per se, since the noise detected on the victim signal is not relevant to the functioning of the system, but to further corroborate the consistency of our simulation method.

Measurements in Fig. 9-a reveal that the coupling is considerably reduced with the revised board, as we can barely distinguish the disturbances caused by the switching node from the random noise. Fig. 9-b shows that similarly, the noise level in the simulation, is reduced in equal proportions.

The results from the revised board therefore validate the accuracy of our simulation flow and confirm the significance of its application in virtual prototyping. In our case study, our simulation method allows us to ensure that the modifications made to the board are indeed efficient, in reducing the coupling on the victim signal. A simulation method is proposed to improve the modelling of the coupling between a power supply and a neighbouring signal in a multi-layer PCB, and depict the crosstalk on the victim signal, by ensuring the current circulation in the ground planes is accurately represented.

A basic analysis doesn't allow the visualization of the coupling noise through ground planes on the victim signal because the extraction tool generates an S-parameter block whose reference is chosen amidst the ground planes assigned to the analysis. Consequently, the GND planes' behaviour described in the S-parameters model is not realistic and imparts an optimistic performance of the case study.

Our method involves adding a virtual ground plane to the design, set it as reference in the modelling tool, and consider the actual ground planes of the board as power planes, thus allowing the crosstalk noise on the victim signal to be fully observed. This simulation process is verified and validated by comparison between the simulation results and the measurements of a case study board bearing a coupling issue, and its revised version.

Our simulation method being validated on a small and rather simple case study, we aspire now to generalize it so we can apply it to all types of boards that Atos develop, which can be much more intricate than the CZBS. We are currently applying it to a larger, denser, and more complex compute board: a 532mm x 184mm x 2.38mm 18-layer PCB, containing 2 CPUs which require very high currents and very fast low voltage signals. The SVID (Serial Voltage Identification Debug) signals on this board present some coupling noise due to their nearness to the VRMs (Voltage Regulator Model) that power the CPUs. These properties make this design a perfect case for testing our simulation method and evaluating the resources it requires, in complex and ultra-dense boards.

Fig. 2 .

 2 Fig. 2. Cross section of the case study

Fig. 1 .

 1 Fig. 1. Top view of the case study

Fig. 4 .Fig. 5 .

 45 Fig. 4. Buck power supply circuit

Fig. 6 .Fig. 7 .

 67 Fig. 6. a-Basic S-parameters structure b-Required S-parameters structure

Fig. 9 .Fig. 8 .

 98 Fig. 9. Time-domain waveforms of aggressor and victim signals in revised case study board


Bert Hardy's Chinatown

Gregory Lee

• The photographs of Chinatown by Bert Hardy were taken in May 1942 during the Second World War. • They were commissioned by the illustrated news magazine Picture Post but were never published.

The collection was subsequently bought by the Getty Foundation. The custody of the photographs has now been entrusted to the Open Eye Gallery in Liverpool. • This lecture was delivered to a mixed audience of academics and members of the local community at the invitatiion of the Open Eye Gallery.

• The men we see are just a few of the 10,000 Chinese seamen who manned the Merchant Navy ships, for Britain during the war, mainly Blue Funnel, Ben Line, and Shell. Chinese citizens also served as cooks and laundrymen for the Royal Navy during the war and even right up until the 1990s. The pictures show them in lodging houses and makeshift clubs while they wait for the next ship.

• This was not the first time the Chinese served the British war effort. The First World War trenches were dug not by British Tommies but by 20,000 or so Chinese men who were enrolled in the British Army Chinese Labour Corps. The Chinese government had supplied these men hoping that in return at the conclusion of the war the German colony of Qingdao (Tsingtao 青岛 ) would be returned to Chinese soveregnty. But the strategy did not work. The Versailles Peace Conference in 1919, decided to transfer the former German territories to Japan, also an ally of France and Britain in the First World War.

Tombstones of three members of First World War British Army Chinese Labour Corps "repatriated" from France/Belgium and buried in the Chinese sector of Liverpool's Anfield cemetery Many of the the Chinese Labour Corps troops were from Shandong (Shantung) 山东 Province one of the objectives of the Chinese authorities in supplying men to serve on the European war front was to recuperate the German colony of Qingdao in Shandong. But the Versailles Peace Treaty transferred the German territorial rights to Japan.

• While many of those serving in the Labour Corps during the Great War came from the northern province of Shandong, during the Second World War the Chinese sailors hailed mostly from the seaport of Shanghai. • The Chinese seamen had no choice to but to serve on the transatlantic supply routes. They were stranded, they could not go back to China where the war was raging against the Japanese occupying forces. The British and American governments even passed laws obliging the Chinese seamen to serve in the merchant navies. They were paid far less than British sailors. And at one point the Shanghai seamen protested and took industrial action.. This act of defiance was not forgotten by the British authorities and the Blue Funnel shipping company.

• After the war, as was revealed just about five years ago, thousands of Chinese men who had settled in Liverpool legally, many of them married to Liverpool women and who had fathered children, were deported. The Home Office and the Holt shipping company had not forgotten the demands for fair pay made by these seamen during the war. • Their landing papers were revoked, they were rounded up at night by police and Home Office officials, and put into cargo ships roughly converted with bunk beds and sent back to a China in turmoil and to an uncertain fate. • Their Liverpool properties were appropriated by the city council, and their families maintained a silence for over half a century.

• Thanks to the efforts of Yvonne Foley the daughter of one of these sailorfathers there is a now a plaque at the Pierhead, Liverpool's landmark waterfront, acknowledging the contribution that was made by the men and the wrong that was done to them.

navvy during the two world wars and who were subsequently deported from the UK in 19461947

Quite possibly some of the men we see in Hardy's photographs fell victim to the furtive postwar deportations.

Liverpool Chinatown, June 1993, before the construction of the arch or pailou 牌楼

• From the 1960s onwards the nature of Chinatown changed.

The new immigrants, who make up the majority in Liverpool today, came from Hong Kong and the Canton region. Whole families came and the former "mixed marriage", bicultural Chinatown gave way to a new kind of Chinatown community.

A Chinatown Museum?