Guglielmo S Aglietti
email: g.aglietti@surrey.ac.uk

Ben Taylor

Simon Fellowes

Thierry Salmon

Ingo Retat

Alexander Hall

Thomas Chabot

Aurélien Pisseloup

A, Zarkesh C Cox

A Mafficini

N Vinkoff

K Bashford

Cesar Bernal

François Chaumette

Alexandre Pollini

Willem H Steyn

The active space debris removal mission RemoveDebris. Part 2: in orbit operations

Keywords: space debris, debris removal, ADR, deorbiting, net, harpoon, vision-based navigation, dragsail

This is the second of two companion papers that describe the development of the RemoveDEBRIS space mission. This second article describes the in-orbit operations that were performed to demonstrate technologies to be used for the active removal of space debris, whereas the first paper described the development of the satellite's hardware.

The RemoveDebris mission has been the world's first Active Debris Removal (ADR) mission to successfully demonstrate, in orbit, some cost effective technologies, including net and harpoon capture; and elements of the whole sequence of operations, like the visionbased navigation.

The satellite was launched the 2 nd of April 2018, to the International Space Station (ISS) and from there, on the 20 th of June 2018, was deployed via the NanoRacks Kaber system into an orbit of 405 km altitude.

During the mission, two 2U CubeSats have been released by the mothercraft platform as artificial debris targets, to demonstrate net capture and cameras to be used for vision based navigation. Harpoon capture has been demonstrated by deploying a target and then firing at it a harpoon tethered to the platform. The various phases of the missions have been monitored using relevant telemetry and video cameras, and this paper reports the results of the various demonstrations.

I. INTRODUCTION

The background to the development of the RemoveDEBRIS mission, is discussed in detail in [START_REF] Forshaw | RemoveDEBRIS: An in-orbit active debris removal demonstration mission[END_REF], [START_REF] Forshaw | Final payload test re-sults for the RemoveDebris active debris removal mission[END_REF] and the companion article [START_REF] Forshaw | The active space debris removal mission RemoveDebris. Part 1: from concept to launch[END_REF] to this paper, where the design and testing of the spacecraft hardware, is also described. Various updates on the development of the mission were given periodically at the leading international conferences [START_REF] Saunders | Mission and systems design for the debris removal of massive satellites[END_REF], [START_REF] Forshaw | Surrey Space Centre: A Survey of Debris Removal Research Activities[END_REF], [START_REF] Forshaw | Review of final payload test results for the RemoveDebris active debris removal mission[END_REF], [START_REF] Forshaw | The RemoveDebris ADR Mission: Launch from the ISS, Operations and Experimental Timelines[END_REF] and [START_REF] Taylor | RemoveDebris Preliminary Mission Results[END_REF]. To avoid repetitions, this article focuses on the in-orbit operations, and this brief introduction puts the work into context.

With future mega-satellite constellations being deployed, where hundreds to thousands of satellites are going to be launched into orbit (e.g. Oneweb or SpaceX' Starlink), to maintain under control the quantity of spacedebris, has become an issue of paramount importance [START_REF] Virgili | Risk to space sustainability from large constellations of satellites[END_REF].

A coherent strategy, along with technological and platform developments, is needed for the de-orbiting, re-orbiting, or servicing satellites, in particular considering the advent of such constellations. The issue of the space debris and related mitigation strategies is receiving growing attention [START_REF] Bonnal | Active debris removal: Recent progress and current trend[END_REF], [START_REF] Innocenti | Clean space -an overview[END_REF], and a variety of concepts have been proposed, see for example [START_REF] Shan | Review and comparison of active space debris capturing and removal methods[END_REF] and [START_REF] Mark | Review of Active Space Debris Removal Methods[END_REF]. In this context, RemoveDebris has been the first space mission to demonstrate technologies for the active removal of space debris.

Cost will be decisive factor in determining the future of ADR (e.g. [START_REF] Federicoj | Conceptualizing an economically, legally, and politically viable active debris removal option[END_REF] [START_REF] Yamamoto | Cost analysis of active debris removal scenarios and system architectures, 7th European Conference on Space Debris[END_REF]) Hence, developing the mission concept, particular emphasis has been placed on minimizing the cost, and this has been a crucial element in the selection of the technologies that have been demonstrated and in the decision to launch the craft via the ISS.

The development of the hardware has met various challenges which related to its launch ( [START_REF] Forshaw | The RemoveDebris ADR Mission: Preparing for an International Space Station Launch[END_REF] [START_REF]Safety Review Process-International Space Station Program[END_REF] [START_REF]Space station CubeSat deployment services[END_REF], from technical issues, for example, related to the exact definition of the launch vibration environment and appropriate levels for testing, to the need to comply with the ISS safety regulations, undergoing the various levels of the NASA safety process. For example, items such as the Cold Gas Generators (CGG) and the platform battery presented some concerns, the first related to the device incorporating a pressurized container and the nature of the chemicals in the device, the latter due to its size and capacity. CGG were used in all the inflatable structures (DSAT#1 and DragSail)) and for the harpoon, and were an essential part of the design. Beside descriptions of all the chemicals and their quantity, the team had to demonstrate that that there were three electrical inhibits for the activation of the devices, which made accidental activation not credible. For the battery, because for the exact model that was used on the RemoDEBRIS platform, no test results were available to reassure against the possibility of a thermal runway, a new test battery had to be produced to demonstrate battery safety under thermal runaway conditions. These are just examples that suggest that early engagement with the launch authorities facilitates the process as the appropriate requirements can be built in the design, rather than retrofitted, as this produces further costs and potential delays.

Nanorack provided technical guidance and logistical support for all the operation related to transportation to the launch site, up to release in orbit using the ISS robotic arm equipped with the KABER interface. .

Once in orbit, mission control was performed by SSTL, using their ground station in Guildford, and for each of the demonstrations, the institution in charge of the demonstration provided support and guidance, as reported in Table 1.

I.1. Paper Structure

Sections 2 describes the operations from the packing of the satellite at the end of the AIT to its release in orbit and satellite commissioning. Sections from 3 and 6 describe the four demonstrations, Net Capture, VBN, Harpoon capture and DragSail respectively. For each demonstration the salient events are described and overall outcomes discussed. Finally, Section 7 concludes the paper and outlines key contributions to the field.

II. FROM LAUNCH TO SATELLITE COMMISSIONING

II.1. Launch

RemoveDEBRIS, complete with its protective panels (see Figure 1), was accommodated inside a foam clamshell (Figure 2), and for launch the craft in its clamshell were accommodated within soft a Cargo Transfer Bag (CTB), strapped to the rack sidewall of Dragon capsule (see Figure 3). The capsule is pressurized and this launch configuration provides a soft ride since the platform is not hard mounted to the launcher. The launch was part of the 14 th Space X CRS (Commercial Resupply Service) to the ISS, and the capsule was launched on the 2 April 2018 at 20:30 UTC on a Falcon 9 Full Thrust rocket from Cape Canaveral Air Force Station Space Launch Complex. The decision to launch the satellite using a commercial resupply mission to the ISS and released in orbit by the ISS robotic arm, rather than a traditional launch, directly in orbit using any of the available launch opportunities (e.g. as a piggy back payload) was taken for two reasons. One concerned the safety of the mission, as the low altitude of the injection point in orbit guaranteed a relatively rapid deorbiting of all the objects. In particular the predictions reported in Figure 4 show that also in case of a malfunction of the de-orbit sail the craft should reenter the atmosphere in ~120 weeks, (i.e. well within the 25 year guideline).

The second reason was commercial, as the service offered by Nanorack via the ISS was at a price lower than the other opportunities available at the time.

The Dragon capsule reached the ISS two days after launch, the 4th of April and was captured by Canadarm2 at 10:40 UTC (see Figure 4) and was berthed to the Harmony module at 13:00 UTC.

II.2. ISS operations

After being taken out of the CTB, the craft, within its foam case was then transferred to the Japanese Experiment Module (JEM). Drew Feustel (ISS Commander) assisted by Ricky Arnold, performed the unpacking maintaining contact with the RomveDEBRIS team via video link. The craft with its protective panels, was extracted from the foam clamshell case (see Figure 5) and the protective panels removed. The craft was then installed on the sliding table (see Figure 6) of the JEM airlock, integrating it to the NanoRacks Kaber Microsat Deployer. Apart from handling tasks, no specific operation like battery charging, electrical checks, etc…, were requested to the astronauts. Once the table was slid into the airlock and the door closed, the airlock was depressurized. The space side door of the airlock was then opened, and the robotic arm with the Special Purpose Dexterous Manipulator (SPDM) grappled RemoveDEBRIS by the KABER interface. The craft is then moved by the robotic arm to a position and orientation from which it can be safely released in orbit (Figure 7), drifting away from the ISS. During all this period, all the platform subsystems remained turned off (Figure 8), as ISS safety imposed that the spacecraft needed to be fully off from 30minutes minimum after deployment from ISS. At the time of the RemoveDebris separation, the ISS was on a ~405km altitude orbit, and timers were triggered by 3 redundant mechanical switches implemented on the NRSS NanoRacks Satellite Separation system (PSC Mark II light band).

II.3. LEOP and Satellite commissioning

After 30mn, timers activated the platform power subsystem, S-band receivers and survival heaters. The attitude of platform was free and not controlled. During the first ground pass over the groundstation in Guildford, a command was sent by the ground segment to activate platform OBC and S-band transmitters. Preliminary communication and checkout tests were performed on the platform. Basic telemetry (e.g. battery charge) showed the craft was performing nominally and the attitude was acquired. This was followed by a slow de-tumbling of the spacecraft that was stabilized into a controlled safe attitude. AOCS commissioning then progressed until the platform was brought under control in a coarse Nadir pointing mode.

The platform commissioning continued checking all subsystem to ensure they survived launch and did not present malfunctions. For example, prime and redundant RF receivers, low rate transmitters and low level command links were checked.

The spacecraft then performed a series of AOCS manoeuvres to verify performance against that required for executing payload experiments.

During commissioning, an issue with the platform Gyros was detected, as the values measured were not consistent. The expected time for recovery was estimated at several weeks, even months. Consequently, it was decided to develop an alternative mode for platform pointing without Gyro and named "nominal gyroless mode". This mode was successfully tested during commissioning and therefore provided the possibility to continue with the mission also without the Gyros. The final phase was the payload commissioning, which included calibration and characterization. The Supervision cameras and VBN camera were tested over a range of exposures and frame rates (see Figure 9) which were planned for use for the experimental demonstrations Overall, the commissioning phase lasted about 8 weeks, and this duration was consistent with the NASA request that no demonstrations should have been performed during a period of 2 month after separation to avoid risk of collisions with ISS. The timing of the various events is shown in Figure 10.

VBN camera

CubeSats Net experiment cover Harpoon's target

III. NET DEMONSTRATION

The Net capture was the first demonstration to be carried out, and during this experiment the CubeSat DSAT#1 was released at low speed from the ISIPOD#1 installed on the mothercraft. DSAT#1 then inflated a structure to increase its size thus becoming more representative of a large space debris, and was captured by the Net launched by the mothercraft. The specific details are reported in the following and in Table 2. First, ISIPOD#1 door was activated by the platform (T0) and DSAT#1 slid out, restrained by the Cubesat Release System (CRS#1) Second, the CRS#1 is activated by DSAT#1 itself at T0+60s to eject DSAT#1 using the energy stored in low stiffness leaf springs that pushed the cubesat gently away.

The DSAT#1 & Net synchronization was done by timers starting the ISIPOD activation (T0). Magnetic sensors were used to confirm DSAT#1 departure. Platform pointing to Net deployment direction, platform compensation after DSAT#1 departure and platform stabilization after Net deployment were performed nominally in the "gyro less nominal mode". All other sub systems (power, coms) behaved as expected. All OBC logs and TLM data was transmitted to ground.

Four different cameras supervised the experiment, Supervision camera 1 & 2 that captured 2305 images each (over a period of 440s), VBN 2D camera, that captured 17 images and finally the VBN 3D camera that captured 29 images. After the release of DSAT#1, a timer commanded the two cold gas generators to inflate the five booms of the assembly (design details in [START_REF] Forshaw | The active space debris removal mission RemoveDebris. Part 1: from concept to launch[END_REF]). There is video evidence of the inflation/deployment of two of the 4 lateral booms and the longitudinal boom, but the other two lateral booms do not appear to have deployed (See Figure 11). The deployed longitudinal boom can be seen in Figure 12, where it is also visible the triangular sail that would not have been deployed and visible unless the longitudinal boom had been deployed. Net deployment was performed 144s after DSAT#1 departure (T0+204s), theoretically leading to a capture at 7m distance. However during the demonstration, DSAT#1 velocity has been assessed based on the analyses of the timing of the images and known size of the objects and camera Filed of View. Results seems indicate that DSAT#1 velocity was slightly higher than expected (about 50% higher than then planned 5cm/s at 7.5cm/s).

The video stream shows that after deployment of its inflatable booms, DSAT#1 started spinning which is consistent with a gas leak from the lateral booms which did not inflate/deploy correctly. It is likely that this has been the result of the lateral door of the boom compartment not opening, which would have prevented the boom from expanding under the pressure of the inflation. Constraining the boom would have generated a further increase of pressure that is likely to have contributed to a rupture of the boom producing a small air jet that acting like a thruster, would have imparted an angular speed to the CubeSat. This might have also contributed to the lateral drift of the cubesat, as it was moving away from the mothercraft. Figure 12 show the net opening up as it travels toward DSAT#1, pulled by the 6 masses at its vertexes, which travelled along diverging trajectories, so that the net unfolds in the shape of a six point "star" where each of the vertexes is pulled by one of the masses. The deployment direction of the target (DSAT#2) for this experiment was driven by two criteria: i) the risk of collision between platform and DSAT#2 after deployment must be minimized, ii) DSAT#2 must enable the acquisition of images with range from [0;1000m], various backgrounds (Earth and black sky), and lighting conditions (daylight and eclipse).

The plan was to observe DSAT#2 drifting away at a speed of 2cm/sec using the VBN: 2D Camera (FoV~18deg ) and the 3D Flash LiDAR (FoV<10deg). The end of the demonstration was planned for 4h30 (~3 orbits) after the deployment of DSAT#2, which by then was supposed to be at 750m from the mothercraft. Also the supervision camera SV0 on the mothercraft observed the DSAT#2 drifting away. SV0 took in total 3691 images over 4h02mn50s: (1947 images during its 1 st orbit, 1161 during the 2nd orbit and 583 during the 3rd orbit ), the 2D and 3D VBN cameras took a total of 361 images each.

For the LiDAR camera, the angle of the surfaces that are observed is crucial to determine their performance, as this governs the ratio between the diffused-reflected sun light and the laser light reflected by the target that hit the camera detector. The best case for the LiDAR (and the worst for the cameras) is when specular light from laser and no, or only diffused light from sun are collected.

During the experiments, the angle with the sun light was to favor the cameras. Hence, images were captured in the less favorable situation for the LiDAR, with specular sunlight reflection into detector. This is illustrated by the pictures in Figure 13, taken at a few seconds interval by the VBN color camera (left) and by the VBN LiDAR (right). The field of view (FOV) of the LiDAR and camera are respectively of 8°x6° and 21°x16°. This design choice resulted from the tradeoff of counter requirements for a large field to increase the chance of viewing the target despite the uncertainty on its relative trajectory uncertainty and for sufficient spatial resolution to resolve a target of 20cm as far as possible (see Figure 14). The ground truth available for the VBN experiment is provided by the GPS integrated in the DSAT#2. With the GPS positions transmitted via the wireless Inter Satellite Link between the target (DSAT#2) and the chaser (mothercraft), we have the possibility to know roughly the relative distance between the two satellites.

In Figure 15, some of the distances measurement by the LiDAR have been superimposed on the plot of the distances generated from the GPS data, showing a good match between the two sets of measurements From Figure 15 it is possible to notice that there is no valid LiDAR measurement between 2100 and 4500 [s].

During this period of time Earth was in the background and distance measurement cannot be derived from the raw LiDAR images. This is illustrated in Figure 16, where the target signal is completely buried in the background noise. The first fact to be highlighted, and shown in Figure 18, is that DSAT#2 remained in the cameras field of view over the entire sequence, meaning that the attitude guidance profile was correctly generated and followed by RemoveDEBRIS, and that dispersions at ejection remained at an acceptable level Different VBN algorithms provide 3DoF (position only) and 6DoF (position and attitude) relative navigation depending on the rendezvous phase. At far and medium range, only 3DoF (position only) relative navigation is needed and performed. At short range, 6DoF (position and attitude) relative navigation is required for successful capture. An a priori knowledge of a 3D model of the target (see Figure 17) is matched to an image in an iterative fashion, leading to an estimate of the relative position and attitude between the camera and the target.

The proposed approach relies on a frame to frame Model based tracking in order to obtain the complete pose of the CubeSat (position and attitude) with respect to the vision sensor. It minimizes the error between visual measurements in the image and the projection of the 3D model of the CubeSat. Tracking and pose estimation are thus simultaneous. It can also be combined with an extended Kalman Filter to improve the navigation precision.

Images (1280x1024 pixels, 12mm lens) were acquired at 3Hz with successively 3 different aperture times, meaning that successive images at same aperture were acquired at 1 Hz. In particular, the CubeSat (DSAT#2) has been successively tracked and localized (see Figure 20). Note that the processing time for each frame is less than 0.1 sec.

The results from the image processing algorithms based on Supervision camera images and VBN camera images are shown hereafter in Figure 19, together with LiDAR data, GPS measurements and the expected, nominal DSAT2 trajectory. Based on the processing of the data, the following statements may be made: SV0 and VBN camera measurements are consistent (2.3 cm/s ejection), LiDAR measurements are close to expected trajectory (2 cm/s ejection), GPS measurements do not seem reliable for short range (<10m). Further analysis has been performed in order to process the GPS data provided by the DSAT2 and mothercraft, with the objective of achieving a reference relative trajectory, to be used as a ground truth against which the VBN performances could be assessed, however these results show some irregularities and more in-depth analysis will be required.

V. HARPOON DEMONSTRATION

This was the third demonstration planned, articulated in two mains steps: the deployment of the fixed target performed on January 26 th 2019, and harpoon firing performed on February 8 th 2019. The Harpoon capture demonstration included three Go/NoGo reviews: 1st to confirm the opportunity window before pretesting platform, 2nd to confirm the correct platform status and allow target deployment, and 3rd to confirm the correct deployment and stability of the target to collegially agree for a Go for harpoon firing.

The main phases of the demonstration were: i) Predemonstration tests (cameras testing), ii) software uploading for HTA target deployment, iii) platform pointing (Earth pointing, as it was requested an alignment with Earth background to better see the harpoon and tether during firing), iv) software uploading for target deployment, v) target deployment, vi) images acquisition to verify correct target deployment, vii Software upload for actual demonstration, viii) harpoon firing.

The decision to fire the harpoon against a fixed target, rather than a free flying one, reduced the unnecessary complexity as only the firing mechanism and harpoon flight/impact needed to be tested. In turn this also increased the payload safety. The dimensions of the elements is reported in Table, and the material (Aluminum honeycomb) was selected specifically to be representative of the composite structures used on large satellites that are potential targets for removal. Two supervision cameras with different parameters were set to record the experiment:100fps / Narrow FOV (40%), 40 fps / Large FOV (100%).

After the demonstration, the retraction of the boom would have been considered as this could have interfered with S-Band HR communication. Two weeks were initially allocated for data downloading. As planned, after the target had been deployed, images were downloaded to the ground station to check the correct positioning and stability of the target. The images showed that during the boom deployment, there was an increasing rotational oscillation of the target (see Figure 21 and Figure 22). This was initially attributed to the the vibration produced by the motor uncoiling the boom, and the low rotational stiffness of the boom; however, even after the deployment was completed, the oscillations continued. Subsequent imagery still showed significant oscillations of the target (up to +/-10 degrees), as shown in Figure 21. The cause of the oscillation was deemed to be the action of the AOCS of the satellite, amplified by the low resonance frequency (and low damping) of twisting mode of the boom-target assembly The mid position of the target was consistent with the images taken during alignment tests on the ground, and a more detailed analysis of the movement showed that even at the maximum of the amplitude of the oscillation a portion of the target was in the line of sight of the harpoon. Therefore provided that the harpoon hit the target within a 10mm radius from its line of sight, the experiment could have been performed. However, it was decided to attempt to stabilize the target reducing the mechanical inputs produced by the AOCS. This was reconfigured in order to reduce the actions of the platform actuators to the minimum. After a period to allow the natural damping of the system to reduce the oscillations, the target appeared stable with only very minor movements (<0.5deg) and therefore it was possible to proceed with the demonstration with minimum risk. Onboard the harpoon deployment chamber were two cold gas generators' (CGG) heaters, these were thermostatically controlled and kept the CGGs at working temperature (between 10°C 40° C).

A preprogrammed command sequence was uploaded to open the harpoon door and fire the projectile. The harpoon was fired as planned and successfully captured the target, achieving a firing speed of 19 m/s. Onboard camera footage recorded the firing performance, with the projectile attitude, accuracy and speed all shown to be consistent with the performance data acquired on ground. A snapshot of the firing sequence can be seen in During tests on the ground it appeared that the mechanical shock produced when the harpoon hit the target (mainly due to the high velocity of the harpoon) was so high that it snapped the target off of the tip of the boom. To address this issue a Delrin clock spring between the tip of the boom and the target was considered, in order to reduce the peak stress on the materials upon impact. This device was successfully ground tested as shown in Figure 25. However as the harpoon is imbedded in the target and tethered to the mothercraft, the separation of the target (braking off from its supporting structure) does not pose particular risks. In addition the spring element added complexity and flexibility to the assembly and therefore it was decided to proceed with a flight model that did not include such a feature. When the harpoon hit the target, the end of the deployable boom snapped, and the target coming off is visible in Figure 26, where it is also possible to see that the harpoon is tethered and therefore the target has been safely captured.

Eventually, after floating for a while in space tethered to the mothercraft, the target and harpoon ended up wrapped around the boom as visible in

VI. DRAGSAIL

Like for any active debris removal mission, the last phase of the mission would be to dispose of the debris, in this case the fourth and final demonstration was the dragsail. This is a drag augmentation device that, thanks to the residual atmosphere in low earth orbit, slows down the spacecraft, progressively reducing the altitude of its orbit until the craft burns in the high atmosphere.

In order to avoid interference with the other payloads, and in particular with the deployable boom of the HTA, (as a malfunction might have left the boom deployed thus preventing the full deployment of the dragsail) the dragsail was mounted on the other side of the spacecraft. The decision to mount the dragsail so that its deployment was in the opposite direction to all the other payloads meant that none of the cameras could actually observe the deployment of the device. Therefore the successful deployment of the dragsail had to be indirectly confirmed by other means, namely, i) examination of the power profile provided to the dragsail during activation, ii) reduction in power generated by the solar panels due to shadowing of the sails, changes in telemetry from the sun sensor, iii) changes in the two-line element set (change in the slope of the altitude versus time curve), iv) observation from the ground (significant changes in brightness).

The command to deploy the dragsail was issued on the 4th March 2019. The first telemetry available was the power profile from the craft, and this was consistent with the electronics operating as expected up to the point of sail deployment. Sufficient and expected power was provided to the burn wire system. Sufficient and expected power was consumed by the cold gas generators used to inflate the mast, however the power profile for the motor during deployment was consistent with a stall condition and a repeated retry. Some observations of the brightness of the object from ground telescope/radar showed some increase of brightness (the brightness fluctuates as it depends on the orientation of the object as it reflects the light, so accurate measurements are often difficult) but the increase in magnitude of the brightness was not consistent with a full deployment of the sail.

No significant changes in the power of the solar arrays or solar sensor signals were detected.

Finally the analysis of the altitude of the orbit of the object showed that there were no significant changes.

Figure 28 shows the 405km altitude orbit of the ISS, the orbit of the mothercraft (slowly decaying and going from the 405km at the moment of deployment from the ISS (20 th June 2018) to 397km at the end of March 2019), and the decaying orbits of the two cubesats. The orbit of DSAT#1, separating from the mothercraft on the 18 th of September 2018, is coming down very quickly as the object produces significant drag (due to the deployable structure and being entangled by the net) so that the object re-entered in the atmosphere the 2 nd of March 2019. The initial prediction was the 4 th of March. The orbit of DSAT#2 is also visible separating from that of the mothercraft on the 28 th of October 2018 and drifting down at a faster rate than the mothercraft, but not as quickly as DSAT#1. The expectation was that deploying the dragsail on the 4 th of March 2018 from here the orbit of the mothercraft would started to decay more rapidly. However, no significant change is visible in the trajectory of the craft.

All the data available is consistent with a partial deployment or the sail not coming out from the container (see Figure 29). The device should have deployed like that shown in Figure 30.

At this stage it is not possible to say whether the issue was caused by the inflatable boom not deploying, and because the sail was still in its container the CFRP boom could not deploy (as constrained by the side walls of the container), or if the sail got tangled and therefore the CFRP booms could not extend as they were constrained by the sail.

VI.1. DragSail further developments

The lesson learned from the design, manufacturing, assembly and integration of the DragSAIL for RemoveDEBRIS were put into practice in the development of the DragSail for the InflateSAIL cubesat ([20] [START_REF] Underwood | InflateSail De-Orbit Flight Demonstration Results and Follow-On Drag-Sail Applications[END_REF]) and two new dragsails that were designed and built for the Space Flight Industries for their SSO-A Mission [START_REF]SSO-A mission update: upper and lower free flyers deorbiting as planned -press release[END_REF].

Although the Dragsail for InflateSail was (successfully) activated at the end of June 2017 (deorbiting the cubsat from an orbit of approximately 500km in approximately 3 months), which is well before the activation of the RemoveDEBRIS' dragsail, the mission benefitted from the lesson learned from the RemoveDEBRIS Dragsail development as well as improvements in the testing processes. The two further Dragsails for the SSO-A mission also benefitted from the RemoveDEBRIS experience. In this case the assembly (shown in Figure 31) did not contain the inflatable boom, and proof of the successful deployment has been given by a significant increase of brightness of both craft that were equipped with the new DragSails, radar cross section, and by an enhanced drag effect. Note that in this mission, the mass of the craft and the altitude of the orbit are both much higher than those of RemoveDEBRIS, and therefore the effect on the orbit is not as pronounced as it would have been for RemoveDEBRIS.

Re-entry of the first of the two craft expected in approximately three years. The first demonstration, the Net capture was successful. The Net tested, fully representative of the operational design, succeeded to capture the target cubsat that had inflated some structures to provide a more realistic representation (in terms of size) of a potential target space debris. There were some deviations from the expected behavior of some elements of the demonstrations, i.e. relatively fast spinning of the target and distance of the capture, however the images acquired have confirmed the proper working of the device. Indeed before using this technology for the capture of real large space debris the hardware will need scaling up and potentially modifications of mesh size or number of masses will be considered. Indeed once captured the debris will need to be secured to the mothercraft in order to be disposed of. This could be achieved by having the net tethered to the mothercraft in order for the debris to be actively towed (thus giving some control in the disposal), or slowed down by the action of a Dragsail on the mothercraft (without active control of the disposal). Further work will be necessary to verify the viability of these (or other) options for the disposal. However, what is important, is that the overall success of the demonstration, established this technology as a variable candidate for capturing large space debris.

The VBN cameras and algorithm performed very well, being the most computationally complex of the four demonstrations, the cameras acquired all the images expected of the target with different lighting conditions, background and range. The database of images collected and ground truth associated is unique, and of prime importance. It enabled us to assess the sensors robustness (camera & Lidar) in real conditions and to assess algorithms performance and robustness for future non collaborative rendezvous on real images.

The harpoon capture experiment was also successful, with the device firing the harpoon with appropriate speed and this centering and capturing the target as planned.

The target was snapped off of the supporting boom, floated away temporarily, retained by the tether line thus demonstrating the correct working of the harpoon & tether concept. Also, this technology has been proven as viable for the capture of large space debris, and, as for the net, further investigations will be necessary to assess the best methodology for the following disposal of the debris that has been captured.

For the final experiment, the DragSail, the absence of monitoring cameras or detailed telemetry prevented complete conclusions on the cause and extent of the malfunction of the device. At the time, this was disappointing, however the lesson learned in developing the hardware has already been put into practice in another three devices (the InflateSail cubesat, and two DragSails provided to a commercial organization) that have been successfully demonstrated in orbit. Hence, as the ultimate purpose of the RemoveDEBRIS in flight demonstrations was to pave the way for industrial/commercial development and exploitation, also the last experiment can be considered successful as it enabled the development of commercial devices.

Besides the technical achievement RemoveDEBRIS has attracted significant media attention, and its various successes have been reported by all the major world news outlets. This has contributed to raise awareness of the issues in the public and policy makers, creating a momentum that hopefully it will help convince the relevant authorities to support the delivery of missions to remove the some of the large debris currently orbiting around the Earth.

Figure 1 :

 1 Figure 1: RemoveDEBRIS Satellite. Left: Satellite in the SSTL cleanroom, 14 December 2017. Right: Satellite with protective panels installed.

Figure 2 :

 2 Figure 2 : Left: the clam shell (foam) that encompasses a RemoveDebris structural model, 2017. Right: loading of the Dragon capsule with a cargo bag. Credit: NASA, SpaceX from [44].

Figure 3 Figure 4 :

 34 Figure 3 Dragon capsule. Left: view of the capsule approaching the ISS and captured by the ISS robotic arm. Right: View from the top hatch of the inside of the cargo bay of the Dragon capsule, showing CTB on the sides of the capsule. Credit: NASA, SpaceX

Figure 5 :Figure 6 :

 56 Figure 5: RmoveDEBRIS with its protective panels, is removed fomr the foam clamshell

Figure 7 :

 7 Figure 7: Release in orbit of RemoveDEBRIS

Figure 8 :

 8 Figure 8: RemoveDEBRIS free flying, image taken by the ISS

Figure 9 :

 9 Figure 9: Examples of pictures used for the calibration of supervision camera and VBN.

Figure 10 :

 10 Figure 10: Launch Sequence. This figure shows the launch sequences for the mission to the International Space Station (ISS). Credit: SpaceX, NanoRacks, NASA [40].

Figure 11 :

 11 Figure 11: DSAT#1 Left DSAT#1 in stowed configuration drifting away from mothercraft. Right: DSAT#1 with lateral inflatable booms deployed

/Figure 12 :

 12 Figure 12: Launch of the net towards DSAT#1 Top: Net opening pulled by the 6 throw masses, Muddle: the net continues to stretch and beyond it DSAT#1 is visble, Bottom: moment of the capture, DSAT#1 entangled in the net.

Figure 13 :

 13 Figure 13: Left, VBN color camera Right LIDAR images taken at 4 seconds interval during VBN experiment Each color on the LiDAR image represents a distance. On the panel A of the DSAT#2, the sunlight reflected into the LiDAR detector is such that the distance measurement is noisy (many different colors on close neighbor pixels). On the panel B, which is protected against specular sunlight reflection by panel C, the colors are much more uniform.

Figure 14 :

 14 Figure 14: Left DSAT#2 in the FoV of both standard camera and LIDAR, right DSAT#2 still in the FoV of standard camera, but out of the LIDAR FoV.

Figure 15 :

 15 Figure 15: Comparison of distance as measured from GPS (continuous line) and LIDAR discrete circles.

Figure 16 :

 16 Figure 16: Superposition of color camera and LiDAR images with earth in the background

Figure 17 :

 17 Figure 17: DSAT#1 Left: As seen by the camera, Right CAD modelThe Airbus in house developed tracking algorithm has been tested on the VBN camera images. This solution had to be improved in order to deal with the low acquisition frequency of the VBN camera (0.1 Hz), showing good results as long as the background is dark. This is illustrated by the green contours illustrate the 3D model of the target which has been matched by the algorithm to the actual DSAT#2 image.

Figure 18 :

 18 Figure 18: DSAT2 within VBN camera field of view over the reference trajectory

Figure 19 :

 19 Figure 19: Comparison of different relative distance estimations between DSTA#2 and RemoveDEBRIS mothercraft.

Figure 20 :

 20 Figure 20: VBN images of DSAT#2 for the removedebris sequence (Inria Algorithm). The cubesat has been successively tracked and localized on 210 successive images (that is, during 210 s))

Figure 21 :

 21 Figure 21: rotational oscillations of the target a) during boom deployment, b) and c) steady state oscillations with the boom completely deployed

Figure 23 :

 23 Figure 23: Harpoon impacting the target.

Figure 23 ,

 23 Figure23, and Figure24shows that the harpoon has actually hit the target as expected, in the centre.

Figure 24 :

 24 Figure 24: Harpoon imbedded in the target

Figure 25 :

 25 Figure 25: Spring/shock absorber, sequence showing the target being hit by the harpoon, rotating thus reducing the shock on the boom and finally returning to its original position

Figure 26 :

 26 Figure 26: Boom snapping off and target coming off, tethered to the mothercraft

Figure 28 :

 28 Figure 28: Top: orbital decay of the various elements of the RemoveDEBRIS mission. Bottom: detail of the curves highlighting ISS orbit raise maneuvers and decay of the CubeSat's and mothership

Figure 29 :

 29 Figure 29: DragSail and schematic of the device in its container mounted on RemoveDEBRIS

Figure 31 :

 31 Figure 31: Deployment test of the DragSail for the SSO-A mission

Table 1 : RemoveDebris On-orbit operations responsibilities.

 1

	Partner	Responsibility
	SSTL	Mission Command & Control,
		Commissioning -LEOP
	Ariane Group	System Engineering
	Airbus Germany	Net
	(support from SSC)
	Airbus France (support	Overall VBN † experiment
	from CSEM Inria)	and algorithms
	Airbus UK (support	Harpoon
	from SSC)
	SSC	DragSail

† vision-based navigation

Table 2 :

 2 Net demonstration

	Parameter	Value
	Solar angle for Net	|| > 40°
	opportunity
	NETCAM mass	6kg
	Net size	5m diameter
	Platform Pointing	Open loop, stability +/5
	Target:	2U Cubesat (DSAT#1) Semi passive
		with inflatable booms to become 1m
		size
	DSAT#1 deployment	V=5cm/s +/20%
	speed
	Target range for	7m +/-1m
	capture
	Complete	440s (DSAT#1 at ~300m distance)
	experiment duration
	Image capture:	B&W, FOV=60°,
	Supervision cameras

Table 2 :

 2 Harpoon experiment data

	Parameter	Value
	Harpoon Target	Mass: 4.3 Kg
	Assembly
	Harpoon (projectile) mass 0.115kg
	Target size &	approx. 10 by 10cm, Aluminium
	material	Honeycomb panel
	Target	fixed at the end of 1.5m (+/-10%)
		deployable boom
	Harpoon velocity	Nominal 20m/s
	Planned image	51s with supV0 (100fps),
	capture	71s with SupV1 (40fps) and
		5820s with VBN (2D&3D)

th March

ACKNOWLEDGEMENTS

This research is supported by the European Commission FP7-SPACE-2013-1 (project 607099) 'RemoveDebris -A Low Cost Active Debris Removal Demonstration Mission', a consortium partnership project consisting of: Surrey Space Centre (University of Surrey), SSTL, Airbus GmbH, Airbus SAS, Airbus Ltd, Ariane Group, Innovative Solutions in Space (ISIS), CSEM, Inria, Stellenbosch University.

The consortium would like also to thank NanoRack for their support to the launch operations.


Sondes Metoui
email: sondes.metoui@ensam.eu

Etienne Pruliere
email: etienne.pruliere@ensam.eu

Amine Ammar
email: amine.ammar@ensam.eu

Frédéric Dau
email: frederic.dau@ensam.eu

Ivan Iordanoff
email: ivan.iordanoff@ensam.eu

A multiscale separated representation to compute the mechanical behavior of composites with periodic microstructure

Keywords: Model reduction, Multiscale simulations, Proper Generalized Decomposition, Composite structures

come

A multiscale separated representation to compute the mechanical behavior of composites with periodic microstructure

Introduction

One of the main challenges in mechanics and engineering is to account for physical phenomena that occur at different scales. A coupling between scales is often observed, generating a real need for multiscale models in many applications. In composite materials, for example, there are at least two or three characteristic scales: the fiber scale, the ply scale and the laminate scale. A major difficulty related to multiscale modeling is the need of multiscale solvers that require a lot of computational resources. Therefore, there is a real need for computational methods able to reduce the cost of such simulations.

In this paper, we focus on problems related to structural mechanics of composite materials, that results in elliptic partial differential equations. The multiscale solvers available in literature can be divided into two main classes:

• The first class considers a set of microscopic volumes, generally representative volume element (RVE). For example, one RVE can be defined in each Gaussian point of the macroscopic finite element mesh. The size of the RVE is often much smaller than the size of the macroscopic finite elements i.e. the fine-scale mesh covers only a very small part of the macroscopic domain. Some well known methods are the Heterogeneous Multiscale Methods (HMM) [START_REF] Engquist | The heterogeneous multiscale methods[END_REF] or the FE 2 methods [START_REF] Feyel | Multiscale FE2 elastoviscoplastic analysis of composite structures[END_REF] in engineering literature. This kind of methods is generally based on an extension of the homogenization theory for non linear materials.

• The second class take into account the microscale (solution) in each point of the macroscopic domain. The fine-scale mesh covers the whole macroscopic computational domain. In this case each fine-scale degree of freedom (DOF) has to be treated at least one time in the algorithm. Therefore, these methods are generally more expansive than the method belonging to the first class. In this class, we find in particular the multigrid methods and the Domain Decomposition Methods (DDM) [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF][START_REF] Lions | On the Schwarz alternating method. I[END_REF][START_REF] Lions | On the Schwarz alternating method[END_REF][START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF].

The approaches belonging to the first classes are efficient in a computational point of view. However, they assume a well separated definition of scales that is not always relevant, in particular when considering the characteristic length of composite materials. Furthermore, this kind of methods requires to introduce some additional hypothesis to make the link between the microscale and the macroscale (generally boundary conditions applied to the microscopic domain [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]). For periodic microstructures, the natural choice is to assume periodic boundary conditions on the RVE but this lead to a loss of precision near the boundaries of the domain.

The approach presented in this paper belongs mainly to the second class in the sense that the whole problem is described. However, the objective is to reduce the number of degrees of freedom as in the methods belonging to the first class, by the use of some additional hypothesis.

The main hypothesis is that only materials with periodic structures are considered with an obvious application in composite materials. The periodicity of the microstructure enables to separate two scales: the scale of the periodic cell and a macroscopic scale. The idea is to use this hypothesis to build the full solution assuming a separated form in the context of the Proper Generalized Decomposition (PGD). With the PGD solver, the number of degrees of freedom to treat can theoretically be reduced by several orders of magnitude.

The PGD philosophy was originally proposed by Ladeveze with space-time problems in the context of the LATIN (LArge Time INcrement) method, which was called radial decomposition [START_REF] Ladeveze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF]. Ammar et al. have generalized this method for multi-dimensional problems in the context of the kinetic theory description of complex fluids [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part II: Transient simulation using space-time separated representations[END_REF]. This method enables the reduction in size of multidimensional and parametric problems [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Nouy | Generalized spectral decomposition for stochastic nonlinear problems[END_REF][START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF]. The idea of the PGD is to build an approximation of the solution based on a separated representation. Assume that we are looking for a field u : Ω → R that depends on some coordinates x j for j = 1, . . . , D where the domain Ω ⊂ R D is defined by:

Ω = Ω x 1 × Ω x 2 × • • • × Ω x D with Ω x i ⊂ R for i = 1 . . . D.
A separated representation of u can be defined as a sum of N products of some functions

F j i (x j ) : Ω x j → R such as: u(x 1 , . . . , x D ) ≈ N ∑ i=1 F 1 i (x 1 ) × • • • × F D i (x D ). (1)
All these functions F j i are a priori unknown. This representation is injected into the weak formulation of the partial differential equations depending on the considered problem. This leads to non-linear equations that are solved iteratively. Separated representations have been used for a wide range of problems. In particular, the PGD has been applied to compute full 3D solutions using an in-plane/out-of-plane separated representation in composite laminates [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF] and for the modeling of laminated and sandwich composite plates [START_REF] Vidal | Composite beam finite element based on the proper generalized decomposition[END_REF][START_REF] Vidal | Proper Generalized Decomposition and layer-wise approach for the modeling of composite plate structures[END_REF]. In our previous works, we have employed this method to simulate delamination of composite laminates [START_REF] Metoui | The proper generalized decomposition for the simulation of delamination using cohesive zone model[END_REF]. For the sake of simplicity, the method will be described firstly for a simple problem defined in a 1D domain Ω = [0, L]. The 2D and 3D cases will be addressed in the next section.

Then, an 1D problem with a periodic geometry and periodic material properties is considered in this section. The weak formulation of the static equilibrium equation for a beam in traction/compression with an elastic material is:

∫ L 0 AE dU ⋆ (X ) dX dU (X ) dX dX = [ AE U ⋆ (X ) dU (X ) dX ] L 0 ∀U ⋆ ∈ V (2)
E is the elastic modulus, A is the area of the beam, X is the coordinate along the beam axis, U ∈ H 1 (Ω ) is the longitudinal displacement, and

U ⋆ ∈ V = { f ∈ H 1 (Ω )} is a test function.
L denotes the length of the beam. We also define: [ f (X )] L 0 = f (L) -f (0) for any function f . The right hand term is related to the boundary conditions on the left (X = 0) and on the right (X = L) of the domain. To enforce the boundary conditions, a penalty method will be used. This method will be described in Section 2.1.4. For now, only the left hand term of the equation will be considered.

In the finite element method, this weak form is approximated using some shape functions over each element of a mesh. The proposed strategy takes advantage of the periodicity of the geometry and material properties to build a periodic mesh. The periodicity assumption, in this context is that the domain may be sliced into identical elementary parts with: E(X + ∆L) = E(X ) and A(X + ∆L) = A(X ). ∆L is the length of each part. If loads are added to this problem, the periodicity is not required.

Then, the mesh can be built as the sum of identical mesh with common nodes on the edge as depicted in Fig. 1 for the 1D mesh. Each occurrence of the periodic cell is associated to an integer denoted k as shown in Fig. 1. k is an integer such as k ∈ {1, . . . , N k }, and N k = L ∆L denotes the number of cells. N x denotes the number of nodes in the periodic mesh.

The coordinate X defining the horizontal position in the beam is written using k:

X = (k -1) ∆L + x (3)
x ∈ [0, ∆L] is the position in each part. The continuous form of the displacement field found using the separated approximation of the solution is given by:

U (X ) = U ((k -1) ∆L + x) = n ∑ i=1 F i (x)G i (k) (4)
where F i (x) : [0, ∆L] → R and G i (k) : {1, . . . , N k } → R are some unknown functions. In the following, the notation U (k, x) = U ((k -1) ∆L + x) will be used.

Using k and x instead of X leads to a double definition of the position at the edge of the elementary parts: U (k, ∆L) = U (k + 1, 0). This is problematic from a numerical point of view because it leads to a multiplication of degrees of freedom on these points. There are two possibilities to treat the edge of the elementary part:

1. Solving the problem using more degrees of freedom than necessary and using boundary conditions to enforce the continuity: U (k, ∆L) = U (k + 1, 0). 2. Deleting some nodes when discretizing the domain in order to suppress the non-necessary degrees of freedom.

The first method requires to enforce some continuity conditions at the interface of each cell. With the separated representation of the solution Eq. ( 4), the continuity conditions give:

n ∑ i=1 F i (∆L)G i (k) = n ∑ i=1 F i (0)G i (k + 1). ( 5
)
As the separated representation is generally built term by term, the continuity must be satisfied for each term of the sum. This leads to:

F i (∆L)G i (k) = F i (0)G i (k + 1) ∀(i, k) ∈ {1, . . . , N k } × {1, . . . , n}. (6)
Thus, if F i (0) ̸ = 0 the function G i is completely defined from the values of F i (0), F i (∆L) and G i (1):

G i (k) = ( F i (∆L) F i (0) ) k-1 G i (1).
This constraint on G i decreases dramatically the convergence rate of the PGD. To circumvent this problem, a possibility is to compute many terms of the separated approximation at the same time instead of a classical term by term algorithm.

In the following, only the second strategy is described.

Let us consider a finite element discretization of the domain with equally distributed nodes as shown in Fig. 1. This distribution of nodes is only for sake of clarity without loss of generality. The size of an element is denoted ∆x. Then, the weak formulation Eq. ( 2) can be rewritten:

N k ∑ k=1 (∫ k∆L-∆x (k-1)∆L AE dU ⋆ (X ) dX dU (X ) dX dX ) + N k -1 ∑ k=1 (∫ k∆L k∆L-∆x AE dU ⋆ (X ) dX dU (X ) dX dX ) + ∫ N k ∆L N k ∆L-∆x AE dU ⋆ (X ) dX dU (X ) dX dX = [ AE U ⋆ (X ) dU (X ) dX ] L 0 ∀U ⋆ ∈ V. (7)
The first sum

∑ N k k=1 ( ∫ k∆L-∆x (k-1)∆L AE dU ⋆ (X ) dX dU (X ) dX dX
) contains integrals defined on the elements that are inside the elementary parts. The second sum

∑ N k -1 k=1 ( ∫ k∆L k∆L-∆x AE dU ⋆ (X ) dX dU (X ) dX dX
) contains integrals defined over the interface elements. The element located on the right boundary needs a special treatment that will be detailed in Section 2.1.4.

With a change of variable between X and x it comes:

I 1    N k ∑ k=1 (∫ ∆L-∆x 0 AE dU ⋆ (k, x) dx dU (k, x) dx dx ) + I 2    N k -1 ∑ k=1 (∫ ∆L ∆L-∆x AE dU ⋆ (k, x) dx dU (k, x) dx dx ) + ∫ N k ∆L N k ∆L-∆x AE dU ⋆ (X ) dX dU (X ) dX dX = [ AE U ⋆ (X ) dU (X ) dX ] L 0 ∀U ⋆ ∈ V. (8)

First iteration

The approximation defined by Eq. ( 4) is built term by term. For now, we focus on the first iteration. The displacement is then approximated by:

U (X ) = F 1 (x)G 1 (k) = R(x)S(k). (9)
To avoid the redundant use of subscripts, we denote R = F 1 and S = G 1 . The determination of R and S involves a non-linear problem that is solved using the classical alternate direction strategy [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF]. At the beginning, R is computed assuming a random value for S. Then S is computed knowing R. And again R is computed knowing S and so on until R and S have converged.

Then, two different problems must be treated:

1. Computing R knowing S 2. Computing S knowing R. The test function is U ⋆ (k, x) = R ⋆ (x)S(k) for the first problem and U ⋆ (k, x) = R(x)S ⋆ (k) for the second problem where R ⋆ (x) ∈ H 1 ([0, ∆L]) = V R and S ⋆ (k) : {1, . . . , N k } → R are some associated test functions.
The first problem is considered in the following. The first integral of Eq. ( 8) can be simply rewritten using the separated approximation:

I 1 = N k ∑ k=1 ∫ ∆L-∆x 0 AE dR ⋆ (x) dx dR(x) dx S 2 (k) dx = (∫ ∆L-∆x 0 AE dR ⋆ (x) dx dR(x) dx dx ) N k ∑ k=1 S 2 (k). ( 10
)
After a finite element discretization it remains:

I 1 = ( R ⋆T KR ) × ( S T I N k S ) (11)
R is the column vector containing the nodal values of R:

R = [ R 1 R 2 • • • R N x ] T (12)
S is the column vector containing the values of S:

S = [ S 1 S 2 • • • S N k ] T (13)
R * is a vector containing the nodal values of the test function R * . K is the stiffness matrix related to the periodic cell and

I N k is the N k × N k identity matrix.
The second integral in Eq. ( 8) requires a little more development since the degree of freedom corresponding to U k (∆L) = U k+1 (0) is defined only on the part k + 1. This integral is defined on the interface elements between two parts. For 1D linear elements the shape function vector over the interface element is:

Φ(x) = [ φ N x (x) φ 1 (x) ]
(see Fig. 1) where φ i (x) denote the shape function associated to the node i. For 3 nodes quadratic elements the shape function vector becomes:

Φ(x) = [ φ N x-1 (x) φ N x (x) φ 1 (x)
] . The matrix of DOF is for linear elements:

Q k = [ R N x S k R 1 S k+1 ] (14)
or for quadratic elements:

Q k = ⎡ ⎣ R N x -1 S k R N x S k R 1 S k+1 ⎤ ⎦ . ( 15
)
Using this finite element approximation, I 2 can be written as:

I 2 = N k -1 ∑ k=1 ∫ ∆L ∆L-∆x AE dU ⋆ (k, x) dx dU (k, x) dx dx = N k -1 ∑ k=1 Q ⋆T ∫ ∆L ∆L-∆x AE dΦ(x) dx T dΦ(x) dx dxQ. ( 16
)
Defining

α i j = ∫ ∆L ∆L-∆x AE dφ i dx dφ j
dx dx and developing the previous equation in the case of linear elements, I 2 becomes:

I 2 = R ⋆ N x α N x N x R N x N k -1 ∑ k=1 S 2 k + R ⋆ N x α 1N x R 1 N k -1 ∑ k=1 S k S k+1 + R ⋆ 1 α 1N x R N x N k -1 ∑ k=1 S k+1 S k + R ⋆ 1 α 11 R 1 N k -1 ∑ k=1 S 2 k+1 . (17)
In general, four operators may be defined as follows:

I 2 = ( R ⋆T M N x ,N x R ) × ( S T D k,k S ) + ( R ⋆T M N x 1 R ) × ( S T D k,k+1 S ) + ( R ⋆T M 1,N x R ) × ( S T D k+1,k S ) + ( R ⋆T M 1,1 R ) × ( S T D k+1,k+1 S ) (18)
where D k,k , D k+1,k+1 , D k+1,k and D k,k+1 are the following N k × N k square matrices:

D k,k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 0 • • • 0 0 1 0 • • • 0 . . . . . . . . . . . . . . . 0 • • • 0 1 0 0 • • • 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ D k+1,k+1 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 • • • 0 0 1 0 • • • 0 . . . . . . . . . . . . . . . 0 • • • 0 1 0 0 • • • 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (19)
D k,k+1 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 0 • • • 0 0 0 1 . . . . . . . . . . . . . . . . . . 0 0 • • • 0 0 1 0 • • • 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ D k+1,k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 • • • 0 1 0 0 • • • 0 0 1 0 . . . . . . . . . . . . . . . . . . 0 0 • • • 0 1 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (20)
and

M N x ,N x , M 1,1 , M N x ,1
and M 1,N x are some N x × N x square matrices coming from the development of Eq. ( 16). In the case of linear elements, these matrices are identified from Eq. ( 17):

M N x ,N x = ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 • • • 0 . . . . . . . . . . . . 0 • • • 0 0 0 • • • 0 α N x N x ⎤ ⎥ ⎥ ⎥ ⎦ M 1,1 = ⎡ ⎢ ⎢ ⎢ ⎣ α 11 0 • • • 0 0 0 • • • 0 . . . . . . . . . . . . 0 • • • 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ (21)
M N x ,1 = ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 • • • 0 . . . . . . . . . . . . 0 0 • • • 0 α 1N x 0 • • • 0 ⎤ ⎥ ⎥ ⎥ ⎦ M 1,N x = ⎡ ⎢ ⎢ ⎢ ⎣ 0 • • • 0 α 1N x 0 • • • 0 0 . . . . . . . . . . . . 0 • • • 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ . ( 22
)
In the separated formulation some operators are non symmetric. However it is interesting to notice that the global problem remains symmetric because: D k+1,k = D T k,k+1 and M 1,N x = M T N x ,1 . The stiffness matrix K defined in Eq. ( 11) is also symmetric.

The weak formulation Eq. ( 8) without boundary gives eventually after eliminating R * :

(KR) × ( S T I N k S ) + ( M N x ,N x R ) × ( S T D k,k S ) + ( M N x 1 R ) × ( S T D k,k+1 S ) + ( M 1,N x R ) × ( S T D k+1,k S ) + ( M 1,1 R ) × ( S T D k+1,k+1 S ) = B 1 (R, S). ( 23
)
And the system for the second problem, that is to compute S knowing R writes:

( R T KR ) × ( I N k S ) + ( R T M N x ,N x R ) × ( D k,k S ) + ( gR T M N x 1 R ) × ( D k,k+1 S ) + ( R T M 1,N x R ) × ( D k+1,k S ) + ( R T M 1,1 R ) × ( D k+1,k+1 S ) = B 2 (R, S) (24)
B 1 (R, S) and B 2 (R, S) are vectors associated to boundary conditions that are defined in Section 2.1.4. The linear systems Eqs. ( 23) and ( 24) are some reduced version of a more global system (i.e, the full FEM system). The nodal values of the global solution U (X ) = R(x)S(k) can be built a posteriori using a simple tensor product U = R ⊗ S. These two systems are then well-posed if the equivalent global problem is also well-posed and if the global solution U computed from S and R satisfies the boundary conditions. Algorithm 1 General PGD algorithm used to make the link between the elementary cell and the global mesh 1: n = 0 2: initialize F n+1 and G n+1 to random values 3: F old n+1 = F n+1 ; G old n+1 = G n+1 4: Compute F n+1 knowing G n+1 (update the function related to the elementary cell) 5: Compute G n+1 knowing F n+1 (update the function related to the global mesh)

6: if max(   F n+1 -F old n+1   2 ,   G n+1 -G old n+1   2 ) > ϵ 1 max(∥F n+1 ∥ 2 , ∥G n+1 ∥ 2 ) then go to 3 7: n = n+1 8: U (x, k) = ∑ n i=1 F i (x)G i (k) (no need to compute U explicitly) 9:
If the normalized Euclidean norm of the residual of the system (Eq. 26) > ϵ 2 then go to 2

Other iterations

For other iterations, F i and G i are assumed known for i = 1, . . . , n. Now, we are looking for R = F n+1 and S = G n+1 such as:

U (x, k) = n ∑ i=1 F i (x)G i (k) + R(x)S(k). ( 25
)
The test function is the same as in the first iteration:

• U ⋆ (k, x) = R ⋆ (x) S(k) if the unknown is R. • U ⋆ (k, x) = R(x) S ⋆ (k) if the unknown is S.
We focus on the case where the unknown is R. Introducing this expression into the weak form, and using a finite elements discretization as described for the first iteration we get:

(KR) × ( S T I m S ) + ( M N x ,N x R ) × ( S T D k,k S ) + ( M N x 1 R ) × ( S T D k,k+1 S ) + ( M 1,N x R ) × ( S T D k+1,k S ) + ( M 1,1 R ) × ( S T D k+1,k+1 S ) = B 1 (R, S) - n ∑ i=1 [ (KF i ) × ( S T I m G i ) + ( M N x ,N x F i ) × ( S T D k,k G i ) + ( M N x 1 F i ) × ( S T D k,k+1 G i ) + ( M 1,N x F i ) × ( S T D k+1,k G i ) + ( M 1,1 F i ) × ( S T D k+1,k+1 G i )] (26)
where F i and G i contain the nodal values of F i (x) and the values of G i (k). A very similar system can be easily written when the unknown is S.

The global PGD algorithm is summed up in algorithm 1. In this algorithm ϵ 1 and ϵ 2 are error tolerances related to the convergence criteria. The convergence of this algorithm has been discussed in [START_REF] Ammar | On the convergence of a greedy rank-one update algorithm for a class of linear systems[END_REF].

Boundary conditions

For now, the method has been described without accounting for boundary conditions. As the method is based on a finite element discretization, stress conditions can naturally be introduced by adding a vector of generalized force to the system to be solved. The only difference with the classical finite element method is that it must be written on a separated form. The direct application of Dirichlet conditions on the separated representation is generally not possible. For instance, if the beam is clamped in x = 0 the boundary condition is:

U (k = 1, x = 0) = 0. ( 27
)
As the solution is built iteratively, the boundary conditions must be satisfied at each iteration, i.e.:

n ∑ i=1 F i (0)G i (1) = 0. ( 28
)
It is easy to prove that this leads to:

F i (0) = 0 or G i (1) = 0 ∀i ∈ {1, . . . , n}. (29)
If we enforce G i (1) = 0 ∀i ∈ {1, . . . , n}, the solution is zero over the first cell and we can get the true solution only in the very restrictive case where there is actually no displacement over this cell. In other hand, if we enforce F i (0) = 0 ∀i ∈ {1, . . . , n}, that leads to:

U (k, 0) = U ((k -1) ∆L) = 0 ∀k ∈ {1, . . . , N k }. (30)
In that case, the displacement vanishes on the left of each cell. In this case again, we cannot get the real displacement or only in some very restrictive cases.

Then, the best solution to apply Dirichlet conditions is to use a penalty method. It consists in adding some new operators that describe the boundary conditions. For instance, to enforce a unitary displacement on the left, the operators related to boundary conditions are:

M bc = β ⎡ ⎢ ⎢ ⎢ ⎣ 1 0 • • • 0 0 0 • • • 0 . . . . . . . . . . . . 0 • • • 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ D bc = ⎡ ⎢ ⎢ ⎢ ⎣ 1 0 • • • 0 0 0 • • • 0 . . . . . . . . . . . . 0 • • • 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ B x bc = β ⎡ ⎢ ⎢ ⎢ ⎣ 1 0 . . . 0 ⎤ ⎥ ⎥ ⎥ ⎦ B k bc = ⎡ ⎢ ⎢ ⎢ ⎣ 1 0 . . . 0 ⎤ ⎥ ⎥ ⎥ ⎦ where β ∈ R + is the purely numerical penalty parameter. M bc is a N x × N x matrix, D bc is a N k × N k matrix, B x bc is a N x × 1 matrix and B k
bc is a N k × 1 matrix. The boundary conditions are specified by defining B 1 and B 2 in Eqs. ( 23) and ( 24):

B 1 (R, S) = B x bc × ( S T B bc k ) -(M bc R) × ( S T D bc S ) B 2 (R, S) = ( R T B bc x ) × B k bc - ( R T M bc R ) × (D bc S) .
The system given in Eq. ( 23) becomes with this boundary condition:

(KR) × ( S T I m S ) + ( M N x ,N x R ) × ( S T D k,k S ) + ( M N x 1 R ) × ( S T D k,k+1 S ) + ( M 1,N x R ) × ( S T D k+1,k S ) + ( M 1,1 R ) × ( S T D k+1,k+1 S ) + (M bc R) × ( S T D bc S ) = B x bc × ( S T B k bc ) . (31)
This define a linear system whose unknowns are the component of R. This system may be solved by the conjugate gradient method or any other classical method.

Another problem has to be discussed relating to boundary conditions. On the right side, there is obviously no interface element. In comparison to other elementary parts, the right part has one node less on the right. This may be a problem. A simple solution to overcome this difficulty is to add a "virtual" elementary part on the right that is only used for one node. The other nodes are only virtual and the element inside this virtual part have to be ignored. In operators D k,k , D k,k+1 , D k+1,k and D k+1,k+1 , the row corresponding to the virtual part (the row N k ) must be filled with 0 in order to prevent the effect of virtual element. There is no undesirable numerical error related to the presence of virtual element even if there is no unique solution (the solution inside the virtual element can be anything). As detailed previously, the enrichment of the separated representation is performed using an alternate direction strategy that is a type of fixed point method. The solution will converge to one among the possible solutions according to the initial values of the vectors R and S (random values in general) [START_REF] Ammar | On the convergence of a greedy rank-one update algorithm for a class of linear systems[END_REF].

Domain decomposition method for periodic 2D and 3D domains

Mechanical model

In this section, a 2D or 3D static problem is considered.

U denotes the displacement vector expressed in the canonical basis:

U (X ) = ⎡ ⎣ U x (X ) U y (X ) U z (X ) ⎤ ⎦ (32)
where X = (X, Y, Z ) T is the coordinate vector. ε denotes the strain tensor using Voigt notation:

ε = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ε x x ε yy ε zz 2ε yz 2ε x z 2ε x y ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 33
)
The stress tensor using the Voigt notation is written as:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ x x σ yy σ zz σ yz σ x z σ x y ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = H ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ε x x ε yy ε zz 2ε yz 2ε x z 2ε x y ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = Hε (34)
H is the matrix describing the constitutive equation. In small displacements, ε is the symmetric gradient of the displacement. With these notations, the weak formulation of the equilibrium equation without dynamic effect and volume forces is expressed as:

∫ Ω ε(U ⋆ ) T (Hε(U )) dΩ = ∫ Γ U ⋆ .(σ.n) dΓ ∀U ⋆ ∈ ( H 1 (Ω ) ) 2 : (35)
Ω is the domain taken by the structure, Γ is the boundary of the domain, U ⋆ ∈ H 1 (Ω ) is a test function and n is the unit vector normal to Γ directed toward the exterior of Ω . The right part of this weak formulation is only used for boundary conditions. It is not included in the following for the sake of readability.

Separated description of 2D problems

As for the 1D case, the domain Ω is decomposed into N k periodic elementary cells Ω k . Ω i denotes the set of elements inside the elementary cell and Ω e is the set of interface elements. The integral over Ω in Eq. ( 35) can be decomposed as a sum of integrals over Ω k = Ω i ∪ Ω e (k). The weak form gives then for all U ⋆ ∈ ( H 1 (Ω ) ) 2 :

N k ∑ k=1 ∫ Ω i ε(U ⋆ ) T (Hε(U )) dΩ    I 1 + N k ∑ k=1 ∫ Ω e (k) ε(U ⋆ ) T (Hε(U )) dΩ    I 2 = ∫ Γ U ⋆ .(σ.n) dΓ . ( 36
)
To avoid the double definition of degrees of freedom at the interface between cells, some nodes need to be deleted at the interface. We choose to delete the nodes on the right and up faces for the 2D case as shown in Fig. 2. The local coordinates in the elementary cells are noted x. approximation associated to a 2D or a 3D problem is then:

U (x) = n ∑ i=1 F i (x)G i (k). ( 37
)
Here, G i is a scalar function and F i is a vector function:

F i (x) = ⎛ ⎝ F x i (x) F y i (x) F z i (x) ⎞ ⎠ . ( 38
)
Fig. 2. Examples of 2D meshes and decomposition in elementary parts.

First iteration and operators assembly

The integral I 1 is treated as for the 1D case. For the first iteration with U (x) = R(x)S(k), I 1 reads:

I 1 = (∫ Ω i ε(R ⋆ ) T (Hε(R)) dΩ ) N k ∑ k=1 S(k) 2 ≈ ( R ⋆T KR ) × ( S T I N k S ) R ⋆ ∈ ( H 1 (Ω i )
) 2 is a test function, K is the stiffness matrix related to Ω i and I N k is the N k × N k identity matrix. R is the column vector containing the nodal values of R x , R y and R z and S is the column vector containing the values of S.

The second integral I 2 needs an assembly of new matrices that are representative of interactions between neighboring cells. The element stiffness matrices of interface elements must be assembled in some global operators written from a tensor product of matrix operators either associated to the nodal values of Ω i or associated to the list periodic cells {1, . . . , N k } .

For instance, in the first 2D problem depicted in Fig. 2, I 2 can be written as:

I 2 = 4 ∑ m=1 ( R ⋆T M r m R ) × ( S T D r m S
) D r m are the matrices defined by: ) T . The second problem depicted in Fig. 2 is more complicated. There are now 3 types of interface elements:

( D r 1 ) i, j = { 1 if i = j
• elements between the cell and its right neighbor (subscript r ),

• elements between the cell and its top neighbor (subscript t),

• an element in the top-right corner that is at the interface of 4 different cells (subscript c).

The assembly results in the following sum:

I 2 = 4 ∑ m=1 ( R ⋆T M r m R ) × ( S T D r m S )    right interface elements + 4 ∑ m=1 ( R ⋆T M t m R ) × ( S T D t m S )    top interface elements + 16 ∑ m=1 ( R ⋆T M c m R ) × ( S T D c m S )    corner interface element .
The right interface elements bring the same operators as in the first case. The operators related to the interface elements are also very similar.

The element on the corner is assembled with 16 operators because it is at the interface of 4 different cells. The assembly of the corner element is detailed as follows.

We assume for illustration that the corner element is the 4 nodes quadrilateral element as shown in Fig. 3. The nodes belonging to the corner element of a cell k ∈ {1, . . . , N k } are numbered from 1 to 4. Each one of these nodes i ∈ 1, 2, 3, 4 belong to a different cell c i (k).

Fig. 3 shows that c 1 = k, c 2 is the cell on the right of k, c 3 is on the top-right of k and c 4 is on the top of k. Using these notations, the matrix containing the degrees of freedom of the corner element is:

Q e k = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ R x 1 S c 1 R x 2 S c 2 R x 3 S c 3 R x 4 S c 4 R y 1 S c 1 R y 2 S c 2 R y 3 S c 3 R y 4 S c 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 39
)
Fig. 4. Example of composite structure with periodic geometry.

Table 1

Material properties and geometric parameters of the virtual cell.

Material properties The element stiffness matrix of the corner element is noted K c . The weak formulation restricted to the corner elements is:

E f (GPa) E m (GPa) ν f ν m V f (%)
∑ k∈C ( Q e⋆T k K c Q e k )
, where C denotes the set of cells having an interface element on the top-right corner. This sum is then decomposed as follows:

∑ k∈C ( Q e⋆T k K c Q e k ) = 4 ∑ i=1 4 ∑ j=1 ( [ R x⋆ i R y⋆ i ] [ K c i,i K c i, j+4 K c i+4, j K c i+4, j+4 ] [ R x j R y j ] × ∑ k∈C ( S c i (k) S c j (k) ) ) . (40)
There are then 16 separated operators. It remains only to assemble these element operators into global operators in order to get M c m and D c m (with m = 1, 2, . . . , 16).

Remark. This method requires a lot of separated operators. There are 25 operators for the 2D case treated here and 129 operators for the full 3D case. Though this is a high number, all the matrices are very sparse (some matrices have a few non-zero elements).

Results

In this section, the method described previously will be considered for solving some test cases in order to evaluate its performance in terms of results precision and computational cost.

Numerical results of 2d problem

Consider a composite structure composed of unidirectional fibers. We assume that this structure has a periodical geometry which can be described by an heterogeneous unit cell, as shown in Fig. 4.

Two modeling scales are then considered: the macroscopic scale which define the structure scale and the microscopic scale represented by the unit cell (a single fiber surrounded by resin).

The material properties and the geometric parameters are summarized in Table 1 [START_REF] Berthelot | Matériaux Composites : Comportement Mécanique et Analyse des Structures[END_REF]. The subscript f is used for the fiber properties and the subscript m for the matrix properties. V f denotes the fiber volume fraction.

In this subsection, the numerical simulations performed using the PGD are compared with results of classical FEM implementation using the same set of parameters. The aim is to validate the PGD approach. The two methods are based on the same mesh of the unit cell and the error related to the finite element discretization is exactly the same in both cases. It is important to notice that the comparison between the PGD and the FEM will only show the error due to the separated approximation. In the best case, the PGD will give the same result as the finite element. Therefore, the finite element method serves as a reference solution in order to assess the performance of the PGD approximation.

The problem was divided into 20 substructures, each substructure corresponding to one cell whose length is ∆L = 9 µm. The loads and boundary conditions applied to the structure are shown in Fig. 6. The linear force is set to -1 N mm -1 and the left and right edges are clamped. The problem is solved in 2D using a plane strain assumption with a unity thickness.

Three meshes are considered to study the mesh convergence. A coarse one composed of 149 nodes, a fine one with 2225 node and an intermediate one with 569 nodes are represented in Fig. 5. The three meshes have been generated by successive refinements. 2D triangular elements with 3-nodes are used for all meshes. The axis x is set to be along the fibers and is then perpendicular to the plane.

The finite element mesh was constructed by assembling the meshes of 20 unit cells. For the fine mesh, a total of 43809 nodes involving 87618 degrees of freedom (two DOF per node) are used.

The displacements of the structure are depicted for the intermediate mesh in Fig. 7. The PGD approach and the FEM give quite similar results for the displacements.

The stresses distributions are also plotted in Fig. 8. The results of the PGD are also in excellent agreement with the 2D FEM solution. The PGD approach is able to capture the stress distribution in the microstructure with almost the same precision as the FEM.

Fig. 9 gives the distribution of the stresses through the thickness (coordinate z) at different positions: the center for σ yy and the left edge for σ yz . These positions are chosen to be the most critical for each stress component. The difference between the two methods is also plotted. Referring to these figures, the PGD solution performs very well with respect to the FEM solution. However, we can see that the error increases slightly near the boundaries.

Table 2 summarizes the value obtained for some stress and displacement components at different positions and for the three meshes. The coordinates are given using the center of the geometry for origin and using the dimension of the unit cell denoted ∆L. We can see that the intermediate mesh gives some results very close to the ones of fine mesh. The coarse mesh is clearly less precise, especially for the shear stress σ yz . The normal stress σ zz on the top should logically converge to the external linear force -1 MPa. This is almost the case for the fine mesh. A Timoshenko beam model using a homogenized transverse Young modulus gives an approximate deflection of 0.0657 mm. This value is quite consistent with the deflection obtained at the center, especially for the two finer meshes. The values obtained by the PGD for around 250 iterations are in good agreement with the one obtained by the FEM. The 250 iterations lead to only n ≈ 30 in the separated approximation because the solution is compressed using a SVD or a simple PGD algorithm with the identity operator. However, we can point out some significant differences. The displacement given by the PGD seems to be more precise than the stress. The value of σ zz on the center/top is -0.748 MPa for the fine mesh, far from the -1 expected. This shows a locally non converged solution. The tolerance in the convergence criterion is difficult to choose correctly. This remains a difficulty in the PGD solver and the development of efficient error estimators is probably an important challenge for some future works.

For now, we have used the classical PGD solver which consist in trying to reach the convergence of the separated approximation Eq. (37) with respect to a defined criterion. Another possible approach is to set the number of terms n in Eq. (37) to a chosen value and to search the solution of the weak equation truncated to this n terms. For that purpose, instead of the Algorithm 1 we can compute the n functions F i (x) (or alternately G i (k)) for i = 1 . . . n at the same time. The interested reader can refer to [START_REF] Pruliere | 3D simulation of laminated shell structures using the proper generalized decomposition[END_REF] for more details. only an approximation of the solution. For instance, with 6 terms, we get mainly an order of magnitude of the solution. With 8 terms, the solution is much closer to the FEM solution and the computational cost remains smaller than the one of the classical PGD solver (71 s instead of 108 s for the intermediate mesh with 20 subdomains). With n = 10, the solution is almost converged and we are even closer to the FEM solution than with the classical PGD algorithm. The advantage of this strategy is that the computational time is weakly correlated to the number of subdomains (expected for very high number of subdomains). On a simple laptop (CPU Intel core i7-4510U at 2.00 GHz, 2 cores, RAM 8 Go DDR3 at 1600MHz) the computational cost required to get the solution with 10 terms is about 115 s for 20 subdomains, 127 s for 980 subdomains, 165 s for 2000 subdomains and 1945 s for 18000 subdomains. With 6 terms, the computational cost required to get the solution is about 34 s for 20 subdomains and 637 s for 18000 subdomains. The computational cost is correlated to the number of degrees of freedom in the separated representation, i.e, n × (N micr o + N macr o ), where N micr o is the number of DOF in the microscopic mesh and N macr o is the number of DOF in the macroscopic mesh. The number of terms n required to catch the physics is dependent on the complexity of the problem and not necessarily on the number of subdomains.

The computational cost is given only for sake of illustration because it depends strongly on the implementation. The codes related to the PGD and classical FEM are both home made codes developed in Python using Numpy [START_REF] Der Walt | The numpy array: A structure for efficient numerical computation[END_REF] and Scipy [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF] for sparse matrices. However, it is important to note that the resolution of the Finite Element is only based on Scipy sparse matrices that include some efficient compiled methods to solve sparse linear systems (umfpack [START_REF] Davis | Algorithm 832: UMFPACK V4.3-an unsymmetric-pattern multifrontal method[END_REF], superlu [START_REF] Xiaoye | An overview of SuperLU: Algorithms, implementation, and user interface[END_REF]). In the case of the PGD, there remains many loops in python (Cpython) that are non optimized and a high factor could be saved with a good implementation using an efficient compiler.

Numerical results of 3d problem

For sake of illustration, a multiscale 3D model is used to simulate the mechanical behavior of woven composite materials. This model is difficult to treat using a classical finite element solver with 3D solid elements. It consists in a two-ply composite made from woven carbon fibers (taffeta) and epoxy resin.

Woven reinforcements are generally periodic media: in fact, they consist of a repeating pattern or a unit cell to reconstruct the complete composite fabric, as shown in Fig. 10. Due to the geometrical complexity of woven composites architecture, the 3D geometric model was built using TexGen software [START_REF]TexGen (version 3.7)[END_REF] as shown in Fig. 11(a). TexGen is a very powerful tool to generate geometric models for composite textiles. The TexGen model is transferred to ABAQUS [START_REF]Abaqus cae, dassault systemes[END_REF]. As could be seen in figure Fig. 11(b), the matrix volume is then created by subtracting the yarns volume. The next step is to construct the finite element mesh of the unit cell. In order to do the simulations, the mesh of the unit cell must be periodic. Fig. 11(c) shown the 3D 4 nodes periodic tetrahedral mesh of unit cell model. The unit cell dimension is 2 × 2 × 0.42 (mm). PGD simulations were carried out using linear elastic behaviors, the matrix and the yarns were assumed to be respectively isotropic and orthotropic. Material properties of the matrix and the effective properties of the yarns are given in Table 4. The subscript 1 corresponds to the longitudinal direction of the yarn, and the subscripts 2 and 3 correspond to the transverse directions. The properties of the yarn have been computed from simple rules of mixtures assuming a fiber volume fraction of 90% (E f = 250 GPa).

The unit cell mesh is composed of 28934 nodes. The PGD problem was divided into 25 substructures which give a total of 706390 nodes involving 2119170 degrees of freedom. The displacement along z of the composite plate is zero on both ends (for x = 0 and x = x max ). A constant pressure of 367 MPa is applied on the top face of the unit cell positioned at the center of the plate (surface force). A visualization of the deformed configuration is presented in Fig. 12 and puts emphasis on the displacements. The maximum displacement value in the z direction logically appears at the middle of the plate where the load is applied.

The stress values for the plate are analyzed and presented in Fig. 13. These stress fields are complex because of the heterogeneity of the model and of the orthotropic behavior of the yarn. With the proposed reduced strategy, these stress fields are accessible with a relatively low computational cost (30 min on a laptop).

Conclusion

This paper presents a new separated representation method for solving efficiently multiscale problems. In this work, the Proper Generalized Decomposition is revisited using a separation of scales in periodic geometries. This separated representation involves both the space coordinates of the microscopic scale as well as the space coordinates of the macroscopic scale. For coupling each subdomain, an efficient algorithm is proposed and validated: the proposed technique has a satisfying precision when the convergence tolerance and the number of terms in the separated representation are chosen adequately. The estimation of convergence error remains a hard point in the algorithm that may be addressed in further studies. When the number of periodic cells is high (several thousands at least), especially in 3D problems, the method can be very impressive in terms of computational cost. Another perspective that has not to be investigated in this article is that the method is naturally capable of treating non-linear behavior since a non linear solver is already used to reach the convergence. However, in the case of strong non-linearity, we can expect that the number of terms required to catch the physics will increase and therefore the computational cost will also increase. In this context, the amelioration of the algorithm that build the Proper Generalized Decomposition for problems with a high number of operators is a key point that also needs further developments.

Fig. 1 .

 1 Fig. 1. 1D mesh and decomposition in elementary parts.

D r 2 )Fig. 3 .D r 4 )

 234 Fig. 3. Corner element: nodes and neighboring cells definition.

Fig. 5 .

 5 Fig. 5. The representative unit cell at the microscale.

Fig. 6 .

 6 Fig. 6. Domain decomposed into subdomains and boundary conditions.

  (a) u y . (b) u z .

Fig. 7 .

 7 Fig. 7. Displacement fields (mm) : FEM solution (left)-PGD solution (right).

Fig. 8 .

 8 Fig. 8. Stresses distributions (MPa) : FEM solution (left)-PGD solution (right).

  Unit cell mesh.

Fig. 11 .

 11 Fig. 11. Periodic tetrahedral mesh of unit cell model.

  (a) u x . (b) u y . (c) u z .

Fig. 12 .

 12 Fig. 12. Displacement fields and deformation of a plain weave based composite: Full model is shown (left)-only yarn type materials are shown (right).

  (a) σ x x . (b) σ yy . (c) σ x y .

Fig. 13 .

 13 Fig. 13. Stress distributions of a plain weave based composite: only yarn type materials are shown. Upper surface (left) -lower surface (right).

Table 3

 3 Results of the PGD with predefined number of terms at different positions with the intermediate mesh: Stress components (MPa) and Displacement components (mm). PGD-fixed number of terms in the separated decomposition

		σ yy	σ yy	σ yz	σ zz	u y	u z
	(x, y)	(0, ∆L 2 )	(0, -∆L 2 )	(-4∆L , 0)	(0, ∆L)	(-4∆L , ∆L)	(0, 0)
	n = 6	-4.83	4.741	-4.092	-0.537	0.0092	-0.0567
	n = 8	-5.084	4.869	-3.81	-0.991	0.00963	-0.0611
	n = 10	-5.406	5.188	-3.736	-1.0176	0.0105	-0.0619

Fig. 10. Woven carbon fiber.

Table 4

 4 Material property.

	Matrix		Effective properties of yarn
	E (MPa)	ν	E 1 (GPa)	E 2 (GPa)	ν 12 /ν 13	ν 23	G 12 /G 13 (GPa)	G 23 (GPa)
	3500	0.3	234.345	45.743	0.327	0.375	20.056	16.628

Results of the FEM and PGD at different positions with 3 meshes: Stress components (MPa) and Displacement components (mm). Table 3 gives some results given by the PGD using this strategy for n = 6, n = 8 and n = 10 with the intermediate mesh. Using a small number of terms permits to reduce significantly the computational cost of the simulation but gives