Nian Zhang

Yujiao Fan

Hui Chen

Sylvain Trépout

Annie Brûlet

Min-Hui Li

Polymersomes with smectic liquid crystal structure and AIE fluorescence

A series of amphiphilic diblock copolymers PEG-TPE-b-PAChol containing a cholesterol-based smectic liquid crystal block PAChol and a tetraphenylethene AIEgen situated between hydrophilic PEG and hydrophobic PAChol blocks were prepared by RAFT polymerization. The liquid crystal properties of the block copolymers were first characterized. One of the copolymers, PEG45-TPE-b-PAChol16, with a hydrophilic weight ratio of 18% was chosen for the self-assembly study using the nanoprecipitation method in dioxane/water and THF/water systems. The co-solvent (THF or dioxane) had significant effect on its self-assembly behaviors. The characterization by cryo-EM showed that elongated ellipsoidal or tubular polymersomes with smectic ordering in the membrane were formed in the dioxane/water system. In contrast, spherical solid nanoparticles were obtained in THF/water system as observed by SEM. All polymersomes and nanoparticles showed AIE fluorescence. The present work shows the possibility to obtain elongated polymersomes with both smectic LC orders and AIE fluorescence. By replacing the prototypical TPE with AIEgens that can be excited in longer wavelength than UV, these cholesterol-based polymersomes and nanoparticles can be designed as fluorescent carriers of drugs for possible theranostic applications.

Introduction

Liquid crystal (LC) polymersomes are a special type of polymer vesicles with thermotropic LC phases displayed in the lyotropic polymer-based bilayer membrane. 1 LC polymersomes exhibit either interesting non-spherical morphology because of the LC phase like smectic one, 2 or stimuli-responsive features due to the sensitivity of chain conformation to external stimuli such as light. 3 LC polymersomes are usually prepared by the selfassembly of amphiphilic block copolymers where the hydrophobic block is a thermotropic nematic or smectic LC polymer. Our group and others have reported several LC polymersomes based on side-chain LC polymers. [4][5][6][7][8][9][10][11][12] The selfassembly of amphiphilic LC block copolymers is affected not only by the relative sizes of two constituting blocks, but also by the LC order and the chain conformation in the LC block, which results in more complicated self-assembly behaviours than those of classical coil-coil block copolymers. 2 Spherical, 4, 5 ellipsoidal 7 and faceted 6 LC polymersomes have been reported, which disclosed the important relationship between the orientational and positional LC order and the curved geometry of a two-dimensional membrane. 2 Among them, polymersomes made from poly(ethylene glycol)-block-poly(cholesteryl-based acrylate) (PEG-b-PAChol) 7 shown in Figure 1a is of special interest for the following reasons. Firstly, the cholesterol is a biosourced and biocompatible LC moiety; it has been exploited as building block of polymer micelles for the delivery of with anticancer drugs due to its favourable hydrophobic interaction with drugs. 13 Secondly, PAChol exhibits a smectic LC phase extending up to at least 100℃ and a glass transition higher than room temperature (precise values depending on the molecular weight). 7 Thirdly, PEG-b-PAChol forms ellipsoidal or elongated polymersomes in water, where the smectic phase or glassy smectic phase play the role of physical cross-linking and ensure the stability and robustness. 7,9 The elongated vesicular shape may be an advantage when they are used as carriers of imaging agents or drugs, because it has been shown that the elongated shape of the carrier does affect the pharmacokinetics and in vivo distribution, 14 for example, by persisting longer after intravenous injection and possessing higher drug loading than their spherical counterparts. 15,16 Therefore, it is interesting to investigate this cholesterol-based LC polymersome for the bioimaging and controlled release. In our previous work, a disulfide bridge (-S-S-) was introduced between the hydrophilic PEG block and the hydrophobic PAChol block (Figure 1b) to get the reduction-responsive LC polymer PEG-SS-b-PAChol. A robust reduction-responsive polymersome based on PEG-SS-b-PAChol was achieved, which has been used to encapsulate the calcein as a model molecule. The release of the calcein from the polymersome was triggered by glutathione (GSH) both in PBS suspension and in vitro inside the macrophage cells with artificially enhanced GSH concentration (≥35 mM). 17 Moreover, cell viability tests have shown that nanoparticles formed by both cholesterol-based polymers (PEG-SS-b-PAChol and PEG-b-PAChol) were not cytotoxic. 17 In this work, we want to explore the cholesterol-based LC polymersome as imaging nanoparticle by introducing a fluorescent moiety (tetraphenylethene, TPE) between PEG and PAChol blocks (PEG-TPE-b-PAChol in Figure 1c). Fluorescent materials, in particular organic dyes, have been widely used as bioimaging agents. When organic dyes are integrated into polymer-based nanostructures, [18][19][20] the optical signals could be amplified due to the incorporation of several fluorophore molecules in the same nanoparticle. In addition, the polymer structures, in some cases, contribute to the uniform distribution of the involved fluorophores, which could avoid aggregation-caused quenching (ACQ) or photo-bleaching to a certain extent. More excitedly, luminogens with aggregation-induced emission (AIE) characteristics, that circumvents the ubiquitous ACQ phenomenon, have been largely developed recently as efficient fluorescent visualizers for organelles imaging and drug delivery monitoring. [21][22][23][24][25][26][27][28] The combination of polymer nanostructures with AIE luminogens (AIEgens) provide then innovative approaches to cell/tissue imaging and to in vivo study of drug bio-distribution. 27,28 For this purpose, several systems of polymersomes with AIE properties have been published. [29][30][31][32][33][34][35][36][37][38][39] In this work, we introduce a typical AIEgen, tetraphenylethylene (TPE) in the polymer structure as shown in Figure 1c. The TPE group has been inserted between the hydrophilic block PEG and the hydrophobic block PAChol, in order to preserve the smectic organisation of the hydrophobic block. A series of block copolymers PEG45-TPE-b-PACholn with different hydrophilic ratios (fPEG = 8 -24 wt%) were synthesized and characterized. One of them, PEG45-TPE-b-PAChol16 (fPEG = 18 wt%) was taken as an example to study the self-assembly by nanoprecipitation in dioxane/water and THF/water systems. Elongated smectic polymersomes were obtained in dioxane/water system, while spherical nanoparticles were obtained in THF/water system. Both polymersomes and spherical nanoparticles exhibited AIE fluorescence under UV light.

Experimental Materials

Bromobenzene (C6H5Br, 99%, Acros), magnesium (Mg, >99%, Sigma-Aldrich), Iodine (I2, 99.5%, Alfa Aesar), carbon disulfide (CS2, Fisher Scientific), hydrochloric acid (HCl, 37%, VWR), maleic anhydride (C4H2O3, 98%, Alfa Aesar), carbon tetrachloride (CCl4, Sigma-Aldrich), N, N'dicyclohexylcarbodiimide (C13H22N2, DCC, Fluka), 4dimethylaminopyridine (C7H10N2, DMAP, Sigma-Aldrich), 2-Hydroxyethyl acrylate (97%, Alfa Aesar), cholesteryl chloroformate (95%, Sigma-Aldrich) and pyridine (99.8%, Sigma-Aldrich) were used as received.

2,2'azobis(isobutyronitrile) (AIBN, Sigma-Aldrich) was recrystallized from ethanol for three times before use. Poly (ethylene glycol) monomethyl ether (mPEG45-OH, Mn = 2000 Da, Fluka) was purified by precipitation in diethyl ether for three times before use. Dry dichloromethane and THF was taken from the solvent purification system (SPS). 1,4-dioxane used for polymerization was dried by refluxing with calcium hydride. The monomer AChol (cholesteryl acryloyloxy ethyl carbonate) was synthesized by a one-step reaction between 2-hydroxyethyl acrylate and cholesteryl chloroformate as detailed in Supporting Information (SI, Figure S1).

Methods

Synthesis of amphiphilic LC diblock copolymer PEG-TPE-b-

PAChol. The diblock copolymer PEG-TPE-b-PAChol was synthesized by RAFT polymerization using TPE-functionalized macro chain transfer agent (CTA), mPEG-TPE-CTA, as shown in Figure 2a. For mPEG-TPE-CTA preparation, benzodithioic acid and 6-hydroxyhexyloxy-substituted TPE (TPE-OH) were firstly synthesized as described previously. 40,41 TPE-OH was further transformed into a TPE derivative bearing an alkene group and a carboxylic acid group (TPE-CH=CH-COOH) by reacting with maleic anhydride. Then TPE-CH=CH-COOH was reacted with benzodithioic acid by thiol-ene coupling reaction to obtain a new TPE derivative bearing a benzodithioate group and carboxylic acid group. Finally, the esterification of this functional TPE derivative by mPEG45-OH in the presence of DCC and DMAP was performed to get mPEG45-TPE-CTA. The product was purified by column chromatography with methanol and dichloromethane mixture (volume ratio 1:30) as the eluents to remove unreacted mPEG45-OH. The obtained product was characterized by 1 H NMR and SEC, as shown in Figure S3-4. The Mn of the prepared mPEG-TPE-CTA was 2680 Da as calculated by 1 H NMR and its polydispersity index (Mw/Mn) was about 1.09 evaluated by SEC results. Taking the synthesis of PEG45-TPE-b-PAChol16 as an example, the typical RAFT polymerization process was as follows: AChol (226 mg, 0.429 mmol), mPEG-TPE-CTA (103 mg, 0.0429 mmol) and AIBN (1.41 mg, 0.0086 mmol) were added into a 15 mL Schlenk tube equipped with a Teflon coated stirring bar. Then 1.5 mL dry 1,4-dioxane was added. After the solution became clear, the Schlenk tube was degassed by three freeze-pumpthaw cycles and then immersed into an oil bath of 80 °C. The mixture was stirred at 80 °C for 24 h, after which the polymerization was terminated by freezing the mixture in liquid nitrogen. The mixture was then added drop-by-drop into cold methanol to precipitate the crude PEG-TPE-b-PAChol diblock copolymer, which was further purified twice using dichloromethane/methanol as solvent/precipitant pair. The dried pure polymer was obtained as a pink solid (220 mg, yield: 67%). The synthesis of other PEG-TPE-b-PAChol with different PAChol block lengths was similar except that the different 1) with THF as the eluent.

amounts of AChol, mPEG-TPE-CTA and 1,4-dioxane were used on keeping the concentration of mPEG-TPE-CTA at 0.028 M. The details are given in Table S1. The molecular weights and molecular weight distributions characterized by 1 H NMR and SEC of all copolymers are listed in Table 1.

Liquid crystalline properties of diblock copolymer PEG-TPE-b-

PAChol. The mesomorphic properties of block copolymers in pure stat were studied by thermal polarizing optical microscopy (POM), differential scanning calorimetry (DSC), and small angle X-ray scattering (SAXS) (see SI for instrumental details). The SAXS sample is a fiber drawn by hand with a pair of tweezers from a molten copolymer at a temperature just below the isotropic-liquid crystal phase. Self-assembly of PEG45-TPE-b-PAChol16. The self-assembly was performed using a classical nanoprecipitation method. Typically, the copolymer was firstly dissolved in THF or dioxane, which was a good solvent for both PEG and PAChol blocks. The initial concentration of the copolymer in organic solvent was 2 mg/mL. Then deionized water was added slowly to the THF or dioxane solution (around 3 uL/min) with gently shaking until the water content reached 75 wt%. THF or dioxane was removed by dialysis against deionized water for 3 days in a 3500 Da cut off cellulose dialysis bag. Finally, the aqueous solution of the PEG45-TPE-b-PAChol16 self-assemblies with a concentration of about 0.5 mg/mL (polymer concentration) was obtained. Studies of AIE properties. The fluorescence quantum yield of the AIE polymersomers or nanoparticles in aqueous solution were obtained by comparing the integrated fluorescence spectra of the self-assembly solution with the fluorescence spectrum of a standard fluorogen, 9,10-diphenylanthracene, in ethanol (excitation wavelength: 355 nm) as described in the literature (see SI for detail). 42 For the detection of AIE phenomena in organic solvent/water mixtures with different water fractions, a classical process was used. Typically, 25.6mg copolymer was dissolved in 1mL THF (or dioxane), and an aliquote of 131 L of the copolymer solution was added into each of the 7 glass vials. Then the mixture of water and THF (or dioxane) with different volume fraction of water (fw = 0%, 20%, 40%, 60%, 80%, 90% and 99%) was added, respectively, into each of the glass vials containing the tinny polymer solution to reach the copolymer concentration of 10 -4 mol/L (3.36 mg/mL). All mixtures in vials were shaken overnight.

Results and Discussion

Synthesis of amphiphilic LC diblock copolymer PEG-TPE-b-PAChol

Five LC amphiphilic diblock copolymers PEG45-TPE-b-PACholn with different PAChol block lengths and hydrophilic ratios (Table 1) were prepared according to processes shown in Figure 2a (see Methods and polymerization using different feeding ratios of monomer/CTA (AChol/mPEG-TPE-CTA), as shown in Table S1. The chemical structures of the copolymers were confirmed by 1 H NMR (Figure S5). The degree of polymerization (DP) and Mn of PAChol was calculated by comparing the integrated areas of the proton peak from the cholesteryl groups (Ig of peak g) and that from the methyl group at the end of PEG chain (Ia of peak a) in the 1 H NMR spectrum (see Figure S5), i.e., DPPAChol = 3*(Ig/Ia). The SEC chromatograms of the five copolymers (Figure 2b) show unimodal peaks with narrow distributions (Đ = 1.17-1.24). So, these amphiphilic LC block copolymers are suitable for the study of their self-assembly behaviours and AIE properties.

Liquid crystal properties of amphiphilic LC diblock copolymer PEG-TPE-b-PAChol

Recall that the homopolymer PAChol10 is a smectic polymer with the phase sequence: g -68°C -SmAd -156°C -I, where the interdigitated smectic A phase (SmAd) has a smectic layer spacing of P = 4.29 nm, which corresponds to a value between l and 2l, l = 2.61 nm being the wholly extended length of the Chol mesogen estimated by Dreiding models. 7 Here, three copolymer samples PEG45-TPE-b-PACholn (n = 11, 14, 16) were studied by DSC, POM and X-ray scattering to elucidate their liquid crystal phases in pure state. Take PEG45-TPE-b-PAChol16 as an example for DSC study (scanning rate ±10°C/min, see Figure S6 and Table S2 in SI). Upon first heating the pristine sample obtained from purification by precipitation, a melting point (Tm) of PEG at 45.7°C, a glass transition Tg of PAChol16 at 66.9°C and several crystals melting points of PAChol16 around 143°C and 152°C were observed until isotropic phase. On cooling from 160°C, the transition from isotropic to smectic (SmAd) was observed at 139.8°C (The assignment was made by POM and SAXS, see below). The smectic phase of PAChol16 was kept until to -50°C without recrystallization, while PEG block recrystallized at Tcr = -34.7°C. On the second heating/cooling scan, SmAd-I and I-SmAd of PAChol16, as well as Tm and Tcr of PEG were again observed. Using POM on cooling slowly (-1°C/min) from isotropic phase, the I-LC phase transition with birefringence appearance was observed at 138.0°C with the development of the typical fan-shape textures of smectic phase. Figure 3a shows textures at 128°C and no significant change of the type of textures was observed until 25°C (see Figure S7 in SI). The fiber sample of PEG45-TPE-b-PAChol16 drawn at 135°C and cooled down to 25°C was studied by SAXS (Figure 3b). Two orders of diffraction was detected at q1 = 1.447 nm -1 and q2 = 2.878 nm -1 (q2/q1 = 2), that are associated with a repeating distance of 4.34 nm of the smectic phase. Its 2D diffraction pattern (Figure 3b) also indicates the type A of the smectic phase with its layer normal aligned perpendicular to the fiber long axis. By comparing all these results to those of the homopolymer PAChol, we can conclude that the liquid crystal phase in the hydrophobic block of PEG45-TPE-b-PAChol16 is the same as that in the homopolymer, i.e., SmAd with layer spacing d of 4.3 nm. Neither the PEG block nor the TPE moiety have disturbed the smectic phase of PAChol block. Note that in a luminescent LC polymer published recently where both TPE and Chol moieties were connected in the same side group, the presence of TPE did influence the LC phases and circularly polarized luminescence properties were achieved. 43 The block copolymers PEG45-TPE-b-PACholn (n = 11, 14) showed similar phase behaviours, only transition temperatures of SmAd-I being slightly lower (TI-SmA = 128.5°C for n = 11, and TI-SmA = 132.8°C for n = 14 on cooling at 1°C/min by POM), because of the relative high ratios of the PEG part (non-LC) in the copolymers. Their POM textures and SAXS patterns are given in SI (Figure S8-9).

Self-assembly of PEG45-TPE-b-PAChol16 by nanoprecipitation in dioxane/water and in THF/water systems

Our previous studies showed cholesterol-containing amphiphilic copolymers with PEG block ratio fPEG > 15% are favourable for vesicle formation. 7,9 Therefore, PEG45-TPE-b-PAChol16 (fPEG = 18%) has been chosen here as a representative polymer to make the self-assembly. The self-assembly of PEG45-TPE-b-PAChol16 was performed by nanoprecipitation using either dioxane or THF as co-solvent. The obtained selfassemblies were characterized by DLS, cryo-EM and SEM. The cryo-EM images in Figure 4 show the morphologies of selfassemblies prepared in dioxane/water system. PEG45-TPE-b-PAChol16 formed elongated vesicles with ellipsoidal or tubular shapes (Figure 4a Another interesting observation of the LC polymersomes obtained in dioxane/water system was that the membrane of ellipsoidal or tubular vesicles presented the periodic stripes perpendicular to their long axis (Figure 4a-c). The inset in Figure 4a shows the fast Fourier transform (FFT) giving the period of the stripe P = 4.2 nm. All these structural characteristics are similar to those observed for the polymersomes of PEG45-b-PAChol16 prepared in dioxane/water system. 7,9 Note that the difference between PEG45-TPE-b-PAChol16 and PEG45-b-PAChol16 is that in the former polymer an acrylate unit with a TPE side group was inserted between PEG and PAChol blocks. As shown in Figure 4d, bending of the membrane perpendicular to the layers makes the layer spacing on both sides of the membrane unequal and therefore costs extra elastic free energy. In contrast, bending the membrane parallel to the smectic layer does not change the layer spacing, and therefore is energetically favoured. Consequently, the smectic layers (stripes) are always perpendicular to the major axis of vesicles.

In brief, polymersomes with smectic order in the membrane and ellipsoidal or tubular shape were also obtained from PEG45-TPE-b-PAChol16.

When the organic solvent in the self-assembly process was changed from dioxane to THF, spherical nanoparticles instead of elongated vesicles were obtained. Figure 5a-b show the morphologies and the size distributions of nanoparticles prepared in THF/water system in the same conditions as in dioxane/water system (the initial concentration of the copolymer in organic solvent was 2 mg/mL). Changing the initial concentration of the copolymer from 2 mg/mL to 1 mg/mL and to 10 mg/mL, respectively, the morphologies of nanoparticles kept unchanged, as shown in Figure 5c, Figure S11-12. These results are similar to those observed in PEG45-b-PAChol16 that also self-assembled into spherical aggregates in the THF/water system. 9 Please do not adjust margins

Please do not adjust margins

The reason why different morphologies of the self-assemblies were formed in THF/water system and dioxane/water system may reside in the different solubility of PAChol block in THF and dioxane. 9 The values of solubility parameter (δ/MPa 1/2 ) of solvents are δTHF = 19.4 and δdioxane = 20.5. [START_REF] Mark | Physical Properties of Polymer Handbook[END_REF] Those of repeating units of polymer blocks are δPEG = 20.2 and δPAChol = 16.1, as calculated by group contribution method. [START_REF] Brandrup | Polymer handbook[END_REF] According to these solubility parameters, THF and dioxane are nearly equivalently good solvent for PEG. However, THF is a better solvent for PAChol than dioxane. This difference in solvent-PAChol interaction was also proved by the higher critical water concentration (CWC) for aggregation formation during the nanoprecipitation in THF/water system than in dioxane/water one. 9 We extrapolated that the PEG45-TPE-b-PACholn series should display similar behaviour, and this CWC difference was effectively demonstrated by AIE properties of PEG45-TPE-b-PAChol16 in dioxane/water and THF/water mixtures that will be discussed in the next section (see Figure 7 below). As the interaction between dioxane and PAChol is not so favourable, the interaction between PAChol and PAChol (mesogenmesogen interaction) may be favoured upon water addition and play then an important role in the aggregate formation. We demonstrated previously that the smectic organization was a kinetic driving force for the smectic polymer aggregates in dioxane/water system. [START_REF] Jia | [END_REF] Small micellar nuclei with smectic order should be formed already at the very beginning of water addition in dioxane solution. A ''liquid crystallization''-driven process leads to the growth of these nuclei and the formation of membranes with smectic stripes, that finally form elongated shapes of vesicles. In contrast, in THF/water system, THF might be expelled from both the PEG and PAChol blocks much latter upon water addition. The nature of kinetical trapping in the polymer self-assembling 47,48 resulted in the formation of solid spherical nanoparticles with both PEG and PAChol blocks frozen inside their cores.

AIE properties of elongated polymersomes and spherical nanoparticles of PEG45-TPE-b-PAChol16

The AIE property of the self-assemblies formed by PEG45-TPE-b-PAChol16 copolymer were also characterized. Figure 6 shows the water dispersions of self-assemblies obtained by nanoprecipitation in dioxane/water and in THF/water systems observed under UV light (λ = 365 nm). The dispersions of polymersomes and spherical nanoparticles emit cyan fluorescence under UV light, while no fluorescence are observed in the polymer solutions in dioxane and in THF (photos not shown). In fact, PEG45-TPE-b-PAChol16 containing a TPE group between PEG and PAChol blocks is soluble in dioxane and in THF. In these solutions, the TPE is free to the intramolecular rotations that cause the non-radiative decay after UV illumination. When PEG45-TPE-b-PAChol16 is self-assembled into polymer bilayer of vesicular membrane or aggregated into spherical nanoparticles, the intramolecular rotation of TPE moieties is restricted. Consequently, under UV illumination the radiative return through fluorescence emission takes place, and the systems display a typical AIE phenomenon. The fluorescence quantum yield of polymersomes was measured as about 0.76% by an indirect method using a standard fluorogen as reported in the literature (see SI). 42 The fluorescence quantum yield of spherical nanoparticles was measured as about 0.74%. The relative low values of quantum yield may be explained by the fact that the only one TPE group in each copolymer chain is a dangling group and located in the intersection of hydrophobic and hydrophilic domains. Consequently, the degree of restriction of TPE intramolecular rotations is relatively low. To increase the quantum yield, triblock copolymers bearing several acrylate unities with TPE side group between PEG and PAChol or after PAChol block can be designed in the future. [33][34][35][36] In addition, the AIE properties of PEG45-TPE-b-PAChol16 in solvent/water mixtures with different water fractions (fw, volume) were also studied. As shown in Figure 7a, at fw = 0%, the polymer solution in dioxane is nonemissive; at fw = 20%, the fluorescence starts to appear. That means the aggregation occurs between 0 and 20 vol % and CWC ≤ 20 vol% in the dioxane/water system. The fluorescence becomes stronger with the increase of fw. The characterization by fluorescence spectroscopy revealed the same phenomena as shown in Figure 7c (λex = 349 nm). The fluorescence increase in the mixtures of THF/water (Figure 7b,d) shows the similar tendency as in the mixtures of dioxane/water, except that the AIE fluorescence appears later, i.e., between 20 and 40 vol % with 20 vol% < CWC ≤ 40 vol% for aggregation. These observations are coherent with the fact that THF is a better solvent for the hydrophobic block than dioxane as above discussed.

Conclusions

A series of LC amphiphilic PEG-TPE-b-PAChol diblock copolymers containing a smectic LC block PAChol and an AIEgen (TPE) between hydrophilic PEG and hydrophobic PAChol blocks has been prepared by RAFT polymerization of AChol monomers using a TPE-functionalized mPEG45 as the macro chain transfer agent (mPEG-TPE-CTA). The liquid crystal properties of the block copolymers were first characterized. The liquid crystal phase in the hydrophobic block of PEG45-TPE-b-PAChol16 is the same as that in the homopolymer PAChol, which is a SmAd phase with layer spacing of 4.3 nm. Neither the PEG block nor the TPE moiety have disturbed the smectic phase of PAChol block. The copolymer PEG45-TPE-b-PAChol16 with a hydrophilic weight ratio of 18% was chosen for the self-assembly study using the nanoprecipitation method. The initial co-solvent (THF or dioxane) has a significant effect on the self-assembly behaviour of PEG45-TPE-b-PAChol16. Elongated ellipsoidal or tubular polymersomes with smectic ordering in the membrane were formed in the dioxane/water system, while spherical solid nanoparticles were obtained in THF/water system. This difference was explained by the different solubility of hydrophobic block PAChol in dioxane and in THF. All polymersomes and nanoparticles showed AIE fluorescence. Their quantum yield was around 0.7-0.8% as measured by an indirect method using a standard fluorogen. The quantum yield and brightness may be improved by introducing more AIEgen units between the hydrophilic and hydrophobic blocks. The present work shows the possibility to obtain elongated polymerosmes with both smectic LC orders and AIE fluorescence. By replacing the prototypical TPE with AIEgens that can be excited in longer wavelength than UV, these LC polymersomes and cholesterol-based nanoparticles can be designed in the future as fluorescent carriers of drugs. 13 This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins

Figure 1 .

 1 Figure 1. Amphiphilic block copolymers with cholesteryl-based LC polymer, PEG-b-PAChol, 7 PEG-SS-PAChol, 17 and PEG-TPE-b-PAChol.

Figure 2 .

 2 Figure 2. (a) Synthetic route to the TPE-functionalized macro chain transfer agent mPEG-TPE-CTA and the LC amphiphilic block copolymer PEG-TPE-b-PAChol. (b) SEC chromatogram of PEG-TPE-b-PAChol copolymers (Table1) with THF as the eluent.

  Figure S1-4 for details). The TPE-functionalized macro chain transfer agent (CTA), mPEG-TPE-CTA, was first synthesized from monomethoxy polyethylene glycol (mPEG45-OH, 2000 Da). The diblock copolymer PEG-TPE-b-PAChol was then synthesized by RAFT This journal is © The Royal Society of Chemistry 20xx Please do not adjust margins Please do not adjust margins

Table 1 .a

 1 LC amphiphilic block copolymers PEG-TPE-b-PAChol. Number average degrees of polymerization (DPPAChol) and molecular weights (Mn, PAChol) of PAChol block evaluated by 1 H NMR. b Number average molecular weights of block copolymers evaluated by 1 H NMR. Mn = Mn, PAChol + 2680 (Da) (2680 Da being the molecular weight of mPEG-TPE-CTA). c Polydispersity index (Đ = Mw/Mn) of copolymers calculated from weight and number average molecular weight (Mw and Mn) evaluated by size exclusion chromatography (SEC) with tetrahydrofuran as eluent and polystyrene as standard. d Weight ratio of the hydrophilic part (fPEG) calculated by the molecular weight of PEG (2000 Da) over Mn of PEG-TPE-b-PAChol.

Figure 3 .

 3 Figure 3. (a) POM image of PEG45-TPE-b-PAChol16 at 128°C showing typical fan-shape textures of smectic A phase. (b) SAXS pattern of PEG45-TPE-b-PAChol16 fiber at room temperature, which was drawn at 135°C and cooled to 25°C. The red arrow indicates the drawing direction of fiber. A layer spacing of P = 4.34 nm was measured for the smectic A phase with its layer normal perpendicular to the fiber long axis as shown in the inset of (b) (Inset: red ellipses represent mesogens and black lines polyacrylate backbones).

  , b), which are similar to other cholesterolbased LC polymersomes. The diameter along the short axis is about 100 nm, but the size along the long axis extends to 500 nm or several microns. Some long vesicles are cut, only because of the cryo-EM sample preparation (the sample was flush frozen Please do not adjust margins Please do not adjust margins in a very thin vitreous ice layer of ~ 100 nm thickness); they should have integral shapes when dispersed in water. Since DLS was not a suitable technique to characterize the non-spherical particles, SEM was used to observe the whole size of these elongated vesicles. As shown in Figure4c, micron-sized crashed ellipsoidal or elongated objects corresponding to the dehydrated vesicles were visible. Despite the heterogeneity in the size of polymersomes, their membrane thickness was homogeneous with a value of 10 ± 1 nm as measured by statistical analysis of the cryo-EM images of 30 different vesicles. As the molecular length of one AChol unit in the extended state on the PEG45-TPE-b-PAChol16 copolymer chain is about 0.3 nm (calculated by Chem3D software), the maximum molecular length of the PAChol block in PEG45-TPE-b-PAChol16 should be 4.8 nm. This value was approximately half of the membrane thickness. It suggests that the LC polymersomes of PEG45-TPE-b-PAChol16 are unilamellar and the vesicle membrane has a bilayer structure in which the stretched PAChol polymers are along the normal of the membrane.

Figure 4 .

 4 Figure 4. Cryo-EM (a-b) and SEM (c) images of the self-assemblies formed by PEG45-TPEb-PAChol16 copolymer in dioxane/water system. The inset in image (a) is the fast Fourier transform (FFT) of the smectic structure present in the red rectangle. d) Schematic illustration of the aniosotropic bending of polymer bilayer with smectic stripes in an ellipsoidal polymersome (adopted from Ref. 7). Red ellipses represent mesogens, black lines polyacrylate backbones, and blue wavy lines PEG chains. Z axis is along the normal of the membrane. e is the membrane thickness and P is the smectic period.

Figure 5 .

 5 Figure 5. SEM images and size distributions by DLS (number profiles) of the selfassemblies formed by PEG45-TPE-b-PAChol16 in THF/water system with the initial polymer concentration of 2 mg/mL (a, b) and of 1 mg/mL (c, d), respectively.

Figure 6 .

 6 Figure 6. Photographs of PEG45-TPE-b-PAChol16 self-assemblies' dispersions in water observed under UV light (365 nm). Left: polymersomes obtained by nanoprecipitation in dioxane/water system. Right: spherical nanoparticles obtained by nanoprecipitation in THF/water system.

Figure 7 .

 7 Figure 7. (a) and (c) Photographs under UV lamp (λ = 365 nm) of PEG45-TPE-b-PAChol16 in the mixtures of dioxane/water (a) and in the mixtures of THF/water (c) with the same copolymer concentration (10 -4 mol/L) but different water fractions fw. (b) and (d) Photoluminescence (PL) spectra of PEG45-TPE-b-PAChol16 in the mixtures of dioxane/water (b) and in the mixtures of THF/water with the same copolymer concentration (10 -4 mol/L) but different water fractions fw (λex= 349 nm).

Acknowledgements

This project has received financial support from the CNRS through the MITI interdisciplinary programs and the French National Research Agency (ANR-16-CE29-0028). Nian Zhang and Yujiao Fan gratefully acknowledge the China Scholarship Council for funding their PhD scholarship. The Multimodal Imaging Centre of Institut Curie is acknowledged for providing access to cryo-EM facility in Orsay.

Author Contributions

The manuscript was written through contributions of all authors. Y.

Conflicts of interest

There are no conflicts to declare.


Emmanuel Rodriguez

Georges-Pierre Bonneau

Stefanie Hahmann

Mélina Skouras

Computational Design of Laser-Cut Bending-Active Structures

Keywords: active bending, laser-cutting, metamaterial, inverse design, computational fabrication

We propose a method to automatically design bending-active structures, made of wood, whose silhouettes at equilibrium match desired target curves. Our approach is based on the use of a parametric pattern that is regularly laser-cut on the structure and that allows us to locally modulate the bending stiffness of the material. To make the problem tractable, we rely on a two-scale approach where we first compute the mapping between the average mechanical properties of periodically laser-cut samples of mdf wood, treated here as metamaterials, and the stiffness parameters of a reduced 2D model; then, given an input target shape, we automatically select the parameters of this reduced model that give us the desired silhouette profile. We validate our method both numerically and experimentally by fabricating a number of full scale structures of varied target shapes.

Introduction

Bending-active structures are curved structures made of initially planar components that have been elastically bent so that they assume a 3D shape [START_REF] Lienhard | Active Bending, a Review on Structures where Bending is Used as a Self-Formation Process[END_REF]. This construction technique offers several advantages: flat panels are easy to manufacture, store and transport, and can be used to create objects of multiple scales, from models to large-scale structures. These multiple benefits explain the popularity of active bending for design applications as well as in architecture [START_REF] Fleischmann | ICD/ITKE Research Pavilion: A Case Study of Multi-disciplinary Collaborative Computational Design[END_REF][START_REF] Laccone | FlexMaps Pavilion: a twisted arc made of mesostructured flat flexible panels[END_REF][START_REF] Kalantar | Informing Deformable Formworks -Parameterizing Deformation Behavior of a Non-Stretchable Membrane via Kerfing[END_REF].

In this work, we are interested in controlling the shapes of the silhouettes of bent planar strips of plywood. More specifically, we want the profile of the deformed structure, subject to external loads, to match a given target curve when at equilibrium. While this can be achieved by playing with the width of the strips (see e.g. the recent work by Hafner and Bickel [START_REF] Hafner | The Design Space of Plane ElasticCurves[END_REF]) or their thickness (at the expense of a more involved fabrication process), we do not want to alter the outer shape of our ribbons. Instead, drawing inspiration from the design of living hinges and from kerfing techniques used in lute-making, we propose to modulate the curvature of our structures by patterning their interior according to a parameterized rectangular motif that is regularly laser-cut. Inserting such cuts to the initially flat panel through laser-cutting is no more difficult than cutting its external layout. Our objects can thus be fabricated using a widespread, scalable and affordable technique.

On the computational side, relying on a parametric cut pattern allows us to treat the cut fiberboard as a metamaterial with programmable bending capabilities and to turn towards two-scale methods when solving the inverse design problem of finding cut parameters that will produce a bending-active structure with the desired target profile. To our knowledge, such an approach, leveraging explicit bending control and building on recent results in bending homogenization [START_REF] Sperl | Homogenized yarn-level cloth[END_REF], has not yet been proposed in the context of computational fabrication.

While we focus on the particular setting of 2D structures and rectangular cuts, we believe that our technique, that decouples Email addresses: melina.skouras@inria.fr (Mélina Skouras), georges-pierre.bonneau@inria.fr (Mélina Skouras) Figure 1: Plots from top to bottom: an input curve with prescribed boundary conditions and external forces, the optimized stiffness profile of the bendingactive structure, the geometric parameters of the cutting patterns. Bottom: the cut layout. Right: the laser-cut bending-active structure whose silhouette matches the input target curve. macro-scale geometry from material inner structure paves the way to generalization. Indeed, beyond the mechanics of the final object, the shapes of the cuts also impact its appearance. Therefore, playing with the cut patterns also allows to play with the aesthetics of the structure. Our methodology, which does not explicitly rely on a specific pattern geometry should be extendable to other motifs, as long as they form a parametric family and exhibit a wide range of mechanical behaviors.

In summary, the main contributions of our work are the following:

• We propose a novel "living-hinge"-inspired parametric family of laser-cut metamaterials with tunable bending capabilities.

• We characterize the bending properties of these metamaterials using a recent bending homogenization method "with a twist".

• We present a two-scale algorithm for the inverse design of bending-active structures of custom deformed shape.

Related works

Fabrication-aware computational design has gained increasing interest in computer graphics for a decade. For a broad overview of the field, we refer the reader to the surveys by Bickel et al. [START_REF] Bickel | State of the Art on Stylized Fabrication[END_REF] and Bermano et al. [START_REF] Bermano | State of the Art in Methods and Representations for Fabrication-Aware Design[END_REF], which cover a wide variety of techniques. In the following section, we focus the discussion on previous work on inverse design of deployables and slender structures, metamaterials and Kirigami-inspired research.

Deployable surfaces. These structures share with ours that they are fabricated flat and assume a 3D shape once deformed. Existing fabrication techniques include the use of flat beams for weaved items [START_REF] Vekhter | Weaving Geodesic Foliations[END_REF][START_REF] Ren | 3D Weaving with Curved Ribbons[END_REF] and for gridshells [START_REF] Panetta | X-Shells: A New Class of Deployable Beam Structures[END_REF][START_REF] Pillwein | On Elastic Geodesic Grids and Their Planar to Spatial Deployment[END_REF][START_REF] Pillwein | Generalized Deployable Elastic Geodesic Grids[END_REF], well adapted to large-scale structures, e.g., architectural installations, sealed inflatable membranes [START_REF] Siéfert | Programming stiff inflatable shells from planar patterned fabrics[END_REF][START_REF] Panetta | Computational Inverse Design of Surface-based Inflatables[END_REF] and 3D printed structures [START_REF] Malomo | FlexMaps: Computational Design of Flat Flexible Shells for Shaping 3D Objects[END_REF] for smaller scale models, wire meshes [START_REF] Garg | Wire Mesh Design[END_REF], as well as laser-cut rigid panels [START_REF] Konaković | Beyond Developable: Computational Design and Fabrication with Auxetic Materials[END_REF][START_REF] Laccone | FlexMaps Pavilion: a twisted arc made of mesostructured flat flexible panels[END_REF] or softer elastic sheets [START_REF] Celli | Shape-morphing architected sheets with non-periodic cut patterns[END_REF][START_REF] Chen | Bistable Auxetic Surface Structures[END_REF]. Common to all these approaches is that curvature is typically obtained from metric frustration. This is why the actual 3D structures generally exhibit non-zero Gaussian curvature and are doubly curved. By contrast, we are interested in controlling the extrinsic curvature of our strips, which cannot be achieved by modulating the surface metric alone.

Rods' design. Closer to our work are the approaches that address the inverse design of 2D and 3D rods and other slender structures so that they match a given target shape at equilibrium. For example, Derouet-Jourdan et al. [START_REF] Derouet-Jourdan | Floating tangents for approximating spatial curves with G1 piecewise helices[END_REF] optimize the physical parameters of elasticae to avoid sagging under gravity, and Bertails-Descoubes et al. [START_REF] Bertails-Descoubes | Inverse design of an isotropic suspended Kirchhoff rod: theoretical and numerical results on the uniqueness of the natural shape[END_REF] prove the uniqueness of the natural shape of the rods solution to the former problem. In the context of computational fabrication, Pérez et al. [START_REF] Pérez | Design and Fabrication of Flexible Rod Meshes[END_REF] control the bending stiffness of networks of rods by locally adjusting their cross-sections and rest centerlines. More recently, Liu et al. [START_REF] Liu | Tapered elasticae as a route for axisymmetric morphing structures[END_REF] explored the design of axis-symmetric structures formed by tapered ribbons. Their theoretical results were extended and generalized by Hafner and Bickel [START_REF] Hafner | The Design Space of Plane ElasticCurves[END_REF] who investigated the design space of planar elastic rods for form-finding applications. This latter work shares similarities with ours but also notable differences: on the fabrication side, bending control of the optimized strips is achieved by playing with their width, while we modify their inner meso-scale geometric structure; on the numerical side, this previous study is focused on clamped-clamped ribbons whose deformed shape is primarily determined by geometric boundary conditions, while we are also interested in structures able to carry external loading.

Metamaterials and two-scale approaches. In a broader context, a popular approach to design an object with a given target shape consists in modifying the structure of the object at a local scale, so as to tune its local geometric and mechanical properties [START_REF] Panetta | Elastic Textures for Additive Fabrication[END_REF][START_REF] Tozoni | A Low-Parametric Rhombic Microstructure Family for Irregular Lattices[END_REF][START_REF] Schumacher | Microstructures to Control Elasticity in 3D Printing[END_REF]. For deployable surfaces, as the ones mentioned above, the intent is generally to be able to control the in-plane behavior of the surface; many related works thus exploit so-called auxetic metamaterials, i.e. engineered materials that can stretch, or contract, in all tangential directions simultaneously [START_REF] Konaković | Beyond Developable: Computational Design and Fabrication with Auxetic Materials[END_REF][START_REF] Malomo | FlexMaps: Computational Design of Flat Flexible Shells for Shaping 3D Objects[END_REF][START_REF] Chen | Bistable Auxetic Surface Structures[END_REF]. The use of structural patterns that are quasi-periodic has the big advantage of allowing for the use of multi-scale approaches when designing an inverse modeling tool. These methods typically proceed in two steps: the local macroscale properties of the material are first homogenized; then they are used as parameters in an inverse problem [START_REF] Panetta | Elastic Textures for Additive Fabrication[END_REF][START_REF] Zhu | Two-Scale Topology Optimization with Microstructures[END_REF]. In this work, we also rely on this type of approach, but rather than exploiting the in-plane properties of the considered metamaterial, we play with its out-of-plane mechanical behavior. While some reduced models may indirectly capture the flexural resistance of the proposed metamaterials [START_REF] Malomo | FlexMaps: Computational Design of Flat Flexible Shells for Shaping 3D Objects[END_REF][START_REF] Leimer | Reduced-Order Simulation of Flexible Meta-Materials[END_REF], as far as we know, the only graphics works that explicitly looked into the homogenization of the bending behavior of structured material sheets are those by Schumacher et al. [START_REF] Schumacher | Mechanical Characterization of Structured Sheet Materials[END_REF] and Sperl et al. [START_REF] Sperl | Homogenized yarn-level cloth[END_REF]. In our work, we build on the latter to characterize the properties of our laser-cut plates. However, rather than applying the method to speed-up forward simulations, we investigate its use in the context of inverse design.

Kirigami. The fabrication technique that we leverage, based on the cutting of a relatively thin plate, is reminiscent from Kirigami, a Japanese art consisting in cutting and folding paper in order to create 3D surfaces. Kirigami has recently inspired several researchers from the material science community in the context of various applications, from the design of nano-structures of complex unfolded shapes [START_REF] Zhang | A mechanically driven form of Kirigami as a route to 3D mesostructures in micro/nanomembranes[END_REF][START_REF] Liu | Nano-kirigami with giant optical chirality[END_REF] to the one of metamaterials formed by thin sheets of material regularly cut according to a (quasi) periodic pattern [START_REF] Rafsanjani | Buckling-Induced Kirigami[END_REF][START_REF] Rafsanjani | Propagation of pop ups in kirigami shells[END_REF]. However, most previous work in the latter category, closest to our work, investigated the mechanical behavior of these metamaterials under tension, while we are interested in their bending response under non-axial loading.

Overview

Given a planar target curve, external planar forces and boundary conditions (Fig. 1, top-left), our method produces a bendingactive structure made of a rectangular piece of wood panel in which holes are laser-cut (Fig. 1, right). The silhouette of the structure in physical equilibrium is aimed to match the target under the given loading conditions. The actual output of our process are the geometric parameters (Fig. 1, curve-plots) defining the cut layout (Fig. 1, bottom) which locally modulates the bending stiffness of the material.

Working on a full scale model, i.e. optimizing the geometry of the cuts in 3D directly, would be intractable. We therefore rely on an intermediate, reduced, 2D representation, and we solve an inverse problem involving the parameters of this reduced model. In practice, our approach involves three stages:

• building the reduced model (a discrete version of a 2D Kirchhoff rod, in our case) and mapping its parameters, here its mechanical properties, to the parameters of our original problem, i.e. the geometry of the cuts . We will do this by locally treating the cut wood as a parametric metamaterial and relying on an homogenization approach (Section 4);

• casting our original inverse design problem as a constrained optimization problem, whose unknowns are the parameters of the reduced model, i.e. its nodal bending stiffnesses (Section 5);

• mapping back the optimized structure's stiffness profile to cut patterns to infer the final cut layout.

We validate our algorithm by designing a number of bendingactive structures, four of which have also been fabricated (Section 6).

Laser-cut metamaterial

Metamaterials are engineered materials designed to exhibit superior mechanical properties to those of the constituent material. They are structured with periodically arranged or spatially varying building blocks which either define a homogeneous, resp. heterogeneous metamaterial. Tailor-made geometry of the building blocks is the key for gaining extraordinary properties.

Designing a metamaterial, and characterizing its mechanical behavior, is not only a matter of geometry, but also of scale. In our setting, where bending-active structures are made of wood panels in which a geometric pattern of holes is laser-cut to tune the mechanical behavior, the latter is modeled at two scales: by a so-called mesoscopic mechanical model at the scale of material, which corresponds to a solid material with geometrically parameterized holes, and by a so-called macroscopic mechanical model at the scale of the complete object.

Furthermore, the models are simplified at 2 levels:

• At the material level: we approximate a heterogeneous structure made of solid material and geometrically parameterized holes by a homogeneous material having another material law.

• At the level of the representation of the full object: we approximate the geometry of a 3D strip with finite thickness, i.e. a 3D object, by a planar curve. Indeed, we are interested in the silhouette of the ribbon, which has constant width and whose patterns remain periodic along the y axis. We therefore expect its motion to stay planar, which allows us to represent the final geometry by a 2D curve. We will model its mechanical behaviour using a Kirchhoff rod model.

We start by defining the geometric structure of our metamaterial and its parameters (Sect. 4.1) followed by the introduction of the mesoscopic and macroscopic mechanical models (Sect. 4.2 and 4.3). To investigate the macro mechanical behavior and link it to the material parameters of our rod model, we turn to numerical homogenization and explain how we handle periodicity (Sect. 4.4). In Sect. 4.5 we build the map from material parameters to geometric parameters.

Geometric structure

In designing the laser-cut patterns, we were inspired by numerous "living-hinge" examples and design tools [START_REF] Brth | Boxdesigner Galerie[END_REF]. A livinghinge can be loosely defined as a flexible piece of material connecting two planar rigid pieces of the same material along a line of hinge, so that the latter pieces can rotate around this line. Living hinges exploit indeed the natural tendency for material to bend around an array of cuts or notches. We choose one of the most standard cutting pattern used in living-hinges, in which rectangular holes are cut and arranged to form an alternating pattern, as illustrated in Fig. 2. Our bending-active structures act as beams where the longitudinal axis will be referred to as the xaxis and the transverse direction will be denoted y-axis. A structure of longitudinal dimension L struct and transversal dimension W struct is periodically cut. The periodic unit cell (dashed rectangle in Fig. 2) has dimensions L c × 2W c . Note, that the longest dimension L c of the unit cell is in the transverse direction of the whole structure. L c (resp. W c ) is also the offset between rectangular cuts in the transverse (resp. longitudinal) direction. The rectangular cuts have dimensions l c × w c . While L c and W c are constant throughout the whole structure, l c and w c vary in the longitudinal direction to be able to modify locally the bending stiffness. Instead of measuring these variations with the absolute values w c and l c in meters, we choose as design parameters the adimensional ratios α = l c L c (transverse direction) and β = w c W c

(longitudinal direction). α and β are spatially variable parameters, depending on the curvilinear x-abscissa (see Fig. 2). By varying α and β in the longitudinal direction we obtain beams which may be locally rigid (low values of α and β) or flexible (high values of α and β).

L struct W struct W c w c L c l c L c /2 x y z

Mesoscopic mechanical model

In this Section we introduce the physical model which we choose to guide the mechanical behavior at the mesoscopic scale of our metamaterial. At this scale the meta-material is modeled in 3D with the geometry illustrated in Fig. 2, extruded in the zdirection.

Base material modeling. The constituent material we use for our physical realizations is MDF (Medium Density Fiberboard). MDF is a composite of fine, uniform wood fibers, which are bonded together by resins, heat and pressure. It is a heterogeneous material at the microscopic scale, and has a non-linear elastic behavior [START_REF] Sliseris | Numerical prediction of the stiffness and strength of medium density fiberboards[END_REF]. We approximate it as a homogeneous and isotropic material at the mesoscopic scale, using a standard Neo-Hookean material model with constitutive equation [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF]

Ψ Neo = µ 2 (Tr(C) -3) -µ log(J) + λ 2 (log(J)) 2 , ( 1
)
where C is the right Cauchy-Green tensor, J = det(C) and λ, µ are the Lam coefficients, derived from the Young's modulus E and the Poisson's ratio ν using λ = Eν (1+ν)(1-2ν) and µ = E 2(1+ν) . The Poisson's ratio is set to ν = 0.3. Since bending is the most salient characteristic governing the deformation behavior in our application, we experimentally fit the Young's modulus with a cantilever test. In practice, we measured a deflection of 35 mm for a beam of length 79 cm and square section of 6 mm side, fixed at one of its ends, from which we deduced a Young's modulus of E = 9.95 GPa using the theoretical Euler-Bernoulli beam model. Numerical model. When simulating the structure at the mesoscale, we need to numerically integrate the energy density [START_REF] Lienhard | Active Bending, a Review on Structures where Bending is Used as a Self-Formation Process[END_REF]. To this end, we discretize the 3D structure using Gmesh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF] with tetrahedral elements of typical edge size of 0.4 mm to 1.6 mm (exact figures are provided where appropriate) and solve the governing equations using the Finite Element solver FEniCS [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF] with Lagrange elements of degree 2.

Macroscopic model

Because of our particular setting, where we want to reproduce 2D curves in the xz-plane, and the external forces are defined in the same plane, our bending-active structure behaves as a rod which deforms in this plane with no twist. Furthermore, since the primary mode of deformation of our structure should be bending, we assume that stretching can be neglected. Therefore at the macroscopic scale it is convenient to model our structure as a Kirchhoff rod in 2D, whose bending energy density (i.e. energy per unit length) is defined as

δE b = 1 2 Y Iκ 2 ,
where Y is the Young's modulus of the homogenized material, I is the second moment of area of the rod's cross section, and κ is the curvature of the rod [START_REF] Audoly | Elasticity and geometry: from hair curls to the nonlinear response of shells[END_REF].

For a rod with a rectangular cross-section bending with respect to the y axis (using the axis notations of Fig. 3), I = 1 12 t 3 W, where t and W are, respectively, the thickness of the rod and its width.

The total bending energy of a rod of length L is then

E rod b = L 0 δE b dl = L 0 1 2 Y Iκ 2 dl, i.e. E rod b is of the form E rod b = L 0 1 2 k b Wκ 2 dl, with k b = 1 12 Yt 3 . (2)
Discrete model. In practice, we represent our rod using a simple sequence of linear and angular springs connecting the discrete nodes. The rod of length L struct corresponding to our object profile is discretized into n nodes with position xi , 1 ≤ i ≤ n, and n -1 edges

e i = xi+1 -xi , 1 ≤ i ≤ n -1.
There are n -2 bending stiffnesses k i b , 2 ≤ i ≤ n -1, in Joule, assigned to each interior node. Following Bergou et al. [START_REF] Bergou | Discrete Elastic Rods[END_REF], we discretize the rod's bending energy [START_REF] Fleischmann | ICD/ITKE Research Pavilion: A Case Study of Multi-disciplinary Collaborative Computational Design[END_REF] as

E bending = n-1 i=2 1 2 k i b W κb i 2 d rest i ,
where W is the width of the structure,

d rest i = e rest i-1 + e rest i 2
is the Voronoi length of the interior node xi and κb i is the discrete curvature binormal vector given by (refer to [START_REF] Bergou | Discrete Elastic Rods[END_REF] for more details):

κb i = 2e i-1 × e i e i-1 e i + e i-1 • e i .
To model quasi inextensibility, we assign a linear spring to each edge with same elongation stiffness k e , in Newton, set to a high value (k e = 1000 N) and compute the discrete stretching energy of the rod as

E stretching = n-1 i=1 1 2 k e e i -e rest i 2 e rest i .
Finally, a mass m i in kg is attached to each node in order to take into account gravity forces due to the structure itself. The total energy at the macroscopic level is then given by:

E macro = E stretching + E bending + n i=1 gm i z i -W ext. forces ,
where z i is the height of xi , g = 9.81m s -2 , and W ext. forces is the work of the external forces. By modifying the geometric parameters α, β along the longitudinal dimension (Section 4.1), we can adjust the values of the n -2 bending stiffnesses k i b of the macroscopic model, and thus control the shape at equilibrium x = ( x1 , • • • , xn ) which is a solution of:

x = argmin xi E macro (3)
s. t. boundary conditions.

In the next section 4.4 we will explain how to compute the bending stiffnesses k i b from the longitudinally varying geometric parameters α, β.

Forward modeling refers to the process of manually tuning the geometric parameters α, β, and solve (3) until the input shape is reproduced with satisfaction. We have implemented the solution of (3) using the optimization library IPOPT [START_REF] Wchter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF], in order to quickly test different design ideas in terms of boundary conditions and external forces.

In Section 5 we will invert this process, and automatically find the geometric parameters allowing to reproduce the input shape.

Homogenization

In this Section we explain how to compute the mapping from the geometric parameters α, β of the model at the mesoscopic scale, to the bending stiffnesses k b of the model at the macroscopic scale. Such a mapping between parameters at mesoscopic and at macroscopic scale can be derived through the framework of numerical homogenization.

The main idea of numerical homogenization is to determine effective macroscopic mechanical properties, such as the elasticity tensor, based on representative simulations of a periodic unit cell at mesoscopic scale (theoretical foundations in [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Guedes | Preprocessing and postprocessing for materials based on the homogenization method with adaptive fini[END_REF]). Application of numerical homogenization in our case is closely related to that of [START_REF] Schumacher | Mechanical Characterization of Structured Sheet Materials[END_REF][START_REF] Sperl | Homogenized yarn-level cloth[END_REF], where periodic boundary conditions for the case of thin shell bending deformation have been considered. We follow the notations of [START_REF] Sperl | Homogenized yarn-level cloth[END_REF]. X is a position at mesoscopic scale and at rest. x(X) is the position after deformation. The theory of homogenization assumes that at any point of a deformed material observed at a macroscopic scale, there exists a mesoscopic scale representation, valid in a representative volume element (RVE), centered around the point. This representation expresses the points (at mesoscopic scale) x(X) as the sum of a first order approximation of the macroscopic deformation, and a mesoscopic scale displacement fluctuation ũ(X). Two conditions must be fulfilled by ũ(X). The first condition is the averaging constraint RVE ũ = 0, which intuitively states that the local fluctuations average out at the macroscopic scale. The second condition is the periodicity constraint ũ+ = ũwhich states that the local fluctuation takes the same values on both sides of the RVE. In other words, periodic copies of the locally displaced RVE join continuously.

Bending homogenization. In the case of a bending deformation, the two-scale representation can be stated as follows, analogously to [START_REF] Sperl | Homogenized yarn-level cloth[END_REF]:

x(X) = R θ (X) + ũ(X), (4)
where R θ (X) denotes the bending deformation observed at the macroscopic scale. Fig. 4-middle shows a schematic illustration of this decomposition. Fig. 3 shows the RVE at rest (dashed The RVE is composed of adjacent periodic unit cells (dashed rectangle in Fig. 2). Choosing the RVE too small would result in overconstraining the local fluctuations ũ. On the other hand, choosing it too large increases the computational time quickly. Using the method explained below, we computed ũ for increasing size of the RVE and found that beyond a size of 2 unit cells in longitudinal direction, and 1 unit cell in transverse direction, negligible difference is observed between the successive displacement functions. Therefore, we choose a RVE of size L = 4W c in longitudinal direction, L c in the transverse direction, and t in the z-direction (t is the thickness of the wood panel). With these notations, the macroscopic bending deformation R θ is given by:

R θ : X =           X Y Z           → R θ (X) =           ( L θ + Z) sin( Xθ L ) Y ( L θ + Z) cos( Xθ L ) -L θ           Figure 4:
The local displacement fluctuation ũ of the RVE is mapped to a local displacement fluctuation ṽ in a computational domain. This mapping allows to simplify the periodic boundary conditions ( 6) and the averaging constraints [START_REF] Bickel | State of the Art on Stylized Fabrication[END_REF].

As pointed out by [START_REF] Sperl | Homogenized yarn-level cloth[END_REF], the averaging and periodicity constraints on the local fluctuations ũ must be adapted to ensure these fluctuations average out properly when applied to a bent RVE. Therefore [START_REF] Sperl | Homogenized yarn-level cloth[END_REF] introduces co-rotated periodicity conditions and averaging constraints. While following a similar approach, we map the two-scale representation (4) to a computational domain, and express the averaging and periodicity conditions in this domain. Referring to Fig. 4, we define a new displacement field ṽ by: ṽ(X) = R -1 θ ( ũ(X) + R θ (X)) -X The two-scale representation (4) maps to:

R -1 θ (x(X)) = X + ṽ(X). (5)
In other words, ṽ displaces points at rest (X) while ũ displaces R θ (X). In the computational domain the co-rotated periodicity conditions on ũ can be expressed as simple periodicity conditions (6) on ṽ. And the co-rotated averaging constraints on ũ are expressed as simple averaging constraints [START_REF] Bickel | State of the Art on Stylized Fabrication[END_REF] on ṽ. This simplifies the implementation since standard FEM libraries have builtin capacities for handling periodic finite element basis functions, fulfilling [START_REF] Sperl | Homogenized yarn-level cloth[END_REF].

ṽ+ = ṽ- ( 6
) RVE ṽ = 0. ( 7
)
Finally, we can now compute the local fluctuations ṽ, and the potential energy E 3D b (θ) of the deformed RVE, by solving the following static equilibrium problem: at (6) and ( 7) are satisfied.

E 3D b (θ) = min ṽ RVE Ψ Neo (ṽ), (8)
s. t. ( 6) and ( 7) are satisfied. Computing the bending stiffness of the reduced model. As mentioned in Section 4.3, the bending term of our reduced model discretizes that of a Kirchhoff rod that reads

E rod b = L 0 1 2 k b Wκ 2 dl
in the continuous setting. Therefore, the bending energy of a deformed rod with uniform curvature κ = 1 r and constant stiffness k b is

E rod b = 1 2 k b W 1 r 2 L.
Denoting by θ the opening angle of the arc formed by the rod (see Fig. 3), E b can be rewritten as a quadratic function of θ as

E rod b (θ) = 1 2L k b Wθ 2 ,
since, in that case, L = rθ.

Our goal is to estimate, for each considered (α, β) pair, the bending stiffness k b of the reduced model that will best capture the macro-scale mechanical behavior of the corresponding lasercut strip of wooden MDF. We do so by matching, in the least square sense, the energies E rod b and E 3D b of the two models corresponding to bending deformations of increasing curvature, gravity turned off, i.e. we find the optimal kb defined by kb = argmin

k b 1 2 i (E rod b (θ i ) -E 3D b (θ i )) 2 = i a i E 3D b (θ i ) i a 2 i , (9)
where a i = 1 2L Wθ 2 i and the (θ i ) are to the opening angles corresponding to these deformation tests.

In practice, we consider RVEs of two unit cells in length that we bend up to a maximum angle of θ max = 0.2 rad. For unit cells of width W c = 3 mm like ours, this corresponds to a radius of curvature of 3 cm, which is enough for our application. We define (θ i ) as a uniform sampling of [0, θ max ] using 11 values.

As shown in Fig. 6, left, the Kirchhoff rod model captures well the behavior of the cut 3D plates for a wide range of bending deformations.

Building the map from bending stiffnesses to geometric parameters

We initially considered varying both the relative length α and width β of our cut pattern. In practice, to explore the range of achievable bending stiffnesses, we uniformly sampled α and β between 0.3 and 0.9, using increments of 0.05, and we estimated the corresponding k b using a relatively coarse 3D mesh (edge size of 1.6 mm) (see Fig. 6, middle). Noting that the length of the cuts had a much bigger impact on the structure's bending stiffness than the cuts' width, we then decided to keep β fixed and to play with the α value only. We set β to 0.6, the middle of the original sampled range and recomputed the bending stiffnesses using a finer sampling for α (increments of 0.01) and a finer discretization for the 3D meshes (maximum element size of 0.4 mm), for a total of about 50k to 100k elements (the maximum that would fit in the memory of our computer), depending on the α value (Fig. 6, right). In practice, to compute one value of bending stiffness kb we need to solve 11 times problem [START_REF] Bermano | State of the Art in Methods and Representations for Fabrication-Aware Design[END_REF]. This process takes between 1h to 5h with these mesh dimensions (see Section 6.1 for the specification of our computing platform). All values of k b are computed and stored as a pre-process.

Despite a bit of noise for small values of α, i.e. for highly rigid materials (likely due to high sensitivity to the mesh discretization in this regime), we observe that the bending stiffness k b strictly decreases as the length of the cut increases, as one would expect. This allows us to invert the relation and to build the map from stiffness to geometric parameter, α(k b ), by linearly interpolating the computed k b values, after smoothing the latter thanks to a basic Gaussian filter (average of immediate neighbours), applied five times (see Fig. 6, right).

Inverse design algorithm

We now address an inverse design problem for laser-cut MDF panels. Remember our initial problem: given a 2D target profile shape and user-defined external forces and boundary conditions, we are searching for a cut layout of a planar panel of fixed length and width, such that the silhouette of the deformed panel matches the target under physical equilibrium.

The approach for solving this inverse modeling problem consists of two main steps. First, a constrained optimization problem is solved to find the nodal stiffness values and node positions of the reduced model at physical equilibrium matching the profile curve. Second, the optimized stiffness values are mapped back to the cut pattern parameters to infer the final cut layout.

In addition to the target curve, the user determines the boundary conditions. Our current implementation (see Section 4.3) enables to prescribe fixed points, tangent vectors at any point, individual spring stiffnesses and a mass added either at one end point or at one interior point, but can be extended to others.

Optimization of the stiffness profile

The purpose of the first step of our inverse modeling approach is to compute simultaneously the deformed shape and the corresponding stiffness values best approximating the target shape. We formulate this inverse design problem in a constrained optimization framework. The target shape is assumed to be a smooth curvature continuous shape equidistantly sampled as a 2D polygonal curve T whose n vertices are denoted by xt ∈ IR 2n . The variables to be optimized are the vertex positions of the deformed object x and the material stiffness values at macro scale k b . The objective function Φ to be optimized is a combination of three metrics defined in the following: [START_REF] Ren | 3D Weaving with Curved Ribbons[END_REF] where k min = 0.19 J and k max = 60 J in our implementation. The lower bound of k b is provided by the maximal feasible cut size. Indeed, the bigger α is, the closer the cuts are to each other and the more the bending stiffness decreases (see Fig. 6, right). The material risks to break. We set k min = 0.19 J, which corresponds to less than 1 mm material left between the holes, i.e. α = 0.95. The upper bound is set by the maximal bending stiffness reached for the full material.

min k b , x Φ xt (k b , x) s.t. ∇E macro = 0 and k min ≤ k b ≤ k max
Shape similarity. The main role of the objective function Φ is to ensure that the final shape best matches the target shape. We encode fidelity to the target shape by relating the curvature values, instead of the curves themselves. This approach follows the idea that curvature is not only an intrinsic shape descriptor, but is known as an efficient shape similarity measure [START_REF] Arkin | An efficiently computable metric for comparing polygonal shapes[END_REF][START_REF] Veltkamp | Shape matching: similarity measures and algorithms[END_REF][START_REF] Gatzke | Curvature maps for local shape comparison[END_REF], since local 2D curvature is a highly informative visual cue for global shape perception, object recognition, and image interpretation [START_REF] Vergne | Flow-Guided Warping for Image-Based Shape Manipulation[END_REF]. The local curvature similarity between the deformed shape and the target shape is defined as

Φ similarity = c(x) -c(x t ) 2 IR n ,
where c(x) is a linearized discrete curvature vector at x = [x 1 , . . . , xn ] with c(x i ) = (x i+1 -2x i + xi-1 ) / ēi 2 . We divide by the squared length of the curve segment to achieve scale independence. In addition to enforce shape similarity, this penalty function also acts as a regularization term, since the target shape is assumed to be a smooth curvature continuous curve.

Distance to target. Euclidean distance between corresponding curve vertices cannot be ignored completely, because it prevents the curves to shift and rotate from each other. Indeed, two planar curves with same curvature vector are identical up to a translation. We therefore add a distance penalty defined as

Φ distance = x -xt 2 IR n
and assign it a very small weight.

Minimal stiffness variation. In addition to the geometric shape, we are searching for optimal bending stiffness. However, a regularization of the stiffness values is required in order to produce a smoothly varying geometric cut pattern (see Sect. 5.3). The latter is a necessary condition for the homogenization approach to be valid. We therefore introduce the following regularization penalty term Fig. 7 illustrates with an example that regularizing the stiffness even with a small weight is efficient and does not significantly alter the distance between the deformed shape and the target.

Φ regularization = n-1 i=2 | k i-1 b -2k i b + k i+1 b | 2 .
In summary, the objective function writes

Φ = ω 1 Φ similarity + ω 2 Φ distance + ω 3 Φ regularization .
In all our examples we have set (ω 1 , ω 1 , ω 3 ) = (10 4 , 10 -2 , 10 -2 ).

To solve (10) numerically, we use the open source software IPOPT [START_REF] Wchter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF], which implements an interior point optimization method to find a local solution.

The quasi-Newton BFGS implementation of IPOPT is applied. We need to compute the gradient of the energy E macro (equality constraint in ( 10)), and its Hessian (Jacobian of the equality constraint). In our current setting, we only consider point masses as external forces, therefore the gradient and the Hessian of W ext. forces are trivial to compute, but more complex external forces could be taken into account.

The optimization problem ( 10) is non-linear and non-convex, so that the result depends on an initial guess of the variables x and k b .

What works well in the present setting is to start from the target curve as initial guess for x and set the bending stiffness uniformly to a constant value 5.14 J, corresponding to α = 0.5. With this choice the optimization starts at the minimum of the objective function Φ but the important physical equilibrium constraint ∇E macro = 0 is not satisfied at the beginning. Alternatively we tried starting with the same uniform initial stiffness, and with x set as the physical equilibrium (solution of (3)). We also tried to set the initial x as the target curve, and to compute the stiffnesses minimizing E macro (which is a linear problem to solve). But we found the first approach generally finds a better local minimum of Φ.

Mass density adjustment

Cutting the material results in a significant loss of mass

ρ(α, β) = ρ • (1 -α • β),
where ρ is the mass density of the material without cuts and ρ the homogenized mass density of the cut material. As an example, for a cut with α = 0.5 and β = 0.5, a quarter of the material is removed. It is therefore important to solve the inverse design problem [START_REF] Ren | 3D Weaving with Curved Ribbons[END_REF] by taking into account the correct mass density of the structure.

Remember that E macro in [START_REF] Ren | 3D Weaving with Curved Ribbons[END_REF] is the total potential energy of the system, composed of internal energy, work of external forces and gravitational energy, that depends on the mass density of the structure. The correct mass density however depends on the amount of material to be cut-off following the cut parameters, which are only known at output of the system. It is therefore not possible to know exactly in advance the mass density of the output structure. To solve this mass adjustment problem we implement an iterative mass reduction algorithm, which iterates the inverse modeling step until convergence on the nodal masses, see Algorithm 1. In practice, it is observed that the difference in masses quickly becomes smaller than 1% after 2 iterations. The inset shows the result of the forward simulation of the output of inverse modeling after taking into account the mass reduction. For this example our algorithm converges in 3 iterations.

Alternatively, we could parameterize the density by the stiffness parameters since the density depends on the size of the cuts (linear relation) and we know the relation between the size of the cuts and the stiffness parameter (piecewise linear relation -so not C 1 ). However, this would mean that we introduce a non-smooth term in the optimization problem. It is well known that this kind of non-regular minimization problem is much more difficult to solve than a problem for which the objective and the constraint are both differentiable. We therefore prefer to proceed in two steps in alternation, which allows us to use a standard gradientbased optimization method. m initial mass when no cuts

m 0 ← 0 m 1 = list of nodal masses while ||m 0 -m 1 || 2 ≥ ε do k b = InverseModeling(m 1 , T ) (α, β) = Mapping(k b ) m 0 ← m 1 m 1 ← Mass(α, β) end while
In order to compute the masses from the geometric parameters α, β we need to know the mass density of the MDF wood panels we are using. We measured the mass density by cutting a piece of known size and weighing it with a precision balance and obtained a density of 878 kg • m -3 .

Cut layout generation

Up to now, we have computed the discrete stiffness parameters of a reduced 2D model that, at equilibrium, closely matches the structure in target pose (Sect. 5.1). We therefore have a piecewise linear relation from arc length of our 2D model to bending stiffnesses.

The goal of this step is now to match the computed stiffness profile with geometric parameters for the cuts. We do this by directly reading the value of α corresponding to the curvilinear abscissa (x-value in Fig. 2) of the mid-point of each cut, using the map built in Section 4.5. In practice, the length of the edges of our discretized rod coincide with the width 2W c of a unit cell so that the α that we compute are exactly those that correspond to the k b values that we have optimized. See an example in Fig. 8. 10) and the corresponding α cutting parameters cutting parameters for 320 unit cells of width 3 mm (right) taken from the light bulb structure shown in Fig. 10(a).

Results

We have used our design system to create a diverse set of bending-active structures whose silhouettes at equilibrium match different target shapes. We demonstrate the versatility of our system by combining the target shapes with various sets of boundary conditions and external loads. In addition to computing numerous results, we have fabricated four structures in order to validate the presented design system as well as the choices and hypothesis made in Sections 4 and 5.

Bending-active structures

We showcase a variety of bendingactive structures in Figures 9,10 and 11. In addition to various shapes with partly strongly varying curvature profiles, we experimented with various sets of boundary conditions (BC) and external loads. The inset shows the legend for the boundary conditions as indicated in the figures throughout the paper. Table 1 summarizes all model statistics.

Experiments were run on a Dell Latitude 5410 (2020) with Intel Core i7 processor running at 4GHz with 16GB RAM.

In Fig. 9 shows the direct output of the inverse design problem [START_REF] Ren | 3D Weaving with Curved Ribbons[END_REF]. In the left column, the target curves xt (black) are superposed by the optimized curves x (orange). The maximal relative distance to target (see Table 1) ranges from 0.01% to 0.2% for all examples, except for the heart with 1.2%. The right column shows the optimized stiffness profiles k b . The cut-layouts α i for each unit cell are then taken directly from the stiffness profile via the monotone inverse mapping described in Section 5.3. We do show a few layouts for other examples in Fig. 10.

Robustness with respect to the boundary conditions has been studied in row 5 of Fig. 9. A set of 4 similar target curves (cubic B-splines) have been created. The BC for each curve consist of only one fixed point and tangent at the bottom extremity and an extra load of 150 g attached at the other. The algorithm converges well in all case. The asymmetric arch (row 3 in Fig. 9) showcases an example, where the BC are symmetric, but not the shape.

Target shapes. It is important to know from the beginning if a given target curve is feasible. Similar to [START_REF] Hafner | The Design Space of Plane ElasticCurves[END_REF], we say that a target shape is "feasible" if there is a solution to the constrained minimisation problem for which the distance to the target curve is close to 0. In our case however, feasibility does not imply fabricability, meaning that the structure will not break during bending. Although we have specified limits on the stiffnesses to avoid this failure, these limits are not a theoretical guarantee. In order to prevent failure during bending, we therefore derive an empirical link between the target's geometry and the optimized stiffnesses by observing above which curvature different laser-cut panels break during bending. To this end, we have constructed 6 panels of 50 cm in length, each with a uniform cut layout corresponding to α ∈ {0.4, 0.5, 0.6, 0.7, 0.8, 0.9}, and bent them around cylinders with decreasing radius of curvature until the panel broke. The critical curvature values κ critical (α) we found are {5, 6.7, 8, 25, 40, 50} cm -1 . As a conclusion, we tag a target shape as fabricable, if none of its curvatures exceeds the critical value κ critical = 50 cm -1 (i.e. none of its curvature radii is lower than 2 cm = 1/50). The examples presented below all satisfy the maximal curvature constraint, except the heart with max(κ) = 56 cm -1 , see Table 1 and discussions below. This fabricability test can be done in real-time and provides immediate feedback to the designer.

Once the algorithm converged to a solution (x, k b ) with a cut profile α it can further be checked if the curvatures of the individual panel elements (unit cells) with same cut parameters α satisfy κ(α) < κ critical (α). Here again, we observed with all examples, that the curvatures of resulting curves never exceed those of the target, in particular when the regularization term is switched on in the objective function (see pentagon example in Fig. 9).

With the heart shape in Fig. 9, the method converges to a minimum energy curve and optimal stiffnesses within the prescribed bounds, even though the distance to target is not close to zero (1.2% error). Note our examples with an extra load extend the set of equilibrium curves considered in [START_REF] Hafner | The Design Space of Plane ElasticCurves[END_REF]. In particular the heart, where the line spanned by the 2 inflection points (dotted line) does intersect the equilibrium curve in 2 further points (see Fig. 5 in [START_REF] Hafner | The Design Space of Plane ElasticCurves[END_REF]).

Choice of external load. The feature offered by our system to add external forces to a point on the structure extends the shape space of feasible shapes, see for example a mass added at the end of the bass clefs or at the mid point of the heart. In order to choose a mass which is not too large (danger of rupture) and not too small (difficulty to reach the target), we checked by forward simulation that the z-coordinate of the target loaded point is between the z-coordinate of this point if we were using the solid material and the z-coordinate of this point if we were using the softest material.

Fabricated prototypes

To further validate our approach, we constructed some physical prototypes. We use a Trotec Speedy 400 laser cutting machine, which allows for cutting sheets with maximal length of 1 m and maximal thickness of 6-10 mm (for MDF). The four bending-active structures we realized are made of MDF panels of 100 cm length and 6 mm thick. Panels with these dimensions withstand a large interval of cut pattern sizes without falling apart even when a reasonable amount of mass is added to its own mass.

Our idea was to realize lamp stands of original curved shapes. The complete set of data and obtained results are shown in Fig. 10 The bass clef and wall lamp examples are more challenging, because of the floating end point. We observed, that achievement of the target is sensitive to a proper choice of the extra mass. In Table 1, wee see that for both examples, the algorithm converged and computed min(k b ) = 0.19 J, which corresponds to large αvalues, i.e. to coarse cut patterns.

The cut-layouts are shown in Fig. 10. The solid material is blue, the void material, i.e. the cuts, are white. The variations of the stiffness profiles (second row in Fig. 10) and consequently the variations of the geometric parameters (third row) are clearly visible. Note, that the two clamped examples have a solid part at both extremities for fabrication purposes.

The photographs in Fig. 10-bottom finally show the fabricated structures. The viewpoint is here orthogonal to the plane of the silhouettes in order to enable measurements for validation purposes (see details in Sect. 6.3). We constructed supports for all clamped points to enforce the prescribed tangent. More aesthetic photographs of the final lamp designs are shown in Fig. 11.

Validation

A direct comparison to a 3D FEM simulation of the entire resulting perforated structure of 1 m length with all its fine geometric details not larger than 3 mm is not possible since it would exceed available computational resources (memory and CPU) of a standard computer. In particular, we would have to simulate a structure of at least 30 million elements to use the same resolution as in the homogenization step, which is a minimum to have

Conclusion

We presented a computational method to design bendingactive structures whose silhouette at equilibrium matched a desired target curve. To this end, we programmed the curvature of the structure thanks to laser-cut rectangular patterns whose sizes were automatically determined using a two-scale algorithm. Our approach, based on the homogenization of the bending properties of the cut MDF wood allowed us to directly play with macroscale stiffnesses, which in turn allowed us to dramatically reduce the number of design variables and solve our inverse problem using an off-the-shelf optimizer.

While very effective, as demonstrated by our various examples, our tool still presents certain limitations that open the door to exciting future work. For example, we currently fix and keep constant the width of the strips, while the width profile could be given as input by the user to play, e.g., with aesthetics. We also limit the optimization to a single parameter of our cut pattern.

Exploring the others, and more generally other parametric families of cuts [START_REF] Zarrinmehr | Interlocked archimedean spirals for conversion of planar rigid panels into locally flexible panels with stiffness control[END_REF], to improve, e.g., the strength or to play with the appearance of the final structure, in addition to its shape, would be another interesting extension. Finally, it would also be very appealing to generalize our work to 3D, in particular to address the design of 3D developables or surfaces exhibiting low Gaussian curvature, that seem a priori difficult to obtain with methods optimizing ribbons' widths or current state-of-the-art approaches specifically tailored to the design of doubly-curved deployables.

Figure 2 :

 2 Figure 2: A bending-active structure of longitudinal dimension L struct and transversal dimension W struct . Dimensions and offsets of the rectangular cuts are shown in blue. The periodic unit cell is highlighted by a dashed rectangle. The main parameters in our design are the adimensional ratios α = lc Lc and β = wc Wc .

Figure 3 :

 3 Figure 3: The bending transformation R θ applied to the RVE.

Figure 5 :

 5 Figure 5: Bending homogenization. The left figure shows the cylindrical deformation R θ (X) of three RVEs. The right figure shows in orange the central RVE deformed periodically by the mesoscopic displacement fluctuations ũ.

Figure 6 :

 6 Figure 6: Left: Potential energy of the RVE with geometric parameters α ∈ {0.5, 0.6, 0.7}, β = 0.6, for bending angle between 0 and 0.2 radians. Middle: stiffness according to α and β. Right: stiffness according to α with β = 0.6 (logarithmic scale to enhanced small values of stiffnesses).

Figure 7 :

 7 Figure 7: Optimized stiffness values computed w and w/o regularization.

Algorithm 1

 1 iteration on mass adjustment Require: function Mapping : (k b ) → (α, β) k b bending stiffnesses, (α, β) geometric parameters InverseModeling : (m, T ) → (k b ) m : list of nodal masses, T : (target shape, boundary conditions) Mass : (α, β) → m(1 -αβ)

Figure 8 :

 8 Figure 8: Optimized stiffness (left) computed with Eq. (10) and the corresponding α cutting parameters cutting parameters for 320 unit cells of width 3 mm (right) taken from the light bulb structure shown in Fig.10(a).

Figure 9 :

 9 Figure 9: Computed bending-active structures (left) and stiffness profiles (right) generated from input shapes and boundary conditions. Input target shape (black), computed bending-active silhouette curve (orange). Boundary conditions (if applied) include fixed points (red dot), fixed tangent vectors (red arrow), external load (blue arrow) of 0.5 kg for the heart example. The ellipse and the pentagon (lines 1,2) are computed without taking gravity into account. The pentagon's stiffnesses have been computed with and without regularization. The bass clefs (line 5) result from 4 different target curves and an external load of 150 g.

  . The designs from left to right represent (a) a light bulb, (b) an arch, (c) a desk lamp in form of a bass clef, and (d) a wall lamp. The four target shapes and the BC are shown in Fig. 10top. The first 2 examples have symmetric BC at both extremities and a symmetric target shape. With the examples (c) and (d) we created the more challenging case where only one extremity is fixed and the other is only constrained by an extra load of 150 g.

  ), fixed tangents (t), a mass added at one point. Maximal curvature of target shape. Minimal and maximal optimized bending stiffnesses k b as output of the inverse design optimization. #It.mass: number of iterations for the iterative mass adjustment algorithm. #It.: number of iterations and timings in seconds for the inverse design problem to converge. Pointwise Euclidean avg/max distances in mm between target and solution. Pointwise Euclidean avg/max distances in mm between target and fabricated structure. All input curves have a length of about 1 m.

Figure 11 :

 11 Figure 11: Laser-cut bending-active structures made of MDF panels.

Table 1 :

 1 Name Boundary κ max Grvty. k b k b #It. #It. inv. IPOPT Dist.(mm) Dist.(mm) conditions target min max mass design inv. design target-sol. fabr.-sol. Statistics of all models. Boundary conditions include fixed points (p

	avg / max	avg / max

Preprint to appear in Journal of Computer-Aided Design 2022

Acknowledgment

This work was supported by the European Unions Horizon 2020 programme under grant agreement No 862025 (FETOpen ADAM 2 ). The authors are grateful to Stan Borkowsky and the Amiqual4Home platform for giving us access to their laser-cut machine.

 1.

sufficiently accurate results. Instead, we performed a visual and quantitative comparison with four fabricated structures. A visual comparison is given in Fig. 10-bottom, where we superposed by hand the optimized silhouette curve x (output of Eq. ( 10)) with the back-side border of the photographed structure. The backside curve is visible thanks to the viewpoint of the camera. All fabricated structures fit almost perfectly. In order to validate quantitatively both the forward simulation and the mechanical models, we computed the average/max Eu-clidean distance between the silhouette curve of the fabricated models digitized from the photo and the computed solution (scale to same length, equidistant sampling, calibration). The resulting errors are remarkably small, less than 10mm max for one meter length structures, giving a maximum relative error of 10 -2 , see the last two columns in Table 1.


Claire Lauvernet

Céline Helbert

Zhu Xujia

Bruno Sudret

Metamodeling approaches to help designing vegetative filter strips and improve the water quality

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1

INRAE, RIVERLY, Villeurbanne, France (claire.lauvernet@inrae.fr) 2 Univ. Lyon, UMR CNRS 5208 Ecole Centrale de Lyon ICJ, France (celine.helbert@ec-lyon.fr) 3 ETH Chair of Risk, Safety and Uncertainty Quantificatio Zurich, Switzerland Significant amounts of pollutant are measured in surface water, their presence due in part to the use of pesticides in agriculture. One solution to limit pesticide transfer by surface runoff is to implement vegetative filter strips (VFS) along rivers. The sizing of these strips is a major issue, with influencing factors that include local conditions (climate, soil, etc.). The BUVARD modeling toolkit was developed to design VFSs throughout France according to these properties. This toolkit includes the numerical model VFSMOD, which quantifies dynamic effects of VFS site-specific pesticide mitigation efficiency. However, the toolkit is quite complex to use with many input uncertain parameters (quantitative -such as the slope, the Curve Number -or qualitative -such as the soil type or the rainfall event), making it not easy to use for risk management.

In this study, a metamodeling (or model dimension reduction) approach is proposed to ease the use of BUVARD and to help users design VFSs that are adapted to specific contexts. Different reduced models, or surrogates, are compared, based on Bayesian learning approaches or not: Polynomial Chaos Expansions, Mixed-kriging, and Deep-GP. Mixed-kriging is a kriging method that was implemented with a covariance kernel for a mixture of qualitative and quantitative inputs. Kriging and Deep-GP are built by couple of modalities and PCE and Mixed-kriging are built considering mixed quantitative and qualitative variables. As a last step, Finally, we perform a global sensitivity analysis with the help of the two surrogate models with the best accuracy. The results show that they give the same ranking of the importance of the input parameters.

The metamodel is a simple way to provide a relevant first guess to help design the pollution reduction device. In addition, the surrogate model is a relevant uncertainty tool, to visualize the impact that lack of knowledge of some parameters of filter efficiency can have when performing risk analysis and management.

Powered by TCPDF (www.tcpdf.org)