Quentin Cossart
email: quentin.cossart@ensam.eu

Frédéric Colas

X Kestelyn

A priori error estimation of the structure-preserving modal model reduction by state residualization of a grid forming converter for use in 100% power electronics transmission systems

Keywords: Model Order Reduction, Simulation, Power System, Power Electronics

A priori error estimation of the structure-preserving modal model reduction by state residualization of a grid forming converter for use in 100% power electronics transmission systems

Introduction

The development of renewable energies and High Voltage Direct Current (HVDC) links rapidly increase the Power Electronics (PE) penetration in the transmission grids. This makes the study of all-converter-interfaced power systems (see Figure 1) necessary [START_REF] Ramasubramanian | Converter Model for Representing Converter Interfaced Generation in Large Scale Grid Simulations[END_REF]. Because of the systems size and the converters complexity, numerical simulations are needed. Converters are physically different from synchronous generators. Therefore, commonly used tools, like the phasor approximation, might not be relevant. Usually, simulations of PE systems are based on detailed EMT [START_REF] Van Der Meer | Advanced Hybrid Transient Stability and EMT Simulation for VSC-HVDC Systems[END_REF] models but, in the case of large systems, the computation time would be large and the analysis complicated.

Consequently, Model Order Reduction (MOR) methods can be useful [START_REF] Gu | Reduced-Order Models for Representing Converters in Power System Studies[END_REF], to simplify the models while keeping a good accuracy. Several methods exist [START_REF] Antoulas | Approximation of Large-Scale Dynamical Systems[END_REF]. Yet most of them use basis changes and truncations, altering the system's structure and changing its state variables. This is why a MOR by state residualization [START_REF] Antoulas | Approximation of Large-Scale Dynamical Systems[END_REF] is applied in this paper to the converters, in order to preserve their physical structure and state variables. It is indeed necessary in power systems stability studies to keep the structure and variables, in order to know which quantities are critical. This paper presents a method to a priori estimate the error induced by this MOR in order to help the user choose the most suitable models, those that minimize this error, depending on the case study under consideration. The first part of this paper presents the structure of the studied converter and its MOR. The second part deals with the method to a priori estimate the induced error and its application to find the most suitable reduced models in two test cases: one converter connected to the infinite grid and two converters connected to each other, a load and the infinite grid.

Structure-preserving model order reduction of a PE converter

In this part an example of PE converter is presented and reduced using the developed structure-preserving MOR method by state residualization. The first subsection presents the converter's structure and model while the second one introduces the MOR methodology. The differential and algebraic equations describing this system and its control are given in Equations ( 1)-(31), starting from the external control (1)-( 13), then the voltage control ( 14)-( 17), the current control (18)-( 21), the converter ( 22)-( 23), the filter (24)-( 27), the transformer (28)-( 29) and the grid (30)-(31). The state variables are in red, the algebraic ones in blue, the parameters in green and the inputs in purple. The equations and parameters of the model are taken from [START_REF] Qoria | Tuning of cascaded controllers for robust grid-forming Voltage Source Converter[END_REF] and more details can be found in this paper. The variables are expressed in the reference frame of the converter's frequency.

Structure of the studied converter

(1)

" # (6) %& ( 7)

%& ( 8
) )& | | )& )& %& ( 9
) )& %& ( 10
) -. )& ( )& ) (11) / 0 1 -. -. ( 12
) %1 1 -. 2 0 3 (13) 4 " )& # (14) 4 " )& # (15) -)& " )& # … … -. 6 4 (16) -)& " )& # … … -. 6 4 (17) 7 " -)& -# (18) 7 " -)& -# (19) -)& " -)& -# … … -. 8 - 7 (20) -)& " -)& -# … … -. 8 - 7 (21) - -)& (22) - -)& (23) 8 0 - - 9 - -. 8 - (24) 8 0 - - 9 - -. 8 - (25) 6 0 - -. 6 (26) 6 0 -

| | )&

1 pu

)& 0.5 pu 2 2 0 pu )& 1 pu )& 0 pu 2 1 pu 2 2

pu

Table 2: Inputs.

To sum up, the converter is modelled by a 15 th order nonlinear differential-algebraic model. The next subsection investigates the MOR of this model, in order to facilitate the analysis and simulation of large transmission systems with several converters.

Model order reduction of the converter

Several MOR methods exist, like the balanced truncation [START_REF] Ramirez | Application of balanced realizations for model-order reduction of dynamic power system equivalents[END_REF], the proper orthogonal decomposition [START_REF] Antoulas | Approximation of Large-Scale Dynamical Systems[END_REF], the modal truncation [START_REF] Chow | Power system coherency and model reduction[END_REF] or Krylov methods [START_REF] Chaniotis | Model Reduction in Power Systems Using Krylov Subspace Methods[END_REF]. However, these methods use basis changes and truncations (projections on subspaces), which change the physical structure of the system. The state variables are indeed altered by the MOR, which is not wanted here for analyses issues. It is indeed mandatory, for stability studies, to keep the variables in order to identify the critical ones. This is why a method using state residualization [START_REF] Antoulas | Approximation of Large-Scale Dynamical Systems[END_REF] is investigated here. The state residualization consists in freezing the dynamics of some chosen states of the system. These states are thus changed into algebraic variables, which reduces the order of the system. For example, in the studied converter, the residualization of the state ? consists in changing Equation ( 32) into (33).

(32) 0 (33)

One advantage of doing state residualization is that it can easily be implemented in a program/language that solves differential-algebraic equations, like Modelica [START_REF] Fritzson | Principles of object oriented modeling and simulation with Modelica 3.3 : a cyber-physical approach[END_REF] (which is used here). This language consists in writing the differential and algebraic equations almost like in the natural language, like shows Figure 3, which is very convenient. With the used reduction method, the variables are kept the same, and so is the structure of the system. To illustrate this; figures 4 and 5 show the evolution of the current controller when the MOR is applied (i.e. when the state variables of the PI controllers of the current control are residualized). During the process, Equations ( 18) and ( 19) are changed into Equations ( 34) and ( 35).

-)& - (34) -)& - ( 35
)
It can be seen on this figures that the reduced system can still be represented in a physical way, which is not the case with methods such as the POD or the balanced truncation (they indeed project the state variables on a subspace, which changes them). The philosophy of the method is then to find the states that can be residualized without too much decreasing the accuracy of the model. To do so, a modal approach using the participation factors is used. The participation factors give the dependencies between the states and the modes/poles of the linearized system. This way it is known which poles are discarded when some states are residualized. They are numbers between 0 and 1, and there is one for each state/pole couple. For each state, the sum of its participation factors in each pole is equal to 1. And for each pole, the sum of its participation factors in each state is equal to 1. To calculate the participation factors, the model first needs to be linearized around its operating point, which gives Equation (36) (the algebraic equations are injected in the differential ones), with

@ " A B B A B C B CA ? D EF G D EF H I I A J C KL M M A N ? # . ∆P A∆@ B∆S (36)
Then the participation factor of T@ U in V W , an eigenvalue of X, is given by Equation (37).

Y U,W S U,W [ U,W (37)
In this equation [ U,W is the \ ] entry of the B ] normalized right eigenvector of X and S U,W is the \ ] entry of the B ] normalized left eigenvector of X. Figure 6 gives an example for one eigenvalue. It shows that the considered eigenvalue is mainly linked to the dynamics of the variables M , M A , I , I A and J C KL . This process is done for each eigenvalue. The participation factors are ranked from the largest to the smallest and the first ones are chosen to have a total participation of at least 90%. It leads to Table 3: Dependencies between the eigenvalues and the state variables of the system.

Thanks to Table 3, it is possible to know which dynamics to freeze in the nonlinear model in order to discard some poles in the linearized one while maintaining the other untouched. This method ensures the stability as the remaining are not moved (see Figure 7, where the fastest poles (i.e. the one with the largest negative real part) are removed for each reduced model. We can see that the remaining poles are really close to the initial ones.).

Figure 7: Eigenvalues of the full and reduced models.

The question now is how to choose which poles to delete. Usually, the fastest ones (with the highest real part in absolute value) are removed [START_REF] Cossart | Model reduction of converters for the analysis of 100% power electronics transmission systems[END_REF] as they are linked to fast transients that are quickly cleared and they don't have any impact on the stability, which is not the case for the poles close to the imaginary axis. But there is no mathematical proof that this is the best solution in terms of accuracy. This is why the next section presents a method that estimates the error induced by the residualization and thus helps the user choose the most suitable reduced model for its studies (i.e. the best poles to keep/discard).

A priori estimation of the error induced by the reduction

This section presents the method to a priori estimate the error made by the MOR and the deduced optimization problem that helps minimizing it in order to have the most accurate reduced model. Two test cases with respectively one and two converters are then studied.

Mathematical method

To simplify, the linearized model is considered in Equation (38) (to simplify the writing, the T are removed).

P A@ BS (38)
The residualization consists in multiplying the first term of this equation by a diagonal matrix i, that has the same size as X, made of 1 and 0. If i(B, B) 0, it means that the state @ U is residualized, if it is equal to 1, the state dynamic is kept. The trace of i gives the order of the reduced model. Table 3 gives some constraints on i in order to keep the remaining poles unchanged. For example I , I A , J C KL , M , M A and N ? must be kept (or residualized) together in order to properly remove the poles. The reduced model is then given in Equation (39).

i P j

A@ k BS (39)

A Laplace transform give Equations ( 40) and (41). l (sI A) n/ Bo (40)

l k (sE A) n/ Bo (41)
Then a subtraction gives Equation (42). l l k q(sI A) n/ (sE A) n/ rBo so (42) s is the transfer function of the error induced by the MOR. It is a matrix of dimension t u N with t u t the size of X and t u N the size of v.

It is possible to choose the best i to minimize the maximum of the transfer function s(B, \) which represents the error for the B ] state variable when considering the \ ] input. This optimization is summed up in Equations ( 43)-(44). In Equation 44, the constraints imposed by Table 3 are put in equation. This way it is possible to have a specific reduced model for each event to be simulated (each \), depending on the state variable that is looked at (each B). Examples to illustrate this will be shown in the two next parts.

Application to find the best reduced model on a single converter test case

The first example is the converter of Figure 2 connected to an infinite grid. The considered state variable is the current in the converter B C G (the same results apply for B C H ) and the considered input is the grid voltage • G Ž . Table 4 gives, for each specified size, the dynamics to freeze/the poles to discard in order to have the most accurate reduced model. These results are obtained by solving the optimization problem of Equations ( 43 V /d,// , V /b,/c , V /_,/`, V g , V f , V e Table 4: Dynamics to freeze/ poles to discard for each reduced model size.

However, in this case, there is most of the time only one possible reduced model that respects Table 3 for a given size. This is why a more complex case with two converters is investigated in the next section. In this case, for each size, there are several possibilities for the reduced model, and the method helps choosing the best one.

Application to find the best reduced model on a twoconverter test case

The system shown in figure 8 is considered. It is made of two identical converters, connected to each other and to a load and the infinite grid. The additional parameters are given in Table 5.

9 • .& 0.006 pu 8 • .& 0.2 pu 9 • • 2 pu Table 5: Additional parameters.
This case is more interesting than the previous one, as for each total size, different models are possible, as shows Table 6. This is the reason why the developed method is interesting here. As an example, the considered state variable is the current in the first converter B C / (same for B CA/ ) and the considered input is the load ' ' " . For the chosen total size of 21, figure 9 shows the estimated error for each possible combination in Table 6. It shows that to keep a good accuracy, converter 1 should not be reduced while converter 2 can. This helps choosing which converter to reduce and how much depending on the event and variable under consideration. On Figure 10, the evolution of the current in the first converter B C/ when ' ' " goes to 0 pu at 1s (short-circuit) and the fault is cleared after 100ms is shown for 3 cases: two full converters, a 15 th order converter 1 and a 6 th order converter 2, a 6 th order converter 1 and a 15 th order converter 2. This figure shows that the case with a reduced first converter gives bad results during the short-circuit. It misses an overcurrent. On the other hand, it adds one that doesn't exist when the fault is cleared. The case with a reduced second converter gives results that are very close to the case with a full second converter. This is what was expected because the fault appears far from this converter and we are looking at the current in the first one. The developed method proved it. Moreover it can be used in cases where this conclusion is not obvious, where there are many converters in a large network.

Conclusion and perspectives

In this paper the MOR of PE systems is investigated. To keep the physical structure, a method by state residualization is used; and to ensure the stability, a method keeping the remaining poles of the system unchanged has been chosen. The main idea of this article is the development of a method that a priori estimates the error made by the reduced model and helps choosing the optimal one for each converter of the system depending on the event to be simulated and the variable to be looked at. A case with two converters has been tested. The method gives the best model for each converter.

It now needs to be tested on systems with many converters. However, as it needs to enumerate all the possible combinations before finding the best one, the process is long. This is why a method to find the poles to eliminate more quickly is investigated at the moment.

Figure 1 :

 1 Figure 1: Structure of a 100% PE power system.

Figure 2

 2 Figure 2 presents the structure of the studied converter. The control is made of a classical cascaded loops structure with a current loop, a voltage loop and an external loop (here an active power and a reactive power droop controls). The power part is made of the DC/AC converter (the DC part is modelled with an ideal DC voltage source and the model for the converter is an averaged one), an RLC filter and a transformer modelled with an RL line. The aim of this grid forming converter is to generate the voltage at the PCC (Point of Common Coupling) .

Figure 2 :

 2 Figure 2: Structure of the studied converter.

Figure 3 :

 3 Figure 3: Modelling of the RLC filter in the Modelica language.

Figure 4 :

 4 Figure 4: Structure of the current control.

Figure 5 :

 5 Figure 5: Structure of the current control after the MOR.

Figure 6 :

 6 Figure 6: Participation of the states in the eigenvalue V / 1.4.

  minimize z {max ~(s(B, \))• (43) €S•\'ƒ3 3": i B † ‡(ˆU), ˆU ∊ {0; 1• ∀B; i(1,1) i(2,2) ⋯ i(5,5) i(6,6); i(8,8) i(9,9); i(10,10) i(11,11) ⋯ i(14,14) i(15,15) (44)

Figure 8 :

 8 Figure 8: Considered 2-converter system.

Figure 9 :

 9 Figure 9: Estimated error for each combination

Figure 10 :

 10 Figure 10: Evolution of the current in the three cases.

Table 1 :

 1 Parameters.

	-. 6	(27)

Table 3

 3

	Eigenvalues 790 a 3821B 759.6 a 3331B 21 a 133.6B V e 31.4 V f 17.2 V g 16.2 15.5 a 28.3B V _ 1.7 V b,c V /_,/` V /b,/c V /d,// V `,h V / 1.4 1.5 a 0.1B	State variables , A , B , B A , B C , B CA ? D EF G , D EF H I , I

. Note that each eigenvalue is linked to several states and each state participates in several eigenvalues (this is a coupled system). This leads to groups of eigenvalues/state variables. A , J C KL , M , M A , N ?

Table 6 :

 6 Possible reduced models for each converter depending of the size of the whole system.

	Total size	Converter 1/Converter 2 size
	12	6/6
	13	7/6 or 6/7
	14	8/6 or 7/7 or 6/8
	15	9/6 or 8/7 or 7/8 or 6/9
	16	9/7 or 8/8 or 7/9
	17	9/8 or 8/9
	18	12/6 or 9/9 or 6/12
	19	13/6 or 12/7 or 7/12 or 6/13
	20	14/6 or 13/7 or 12/8 or 8/12 or 7/13 or 6/14
	21	15/6 or 14/7 or 13/8 or 12/9 or 9/12 or 8/13 or 7/14 or 6/15
	22	15/7 or 14/8 or 13/9 or 9/13 or 8/14 or 7/15
	23	15/8 or 14/9 or 9/14 or 8/15
	24	15/9 or 12/12 or 9/15
	25	13/12 or 12/13
	26	14/12 or 13/13 or 12/14
	27	15/12 or 14/13 or 13/14 or 12/15
	28	15/13 or 14/14 or 13/15
	29	15/14 or 14/15

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 691800. This paper reflects only the authors' views and the European Commission is not responsible for any use that may be made of the information it contains.


Executive summary

This deliverable specifies the development of the proof of concept (POC) which will allow the linking of publications deposited in HAL (including its HAL-SHS portal), the French open archive (https://hal.archives-ouvertes.fr/) developed by the CCSD (https://www.ccsd.cnrs.fr/), and data or datasets deposited in Nakala (https://nakala.fr/) developed by Huma-num (https://www.huma-num.fr/), the French research data repository for SSH data.

It describes the workflow to be set up between the two repositories, the creation of the bidirectional relationship between publications and data, prospection for the implementation of tools to visualize these links in the two repositories, and finally how to expose these relationships and publish them, so that they are available, findable and usable by research communities (and by all citizens).

Introduction

The new publishing models combine articles and underlying data, deposited in publisherspecific or publisher-recommended repositories, which contributes to the transparency and reproducibility of research, allows for better evaluation of results and promotes data reuse.

But beyond the open access publishing model, open science principles require that publications be deposited in open archives and that data be managed according to FAIR principles, and deposited in data repositories. Today, there are different types of repositories, publication repositories and data repositories, which can be generalist, disciplinary or institutional, independently of each other, and without the deposited objects being linked to each other. The search for, and the exploitation of, these research objects are therefore decorrelated. Linking the different objects together, publications, data, as well as the software codes produced during a research project makes it possible to put the research in context, to increase its visibility and its reusability.

The aim of this POC is thus to increase interoperability between repositories. The objective is to create the relationships between the publications deposited in the HAL open archive, and the data deposited in the Nakala data repository, and to expose the relationships thus created so that they can be easily found and exploited. A second objective is to test tools for displaying and disseminating these links, so that they can be exploited in the European Web of data (e.g. Scholexplorer https://scholexplorer.openaire.eu/, or B2Note https://esdf.github.io/b2note-docs).

Description of the two repositories

HAL

The HAL open archive (https://hal.archives-ouvertes.fr) is the common platform, shared by the French academic community, for the open access dissemination of scientific production.

HAL provides access to the full-text of journal articles, or papers, theses, reports and so on. The document may or may not have been published, so both journal articles and preprints can be found in the archive. Documents from all academic fields can be submitted to HAL. The humanities and social sciences represent 20 percent of the documents (files) submitted to HAL.

HAL is a central repository for infrastructure (technical centralization), but is multifaceted (distributed portals), based on the metadata describing each deposit. One of these portals is HAL-SHS, which shows archives and disseminates scientific literature in all human and social science disciplines. HAL provides:

 an OAI infrastructure ensuring interoperability (OAI-PMH), API and Sword protocol  dissemination of metadata through a Triple Store  assignment of a PID (Persistent IDentifier), making data and metadata citable  stability of identifiers (URL of deposits in particular)  preservation of documents thanks to a partnership with CINES for archiving  scientific quality of documents deposited as well as the details describing them: all deposited documents are therefore checked before being put online  time stamping of deposits guaranteeing the intellectual property rights of the text deposited  authority files available via the AURéHAL platform, which is interoperable with other listings.

HAL has been developing relationships and links with other international repositories, such as ArXiv, PubMedCentral and Software Heritage. When submitting a document or software code in HAL, the researcher can ask to push it towards one of these repositories, and HAL retrieves the repository identifier and displays it in the identifiers metadata of HAL.

In HAL, there are already different ways to create a relationship between a publication and another resource:

 A specific interface for linking documents together. o This functionality is only available for documents within HAL (with possible typing of the relationship). It requires being connected to HAL and having rights on the repository (property of the repository).  Cf. appendix for the types of relations proposed in HAL.

 Specific metadata at the time of the deposit: several fields make it possible to link the deposit to another resource:  Other identifier: corresponds to the identifier of the document in another repository (arxiv, pubmed, ...)  See also: URI to a document associated with the publication  Associated data: DOI to data (or a dataset) deposited in another repository  Linked publications: identifier of the associated publication (field accessible for software and data deposit (MediHal).  Retrieval of bibliographical references associated with the publication using the GROBID tool (https://grobid.readthedocs.io/en/latest/Introduction/).

These various links are possible today, but are not exploited or are underexploited.

Relationships can be added in several places: in the repository (link to a repository outside HAL) and as a feature outside the repository (specific interface to link repositories in HAL together).

Characteristics of the Relationships: the entry is manual in both cases (DOI or URL linked in the publication metadata / establishment of a relationship to link two objects in HAL).

NAKALA

NAKALA is an interoperable and secure service for depositing all types of data (text files, audio, video, images etc.) in order to share them. Based on OpenSource technologies such as Symfony, Elasticsearch, React, etc. this repository mainly provides different types of services:

 assignment of a PID (Persistent IDentifier) making data and metadata citable;  permanent data access through Web frontend, an integrated search engine and APIs;  dissemination of metadata through a Triple Store and OAI-PMH;  organization of data into collections;  presentation of collections in specific websites through Nakala_Press.

The default format used by NAKALA is the DublinCore format qualified for metadata but more broadly any type of field of a specific format can be adopted.

Each data is visible on a dedicated page (landing page) which presents the visualization of the data (different viewers for common format files, IIIF, etc.) as well as the associated metadata. This page can be accessed via the DOI assigned to the data.

User stories

In order to best reflect the needs of researchers (user side) as well as the possibilities of implementing the link between data and publications (infrastructure side), we have developed the following use cases:

 An author wishes to link a publication already deposited in HAL with the data present in NAKALA.  A person wishes to access/view the data associated with a publication.  An author wishes to link data already stored in NAKALA with publications already published in HAL.  A person wishes to access/view the publications associated to a data.  As HAL manager, I would like to retrieve the relation(s) created in NAKALA (and vice versa)  As a repository manager, I wish to make statistics on the citation of the data deposited in my repository.

4 Describing relationships within both repositories The choice of the typology of the relationship will have to be made upstream of the developments, giving priority to interoperability. We will choose a vocabulary that already exists and evolves. The chosen vocabulary will be implemented in the repositories, used in the export formats, and will have to be made compatible with the vocabularies used in the data exposure formats (DC Terms, ScholeXplorer schema for example).

Relationship between objects

The licence granted to the relationship will be a CCO licence https://creativecommons.org/publicdomain/zero/1.0/deed.en, which allows full reuse.

Relationship building workflow

Fig.1 -Schema of the POC's workflow between HAL and NAKALA

The user can create a relationship between data deposited in HAL and data already stored in Nakala:

 When creating a new repository in HAL  On a deposit already made in HAL To be able to create a relationship, the user must have an account in HAL and be the owner of the data deposited in HAL for which he is creating a relationship (or have rights to edit the data). It is not necessary for the user to have rights to the data in NAKALA.

Searching for data in NAKALA:

 The user can indicate the data identifier (Handle or DOI).  The user can search the data in Nakala thanks to a search API made available by NAKALA. The search will be carried out on the different metadata describing the data in NAKALA. However, the types of data used by SSH are by nature very varied, which does not always make them easy to cite: the use of so-called dynamic data from social networks is a good example of this.

In brief, citing data especially in the context of SSH is important for different reasons:

 Providing a way to reproduce research which will in turn enhance the quality and effectiveness of research  Reusing data for different research purposes in other contexts  Giving credit to the creator and the funder of the data  Proving the usefulness of infrastructures  Enhancing links between data and publications An inventory of current citation practices has been carried out within the framework of SSHOC (https://sshopencloud.eu/), a European SSH cluster project bringing together the main SSH actors: Inventory of SSH citation practices, and choice for SSHOC citation formats and implementation planning (https://doi.org/10.5281/zenodo.3595964).

The citation of data is generally based on various recommendations made in the framework of working groups, for example: More recently, the creation of the notion of "Fair Digital Objects" (https://github.com/GEDE-RDA-Europe/GEDE/tree/master/FAIR%20Digi tal%20Objects/FDO-declaration) a synthesis of various previous studies (e.g. "A Framework for Distributed Digital Object Services"http://www.cnri.reston.va.us/home/cstr/arch/k-w.html) will facilitate the citation of data even if the main purpose is to provide a framework for their interoperability. The next step is to make these citations "actionable", i.e. machine-readable and machineprocessable. The W3C has published recommendations for making data catalogues on the World Wide Web interoperable (https://www.w3.org/TR/vocab-dcat/).

 APA & MLA-
Other services are being developed, for example by Crossref and Datacite to link datasets and other objects identified by DOIs. In this framework these institutions have developed a very comprehensive DOI citation formatting service (https://citation.crosscite.org/).

Large companies are also increasingly present on these subjects. One cannot fail to mention the "Google Data Search" service which has recently gone into official production (https://datasetsearch.research.google.com/), and traditional publishing houses are also positioning themselves such as Elsevier with its service for Mendeley data (Mendeley Datahttps://data.mendeley.com).

It is therefore crucial today to have tools for linking data and publications for the trusted repositories offered to communities.

Repositories web interface

In HAL, the relationships created will be displayed on the landing page of a publication.

The citation of the data will be displayed, making it possible to bounce on the landing page of the data in NAKALA in order to have a more precise description of the data.

APIs

Parameter

Identifier (DOI, Handle, HALID)

Result

Citation of the resource

Implementation of an API enabling the retrieval of relations made in the other platform from an identifier. In this way a relation created in HAL can be retrieved in NAKALA and vice versa.

Exports in different formats

 Present in the different exports available in which it can be integrated (TEI, RDF, bibteX)  Exhibited in HAL's Triple Store (RDF)  Indexed in the HAL search engine (SOLR)  Exhibited in HAL's OAI-PMH repository: DC terms 1 To support the seamless operation of services in the future European Open Science Cloud, the EOSC-hub has proposed a number of technical specifications and interoperability guidelines covering both common and federation services

External tools to disseminate relationships

For the moment, we have considered the use of two different tools to disseminate the relationships created between HAL publications and NAKALA data.

B2Note

B2Note (https://www.eudat.eu/catalogue/B2NOTE) is a tool for making annotations on resources. The tool also allows the user to create relationships between resources. The relations created between HAL and NAKALA can be added in B2Note in order to be more widely disseminated.

 Documentation: https://e-sdf.github.io/b2note-docs  Service : https://b2note.bsc.es/  EOSC Technical specifications and interoperability guidelines 1 : annotation service EUDAT-B2NOTE ; https://wiki.eoschub.eu/display/EOSCDOC/Metadata+Management+and+Data+Discovery?previe w=/68223176/68223179/EOSC_Technical_specification_AnnotationService_Prop osal_v1.pdf

SCHOLIX, A Framework for Scholarly Link eXchange

The objective of the Scholix initiative is to establish a high-level interoperability framework for the exchange of information on the links between scientific literature and data.

It might be interesting to trace these relationships also in the Scholexplorer platform of OpenAire: OpenAire, as one of the Scholix hubs, assists in the global aggregation of dataliterature link information.

Conclusion -Prospects

The first step is to build the relationship between the publications deposited in HAL and the data deposited in Nakala, using the APIs available in each of the repositories. The relations thus created will be displayed, exported and harvestable.

A second step will consist, if possible, in allowing the simultaneous deposit of publications and data in the same repository, HAL, and then transferring the data into Nakala, creating the relationship between the two repositories automatically, on the model of what has already been developed in HAL for software codes, in partnership with the Software Heritage repository.

It is also conceivable, by means of a survey of users of the two data repositories, to deepen the bi-directional character of the "simple" link initially created, as well as its visualization.

Appendices

Relationship typing

Depending on the repository they belong to, vocabularies are already available to define a relationship between documents.

Repository Relations

HAL  illustrate  is illustrated  requires  is required by  has part  is part of  references  is referenced by  has format  is format of  conforms to  has



  http://www.scholix.org/  https://scholexplorer.openaire.eu/#/ MS26 Communications PoC between NAKALA and HAL www.eosc-pillar.eu Page 18 of 22



  HasVersion (indicates A has a version B)  IsVersionOf (indicates A is a version of B)  IsNewVersionOf (indicates A is a new edition of B, where the new edition has been modified or updated)  IsPreviousVersionOf (indicates A is a previous edition of B)  IsPartOf (indicates A is a portion of B; may be used for elements of a series)  HasPart (indicates A includes the part B)  IsReferencedBy (indicates A is used as a source of -information by B)  References (indicates B is used as a source of information for A)  IsDocumentedBy (indicates B is documentation about/explaining A)  Documents (indicates A is documentation about/explaining B)  IsCompiledBy (indicates B is used to compile or create A)  Compiles (indicates B is the result of a compile or creation event using A)  IsVariantFormOf (indicates A is a variant or different form of B, e.g. calculated or calibrated form or different packaging)  IsOriginalFormOf (indicates A is the original form of B)  IsIdenticalTo (indicates that A is identical to B, for use when there is a need to register two separate instances of the same resource)  IsReviewedBy (indicates that A is reviewed by B)  Reviews (indicates that A is a review of B)  IsDerivedFrom (indicates B is a source upon which A is based)  IsSourceOf (indicates A is a source upon which B is based)  IsRequiredBy (indicates A is required by B)  Requires (indicates A requires B)

  Contents 1 Introduction ....................................................................................................................7 2 Description of the two repositories ................................................................................. 8 2.1 HAL............................................................................................................................8 2.2 NAKALA.....................................................................................................................9 3 User stories .................................................................................................................. 10 4 Describing relationships within both repositories.........................................................11

4.1 Relationship between objects..................................................................................11 4.2 Relationship building workflow ................................................................................ 12 5 Exposure of relationships.............................................................................................14 5.1 About data citation practices and links to publications in SSH and beyond.............14 5.2 Repositories web interface ......................................................................................15 5.3 APIs .........................................................................................................................15 5.4 Exports in different formats......................................................................................16 6 External tools to disseminate relationships..................................................................17 6.1 B2Note.....................................................................................................................17 6.2 SCHOLIX, A Framework for Scholarly Link eXchange............................................17 7 Conclusion -Prospects ................................................................................................ 18 8 General Bibliography....................................................................................................19 9 Appendices...................................................................................................................20 9.1 Relationship typing..................................................................................................20

Table 1 -Definition of the relationship between objects

 1

	Local PID		PID
	Local repository HAL or NAKALA for POC
	Type relationship	of	Type of relationship between objects. See the Appendix for what is already in place in some repositories
	Remote identifier		Identifier
	Remote repository		HAL or NAKALA for POC
	Creation date		Date of creation of the relationship

Table 2 -Search API Description

 2

	5 Exposure of relationships
	5.1 About data citation practices and links to publications in SSH and beyond Although data citation has been a long-standing concern (see "Out of Cite, Out of Mind" -http://doi.org/10.2481/dsj.OSOM13-043), data citation practices in SHS are quite disparate and more recent even though the Social Sciences have developed considerable experience in this area.

Allows you to search in a repository. The search is carried out on the different metadata of the objects

Parameter

Text searched for (identifier or metadata title, author, ...)

Result

List of objects containing:

 the identifier of the resource  the citation of the resource This mechanism will be set up in both repositories in order to create relationships in HAL and NAKALA.

NB: In a first step, the applicant and users with the administrator role on the data are able to create the relationship in NAKALA.

This is due to the fact that research data were previously considered secondary. However, these practices are evolving rapidly, encouraged in particular by the need to establish data management plans which are required by agencies to obtain funding. In addition, data paper type objects are appearing in SSH (e.g. Cybergeo: European Journal of Geography has been publishing data papers since 2017: https://journals.openedition.org/cybergeo/28545), evidencing the interest in data today.

  type citations that were originally intended for publications but have been adapted for datasets (See "MLA handbook for writers of research papers, New York: Modern Language Association of America, 2009, 7th ed.");  The citation extension for electronic documents "ISO 690";  RDAs Data Citation WG (https://www.rd-alliance.org/groups/data-citation-wg.html);  Data Citation Synthesis Group: Joint Declaration of Data Citation Principles. Martone

	 The DataCite Metadata Working Group (See "DataCite Metadata Schema" https://schema.datacite.org/ and http://schema.datacite.org/meta/kernel-3/doc/DataCite-MetadataKernel_v3.1.pdf)

M. (ed.) San Diego CA: FORCE11; 2014 https://doi.org/10.25490/a97f-egyk``

Table 3 -Search API description Search API Description Allows

 3 you to retrieve the citation of a data or publication.

Table 4 -Relationship recovery API Description Allows

 4 you to retrieve the relationships created in a repository based on an identifier

	Parameter	Identifier (DOI, Handle, HALID)
	Result	List of relations

  IsCitedBy (indicates that B includes A in a citation)  Cites (indicates that A includes B in a citation)  IsSupplementTo (indicates that A is a supplement to B)  IsSupplementedBy (indicates that B is a supplement to A)  IsContinuedBy (indicates A is continued by the work B)  Continues (indicates A is a continuation of the work B)

		 documents this upload  is documented by this upload  is compiled/created by this upload  compiled/created this upload  is the source this upload is derived from  has this upload as its source  is identical to this upload  is an alternate identifier of this upload
	Datacite	https://support.datacite.org/docs/relationtype_for_citation
		version  is version of
	Scholix	https://zenodo.org/record/1120265#.Xy1HlhP7TUI  IsSupplementTo (indicates that A is a supplement to B when both are published together)  IsSupplementedBy (Indicates that B is a supplement to A when both are published together)  References (Indicates B is used as a source of information for A)  IsReferencedBy (Indicates A is used as a source of information by B)  IsRelatedTo (Indicates a generic relation between A and B)
	Zenodo OpenAire	https://developers.zenodo.org/#representation  cites this upload  is cited by this upload  is supplemented by this upload  is a supplement to this upload  is continued by this upload  continues this upload  is new version of this upload  is previous version of this upload  references this upload  is referenced by this upload h t t p s : / / o p e n a i r e -g u i d e l i n e s -f o r -l i t e r a t u r e -r e p o s i t o r y -managers.readthedocs.io/en/v4.0.0/vocab_relationtype.html#vocab-relationtype-relationtype

 has this upload as part  is part of this upload  reviews this upload  is reviewed by this upload  IsCitedBy

 Cites  IsSupplementTo  IsSupplementedBy  IsContinuedBy  Continues  HasMetadata  IsMetadataFor  IsNewVersionOf  IsPreviousVersionOf  IsPartOf  HasPart  IsReferencedBy  References  IsDocumentedBy  Documents  IsCompiledBy  Compiles  IsVariantFormOf  IsOriginalFormOf  IsIdenticalTo  IsReviewedBy  Reviews  IsDerivedFrom  IsSourceOf   IsDescribedBy (

indicates A is described by B)  Describes (indicates A describes B)  HasMetadata (indicates resource A has additional metadata B)  IsMetadataFor (indicates additional metadata A for a resource B) MS26 Communications PoC between NAKALA and HAL www.eosc-pillar.eu Page 22 of 22


domain: spi phys math info

text:
Joan Mas Colomer

Nathalie Bartoli
email: nathalie.bartoli@onera.fr

Joaquim Thierry Lefebvre
email: thierry.lefebvre@onera.fr

R R A Martins

Joseph Morlier
email: joseph.morlier@isae-supaero.fr

Joan Mas Colomer

Joaquim R R A Martins

An MDO-based methodology for static aeroelastic scaling of wings under non-similar flow

Keywords: Aeroelastic Scaling, Multidisciplinary Design Optimization, Aircraft Design Nomenclature

() c

Relative to the cruise condition () m Relative to the scaled model () r Relative to the reference aircraft () T Transpose of a matrix or vector Searching for more efficient airplane designs drives aerospace engineers to investigate innovative solutions, such as the blended wing body, the box wing, or the strut-braced wing aircraft configurations. Whereas the conventional tube-and-wing configuration is well known and studied, the in-flight behavior of these new aircraft types often remains unknown. The development, construction, and experimentation of unmanned aeroelastically-scaled prototypes is a way to study the in-flight behavior of these innovative aircraft configurations. This mitigates the economical and operational risk of building and testing a full-scale prototype. In addition, aeroelastically-scaled prototypes can be used to show the impact of current aircraft layout modifications (e.g., an engine upgrade) on the in-flight response.

Traditional aeroelastic scaling theory of aircraft considers that the nondimensional fluid characteristics of the full-scale aircraft are preserved, or that the effect of their discrepancy (particularly the Mach number) can be neglected. However, attaining the same Mach number is sometimes impossible given the operational boundaries of the scaled unmanned prototypes. [START_REF] Bisplinghoff | Aeroelastic Scaling of a Joined Wing for Nonlinear Geometric Stiffness[END_REF] defined the basics of classical aeroelastic scaling theory in 1955. [START_REF] Scruton | Similarity Requirements for Flutter Model Testing[END_REF] expanded the method to include the implications of thermal transfer and compressibility. To match the static characteristics of finite-element models intended for aeroelastic scaling, [START_REF]An application of structural optimization in wind tunnel model design[END_REF] presented an optimization-based strategy. [START_REF] French | Aeroelastic model design using parameter identification[END_REF] later proposed a two-step strategy to also include dynamic aeroelastic scaling. First, they designed the structure by minimizing the discrepancies in static displacements. Then, they determined the values of non-structural masses through a minimization problem where the difference in mode shapes was minimized and the reduced modal frequencies were set as equality constraints. [START_REF] Pereira | Aeroelastic scaling and optimization of a joined-wing aircraft concept[END_REF] defined a method that optimized the natural frequencies and applied it to the design of a joined-wing wind tunnel scaled prototype to obtain equivalent dynamic and static aeroelastic scaling. Bond et al. (2012) proved that in addition to matching the natural frequencies, fitting the corresponding vibration modes is also required. [START_REF] Richards | Design of a scaled rpv for investigation of gust response of joined-wing sensorcraft[END_REF] established a comparison between a direct modal one-step strategy and a two-step methodology where the modal optimization was performed after establishing the stiffness distribution.

Like Bond et al. (2012), other authors have accounted for static scaling nonlinearities, extending the strategy outlined by [START_REF] French | Aeroelastic model design using parameter identification[END_REF] to include geometric nonlinearities in the aeroelastic scaling method. [START_REF] Ricciardi | Nonlinear Aeroelastic Scaling of a Joined Wing Aircraft[END_REF] adapted the twostep strategy suggested by [START_REF] Richards | Design of a scaled rpv for investigation of gust response of joined-wing sensorcraft[END_REF] to include the equivalent nonlinear static deflections in the first optimization cycle. These two approaches were developed in the context of the experimental testing of a joined-wing SensorCraft prototype. [START_REF] Ricciardi | Nonlinear Aeroelastic-Scaled-Model Optimization Using Equivalent Static Loads[END_REF] later proposed to use a single-step strategy where linear and nonlinear static responses were scaled while satisfying the natural frequency constraints. [START_REF] Wan | Geometrically Nonlinear Aeroelastic Scaling for Very Flexible Aircraft[END_REF] provided a scaling method considering geometrical nonlinearities for very flexible airplanes, extending the linear scaling factors and similarity rules. Recently, [START_REF] Ricciardi | Nonlinear Aeroelastic Scaled-Model Design[END_REF] described a systematic strategy for the definition of scaled prototypes that are equivalent in terms of the aeroelastic properties. They proposed an optimization formulation that sought to match vibration and buckling modes, as well as the linear static displacements. They also highlighted the mode swapping and mixing phenomenon in the optimization process. Mas [START_REF] Colomer | Similarity Maximization of a Scaled Aeroelastic Flight Demonstrator via Multidisciplinary Optimization[END_REF] discussed the issue of aeroelastic scaling when flow similarity is not matched exactly. They also discussed the issue of mode crossing through the implementation of a mode tracking approach in the framework of classic aeroelastic scaling (i.e., with presumed flow similarity). [START_REF] Cavallaro | Challenges, Ideas, and Innovations of Joined-Wing Configurations: A Concept from the Past, an Opportunity for the Future[END_REF] offered a thorough literature study on the works related to the joined-wing Sensor-Craft prototypes for aeroelastic scaling. Additionally, [START_REF] Afonso | A review on non-linear aeroelasticity of high aspect-ratio wings[END_REF] assessed the state-of the-art on the topic of nonlinear aeroelastic scaling of wings presenting a high aspect ratio. [START_REF] Spada | Nonlinear aeroelastic scaling of high aspectratio wings[END_REF] illustrated the implementation of the two-step technique to high aspect ratio wings. More recently, [START_REF] Pontillo | Flexible High Aspect Ratio Wing: Low Cost Experimental Model and Computational Framework[END_REF] introduced a suite of techniques for the design of a cantilevered wing model that replicates the static and dynamic response of a conceptual prototype. Regarding the dynamic aeroelasticity and multidisciplinary optimization, [START_REF] Jonsson | Flutter and postflutter constraints in aircraft design optimization[END_REF] analyzed how flutter constraints can be integrated into multidisciplinary optimization of aircraft earlier in the design phase.

All the mentioned references assume that flow similarity is ensured between the two scales, or that at least the discrepancies in nondimensional airflow parameters (e.g., Mach and Reynolds numbers) can be neglected. This may not, however, be applicable to some flying prototypes with a limited flight envelope (in terms of altitude and speed). Consequently, it is necessary to establish a method for the design of aeroelastically scaled models that exhibit the same scaled aeroelastic behavior.

In general, classical aeroelastic scaling theory is based on the fact that flow similarity exists between the two scales. The existing works cited in this article assume flow similarity or expect that the existing differences are small enough to be neglected.

In this paper, we propose a methodology to establish the design of scaled aircraft wings in the case where the differences in the nondimensional airflow conditions result in substantial discrepancies in the static response (i.e., in the total aerodynamic load and structural wing deflection). To do so, we maximize the static aeroelastic similarity of the wing as a whole-in terms of the scaled total force and wing deflection. Once the airflow is not similar, the local data of the flow field and the structural response cannot be used to extrapolate the local response of the full-size aircraft.

To address this issue, we propose a multidisciplinary design optimization (MDO) approach to modify the geometry of the wing (while preserving the overall scaled span) and the structural properties simultaneously to obtain an equivalent scaled response (in terms of the total load and overall wing deflection).

Typically, MDO in the context of aircraft design is applied to maximize performance, such as aircraft mass or fuel burn. For example, [START_REF] Kenway | Multipoint High-Fidelity Aerostructural Optimization of a Transport Aircraft Configuration[END_REF] apply high-fidelity aerostructural MDO to solve takeoff mass and fuel burn minimization problems. However, unlike typical applications, in this paper we use MDO to minimize the differences in the aeroelastic response between a scaled model and its full-size counterpart.

Given an MDO problem, there are many architectures that can be used to solve it. The work by [START_REF] Martins | Multidisciplinary design optimization: a survey of architectures[END_REF] presents a survey of all the available architectures for MDO problems. As we will detail on Section 2.3, in this paper we use the multidisciplinary feasible (MDF) approach.

The method that we propose is intended to be used to assess the static aeroelastic qualities of a given aircraft configuration, by matching the overall wing load and deflection, if flow similarity cannot be achieved.

Optimization Problem for Static Aeroelastic Similarity

In this paper, we focus on the static aeroelastic scaling in the case where flow similarity cannot be achieved. This is motivated by the limitations in altitude and speed of scaled flight demonstrators, which make it impossible to match both Mach number and Reynolds number. As shown by [START_REF] Bisplinghoff | Aeroelastic Scaling of a Joined Wing for Nonlinear Geometric Stiffness[END_REF], preserving both shape and flow parameters is a requirement for complete aeroelastic similarity. Therefore, if there is no flow similarity, the final deformed state cannot be matched if the scaled structural stiffness and the aerodynamic shape are kept constant.

The goal of this paper is to develop a method for designing a scaled wing model whose overall deflection is as close as possible to the scaled response of the fullscale aircraft for the same angles of attack, even in the presence of a non-similar flow in the scaled model. To achieve this goal, we relax some geometrical design variables, namely the chords that define the wing planform, as well as the sweep angle and the mounting angle of the wing. The geometrical design variables have an impact on the stiffness of the wing because they determine the shape of the wingbox. Apart from these geometrical variables, we also set structure-specific variables for the wingbox properties (e.g., panel thicknesses and stringer sections).

Effect of Flow Similarity on Static Aeroelastic Scaling

In the case where there is flow similarity (shown in Figure 1), no aerodynamic analysis is required to ensure complete aeroelastic similarity between both scales because the load equivalence is guaranteed through a nondimensional analysis of the aerodynamic equations.

In practice, it might not be possible to match the Mach number due to the operational limitations in airspeed and altitude of the scaled demonstrator. However, we can match the scaled outer mold line (OML), angle of attack, and scaled stiffness. In that case, the final in-flight shape of the wing of the scaled demonstrator, as well as the load distribution, are not scaled versions of the ones on the reference wing at the reference conditions. This is because for a given geometry at a certain angle of attack, the aerodynamic loads are not equivalent unless the nondimensional parameters are the same.

To understand the effect of Mach number, we can use the Prandtl- [START_REF] Glauert | The effect of compressibility on the lift of an aerofoil[END_REF] correction to provide an estimate of the aerodynamic forces at high subsonic Mach numbers based on the forces computed in the incompressible regime (M < 0.3). The Prandtl-Glauert correction for the lift coefficient of a two-dimensional airfoil is given by

C l = C l0 √ 1 -M 2 , ( 1
)
where C l0 is the lift coefficient at incompressible conditions, which in turn is given by

C l0 = l 0 1 2 ρV 2 c , (2)
where l 0 is incompressible lift per unit span. If we tried to reproduce the static aeroelastic behavior of a typical airliner (M ∼ 0.85) using a model with scaled stiffness, shape, and angle of attack flying at a lower Mach number (e.g., M ∼ 0.3), we would find that the aerodynamic loads, and hence the displacements, would be noticeably lower, according to the Prandtl-Glauert (1). This is illustrated by Figure 2 and shown with a particular example of application on Section 5.

If we now consider the case where we attempt to optimize the wing geometry and structural properties at the same time to match the static aeroelastic deflections despite significant differences in the flow nondimensional parameters, we cannot use the nondimensional analysis of the aerodynamic equations, as in this case the scaled OML is not preserved. In that case, we need the in-flight deflections of the reference aircraft for a particular angle of attack. To evaluate the fitness of a particular design, we need to perform a complete static aeroelastic analysis for the given design variables and the Mach number of the scaled demonstrator. Then we can compare the in-flight shape to the scaled version of the reference one and evaluate the lift and stress constraints.

Aeroelastic Multidisciplinary Analysis

In this paper, we use Panair (Magnus and Epton, 1980), a panel aerodynamics code implementing the Prandtl-Glauert compressibility corrections coupled to Nas-tran95 (MacCormick, 1972), a finite element solver. An independent aerodynamic analysis is performed for each deformed configuration of the wing. The XDSM (Lambe and [START_REF] C032150 ; Lambe | Extensions to the design structure matrix for the description of multidisciplinary design, analysis, and optimization processes[END_REF] diagram corresponding to this MDA is shown in Figure 3. Within each MDA cycle, the aerodynamic loads are computed and transferred to the FEM model. Then, the structural displacements are solved for and interpolated back on the aerodynamic grid. We have two flight conditions: cruise (to ensure steady flight), and limit load (to ensure structural integrity). Therefore, a separate MDA is solved for each one.

Optimization Formulation

Once we have established the MDA for this problem, we present the formulation of the optimization problem. To solve the optimization problem, we adopt the multidisciplinary feasible (MDF) approach, which is described by [START_REF] Martins | Multidisciplinary design optimization: a survey of architectures[END_REF]. Using the MDF approach, the complete MDA is solved completely for each iteration of the design variables. This approach is both easy to implement and computationally efficient. Figure 4 shows the XDSM for the aeroelastic optimization. Because each optimization loop includes an MDA for each angle of attack α i , they are represented in a stacked manner as they can be performed independently. One of the goals of this study is to ensure that the scaled aerodynamic loads and wing deflections of the scaled model are equivalent to the ones of the reference aircraft for the same angles of attack. In that manner, the response of the scaled wing loads and displacements with respect to the angle of attack is similar. The values of α i correspond to fixed values of the angle of attack for which the lift force on the reference aircraft is known.

X 0 a , X s y t α, V, ρ a , S w , M, X 0 a E, ν, X s ,
The goal of this optimization problem is to find the wing design parameters that minimize the error in the in-flight shape with respect to the scaled wingtip deflection of the reference aircraft while keeping the similarity at certain angles of attack. For that purpose, we require the scaled lift to be matched for two angles of attack. By doing so, we ensure that the derivative of the aerodynamic load with respect to the angle of attack C Lα is preserved. If similarity of airflow, aerodynamic shape, and wing elasticity existed between the reference aircraft and the scaled model, the identity of C Lα would be automatically satisfied. Since in our case similarity in airflow and external shape is no longer preserved, we ensure the equality of C Lα using optimization constraints.

To ensure the structural integrity for the most critical flight condition, we require the von Mises stress in any point of the wingbox structure to be lower than the allowable stress of the structural material for the +2.5g maneuver condition.

The optimization problem formulation that we propose is

minimize f = u m tip -λ l u r tip 2 2 w.r.t. x ∈ R n subject to: C m,c L -C r,c L C r,c L = 0, C m,2.5g L -C r,2.5g L C r,2.5g L = 0, max(σ V M ) -σ y σ y ≤ 0, x m,r c 2 -λ l x r,r c 2 = 0, x m,b c 2 -λ l x r,b c 2 = 0, x m,t c 2 -λ l x r,t c 2 = 0, (3)
where u m tip is the wingtip displacement vector of the scaled model wing during cruise flight, u r tip is the corresponding quantity for the reference aircraft, and λ l is the overall geometrical scale factor (i.e., wing span ratio). The objective is to minimize the squared norm of the difference between the scaled reference displacement vector of the wingtip and the one of the scaled model. Thus, we find a wing with a similar overall flexibility while ensuring that the total scale wing load is equivalent.

[t] 0 , [s] 0 , c 0 r , c 0 b , c 0 t , Λ 0 , θ 0 , k 0 , x 0 r y t(0) αc, V, ρa, Sw, M E, ν, ρs u r tip W, ρa, V, Sw, n σy, F S, n x r r , c r r c r b , c r t Optimizer c 0 r , c 0 b , c 0 t Λ 0 , θ 0 ,
The nondimensional lift is matched through the lift coefficients, which are already nondimensional. The cruise lift coefficient of the reference aircraft is

C r,c L = W r g 1 2 ρ r V 2 r Sr , (4)
where W r is the weight (in units of mass) of the reference aircraft. For the case of the scaled model, the lift coefficient is

C m,c L = L α=αc 1 2 ρ m V 2 m λ 2 l Sr . (5)
In Eq. ( 5), the reference surface is λ 2 l S r and not the actual wing surface of the model wing. This is because we want to match the scaled lift of the reference aircraft (which is determined by λ l ) for the reference cruise angle of attack. For the lift constraint at 2.5g, the definition of the lift coefficient is the same except that the apparent weight is now 2.5W , and the lift force is computed at the angle of attack that produces that lift on the reference aircraft (α = α 2.5g ).

The inequality constraint on Eq. (3) ensures that the von Mises stress for any element on the structure is lower than the allowable stress, at the 2.5g maneuver condition, as defined by the Civil Aviation Regulations FAR-25 (2020). Since the optimizer that we use to solve this problem does not require the derivatives of the objective function and constraints, we use the max() function applied to the von Mises stress of all elements-whose derivatives might are continuous. The last three equality constraints in Eq. ( 3) ensure that the longitudinal positions of the mid-chord points of the three characteristic wing sections of the scaled model remain the same as the scaled longitudinal positions of the reference aircraft.

In the next section, we apply this optimization formulation a scaled version of the undeflected Common Research Model (uCRM) wing [START_REF] Brooks | Benchmark Aerostructural Models for the Study of Transonic Aircraft Wings[END_REF].

Methods and Tools

Method for Displacement and Load Transfer

In this section, we present the method for force and displacement transfer between the aerodynamics and structures disciplines, which is needed for the MDA described in the previous section.

The aerostructural MDA involves four state vectors: the displacements of the structure nodes u s , the displacements of the aerodynamic mesh u a , the forces acting on the aerodynamic grid points f a , and the aerodynamic forces acting on the nodes of the structural model f s .

Given that the grid points of the structural and aerodynamic models do not coincide, it is necessary to transfer the displacements on the structural nodes u s to the aerodynamic mesh to determine its displacements u a . We use the displacement and load transfer method described by [START_REF] Rendall | Unified fluid-structure interpolation and mesh motion using radial basis functions[END_REF]. This method is based on creating a displacement field that satisfies the interpolation condition on the structural nodes, on which the displacements are known from the FEM analysis. Through the use of radial basis functions (RBF), defined by [START_REF] Broomhead | Radial basis functions, multi-variable functional interpolation and adaptive networks[END_REF], each one of the three components of the displacement field (we write the equation for the x component, but the same holds for y and z) is written as

u x = Ns i=1 β x i φ( x -x i ) + γ x 0 + γ x x x + γ x y y + γ x z z, ( 6
)
where x is the vector of space coordinates, x i represents the nodes at the centers of the radial basis functions, and φ(r) is the RBF. We chose to use the thin plate spline (TPS) RBF, which can be written as φ(r) = r 2 ln r. This choice was based on the recommendations given by [START_REF] Lombardi | Radial basis functions for inter-grid interpolation and mesh motion in FSI problems[END_REF]. After comparing many available functions, they concluded that TPS functions yielded the best, most accurate, and most robust results for interpolation because they do not require any parameter tuning. The terms β x i are the RBF coefficients and γ are the coefficients of the linear polynomial part. The interpolated displacement field evaluated on the aerodynamic grid points is obtained through a matrix-vector product from the vector of known structural displacements u s and a matrix [H]. This matrix [H] is obtained by imposing the interpolation condition and by evaluating each RBF on the aerodynamic points. The displacement transfer is then expressed as where [H] depends only on the positions of the structural nodes and the grid points of the aerodynamic mesh, as well as on the RBF type (such as Gaussian or TPS). Figure 5 illustrates how the displacement field, known at the structure grid (in orange), is transferred to the aerodynamic grid (in green).

u a = [H]u s , (7)
Once the displacement transfer matrix is established, [START_REF] Rendall | Unified fluid-structure interpolation and mesh motion using radial basis functions[END_REF] determined the matrix for the transfer of the aerodynamic forces from the aerodynamic grid f a to the structural nodes f s . This was done by using the principle of virtual work for the conservation of the energy. The virtual work is

δW = δu s T f s = δu a T f a , (8)
where δu s represent the virtual displacements of the structural nodes used for the displacement interpolation and δu a are the virtual displacements of the aerodynamic grid points. By using the displacement matrix [H], the virtual displacements of the aerodynamic grid are

δu a = [H]δu s . (9)
After substituting Eq. ( 9) into Eq. ( 8), the aerodynamic forces on the structural nodes are

f s = [H] T f a . ( 10
)
The construction of the interpolation matrix for each geometry of the aerodynamic and structural meshes allows us to exchange displacements and forces through a matrix-vector product.

Geometry and Mesh Generation

For the generation of the OML from the planform design variables, we use the baseline uCRM wing geometry for the twist distribution and the vertical coordinates of the wing sections. The airfoil shapes are those of the original CRM wing, including the relative Y X Z Y X Z Fig. 6: The structural mesh (orange) is initially known only in its baseline configuration (left). Through the new shape of the OML (in green), it is automatically adapted to new shape using the RBF methodology (on the right). For illustration purposes, the sweep is increased while keeping chords and span constant. thickness, relative camber, and spanwise position of the wing sections. This baseline geometry defines the airfoil shape and twist at eight sections from the root to the wingtip. The planform design variables determine the leading edge position and chord length. The spanwise and vertical position of the sections, as well as their relative thickness and camber, are kept constant. The values for these correspond to the scaled reference wing. From this definition of the wing through its section properties, we define the jig OML corresponding to each set of design variables.

The OML and a structured surface mesh for the aerodynamics panel code are generated using the GMSH software package, developed by [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and postprocessing facilities[END_REF]. For each vector of design variables, the topology of the structured surface mesh remains the same.

The geometric variables described above first determine the shape of the jig OML, which then determines the new shape of the jig wingbox. To find the structural wingbox geometry for a given OML, we use a morphing technique based on the RBF method that we used for the interpolation of displacements, defined by [START_REF] Rendall | Unified fluid-structure interpolation and mesh motion using radial basis functions[END_REF] (see Section 3.1). This is done by treating the changes in the OML as displacements, which are interpolated to find the corresponding changes to the wingbox geometry and mesh, as shown in Figure 6.

Application to the Scaled uCRM Wing

In the previous section, we defined the optimization problem used to match the static aeroelastic behavior of scaled demonstrators when flow similarity cannot be achieved. In this section, we apply the aforementioned methodology to a scaled version of the uCRM wing. The uCRM wing, established by [START_REF] Brooks | Benchmark Aerostructural Models for the Study of Transonic Aircraft Wings[END_REF], is based on the original CRM wing, defined by [START_REF] Vassberg | Development of a common research model for applied CFD validation studies[END_REF], and represents the jig shape of a wing such that when in flight its shape matches the CRM one at its nominal flight condition.

Definition of the Scaling Factors

In this example problem, we consider, as the reference wing, the uCRM wing flying at 35, 000 ft and M = 0.85 with an aircraft weight of 500, 000 lbm (∼ 226796 kg) 1 For the scaled model, whose design we establish through optimization, we consider a reduced model scaled by an overall factor λ l = 1/5. For the flight conditions of the model, we consider an altitude of 2, 000 m (∼ 6, 562 ft). From the Froude number matching, as described by [START_REF] Pires | LINEAR AEROELASTIC SCALING OF A JOINED WING[END_REF], the scaled airspeed is

V m = √ λ l V r = 112
.77 m/s. At the scaled model altitude, and using the international standard atmosphere model, this yields a Mach number of M = 0.34.

As described by [START_REF] Pires | LINEAR AEROELASTIC SCALING OF A JOINED WING[END_REF], the scaling factor and the air density ratio λ ρ = ρ m /ρ r determine the total mass of the scaled model,

W m = λ ρ λ 3 l W r , (11)
which dictates the amount of lift that needs to be produced for the several flight conditions.

Design Variables

The geometric design variables of this problem are the chord lengths at three span sections (root, break, and tip-c r , c b , and c t respectively), the sweep angle Λ, the wing mounting angle (the angle between the root section of the wing and the longitudinal axis of the fuselage), as well as the longitudinal position of the wing leading edge at the root. The thickness-to-chord ratio of the sections is obtained by multiplying the baseline distribution of relative thicknesses along the wing span by a scaling factor k. We also set k as a design variable to control the flexibility of the wing structure through the height of the wingbox. For the aerodynamic analysis, we define the angle of attack as the angle of the longitudinal axis of the fuselage with respect to the incoming airflow. The geometric planform variables are illustrated in Figure 7.

The structural design variables are the plate thicknesses of the shell elements of the wingbox, as well as the stringer section modeled by beam elements. For the definition of the shell thicknesses, the wingbox is divided into a total of 12 different thickness regions: 6 of them define the thickness of the upper, skins lower skins, and ribs, at 6 wing segments along the span. The other 6 regions define the thickness of the elements on the front and rear spars: the two flanges and the web. For the stringer sections, there are a total of six variables that define the cross-section values of all the stringers in the model, which belong to six different groups. These regions of different stringer sections are displayed on Figure 9. With these structural design variables, we can tailor the wingbox stiffness.

In contrast to the usual structural optimization problems that aim to improve the performance of an existing design, here we aim to find the whole design of the scaled wing.

The displacement interpolation matrix [H] (defined in Section 3.1) has to be recalculated for each optimization iteration before performing the MDA. This is because [H] depends on the coordinates of both aerodynamic and structural meshes, as shown in Figure 19.

Definition of the Reference Quantities

According to the optimization problem formulation (3), we need the wingtip displacements of the in-flight reference wing. To obtain them, we perform an MDA loop with the baseline full-scale wing geometry and structure at the nominal flight conditions. To find the correct angle of attack, we perform several iterations until the generated lift matches the reference aircraft weight, W r g. To determine the angle of attack at the limit load condition, we perform the same operation, but with a target lift of 2.5 times the weight of the reference aircraft. This yields a cruise angle of attack of α c = 1.34 • and α 2.5 = 9.66 • for the limit load condition.

The reference document for the CRM finite element model specifies the following elastic properties: Elastic modulus E = 68900 MPa, Poisson's ratio ν = 0.31, and density ρ s = 2795.67 kg/m 3 . Since the allowable stress is required to evaluate the stress constraint and that this value is not available in the CRM report, we use the maximum von Mises stress computed for the baseline wing at 2.5g. This yields a yield stress of about 600 MPa, which is representative of a 7068 aluminum alloy.

If we used an aluminum alloy with a lower allowable stress than the maximum observed on the reference wing at limit load condition, the target displacements could probably not be matched, since the flexibility of the scaled wing would be limited by the stress constraint. In that case, the optimizer would need to increase the structural thicknesses to satisfy the stress constraint, thus reducing the required flexibility to match the in-flight shape.

Optimization Problem Setup

The objective function and all the constraints are evaluated for each optimization iteration at the end of the MDA cycle. Figure 4 shows the XDSM diagram of the optimization problem described in this section. A more detailed version of this diagram is shown in the Appendix 6. The optimization problem is summarized in Table 1, which lists the 25 design variables with the corresponding bounds and 6 constraints. The values of [t 0 ] and [s 0 ] are estimates for the structural thicknesses and stringer cross-section variables, and are obtained in the same manner as for the modal scaling problem, in Section 4.5, as described by Eqs. ( 19) and ( 22). For the values c 0 r , c 0 b , and c 0 t , they are obtained by scaling the baseline section chords by the overall scaling factor λ l .

The optimization algorithm we use is COBYLA [START_REF] Powell | A direct search optimization method that models the objective and constraint functions by linear interpolation[END_REF], a gradient-free optimizer, available within the SciPy library [START_REF] Jones | Scipy: Open source scientific tools for python[END_REF]. We use this approach because neither Panair or Nastran95 provide the derivatives of their outputs.

Establishing the starting point

For the scaling of the natural frequencies and the overall mass, [START_REF] Pires | LINEAR AEROELASTIC SCALING OF A JOINED WING[END_REF] detailed how to obtain the ratios λ between the scaled quantities and the reference ones, which are

λ ω = 1 √ λ l , (12)
and

λ m = λ ρ λ 3 l . ( 13
)
Given that the COBYLA optimization path depends on the starting point, we developed a method for finding a set of design variables that is close enough to the optimum. This method is based on the Euler-Bernoulli beam theory and the bending properties of thin-walled structures to find an approximation for scale-down factors of the thicknesses of the structural panels and stringer sections. We consider a cantilever beam of length L (corresponding to the full scale), bending inertia I, total mass m, and Young's modulus E. We start by computing the bending inertia I , whose length and height are scaled according to the ratio λ l = L /L. The overall length ratio affects the length, height, and width of the beam uniquely. Since we assume that the same material is used for both beams, E = E. The two beams are shown in Figure 10. By considering the equation for the natural frequencies of the bending modes of a cantilever beam

ω i = K i L 2 EIL m , (14)
where K i denotes the constant value that depends on the vibration mode, we obtain the frequency ratio as

λ ω = ω i ω i = Ki L 2 E I L m Ki L 2 EIL m . ( 15
)
Objective Function Dimension In-flight shape difference minimization

f = u m tip -λ l u
C m,i L -C r,i L C r,i L = 0 2 Stress constraint max(σ 2.5g V M )-σ y σ y < 0 1
Longitudinal positions of the chord centers

x m,i c 2 -λ l x r,i c 2 = 0 3
Total constraints: 6 Substituting Eqs. ( 12) and ( 13) into Eq. ( 15), we get

λ I = I I = λ ρ λ 5 l . (16)
After establishing the bending inertia ratio λ I , we determine the thickness and stringer section ratios, λ t and λ s , respectively. By making assuming that the thickness of the beam flanges (t) are much smaller than their distance from the neutral axis (h), we estimate the contribution of a flange of width w (shown in Figure 11) to the bending inertia of the beam I as

∆I t = wth 2 . ( 17
)
We determine the bending inertia ratio by applying Eq. ( 17) to the geometrically scaled beam

λ I = w t h 2 wth 2 = λ 3 l λ t . (18)
Then, we establish the thickness ratio for the initial design as

λ t = λ ρ λ 2 l . (19)
In a similar manner as presented above for the thickness ratio, we compute the contribution of a stringer placed at a distance h from the neutral axis with a cross-section s (see Figure 12) as

∆I s = sh 2 . ( 20
)
Then, the bending inertia ratio of the contribution of the stringer is As in the case for the plate thickness ratio, we find the ratio of stringer cross-sections by substituting Eq. ( 16) into Eq. ( 21), thus giving

λ I = s h 2 sh 2 = λ 2 l λ s . (21)
λ s = λ ρ λ 3 l . (22)
Using this method, we obtain a first guess of thicknesses and cross-sections of the different components of the structure that we aim to scale, to use it as the initial set of design variables for the optimizer. We obtain these design variables by simply multiplying the thicknesses and cross-sections of the full-scale model by λ t and λ s , respectively.

Results and Discussion

As previously mentioned we used the gradient-free optimizer COBYLA [START_REF] Powell | A direct search optimization method that models the objective and constraint functions by linear interpolation[END_REF] to solve the optimization problem (3). The stopping criterion is based on the absolute tolerance on the size of COBYLA's trust region, which we set to 10 -3 .

The optimization history (Figure 15) shows that the optimum satisfies the stress constraint. The objective function (the wingtip displacement error) converged to 1.5%. Both lift constraints-the cruise condition one and the limit load condition-are satisfied within an error lower than 1%. The lift constraint at the cruise condition has an error lower than 0.1%. The three equality constraints-shown on Figure 16-show that the longitudinal positions of all the wing sections are well preserved for the best design found. This can also be visualized on the optimized planform (Figure 13). The optimization history shows that promising values of the objective function that satisfy the stress constraint were already found before iteration 100. However, the lift constraints are much below zero for such cases. These probably correspond to wing designs with a very low flexibility that reached the desired level of wing displacement, but under a wing load lower than the target value. The in-flight shape corresponding to the optimized design is shown on Figure 18.

To see the importance of performing an aerodynamic and structural re-design when the Mach number cannot be matched, we compute the errors in terms of loads and displacement if the aerodynamic geometry and scaled flexibility of the reference aircraft were left constant. To do so, we present the results that we would get if we used the exact scaled aerodynamic geometry Fig. 13: To compensate for the lift loss due to lower Mach, the optimized planform (red) has an increased wing area compared to the scaled reference planform (blue). of the reference aircraft along with a structure having the exact scaled shape and the scaled stiffness and using the flight Mach value corresponding to the scaled demonstrator conditions. In that case, whose in-flight shape is shown on Figure 17, the error in terms of the air loads is 16%, while the error in wingtip displacement is 17.5%. These results show that the proposed method can be used effectively to ensure loads and displacements similarity if the airflow conditions of the scaled model are substantially different to the ones of the reference aircraft.

After having presented the results of the last example of application, in the following section we present the conclusions and perspectives of this paper.

Conclusion

In this paper, we maximize the similarity of the static response of the wing as a whole (i.e., in terms of the total air load and the overall wing deflection). This enables to assess the aeroelastic qualities of a certain wing concept, but not to extrapolate local data of the airflow or the structural response of the full-scale aircraft.

The dominant effect of operating the scaled model at a lower Mach number than its full scale counterpart is the reduction in the lift forces for the same angle of attack. With the baseline geometry and scaled stiffness, the relative error in the total air load in cruise flight is 16%, while the relative error in the wingtip displacement is 17.5%. These substantial errors justify the changes in the wing geometry and the structural properties when the airflow conditions change substantially. With the optimized wing, we observe an increase in the area to produce equivalent lift with the imposed angles of attack of the reference conditions. The results show a good in-flight shape matching, with an error of less than 1% in the lift constraints. Since the shape of the wing is modified, which has an impact on the structural stiffness, the properties of the structural components are modified accordingly to give the wing the appropriate stiffness in spite of the changes in geometry.

The results in this paper are obtained using the panel method with compressibility corrections. However, more physically accurate results could be obtained by using the same optimization formulation and simply replacing the current aerodynamic model by a more realistic CFD analysis.

Fig. 1 :

 1 Fig. 1: Reference wing deflection and load distribution (a) and its scaled counterpart when complete aeroelastic similarity exists (b).

Fig. 2 :

 2 Fig.2: Reference wing deflection and load distribution (a) and its scaled counterpart when there is no flow similarity (b). In that case, the wing deflection and loading are no longer similar.

Fig. 3 :

 3 Fig. 3: XDSM diagram for the static aeroelastic MDA.

Fig. 4 :

 4 Fig. 4: XDSM diagram for the static aeroelastic MDAO. The orange box represents the MDA cycle. Several MDA cycles are stacked (one for each angle of attack) since they are performed independently one from another.

Fig. 5 :

 5 Fig.5: Transfer of the displacements from the structural grid (orange) to the aerodynamic surface mesh (green).

Fig. 7 :

 7 Fig. 7: Planform variables for the static optimization problem: longitudinal position relative to the fuselage x LE r , sweep angle Λ, root chord c r , break chord c b , and tip chord c t .

  These regions of different thickness are shown in Figure 8.

Fig. 8 :Fig. 9 :

 89 Fig. 8: Thickness regions for the FEM model, which consists of 12,475 shell elements and 10,804 nodes.

Fig. 10 :

 10 Fig. 10: Full scale and scaled cantilever beams.

Table 1 :

 1 Optimization problem statement of the uCRM wing. The baseline thickness [t 0 ] and cross-section [s 0 ] vectors are the ones defined in Section 4.5, as the design variables and the scale are the same. For the baseline chords, the values are: c 0 r = 2.72 m, c 0 b = 1.45 m, and c 0 t = 0.55 m.

Fig. 11 :

 11 Fig. 11: Plate thicknesses on different scaled beams.

Fig. 12 :

 12 Fig. 12: Stringer sections on different scaled beams.

Fig. 14 :

 14 Fig. 14: Stress distribution (von Mises) for the limit load condition of the best design found. Values in Pa.

Fig. 15 :Fig. 16 :

 1516 Fig. 15: Evolution of the objective function, the maximum stress constraint, and the two lift constraints with the number of iterations.

Fig. 17 :

 17 Fig. 17: Not using any re-design to account for difference in Mach number leads to a discrepancy between the wingtip displacement of the reference wing under reference conditions (blue) and the same scaled wing in test conditions (green). Front view (a), side view (b), top view (c), and detail of the wingtip (d).

Fig. 18 :

 18 Fig. 18: The wingtip displacement of the optimized design (red) matches the one of the target shape (blue). Front view (a), side view (b), top view (c), and detail of the wingtip (d).

https://commonresearchmodel.larc.nasa.gov/files/

2014/02/CRM_wingboxFEM_description_1.pdf

Acknowledgements This work is part of the activities of ONERA-ISAE-ENAC joint research group. This work was funded by the Airbus CEDAR Chair and the EU Project 658570-NextGen Airliners.

Fig. 19: XDSM diagram for the static aeroelastic MDAO. The first MDA loop determines the cruise angle of attack to evaluate the objective function as well as the cruise lift constraint. The second MDA loop determines the angle of attack of the limit load condition to evaluate the lift constraint at the limit load as well as the stress constraint. Both MDAs are independent from each other and can be executed in parallel.

Conflict of interest

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

Replication of results

The results presented in this paper can be replicated using the publicly available data of the CRM wing [START_REF] Vassberg | Development of a common research model for applied CFD validation studies[END_REF], the uCRM wing [START_REF] Brooks | Benchmark Aerostructural Models for the Study of Transonic Aircraft Wings[END_REF], as well as the open-source structural solver NAS-TRAN95 [START_REF] Maccormick | The NASTRAN User's Manual:(level 15)[END_REF], and the open-source aerodynamics solver PANAIR (Magnus and Epton, 1980).