A Dehbaoui
email: amine.dehbaoui@cea.fr

J M Dutertre

B Robisson

P Orsatelli

P Maurine

A Tria

Injection of transient faults using electromagnetic pulses Practical results on a cryptographic system

This article considers the use of magnetic pulses to inject transient faults into the calculations of a RISC micro-controller running the AES algorithm. A magnetic coil is used to generate the pulses. It induces computational faults without any physical contact with the device. The injected faults are proved to be constant (i.e. data independent) under certain experimental conditions. This behaviour, combined with the ability to choose the faulted bytes thanks to timing accuracy in the fault injection process, makes it possible to implement most of the state-ofthe-art fault attack schemes.

Introduction

Since the electronic devices that implement cryptography (such as "smart cards") are key components to secure communications, they are subjected to 'attacks'. Among them, fault attacks are considered as very powerful. They consist first in modifying the behavior of the chip with dedicated experimental setups and then recovering the secret information by using cryptanalysis techniques such as Differential Fault Analysis [START_REF] Biham | Differential fault analysis of secret key cryptosystems[END_REF][START_REF] Piret | A differential fault attack technique against spn structures, with application to the aes and khazad[END_REF][START_REF] Giraud | Dfa on aes[END_REF][START_REF] Roche | Combined fault and side-channel attack on protected implementations of aes[END_REF] (or DFA), safe-error [START_REF] Yen | Checking before output may not be enough against fault-based cryptanalysis[END_REF], fault sensibility analysis [START_REF] Li | Fault sensitivity analysis[END_REF], collision [START_REF] Blömer | Fault based collision attacks on aes[END_REF], round reduction [START_REF] Choukri | Round reduction using faults[END_REF], etc.

Various experimental setups are commonly used to modify the behavior of the chip. Under-powering a device or over-clocking it [START_REF] Selmane | Practical setup time violation attacks on aes[END_REF][START_REF] Barenghi | Low voltage fault attacks on the rsa cryptosystem[END_REF] leads to setup time violations resulting in the injection of errors. Transient deviations under nominal values of the power supply or of the clock period (a.k.a. power or clock glitches) [START_REF] Hamid | The sorcerer's apprentice guide to fault attacks 94[END_REF][START_REF] Agoyan | When clocks fail: On critical paths and clock faults[END_REF][START_REF] Li | Fault sensitivity analysis[END_REF], also offer similar fault injection abilities, besides providing higher timing accuracy. Even a rise in the temperature of operation may be used to inject faults [START_REF] Hamid | The sorcerer's apprentice guide to fault attacks 94[END_REF]. Another important means of fault injection is the use of optical radiations: intense white light (e.g. from a flash bulb) or a laser beam [START_REF] Skorobogatov | Optical fault induction attacks[END_REF]. The latter is widely used while assessing the security of cryptographic systems against fault attacks [START_REF] Van Woudenberg | Practical optical fault injection on secure microcontrollers[END_REF] for certification purposes. It may offer the ability to inject fault affecting a byte or even a single bit of sensitive data [START_REF] Agoyan | How to flip a bit?[END_REF]. At last, a several scientific articles in the field of hardware security refer to the use of electromagnetic (EM) field in order to induce faults into cryptographic systems.

In this article, an experimental setup, dedicated to such EM fault injection is described. The effect of such injection on a not protected implementation of a cryptographic algorithm is proposed. The results obtained lead to a fault model which appears to be different from those classically used in the literature or from those obtained with other injection methods.

This article is organized as follows. A short review of the state-of-the-art EM active attacks is given in section 2. We describe our experimental setup and provide a short review of the effect of the magnetic pulses in section 4. Whereas the section 5 reports the results of the implementation of the method described in [START_REF] Piret | A differential fault attack technique against spn structures, with application to the aes and khazad[END_REF] and the section 6 reports the nature of the injected faults. As a conclusion, section 7 summarizes our findings and draws some prospects.

EM Active Attacks : State-of-the-Art

As already mentioned in the introduction, the EM medium may be used to conduct active attacks. Two kinds of near-field EM perturbations are usually considered: transient pulses and harmonic emissions.

Regarding transient pulses, two articles report actual results of successful fault injection. J.-J. Quisquater and D. Samyde [START_REF] Quisquater | Eddy current for Magnetic Analysis with Active Sensor[END_REF] in 2002 described the use of an active probe to apply an intense and transient magnetic field on a microprocessor. This results in faulting RAM and EPROM memory cells. Moreover, faults on the device's bus were also obtained. More recently, in 2007, J.-M. Schmidt and M. Hutter reported the use of a spark generator to fault a CRT-based RSA algorithm running on an 8-bit micro-controller [START_REF] Schmidt | Optical and em fault-attacks on crt-based rsa: Concrete results[END_REF]. The injected fault leads to a successful attack as it permits them to factorize the RSA modulus.

Concerning harmonic emissions, Alaeldine et al. studied the electromagnetic compatibility (or EMC) of integrated circuits (or IC) to near-field injection with frequencies up to 1 GHz [START_REF] Alaeldine | Assessment of the Immunity of Unshielded Multi-Core Integrated Circuits to Near-Field Injection[END_REF]. They investigated the effects of both electric and magnetic fields along the x-, y-, and z-axes. The immunity criterion they used was a deviation of the amplitude or a jitter of the output signals over 20% and 10% respectively. Their test circuits were found sensitive to both magnetic and, to a greater extent, electric fields. Recently, Poucheret et al. [START_REF] Poucheret | Local and direct em injection of power into cmos integrated circuits[END_REF] considered the effect of an 1 GHz electric field applied to an IC with an embedded ring oscillator (or RO). The main component of that electric field was the transverse one (i.e. parallel to the surface of the chip). The perturbation impacted the output frequency of the RO. Monitoring the effect of that perturbation enabled them to draw a cartography of the sensitive areas of the chip. A cross examination between the layout of the device and the cartography demonstrates that the coupling between the injection probe and the circuit lies mainly in the powerground network (PGN). Besides, this kind of direct power injection through the PGN permits them to increase the output frequency by up to 50%. This result reinforces the threats against the security of secure devices using RO-based true random pattern generator as a source of entropy [START_REF] Markettos | The Frequency Injection Attack on Ring-Oscillator-Based True Random Number Generators[END_REF].

In a complementary way with this previous works, we plan 1. to focus the fault injection on a small part of the targeted device, 2. to focus a very precise instant, 3. to set the energy of the fault injection.

The first and the third features are necessary, for example, to override countermeasures that detect glitches on the power supply or to avoid disturbing the circuit in its entirety (and thus to reduce the risk of crash or system denial).

All the three features are also important to impact only a particular internal computation of the cryptographic algorithm and thus to be able to apply a as wide as possible number of cryptanalysis techniques.

Those features lead us to mount an electromagnetic injection bench able to generate localized short transient pulses with tunable energy.

Experimental setup

In this section, the electromagnetic injection bench and the circuit under attack are described.

Magnetic pulse injection bench

The fault injection bench is built of: a control PC, the targeted device (i.e. the victim), a motorized stage, a smart card reader, a pulse generator, and a magnetic probe (see Figure 1). The target (described in section 3.2) is fixed on the x-y-z motorized stage. Every element of the bench is controlled by the control PC, and the communication with the target is established through the smart card reader.

The pulse generator is used to deliver voltage pulses (with low jitter) to the magnetic coil. It has a constant rise and fall transition time of 5ns. The amplitude range (respectively the width) of the pulses extends from 1V to 100V (respectively from 10ns to 100ns). We use a magnetic probe (composed of 10 turns with a diameter of 1mm Figure 2) in order to disturb only a small part of the targeted device. This spatial accuracy is possible thanks to the directive magnetic flux in the near field domain. This magnetic probe incorporates a ferromagnetic core, which acts as a flux concentrator. The resulting magnetic flux is then enhanced (according to the value of the material's permeability µ r ) by comparison with an air coil [START_REF] Tumanski | Induction coil sensors, a review[END_REF]. However, this enhancement is obtained at the expense of sensor linearity. In our case, this parameter is not impeding.

Target

We use a smart card emulation board built in the laboratory. It is composed of an 8-bit 0.35µm AVR micro-controller with integrated 128KB flash program memory, 4KB EEPROM and 4KB SRAM. This none secure micro-controller has an operating voltage in the range 4.5-5.5V and runs at 3.59 MHz frequency. A smart-card-like OS (called SOSSE [START_REF] Bruestle | sosse -simple operating system for smartcard education[END_REF]) is used for communication purposes. A software implementation of the AES encryption algorithm [START_REF]NIST: Announcing the Advanced Encryption Standard (AES)[END_REF] using a key size of 128 bits (hereafter the AES-128) is embedded. The AES-128 ciphers a 128-bit long data (called plaintext or input ) by using a 128-bit long key to obtain a 128-bit long data (called ciphertext or output). The ciphering consists first in transforming the input data into a two-dimensional array of bytes, called the "State", as illustrated in Figure 3. This array, denoted by the symbol s, is composed of four rows and four lines. Each individual byte s r,c has two indices, with its row number r in the range 0 ≤ r < 4 and its column number c in the range 0 ≤ c < 4.

Next, after a preliminary XOR between the input and the key, the AES-128 executes 10 times a function (called round ) which operates on the "State". The operations used during these rounds (which are identical except the last one) are the following: At the end of the ciphering operations, the "State" is copied to the output as depicted in Figure 3.

Preliminary experiments

In this section, two experiments are conducted on the target described in 3.2. During both of these experiments, the probe was located at a fixed position above the surface of the circuit (500 µm). This corresponds to an area containing the Arithmetic and Logic Unit (ALU) and the SRAM.

Experiment 1

First, we investigated the relationship between the parameters of the voltage pulses delivered by the pulse generator (amplitude and width) and their effect on the target. Note that both the rise and fall times of the pulses were equal to 5ns. We measured the induced transient voltage between the power supply and ground pins of the unpowered chip. The pulse duration and amplitude were set to values ranging from 10ns to 100ns and from 1V to 100V respectively. By doing so, we intended to identify the maximum coupling between the injection probe and the victim without recording the noise caused by the internal operations of the chip. Figure 4 reports the maximal amplitude of the induced voltage as a function the pulse voltage amplitude.

The maximum voltage of the perturbation increases linearly with the pulse amplitude. At 100V, a perturbation between 170mV and 205mV is obtained, depending on the pulse width. The maximum coupling is obtained for the smallest durations. The explanation lies in the closeness of the rise and fall edges of the pulses (the coil acts as a differentiator). Figure 5 presents the shape of the induced perturbation resulting from a pulse with a 20ns width and a 50V amplitude.

These voltage variations may seem quite small to be able to induce faults into the device computations. However, because these measures were done outside the chip package, a large part of the perturbation is filtered out.

Experiment 2

Second, we targeted the execution of the last round of the AES (i.e. the 10th round for the AES-128). We intended to identify the type of faults induced by the EM pulses. Figure 6 shows the power consumption of the target during the In this experiment, the trigger signal is activated at the beginning of the 10th round for ease of synchronization. We swept the instant of the EM pulse injection from this beginning to the end of the 10th round (90µs) by steps of 10ns. For each of these steps, a ciphering was executed with and without fault injection. The results of these two computations were compared and the faulted byte determined. Figure 7 indicates the timing of the EM perturbations resulting in faulting every bytes of the AES state. As reported, we were able (a) to create a fault on only one byte of the computation (thanks to the 8-bit architecture) and (b) to select the faulty byte by choosing the right instant of injection. Then, we tested the correctness of this ability to inject byte-wise faults at every byte location of the AES state during the 9th round (see section 5).

Attack

In this section, we describe the attack path that we have chosen and the experiments that have been done in order to validate this attack path on the target described in 3.2.

Chosen DFA scheme

The experiments described in section 4.2 show that our injection method is:

able to create a fault on only one byte of the computation, able to select the faulty byte by choosing the right instant of injection. These two properties seems perfectly adequate to mount the first attack described in [START_REF] Piret | A differential fault attack technique against spn structures, with application to the aes and khazad[END_REF] (further called the P&Q method). According to this scheme, the attacker has to fault one byte of the AES state between the start of the 9th round and the MixColumns transformation. Due to the diffusion properties of the AES, a ciphertext that contain 4 erroneous bytes is then outputted. These 4 errors leak information on 4 bytes of the 10th round key. Depending on the location of the faulted byte inside the AES state of the 9th round, 4 different error patterns are encountered:

-Error pattern "0" (hereafter denoted p 0 ) is obtained for a byte wise fault injected among bytes s 0,0 , s 1,0 , s 2,0 , or s 3,0 either on state "round [START_REF] Derbez | Meet-in-the-middle and impossible differential fault analysis on aes[END_REF].start" or on state "round [START_REF] Derbez | Meet-in-the-middle and impossible differential fault analysis on aes[END_REF].m col" resulting in faulting bytes s 0,0 , s 3,1 ,s 2,2 , s 1,3 of the ciphertext (i.e. the error pattern),

similarly error patterns p 1 , p 2 , p 3 are depicted on figure 8.

Every pattern also indicates the location of the information leakage among the bytes of the 10th round key. For every pattern, by applying the P&Q method on two pairs of correct/faulty ciphertexts, 4 bytes of the 10th round key are retrieved with a success rate of 97%. As a consequence, to discover the whole key the attack had to be completed for the 4 patterns. It was achieved successfully by changing the fault injection time (and consequently the location of the faulted bytes in the AES state). Note that a more efficient implementation of this attack consist in inducing the fault before the MixColumns transformation of the 8th round.

Fig. 8. Fault propagation depending on the location of the fault on the state "round [START_REF] Derbez | Meet-in-the-middle and impossible differential fault analysis on aes[END_REF].start".

Besides, a given error pattern may be produced for a byte wise fault induced either in the states "start", or "s box", and "s row" of the 9th round.

DFA scheme in practice

As described above, the P&Q method requires a 1-byte fault to be injected before the "MixColumns" transformation of the AES-128 9th round. Note that a trigger signal is activated during the 9th round for ease of synchronization. An error pattern among the 4 defined in 5.1 is selected. Let's consider, for example, the pattern p 0 . The instant of the EM pulse injection is swept from the beginning of the 9th round over a range of 200µs by steps of 10ns. For each of these steps, a ciphering is executed with and without fault injection. The results of these two encryptions are compared and the error pattern is determined. If the retrieved pattern is not p 0 , the injection time is incremented to the next time step. Whereas, if it corresponds to pattern p 0 , the time step is recorded: noted F [p 0 ].f irst. The P&Q method is also applied to recover the corresponding bytes of the 10th round key. The same experiment is performed for the three other error patterns. The gathered data make it possible to recover the full 10th round key. The algorithm which describes more precisely the attack is reported in Annex (Algorithm 1). In this section, we describe a method to compute some statistics about the faults induced by the EM pulses. Then a fault model is proposed.

Algorithm 1 DFA scheme in practice

Characterization method

We intended to study the effect of the handle data on the value of the injected fault. The experiment described in this section was conducted for 4 different injection times corresponding to the 4 error patterns. For the sake of clarity, we will consider in the following the error pattern p 0 . The magnetic pulse was induced at the corresponding time step F [P 0 ].f irst for 1000 encryptions made with random plaintexts (however the key was kept constant). Then statistics were built. We determined the percentage of erroneous ciphertexts obtained during these encryptions. Moreover, since we knew the key we were able to run backward the AES encryption steps from the erroneous ciphertext to retrieve the byte value of the injected faults. However, from that analysis a doubt remains on the exact step of the fault injection among the transformations of the AES. We have not investigated further. The fault value we consider is calculated at the step "start". The algorithm which describes more precisely the characterization method is reported in Annex (Algorithm 2).

Results

The results obtained by the experiments described in the algorithm 2 are summarized in the Table 9. The results show that the proposed EM injection:

enables to create a fault on only one byte of the state during round 9, enables to fault the 16 bytes of the state during round 9 independently (by selecting an adequate injection time step), induces a fault that is independent of the value of the faulted byte.

Fault model

From the results reported above, it appears that a magnetic pulse injected during the 9th round is able to force a byte of the AES state to a constant value (i.e. independent from the plaintext). However, this value differs from one byte to the In this article, a fault injection technique based on the generation of an electromagnetic field close to a circuit, has been presented. This technique enables to fault every byte of the AES state on a non protected software implementation of an AES, running on an 8-bit micro-controller. We have shown that this feature enables to mount a classical DFA scheme on our target but also should be suitable to implement any DFA method that requires the injection of the byte-wise fault. Moreover, a careful review of the nature of the faults created on micro-controller showed that they were constant (i.e. data independent) for each erroneous byte considered individually. Yet, these results were found to be reproducible with a success rate close to a hundred percent. Thus, it may be used to implement fault attacks that are based on a constant fault model. At last, in comparison with other optical fault injection means, the proposed EM injection experiments were carried out both on decapsulated (i.e. frontside opening of the package) and on untampered ICs without any noticeable differences in their responses to the perturbations. We now plan to create fault with the EM injection technique described in this paper on an hardware implementation of cryptographic algorithms and to apply this technique against protected implementations.

Fig. 1 .

 1 Fig. 1. EM pulse injection bench.

Fig. 2 .

 2 Fig. 2. Magnetic injection probe over an open device (left) and an untampered device (right).

  Fig. 3. AES Notations

Fig. 4 .

 4 Fig. 4. Amplitude of the measured perturbation versus the amplitude of the voltage transient issued from the pulse generator(given for pulse width ranging from 10ns to 100ns).

Fig. 5 .

 5 Fig. 5. Pertubation induced by a magnetic pulse on the targeted IC.

Fig. 6 .

 6 Fig. 6. Power consumption trace during magnetic fault injection.

Fig. 7 .

 7 Fig. 7. Timing map of the obtained erroneous bytes.

  Require: Input K : Fixed unknown key Require: Input S : Time steps (vector of size M = 9 * 10 3 ) Require: Input P : Types of errors (vector of size 4 and P = [p0, p1, p2, p3] Ensure: Output: Kr retrieved value of the key K Ensure: Output F : Statistics related with the 4 types of errors Ensure: With F [pi].f irst, the first time step when an error of type pi is detected (scalar) for All type pi ∈ P of error do Set s = 0 Set Stop=FALSE while While Stop==FALSE do Select a text t at random Compute the associated ciphertext c = AES(t, K) without fault injection Inject a pulse at step S[s] during the ciphering of t with K. Let c being the result of this computation. if If c ⊕ c is an error of type pi then Set Stop=TRUE Set F [pi].f irst = S[s] Perform the attack described in [15] to retrieve 4 bytes of the key Kr by injecting fault at time step S[s] during the ciphering of approximately 2 plain texts chosen randomly end if Increment s end while end for 6 Towards a fault model

Fig. 9 .

 9 Fig. 9. Results of the fault characterization step

Compute offline all the intermediate variables of AES -1 (c, K). Let Round [START_REF] Agoyan | When clocks fail: On critical paths and clock faults[END_REF].isbox(c, k) be the value of the state "isbox" during the second round of the deciphering process(equivalent of the Round [START_REF] Derbez | Meet-in-the-middle and impossible differential fault analysis on aes[END_REF].start(t, k) in the ciphering process) Compute offline all the intermediate variables of AES -1 (c , K). Let Round [START_REF] Agoyan | When clocks fail: On critical paths and clock faults[END_REF].isbox(c , k) the value of the state "isbox" during the second round of the deciphering process. Store the "reconstructed" value of the fault in F [pi].f aults[m] = Round [START_REF] Agoyan | When clocks fail: On critical paths and clock faults[END_REF].isbox(c, k) ⊕ Round [START_REF] Agoyan | When clocks fail: On critical paths and clock faults[END_REF].isbox(c , k) end for end for other. The value depends also on the injection time, as the way to fault different bytes is to sweep the instant of injection. In other words, the function which associates the value x of a byte and the faulted value when the EM injection is performed at time t i is the function f ti : x ← D with D ∈ 00 H , . . . , F F H for all x ∈ 00 H , . . . , F F H . Such a type of fault could be considered as a combination of set and reset faults which modify all the bits of the byte. Note that the value of D may be known by the attacker if, like in the side-channel attack called "template", he is able to performed previous measurements on the same circuit than the target but with the access on the value of the key.

Other possible attack schemes

In the previous Sections, we have shown that our EM fault attack setup allow us to inject a fault on one byte, with a high accuracy, and we have experimentally verified that, assuming a precise temporal trigger, the probability of the fault repeatability is close to 100%.

Template EM attack Let's consider that the attacker is able to perform measurements on a circuit and on a software implementation of the AES which are strictly identical to the target but with a known key. During a first step (called characterization step), the attacker scans the round 10 in the time domain to detect the instant of the EM injection which modifies a byte of the output (exactly as it has been done in 4.2). Because the value of the key is known, he is able to compute a fault on the state "round [START_REF] Giraud | Dfa on aes[END_REF].s row" which could have produced such an output. By doing the same for the all the bytes, he's able to associate each byte j of the output with the value of the fault (equivalent to a fault created on the state "round [START_REF] Giraud | Dfa on aes[END_REF].s row" and noted D i ) and with an instant of injection (noted t j ). During the attack step, the attacker performs the fault injection at one of those different instants t j . He verifies that the byte j is impacted at the output. If it is the case, the value of this faulted byte out j is stored. The byte j of the round key 10 is computed by using the formula D j ⊕ out j . He performs the same for the others bytes.

Combined Attack to defeat protected implementation A recent work of Roche et al. [START_REF] Roche | Combined fault and side-channel attack on protected implementations of aes[END_REF] propose a combined fault and side-channel attack, allowing to defeat protected implementations of AES. More precisely, the attack is able to defeat a masked and fault protected implementation, by injecting a fault during the penultimate key schedule, and by measuring the leakage of the faulty ciphertext. Nevertheless, they need to be able to inject a fault with a good repeatability, in other ways it means that the effect of the fault on the faulted byte can be rewritten as a XOR with a constant value. Considering this requirement, our experimental results seem showing that the good repeatability of the fault we observe could allow us to practically perform this attack.

Meet-in-the-Middle and impossible DFA on AES Other fault attack schemes seem to be practically applied in the same way. For instance, the recently published work of Derbez et al. [START_REF] Derbez | Meet-in-the-middle and impossible differential fault analysis on aes[END_REF] propose several new Differential Fault Attacks on AES, allowing to retrieve the key by injecting random faults on one byte between the 6th and 7th MixColumns ( in the case of AES128). Our practical results show that the fault model we get could allow us to perform these attacks. Considering the constant fault model, the attack proposed needs only 5 faults and its complexity in memory is reduced to 2 24 while other fault models require either 1000 or 45 faults depending on the fault model and recover the secret key with a time and memory complexity around 2 40 .


Ignacio Fuenzalida-Henriquez

Paul Oumaziz
email: paul.oumaziz@insa-toulouse.fr

Emilio Castillo-Ibarra

Jorge Hinojosa
email: jhinojosa@utalca.cl

Global-Local Non Intrusive Analysis with Robin Parameters: Application to Plastic Hardening behavior and Crack Propagation in 2D and 3D Structures

Keywords:

The global/local analysis allows to embed a specific local zone of interest with a different behaviour in a global coarse model. In this local model, fine meshes are usually used to model some structural details and potentially non-linear behaviours, such as plastic hardening and crack propagation. The standard global/local approach can be observed as a Dirichlet-Neumann iterative algorithm where a Dirichlet problem on the local model and a Neumann problem on the global one are solved successively. This paper proposes a new approach for the global/local framework as Robin parameters are considered on both local and global models to obtain more flexibility and improvement for convergence. Particularly, Robin parameters are obtained using a pre-defined strip of elements and the results are later improved by means of single-objective optimization, minimizing the number of iterations to achieve convergence. This improvement is illustrated for cracked domains and plastic hardening in 2D problems and 3D elements within a non-intrusive framework, allowing the usage of commercial finite element software along with open-source research finite element software.

Introduction

The objective of non-intrusive frameworks is to carry out advanced numerical methods as a result of efficient linear and non-linear solvers applied in the commercial software and used as "black boxes". An example of the above is that it is feasible to operate the free open-source industrial finite elements software code-aster [START_REF] Edf | Code Aster, Analysis of Structures and Thermomechanics for Studies and Research[END_REF] with the Python interface. Advanced algorithms can be designed within this Python interface and call non-linear methods of codeaster to consider complex phenomena. In addition, depending on the perspective, the non-intrusivity might be the alternative to "enrich" a global model with those that are more detailed and without modifying the global one. For example, the global model could be the product of a long industrial process, which might not be altered easily. In this regard, methods considering local models without altering the global one can be perceived as non-intrusive.

The global/local approach [START_REF] Whitcomb | Iterative global/local finite element analysis[END_REF] is a very good option for implementing the non-intrusive method [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF]. In essence, the idea is to introduce some more localized details into a global and coarse model and potentially add non-linear behaviours in specific areas without altering the global model. Summarizing, a global and a local model connected through an interface coexist. The global/local approach is about an iterative Dirichlet-Neumann algorithm and an iteration consisting of two steps: [START_REF] Aldakheel | A global-local approach for hydraulic phase-field fracture in poroelastic media[END_REF] a problem in the local model with Dirichlet boundary conditions on the interface, and (2) a problem in the global model with Neumann boundary conditions on the interface. Connections with optimized Schwarz domain decomposition methods can be observed in [START_REF] Gosselet | Non-invasive global-local coupling as a Schwarz domain decomposition method: acceleration and generalization[END_REF]. The global/local structure has been expanded to a domain decomposition method [START_REF] Duval | Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition[END_REF] and has been implemented to different types of non-linearities such as crack propagation [START_REF] Passieux | Local / global non-intrusive crack propagation simulation using a multigrid X-FEM solver[END_REF][START_REF] Duval | Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition[END_REF][START_REF] Noii | An adaptive global-local approach for phase-field modeling of anisotropic brittle fracture[END_REF], structural joints and assemblies [START_REF] Guguin | On the computation of plate assemblies using realistic 3D joint model: a non-intrusive approach[END_REF], local plasticity [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF] and cycling viscoelastic behaviour [START_REF] Blanchard | Space/time global/local noninvasive coupling strategy: Application to viscoplastic structures[END_REF]. It is possible to find other studies regarding different applications and improvements for non-intrusive frameworks and the global/local method in [START_REF] Liu | A non-intrusive global/local algorithm with non-matching interface: Derivation and numerical validation[END_REF][START_REF] Gendre | A two-scale approximation of the Schur complement and its use for non-intrusive coupling[END_REF][START_REF] Gerasimov | A non-intrusive global/local approach applied to phasefield modeling of brittle fracture[END_REF][START_REF] Gosselet | Non-invasive global-local coupling as a Schwarz domain decomposition method: acceleration and generalization[END_REF][START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: Geometric details, local refinement, inclusion, fracture[END_REF][START_REF] Aldakheel | A global-local approach for hydraulic phase-field fracture in poroelastic media[END_REF].

Considering the St. Venant's principle, the interface between the local and global models should be nowhere near the local details to avoid the inaccuracy of the Dirichlet boundary conditions from the linear global model close to the area of potential non-linearity. This problem could be solved by using Robin parameters on the interface. The first proposal to apply Robin parameters on the interface for a local plastic model can be found in [START_REF] Gendre | A two-scale approximation of the Schur complement and its use for non-intrusive coupling[END_REF]. Their idea is to design a "quasi-optimal" Robin parameter avoiding the full computation of the Schur complement of the complementary models. They prefer two-scale Robin parameters that include both short and long-range effects. Very recent works [START_REF] Noii | An adaptive global-local approach for phase-field modeling of anisotropic brittle fracture[END_REF][START_REF] Aldakheel | A global-local approach for hydraulic phase-field fracture in poroelastic media[END_REF] applied a global/local approach with Robin conditions for phase field modelling for mechanical and hydraulic fractures in porous media. They chose the Schur complement of the complementary domain and observed very good performance, improving drastically the convergence of the global/local algorithm.

Finally, a formal mathematical approach to Mixed Domain Decomposition Methods and the application for non-intrusive analysis with Robin Parameters can be found in [START_REF] Negrello | A new impedance accounting for short-and long-range effects in mixed substructured formulations of nonlinear problems[END_REF][START_REF] Oumaziz | Analysis, improvement and limits of the multiscale Latin method[END_REF][START_REF] Noii | An adaptive global-local approach for phase-field modeling of anisotropic brittle fracture[END_REF]. This paper focuses on the global/local method with a non-intrusive implementation of Robin parameters on the interface. These Robin parameters are developed with a strip of elements linked to the interface similarly as [START_REF] Gendre | A two-scale approximation of the Schur complement and its use for non-intrusive coupling[END_REF]. The choice of these Robin parameters is crucial for the convergence of the method. The optimal choice is classically known but with an unaffordable cost as it would involve the computation of a Schur complement of a large model. Therefore, approximations are sought and we propose to carry out an optimization of the parameters to obtain a better convergence of the mixed method for the global/local analysis and study the impact of the Robin parameter as well.

The article is structured as follows: Section 2 formulates the global/local method as a Dirichlet-Neumann algorithm; Section 3 pursues the global/local method with the Robin parameters; Section 4 and subsequent sections illustrate the improvements of our implementation with numerical examples and different non-linear problems in 2D and 3D models, as well as the optimization of the Robin parameters.

Global/Local Problem Formulation

Reference problem

A mechanical model of a structure determined on a domain Ω R is considered. This domain consists of two nonoverlapping domains Ω C , the complementary domain, and Ω L , the local domain. Ω C considers elastic linear isotropic assumptions with small perturbations, meanwhile in Ω L , the mechanical model may be non-linear, such as a crack propagation model or a plastic material. This local model allows embedding localized richer contents in the structure's simple global model. Γ is noted as the interface between the complementary and local models.

Defining the admissible space of displacements as:

V (Ω R ) = {v ∈ H 1 (Ω R ) v = u d on ∂ u Ω R } (1)
.

The mechanical problem is equivalent to:

Find u ∈ V (Ω R ) , a R (u, v) = l R (v), ∀v ∈ V 0 (Ω R ) (2)
with a R the bi-linear form representing the structure's equilibrium and l R the linear form symbolizing the Neumann boundary conditions. The Dirichlet boundary conditions are taken into account in the affine space V (Ω R ).

A discretization of standard finite elements with Lagrange shape functions is considered in order to obtain discrete models for which a conforming mesh in the interface between the complementary and local models is adopted.

A full elastic linear isotropic material over the whole domain Ω R = Ω C ∪ Ω L is considered to detail the methods and the equation. Under this type of hypothesis, the discrete problem becomes as follows:

Find u R solution of K R u R = f R d (3)
with K R the stiffness matrix, u R the discrete unknown of displacements and f R d the right-hand side conforming to the boundary conditions. The null Dirichlet conditions are expected to be eliminated. This reference problem can be seen as the coupling between the complementary model in Ω C and the local model in Ω L , as presented in Figure 1. We note that K C and K L are the stiffness matrices of the complementary and local models respectively, u C and u L are the displacements in Ω C and Ω L , f C d and f L d are the load vectors in Ω C and Ω L . The perfect coupling of the two models to impose the continuity of the displacement is enforced with a Lagrange multiplier λ into the following Lagrangian:

L (u C , u L , λ) = 1 2 u T C K C u C + 1 2 u T L K L u L -f C d T u C -f L d T u L + λ T (C C u C -C L u L ) (4)
Where C C and C L are coupling operators. For conforming meshes, these are trace operators that extract displacements from Ω C and Ω L on the interface Γ. They are sparse matrices filled with number ones in the indexes associated to the interface degrees of freedom.

The minimization of the Lagrangian leads to:

• Equilibrium of the coarse model:

∂L ∂u C = K C u C + C T C λ -f C d = 0 (5)
• Equilibrium of the local model:

∂L ∂u L = K L u L -C T L λ -f L d = 0 (6)
• Continuity of the interface displacements:

C C u C -C L u L = 0 (7)
Therefore we obtain the standard coupling system:

Find (u C , u L , λ) ,   K C 0 -C T C 0 K L C T L -C C C L 0     u C u L λ   =   f C d f L d 0   (8)
This problem can be rewritten as two sub-problems: one over Ω C with Neumann conditions of the interface Γ and the other problem over Ω L with Dirichlet conditions on the interface Γ:

Find u C , K C u C = f C d + C T C λ Find (u L , λ) , K L C T L C L 0 u L λ = f L d -C C u C (9)
A Dirichlet-Neumann iterative algorithm is used to solve this problem. The iteration consists in the following two steps:

1. Knowing a (u n L , λ n ) solution in Ω L to solve a Neumann problem on Ω C : Find u n+1 C , K C u n+1 C = f C d + C T C λ n (10)
2. Knowing a u n+1 C solution in Ω C to solve a Dirichlet problem on Ω L :

Find u n+1 L , λ n+1 , K L C T L C L 0 u n+1 L λ n+1 = f L d C C u n+1 C ( 11
)
and obtain λ n+1 , the opposite of the reaction forces λ n+1 L on the interface of the local model.

Remark 1. This type of algorithm could be considered with multiple local models and conceived as a non-overlapping domain decomposition [START_REF] Duval | Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition[END_REF].

Primal-Dual Global/Local Analysis

The subtlety of the global/local method is to slightly modify the problem in the complementary domain with the aim of not solving it in Ω C , but in a "global" domain Ω G , which corresponds to the union of Ω C and an auxiliary domain Ω A . Globally, this auxiliary domain corresponds geometrically to the local domain Ω L but coarsely meshed as the complementary model. As shown in Figure 2, the local model is substituted by an auxiliary model to obtain a global model. This can also be interpreted differently. Starting from a global model with a coarse mesh and an elastic linear isotropic material, this model can be enriched with localized details such as structural details (holes or a particular geometry) or complex non-linear behaviour (plasticity or cracks). This is especially suitable in an industrial design process where the global model may not be easily changed to incorporate localized details.

The complementary problem to be solved is defined as:

Find u n+1 C , K C u n+1 C = f C d + C T C λ n (12)
The contribution of the auxiliary model is added to the two sides of the equation where the exponent n + 1 and n are omitted in order to simplify the notations:

K C u C + K A u A = f C d + f A d f G d + C T G λ + K A u A -f A d C T G λ A (13)
where the notation X means that the variable or operator specified in Ω X is extended with zeros to the domain Ω G for the rest of the degrees of freedom. In addition, f . For more details see [START_REF] Gosselet | Non-invasive global-local coupling as a Schwarz domain decomposition method: acceleration and generalization[END_REF]. Therefore, as u C = u A on the interface Γ, the equation becomes as follows:

K G u n+1 G =f G d + C T G λ n + C T G λ n A =f G d + C T G P n (14)
with u n+1 G the displacement specified on the global model. The reaction forces λ n A are computed with displacement u n G from the previous iteration of the algorithm. P n = λ n A -λ n L represents the correction of forces that will be applied to the global model in order to consider the local model. It can also be observed that P n verifies the following equations assuming exact solutions for the global and local problems [START_REF] Gosselet | Non-invasive global-local coupling as a Schwarz domain decomposition method: acceleration and generalization[END_REF]:

P n = λ n+1 A + λ n+1 C P n+1 = P n + r n+1 with r n+1 = -λ n+1 L + λ n+1 C ( 15
)
Under this notation, r n+1 represents the disequilibrium of forces between the complementary model and the local model. When convergence is fulfilled, r n+1 = 0 and the forces are balanced: the solution on Ω C ∩ Ω L reaches the solution of the reference model.

The error indicator for the convergence of the method is the residual norm for each iteration normalized regarding the first residual obtained from global/local iterations, as presented in Eq. ( 16).

η = r n+1 2 r 0 2 (16)

Global/Local algorithm

The iterations described above lead to the algorithm 1 of the global/local method. Solve the global problem:

Find u n+1 G , K G u n+1 G = f G d + C T G P n 6:
Solve the auxiliary problem to compute λ n+1 A :

Find u A , λ n+1 A , K A C T A C A 0 u A -λ n+1 A = f A d C G u n+1 G 7:
Solve the local problem:

Find u n+1 L , λ n+1 , K L C T L C L 0 u n+1 L λ n+1 = f L d C G u n+1 G 8:
Compute the new P n+1 :

P n+1 = λ n+1 A + λ n+1 9:
Calculate the residual of the iteration r n+1 = P n+1 -P n 10:

Update the new P n+1 including relaxation factor µ:

P n+1 = µP n+1 + (1 -µ)P n 11:
Estimate the error of the iteration step as η = r n+1 2 / r 0 Neumann fixed-point algorithm. However, in the same way as Schwarz domain decomposition methods [START_REF] Gander | Schwarz methods over the course of time[END_REF], Robin conditions can be introduced to enhance the convergence and obtain more flexibility. In fact, the kinematic compatibility between the complementary and local model is broken and replaced by Robin conditions written on the interface Γ. The first studies in the global/local method with Robin conditions are introduced in [START_REF] Gendre | A two-scale approximation of the Schur complement and its use for non-intrusive coupling[END_REF] and reminded in [START_REF] Gosselet | Non-invasive global-local coupling as a Schwarz domain decomposition method: acceleration and generalization[END_REF]. The proposal in this paper is to go further and derive a complete mixed global/local method.

Derivation of Global/Local Analysis with Robin Parameters

The mechanical problem is presented differently in this section as it follows a mixed domain decomposition framework. Instead of using a Lagrange's multiplier specified on the interface Γ, reaction forces λ L and λ C are considered as full unknowns on the interface. The conditions to be enforced at the interface are provided by writing the problem as a domain decomposition:

λ L + λ C = 0 Equilibrium of forces C L u L -C C u C = 0 Continuity of displacements ( 17
)
These two equations interpreting the interface behaviour are written with Robin conditions:

λ L + λ C -k C (C L u L -C C u C ) = 0 λ L + λ C + k L (C L u L -C C u C ) = 0 (18)
k C and k L are Robin parameters, which are stiffness operators. Symmetric definite positive operators are selected to ensure the equivalence with Eq. ( 17) of the interface. These Robin conditions are determined and written on the interfaces, connecting reaction forces λ C and λ L to the interface displacements C L u L and C C u C . Thus, the new system to be solved is:

K C u C = f C d + C T C λ C λ L + λ C -k C (C L u L -C C u C ) = 0 K L u L = f L d + C T L λ L λ L + λ C + k L (C L u L -C C u C ) = 0 (19)
Similarly, as the standard global/local method, the problem in the complementary model is generalized to the global model. Therefore, the first equation of the system becomes as follows:

K G u G = f G d + C T G λ C + C T G λ A (20)
In addition, as u G and u C are equal on the interface Γ, Robin conditions can also be written as:

λ L + λ C -k G (C L u L -C G u G ) = 0 λ L + λ C + k L (C L u L -C G u G ) = 0 (21) with k G = k C .
Remark 2. λ C is the reaction force of a problem in the complementary domain. Therefore, it is preferable not to use the notation λ G in order to avoid a confusion on the nature of λ G .

The first equation of the Robin conditions allows to obtain the expression of λ C to incorporate it into the equilibrium of the global model, which led to:

K G u G = f G d + C T G λ A -C T G (λ L -k G C L u L -k G C G u G ) (22)
The second equation of the Robin conditions allows to obtain the expression of λ L to incorporate it into the equilibrium of the local model, which led to the following equation:

K L + C T L k L C L u L = f L d + C T L (k L C G u G -λ C ) (23)
Similarly as the global/local method, a fixed-point algorithm can be derived where an iteration consists of the following steps:

1. Global problem: knowing (u n L , λ n L , λ n A ), find u n+1 G solution of:

K G u n+1 G = f G d + C T G (λ n A -(λ n L -k G C L u n L -k G C G u n G )) P n (24)
2. Auxiliary problems: computing the reaction forces on Ω A and Ω C on the interface Γ:

λ n+1 A = C A K A u G|Ω A -f A d λ n+1 C = C C K C u G|Ω C -f C d (25) 3. Local problem: knowing u n+1 G , λ n+1 C
, find u n+1 L solution of:

K L + C T L k L C L u n+1 L = f L d + C T L k L C G u n+1 G -λ n+1 C ( 26
)
Remark 3. As to preserve the global model without any change, the contribution of u G is considered on the righthand side of the previous equation with the term -C T G k G C G u n G , due to the Robin parameter k G . In fact, the contribution of Robin conditions leads to local changes of global stiffness operators regarding the interface degrees of freedom. For another framework of dealing with Robin conditions see [START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF] As mentioned in the previous section, auxiliary problems are not needed if finite element codes can get reaction forces on an immersed surface. Nevertheless, if this is not the case, the auxiliary problem in Ω C is still not necessary. Indeed, it can be computed from known and easy computations without assembling the operator K C :

λ n+1 C = C C K C u n+1 G |Ω C -f C d = C G K G u n+1 G -f G d -K A u n+1 G -f A d = C G C T G P n -λ n+1 A λ n+1 C = P n -λ n+1 A ( 27
)
The same relation as the standard global/local method is given where

P n = λ n+1 C + λ n+1 A .
In addition, since

P n = λ n A -λ n L + k G C L u n L -k G C G u n
G , P n+1 can be written in function of P n , a new rest r n+1 will appear:

P n+1 = λ n+1 A -λ n+1 L + k G C L u n+1 L -k G C G u n+1 G = P n -λ n+1 C -λ n+1 L + k G C L u n+1 L -k G C G u n+1 G = P n -λ n+1 C + λ n+1 L + k G C G u n+1 G -C L u n+1 L = P n + r n+1 (28)
The new rest r n+1 now incorporates a mixed contribution of the discontinuity of displacements and the disequilibrium of forces:

r n+1 = -λ n+1 C + λ n+1 L Disequilibrium -k G C G u n+1 G -C L u n+1 L discontinuity (29)

Algorithm for Global/Local method with Robin Parameters

The mixed global/local methodology with Robin conditions can be located in Alg. 2. Solve the global problem:

Find u n+1 G , K G u n+1 G = f G d + C T G P n 6:
Solve the auxiliary problem:

Find u A , λ n+1 A , K A C T A C A 0 u A -λ n+1 A = f A d u n+1 G |Γ 7:
Compute the reaction force λ n+1 C :

λ n+1 C = P n -λ n+1 A 8:
Solve the local problem:

Find u n+1 L , K L + C T L k L C L u n+1 L = f L d + C T L k L C G u n+1 G -λ n+1 C 9:
Compute the reaction force λ n+1 L :

λ n+1 L = -λ n+1 C -k L C L u n+1 L -C G u n+1 G 10:
Compute the residual of the iteration

r n+1 = -λ n+1 C + λ n+1 L + k G C G u n+1 G -C L u n+1 L 11:
Compute the new P n+1 :

P n+1 = P n + r n+1 12:
Update the new P n+1 including relaxation factor µ:

P n+1 = µP n+1 + (1 -µ)P n 13:
Estimate the error of the iteration step as η = r n+1 2 / r 0 2 14:

Verify convergence criteria: η < tol with tol = 10 -4

The method error is calculated by using the same expression as for the primal-dual strategy presented in Eq. ( 16). However, each residual term in Eq. ( 16) is calculated differently for the mixed global/local analysis (it is 155 estimated as a combination of interface forces and displacements). It could be interesting to take into consideration the error of the results reached when using the global/local analysis with mixed conditions and the results regarding the monolithic computation. This error is calculated as the rule of displacement difference on the interface between the global and the monolithic model, as presented in Eq. [START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF].

η M ono. = C G u n+1 G -u M ono. 2 u M ono. 2 (30)
Figure 3 shows a flowchart of the iterative method implementation of the global/local with Robin parameters. Remark 4. In the Code Aster implementation, the Robin operator is calculated as follows: the strip of elements is degrees of freedom on the interface. This strip of elements is connected directly to the interface and with fixed support boundary conditions. On the other hand, the strip of elements is used to compute the part on the right-hand side for the global model, which corresponds to the Robin conditions in Eq. (28). It is not assembled with the global stiffness. The strip of elements considers a Poisson ratio of 0.3, a Young Modulus variable (α G E G ) and fixed boundary conditions on the opposite side of the interface, as presented in Figure 5. crack tip; Therefore, based on St. Venant's principle, only the effects of these singularities are transmitted between models. In order to study specific aspects of the non-linear behaviour, the local model should be further analyzed. Finally, the "total" solution of the analysis should be considered as the superposition of the global linear model and the local non-linear model.

Results

5 Optimization of the Robin parameter Global/Local Analysis with Robin Conditions in 2D plane stress

Description of the optimization process

As presented in [START_REF] Lorenzis | Modeling in Engineering Using Innovative Numerical Methods for Solids and Fluids[END_REF], the convergence of the global/local analysis with mixed conditions can be assured for a certain range of values for the Robin parameters. Within this range, an optimum that minimizes the number of iterations can be found. However, as presented in [START_REF] Saavedra | An enhanced nonlinear multi-scale strategy for the simulation of buckling and delamination on 3D composite plates[END_REF][START_REF] Fuenzalida-Henriquez | Optimizing the search directions of a mixed DDM applied on cracks[END_REF], the value depends on the problem type, the applied loads, the boundary conditions and the form of the Robin operator used in the analysis. Therefore, in order to find the Robin operator that guarantees convergence and also minimizes the number of iterations until convergence, an optimization process is performed for the different problems presented. Optimization algorithms have been used widely in engineering problems and other scientific areas, such as the Basin-Hopping Algorithm [START_REF] Locatelli | On the Multilevel Structure of Global Optimization Problems[END_REF][START_REF] Olson | Basin Hopping as a General and Versatile Optimization Framework for the Characterization of Biological Macromolecules[END_REF], already distributed within Scipy library [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF].

Remark 6. The values of the Robin operator assuring the convergence of the non-intrusive method with Robin conditions can be estimated analytically for simple examples with few degrees of freedom. For industrial/real structures, the range of the values are not easily found directly, so these values are found by iterative methods as presented in [START_REF] Hinojosa | Domain decomposition methods with nonlinear localization for the buckling and post-buckling analyses of large structures[END_REF][START_REF] Fuenzalida-Henriquez | Optimizing the search directions of a mixed DDM applied on cracks[END_REF] or by using two scale strategies as shown in [START_REF] Negrello | A new impedance accounting for short-and long-range effects in mixed substructured formulations of nonlinear problems[END_REF][START_REF] Gendre | A two-scale approximation of the Schur complement and its use for non-intrusive coupling[END_REF].

For this article, the Basin-Hopping algorithm is used in the optimization process, allowing the specification of an initial "guess" of the problem within the search space and stochastically start the refinement of the solution until finding the best possible solution. The Basin-Hopping algorithm is based on the transformation of the "potential energy landscape" of an objective function, finding local minima using the Limited-Memory Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS) [START_REF] Byrd | Representations of quasi-Newton matrices and their use in limited memory methods[END_REF]. This potential energy surface is being explored by means of random jumps for each evaluation based on the Boltzmann probability distribution. The number of random jumps is defined by a variable N, in which the probabilities for finding the best local minimum are increased for larger values of N.

Therefore, the Basin-Hopping algorithm has been used successfully to find a minimum using N=110 or larger [START_REF] Wales | Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms[END_REF]. This parameter is also a measurement of the computational cost, because it is a direct measurement of the amount of random jumps between each solution. For our work, N has been defined as the default option provided in Scipy, with a value of N equal to 100. Finally, the L-BFGS method has a computational complexity of O(n 2 ) and must be repeated N times in accordance with the number of random jumps [START_REF] Byrd | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF].

Nevertheless, a global optimum cannot be assured in metaheuristics, due to the non-convex nature of the problem. However, a good approximation of a global minimum solution can be found. This is due to the internal behaviour of the optimization process: it is stopped through an internal convergence criterion and not by founding a global minimum [START_REF] Bianchi | A survey on metaheuristics for stochastic combinatorial optimization[END_REF]. More information and applications on the different optimization algorithms can be seen in [START_REF] Das | Differential evolution: A survey of the state-of-the-art[END_REF][START_REF] Kukkonen | Generalized differential evolution for general non-linear optimization[END_REF][START_REF] Babu | Differential evolution for multi-objective optimization[END_REF][START_REF] Li | Adaptive Differential Evolution with Evolution Memory for Multiobjective Optimization[END_REF][START_REF] Marinho | Application of differential evolution to multi-objective tuning of vibration spectrum analyzers based on microelectromechanical systems[END_REF][START_REF] Burnham | Crystal Structure Prediction via Basin-Hopping Global Optimization Employing Tiny Periodic Simulation Cells, with Application to Water-Ice[END_REF][START_REF] Zhou | Augmenting Basin-Hopping With Techniques From Unsupervised Machine Learning: Applications in Spectroscopy and Ion Mobility[END_REF][START_REF] Schönborn | The performance of minima hopping and evolutionary algorithms for cluster structure prediction[END_REF], among others. Therefore, considering that the convergence of the method depends on the Robin parameters as presented in [START_REF] Lorenzis | Modeling in Engineering Using Innovative Numerical Methods for Solids and Fluids[END_REF], the initial Robin parameter (calculated as the condensed stiffness on the interface estimated from an initial strip of elements and presented in Figure 5) can be modified by means of two factors: α G for the global model and α L for the local fine model, as shown in Eq. ( 31):

k L = α L k L init k G = α G k G init (31)
The optimization algorithm will be used to find an optimum value for α G,L that minimizes the number of iterations until convergence.

It is important to note that the search space for the values of α G,L to be used in the optimization process are constrained between 0 and 100 in order to find an overall optimum that minimizes the number of iterations, starting with the initial search direction k init G,L , i.e., α G = 1 and α L = 1. in Table 5 and presented in Figure 11 for the cracked case and in Figure 12 for the plastic hardening case. The use of the optimized α G,L implies a reduction on the number of iterations from 10 to 9 for the cracked case with 3 propagation steps, meaning a 18% faster respect to the primal-dual and a 10% improvement in relation to the global/local with non-optimized Robin conditions. For the hardening problem, the optimized α G,L reduces the number of iterations from 10 to 8, which corresponds to an improvement of 20% with respect to the non-optimized method with Robin parameters using this type of analysis. The feasibility of coupling the global/local method with phase field approaches to deal with crack propagation has been studied successfully in [START_REF] Gerasimov | A non-intrusive global/local approach applied to phasefield modeling of brittle fracture[END_REF][START_REF] Noii | An adaptive global-local approach for phase-field modeling of anisotropic brittle fracture[END_REF].

The studied strategies for solving the global/local problem, i.e., primal-dual and with Robin conditions, work without problems for incompatible patches. Nevertheless, a sensitivity analysis can be performed to study the effect of incompatible meshes for different structures and behaviours, and the error that can be induced due the interpolation of forces and displacements in the different interfaces.

Basin-Hopping is a meta-heuristic that can be easily implemented, but as presented in [START_REF] Blum | Metaheuristics in combinatorial optimization: Overview and conceptual comparison[END_REF], it may not be the fastest to obtain convergence. Therefore, a more efficient way to optimize the Robin parameters and reduce computational times must be studied and applied, such as the Particle Swarm Optimization [START_REF] Yasear | Review of the multi-objective swarm intelligence optimization algorithms[END_REF] and exploiting the ability to evaluate multiple trials with parallel architecture. In addition, the efficiency of the algorithm on larger cases will be studied. For now, due to the small problems we studied, the iterative algorithm cannot compete with a full-scale computation. We expect that larger problems may reverse this trend.

Finally, the effect of the discretization for the 2D and 3D models, the degree of the shape functions used for the FEM method elements and also the coupling between 2D global models and 3D local non-linear models must be analyzed, in order to study the convergence of these elements and how affects the response.

Figure 1 :

 1 Figure 1: Reference problem.

Figure 2 :

 2 Figure 2: Linear problem.

Ad

  corresponds to the boundary conditions that could be implemented on the auxiliary model in order to represent an approximation of those on the local model. The term K A u A -f A d symbolizes the reaction forces of a Dirichlet problem on the auxiliary model with displacement u G|Γ imposed in Γ and a Neumann boundary condition through f A d

Algorithm 1 1 :

 11 Global/Local Primal-Dual Solution Initialize interface vector P 0 = 0 2: Initialize error η = 1 and tolerance tol 3: Start iterations, n = 0 4: while η > tol do 5:

Algorithm 2 1 : 2 : 3 :

 2123 Global/Local algorithm with Robin Conditions Initialize interface vector P 0 = 0 Initialize error η = 1 and tolerance tol Start iterations, n = 0 4: while η > tol do 5:

  160

Figure 3 :

 3 Figure 3: Flowchart of the algorithm of the Global/Local Analysis with Robin Parameters.

Figure 5 :

 5 Figure 5: Element strip of elements used for the analysis.

Figure 6 :

 6 Figure 6: Method error for Global/Local Analysis with Robin Parameters and cracked domain with 4 propagation steps.

Figure 8 :

 8 Figure 8: Method error for Global/Local Analysis with Optimized Robin Parameters and cracked domain.

Figure 9 :

 9 Figure 9: Method error for Global/Local Analysis with Optimized Robin Parameters and plastic hardening domain.

Figure 11 :

 11 Figure 11: Method error for Global/Local Analysis Robin Parameters for 3D cracked domain and optimized α G,L .

12:Verify convergence criteria: η < tol Moreover, the case of matching meshes on the interface between auxiliary and local models is assumed. For the case of non-matching meshes, a projection step needs to be developed to communicate the displacement field from the global to the local model, and the reaction forces from the local to the global auxiliary domain.This algorithm is no more than a fixed point for which the convergence is known to be slow. Full linear models can be considered as preconditioners for a Krylov solver to improve the convergence rate. In the case of a non-linear local model, quasi-Newton or Aitken's relaxation must be taken into account[START_REF] Duval | Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition[END_REF][START_REF] Gosselet | Non-invasive global-local coupling as a Schwarz domain decomposition method: acceleration and generalization[END_REF].

Global/Local method with Robin ParametersAs previously stated, the standard global/local analysis can be observed as a coupling of the complementary and local model by means of a Lagrange's multiplier. This coupling problem can be solved with an iterative Dirichlet-

Acknowledgements

This work was partly funded by the international project ECOS-CONICYT C17E04 between LMT and Universidad de Talca, ANID/CONICYT grant number ECOS C17E04. In addition, this research was funded by ANID, grant number PFCHA / DOCTORADO BECAS CHILE / 2018 -21181707 and BECA ESTUDIO DE DOCTORADO, UNIVERSIDAD DE TALCA

The Robin parameter is chosen as the stiffness of an element strip [START_REF] Roux | Domain decomposition methodology with robin interface matching conditions for solving strongly coupled problems[END_REF], as presented in Figure 5. Robin operators k L and k G are determined independently: k L is the linear elastic stiffness of the strip on the side of It is clear that the global/local Analysis with Robin parameters converges slower with respect to the primal-dual Method (from 8 to 6 iterations compared to the primal-dual for the cracked model). For the hardening problem, it converges from 8 to 6 with respect to the primal-dual and considering the number of iterations. The Aitken δ 2 dynamic relaxation converges both models rapidly, with a faster convergence rate for the primal-dual Method.

The complete global/local analysis results with the initial Robin operator are shown in

Summary of the Results

As presented in Section 4, the mixed method converges slower with respect to the primal-dual method. Regarding computational times, these were 18 seconds for the cracked model and 10 seconds for the non-linear hardening problem. However, these times are referential and can be improved by optimizing the python code. Finally, as the analysis considers an arbitrary strip of elements used to start the Robin operator, an optimization process is carried out with the intention of reducing the number of iterations until convergence and analysing the results for different Robin parameter values. In addition, the strip of elements are very stiff and do not represent the neighbour structure, i.e., the theoretical optimum.

Remark 5. Is important to mention that non-linear plastic behaviour or crack tip propagation are not being transferred between the local and global models. The interfaces are considered far from the source of plasticity or the

The optimization problem is carried out using the initial configuration of the domains and considering the initialized Robin conditions and multiplied by the α factor. The analysis contemplates linear elasticity behaviour for both problems for all optimization trials. Finally, the results for the optimum α G,L will be used for non-linear problems (all propagation steps and post-yielding hardening).

The following procedure is applied to perform the optimization of the Robin parameter, using the Basin-Hopping [START_REF] Olson | Basin Hopping as a General and Versatile Optimization Framework for the Characterization of Biological Macromolecules[END_REF] algorithm, as well as for linear elastic initialized models, i.e., without plasticity or crack propagation:

• Choose global and local models to be analysed considering the elastic linear behaviour, as presented in Figure 4a and Figure 4b.

• Perform the analysis on the model considering the iterative non-intrusive scheme with α L = 1 and α G = 1 (Initial guess or initial trial for the Robin parameter).

1. Choose new values of α L and α G found in the local neighbourhood of the initial guess.

2. Perform the non-intrusive analysis using Robin parameters and recover the number of iterations until the analysis converges for each α L and α G in the neighbourhood of the current trial.

3. Study the convergence values of the previous neighbourhood and select random values in the best direction (gradient) found in the neighbour to proceed with updated α G,L trial.

• When the optimization algorithm converges and the best solution for the linear model is obtained, the best factors α opt L and α opt G are recovered.

• Apply α opt L and α opt G to the non-linear model (hardening or cracked model) and study the convergence.

Results

For the 2D cases, the optimization results are presented summarized in The optimized α found for the linear cases are used directly as an input for the non-linear cases. These results, regarding the evolution of the method error with respect to iterations, are introduced in Figure 8 for the cracked case and in Figure 9 for the plastic hardening case. Figure 8 and Figure 9 show that the number of iterations are reduced by 25% for the cracked case and 12.5% for the plastic hardening case respecting the non-optimized solution (initial guess of the Robin parameter). However, with regard to the primal-dual solution, the cracked optimized model equals the performance in terms of the number of iterations. For the non-linear hardening case, the optimized model with Robin parameters performs slower with respect to the primal-dual method. The results for the cracked model and non-linear hardening, using the optimized α G,L factors, are presented in Table 3.

Summary of the Results

When optimizing the linear elastic analysis of each model (including the initial configuration of the local geometry) using the Basin-Hopping algorithm, the convergence of non-linear models is improved as expected. The optimized α opt L and α opt G used for the non-linear problems obtain better results with respect to the non-optimized problems, similar to the performance of the primal-dual method. This approximation using the linear elastic optimized model works for localized patches with non-linearities. However, the solution found with the Basin-Hopping algorithm for the linear models may not be the actual optimum for the non-linear models. Finally, as the Code Aster code with the python interface was implemented, the optimization times are elevated in comparison with the non-optimized α case. However, these times can be optimized by using better coding practices and more efficiency in computing resources.

Global/Local Analysis with Robin Conditions for 3D structures

As presented in Section 5, the same approach will be applied in two models: a cantilever beam with an initial crack near the fixed support and a beam with two 50 (mm) perforations and non-linear hardening behaviour. The properties are the same as the previous case but including different lengths, loads applied and crack size, as shown in Figure 10. The crack propagates for the 3D case, including 3 propagation steps and 10 (mm) for maximum crack advance per step for the cracked domain. As the discretization carried out in the global and local 3D models are geometrically complex, incompatible meshes were generated. The incompatibility is managed by Code Aster projecting the fields between the models and the interface. The initial Robin parameter will be used in order to find the results for the non-optimized α G,L for both models, which are summarized in Table 4. Results using the optimized Robin parameter, following the same procedure as shown in Sec. 5, are summarized Table 5: Results for the 3D cracked models 3 propagation steps.

Model

The final deformed shape for the cracked model is presented in Figure 13 for the Global/Local 3D model with Robin conditions (optimized) and in Figure 14 for the hardening case.

Summary of the Results

As in Section 5, the optimization improves the convergence with respect to the initialized Robin parameter.

The optimization process, when applied to the 3D cracked model, improves the convergence rate, resulting in better performances, considering the number of iterations with respect to the primal-dual methodology.

Discussion

In this article, the global/local non-intrusive analysis with Robin parameters was implemented in 2D and 3D problems with cracked and plastic hardening behaviour. Subsequently, results were optimized using the Basin-Hopping algorithm and compared to a primal-dual non-intrusive analysis.

The method presented in this study has a good convergence rate for 2D and 3D non-linear plastic hardening behaviour and cracked models. However, as for any mixed method, the selection of the Robin parameters is crucial for convergence. Our initial choice is based on studying the effect of a simple element strip, but it may not be very efficient compared to the primal-dual approach.

The optimization of the mixed global/local method was performed using the Basin-Hopping algorithm, allowing an improvement in the performance of the global/local analysis with Robin parameters. The optimized α is different for each case, which is consistent with results presented in [START_REF] Saavedra | An enhanced nonlinear multi-scale strategy for the simulation of buckling and delamination on 3D composite plates[END_REF], indicating that the values of α are specific to each model (geometry, behaviour, etc.), loads and boundary conditions. Nevertheless, the optimization does not guarantee that the mixed global/local analysis works better with respect to the primal-dual non-intrusive analysis, although a local minimum could improve the mixed analysis results.

Future research

The Robin parameter and the flexibility added to the structure are a great advantage of the mixed model, as presented in [START_REF] Hinojosa | Domain decomposition methods with nonlinear localization for the buckling and post-buckling analyses of large structures[END_REF]. Therefore, the non-intrusive analysis with Robin parameters can be studied for structures with geometric non-linear behaviour and large displacements.

Another topic to be studied further is to extend the analysis to larger structures, associated with multiple domains and with multiple stiffeners or variable cross sections. This is useful to study the convergence of the method in structures with strong modifications in the stiffness. Consequently, the original approximation of the Robin parameter may not be the best for each interface and these factors could be different for each analysis direction, as presented in [START_REF] Fuenzalida-Henriquez | Optimizing the search directions of a mixed DDM applied on cracks[END_REF].

Different crack propagation techniques must be considered in future studies, for example, the Cohesive Zone Model (CZM) and the Virtual Closing-Crack Technique (VCCT) to study the effectiveness of the global/local

Conflict of interest

The authors declare that they have no conflict of interest.