# Digital Thread-Based Optimisation Framework for Aeronautical Structures: A Vertical Tail Plane Use Case

Sébastien de Longueville, Christophe Bouvet, Emmanuel Bénard, Joël Jézégou, Yves Gourinat

► To cite this version:

Sébastien de Longueville, Christophe Bouvet, Emmanuel Bénard, Joël Jézégou, Yves Gourinat. Digital Thread-Based Optimisation Framework for Aeronautical Structures: A Vertical Tail Plane Use Case. Aerospace, 2024, 12 (1), 10.3390/aerospace12010002. hal-04918188v2

HAL Id: hal-04918188

https://hal.science/hal-04918188v2

Submitted on 5 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Digital Thread-Based Optimisation Framework for Aeronautical Structures: A Vertical Tail Plane Use Case

Sébastien de Longueville <fnref n="1" />, Christophe Bouvet <fnref n="2" />,* ID, Emmanuel Bénard <fnref n="1" /> ID, Joël Jézégou <fnref n="1" /> ID and Yves Gourinat <fnref n="2" /> ID

**Abstract:** In this modern era, the constant increase in computational and sensing power has lead to the development of multiple data-driven concepts. Amongst these, the 'digital thread' is an architecture that aims at optimising the knowledge of a system by merging prior knowledge of the product with information from multiple stages of its lifecycle in order to improve the performance of new products to be designed, thanks to the increased accuracy of this updated knowledge. Even though the use of these data-driven architectures is becoming increasingly widespread, most of the corresponding developments remain currently limited to the component level. In this respect, this article extends the application of the digital thread from the component level to the structural assembly level and enriches it with additional multi-physics considerations and non-linear failure constraints. To this end, it details the development of a digital thread dedicated to an aircraft vertical tail plane structure made of carbon fibre reinforced polymer, as well as the tools and models required to implement such an approach using Bayesian inference, multi-physics simulations, and empirical models.

**Keywords:** data-driven design; digital thread; digital twin; data fusion; Bayesian inference; composite structure modelling; ACCEM; aerodynamics; uncertainty analysis; multidisciplinary design; structural optimisation; design under uncertainty

<figure><img src="image_2.png" /></figure>

check for
updates

Received: 1 November 2024

Revised: 16 December 2024

Accepted: 17 December 2024

Published: 24 December 2024

**Citation:** de Longueville, S.; Bouvet, C.; Bénard, E.; Jézégou, J.; Gourinat, Y. Digital Thread-Based Optimisation Framework for Aeronautical Structures: A Vertical Tail Plane Use Case. *Aerospace* 2025, 12, 2. https://doi.org/10.3390/aerospace12010002

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

# 1. Introduction

## 1.1. Context

Over the last decades, data-driven applications have greatly benefited from the rise of powerful smart systems [1], collecting an increasing amount of data throughout their lifecycles. These data are used inter alia to improve the design performance of these systems by more accurately characterising their actual state and environment [2–4]. This has led to the emergence of data-driven strategies such as *Digital Twins*, or *Digital Threads* [5,6].

A digital twin [7] is a virtual replica of a product, with a bidirectional connection between this replica and the real system to which it is paired (often referred to as the *physical asset*). It is adapted to the specificities of a given instance and is updated throughout the life of a product according to the history of punctual events encountered by this system, such as maintenance, repairs and associated feedbacks from the controls, or continuous variations in the product properties due to the normal wear and tear of the system and variations in the operational environment of the product [8].

<footnote n="1">1 Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, 31400 Toulouse, France; sebastien.de-longueville@isae-supaero.fr (S.d.L.); emmanuel.benard@isae-supaero.fr (E.B.); joel.jezegou@isae-supaero.fr (J.J.)</footnote>

<footnote n="2">2 Institut Clément Ader (UMR CNRS 5312)-ISAE-SUPAERO/INSA/IMT/Mines Albi/UPS, Université de Toulouse, 31400 Toulouse, France; yves.gourinat@isae-supaero.fr</footnote>

<footnote>* Correspondence: christophe.bouvet@isae-supaero.fr</footnote>
Digital threads and twins can be used to simulate the performances and be-
haviour of a system in different operational conditions, for example to predict system
failure [9–11], for health monitoring [12], and for condition-based maintenance, providing a
safer and more effective approach than the traditional time-based and usage-based preven-
tative maintenance strategies [13]. They enhance data management in Industry 4.0 [14]. In
smart manufacturing, digital twins play a crucial role in connecting physical facilities and
their digital components [15,16]. They can also be used for design and analysis purposes to
assist dynamic decisions based on up-to-date knowledge of the internal properties of the
system [17].

In this respect, several applications have emerged in the military aeronautics industry, such as for the Lockheed Martin F-35 program [18] or in the US Air Force [19–21], as well as for the definition of standardized architecture and implementation strategies for structural components [22] or for commercial airlines [23] and the proof of concept for a minimum viable product [24]. Eskue [25] highlights the importance of the digital thread in aerospace manufacturing and proposes a roadmap for implementing a fully realised digital thread. In addition, a multi-phase plan [26–30] was developed to implement the digital thread in engineering design through the detailed design phase of a wing rib and the related architectural implementations on the structural component [22,31], aiming at integrating multidisciplinary processes and analysis disciplines to more effectively evaluate product costs and performances [2,32]. The wing rib was also used as a proof of concept of STAnDD implementation [33] in a multi-phase plan. However, the existing applications remain focused on the component or conceptual level, while the industrial projects already developed at the system level, such as [18], remain focused on the integration capabilities provided by the digital thread, without taking advantage of the potential uncertainty reduction capabilities brought by a comprehensive digital thread [34].

Amongst the emerging data-driven architectures, the digital thread [22] opens new perspectives for the design of aeronautical systems under uncertainty [31,34,35] in order to make informed decisions for optimal product design. Data fusion is ensured using Bayesian inference to more accurately estimate the true properties of a system, based on related measurements, and to improve the accuracy of identification using difference sources of data [36]. The data flow presented below, which links each stage of the lifecycle of a ‘physical entity’ from design, to manufacturing, to field use, builds upon the framework proposed by Singh [22], illustrated in Figure 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of the digital thread implemented in [22], where the information gathered throughout the product lifecycle is fused at the design stage, before using it to make new design decisions.</figcaption></figure>
As depicted in Figure 2, a digital thread can therefore be considered a data structure that merges probabilistic knowledge of the system at a given design stage and all the resources and information available at that stage, with the aim of providing decision support tools to assist design processes [37,38] in the context of data-driven applications. To this end, this tool aims to reduce the level of uncertainties of the intrinsic state and characteristics of the system by gathering and filtering data generated during the different stages of the product lifecycle of the product, such as measurements from operational data and offline experiments (including loads, material properties, and manufacturing costs), design decisions, system requirements, including certification constraints and functional requirements, as well as all of the resources used with the given system of interest (software, etc.).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Data-driven design environment [39] gathering and updating knowledge of the system state, usage, and properties in a multidisciplinary environment used as an authoritative database.</figcaption></figure>

Authoritative data sources, data linkage, and model integration are essential in the process of implementing the digital thread framework [40–42]. A common definition and values are given for the digital thread in [34]. The importance of keeping traceability and continuity of these data [43] by connecting all information generated during the product lifecycle [34] is recognized, as well as the need to improve communication and collaboration by ensuring that all stakeholders can transparently and quickly access the necessary information from a single authoritative source of information [17,23,44]. Real-time monitoring of system performance and status improves the decision-making process [6,45] and helps reduce development time [46] while improving consistency [47].

## 1.2. *Challenges*

From the above, numerous challenges for the future development of the digital thread can be identified that essentially relate to scalability [34,48], flexibility [48], standardization [16], model relevance and accuracy, and the probabilistic approach for certification through exploitation of uncertain data [34,49,50].

The paper proposes original answers to some of these issues: scalability, uncertainty, filtering, and optimisation.
### 1.3. Objectives

The objective of this study is to provide a decision aid tool to assist design processes in the context of data-driven applications. For that purpose, this tool intends to reduce uncertainties about the intrinsic state and characteristics of the system by gathering and filtering data generated during the different stages of the lifecycle of the product, such as:

* Measurements (Loads, materials, manufacturing cost)

- Operational Data

- Offline experiments

* Design decisions

* System requirements

- Certification constraints

- Functional requirements

* Resources (software, etc.)

In that respect, this paper investigates how the design of complex aeronautical structures can benefit from a digital thread and, to this end, how such a concept could be enriched by means of multidisciplinary analyses, tools dedicated to aircraft design, and additional loading and failure considerations. As a use case, an aircraft structure that presents a good balance between the additional complexity required by the integration of aero-structural analyses in a digital thread and the computational requirements of the underlying models was selected. With this respect, the vertical tail plane (VTP) presents an optimal scale and level of complexity for such an integration, by being sufficiently complex to facilitate the integration of load cases and multidisciplinary analyses tools in the pre-existing architectural definition of the digital thread, while remaining manageable in terms of calculation time. This makes it an excellent intermediary framework for integrating advanced techniques within sufficiently complex systems involving a variety of physical domains in its dimensioning. This is why this study focuses on the implementation of a digital thread dedicated to the VTP of an A320-like aircraft.

To this end, three main features have been identified for this digital thread. On one hand, it aims at reducing the uncertainty related to the real state of a system of interest. This is achieved by merging the knowledge and information coming from measurements, design decisions or digital models. In addition, this data structure is enriched by interfacing the knowledge about this system with different models providing knowledge-based information regarding the geometry and the structural and aerodynamic properties of this system, and the cost required to manufacture this system [51]. Finally, this digital thread contains an optimiser, whose purpose is to identify the best design choice, based on the uncertainty level related to the knowledge of the real system.

For the first feature, a Bayesian filter is used in order to improve the knowledge of the system and its environment, based on measurements performed during the different phases of the lifecycle of the system (manufacturing, operation...). The choice of a Bayesian filter is supported by the good results shown in [35]. These measurements are then used to infer the real state of the system. To this end, the current use case focuses on the following states: the loads experienced by the structure during operation, the material properties of the system, and, finally, the detailed geometry resulting from the manufacturing process. Since this information is not directly measurable during the operational life of the vertical tail plane, the filter uses indirect methods to infer the system state based on these measurements.

This is achieved thanks to the second feature of the digital thread, which allows it to interact seamlessly with physical and empirical models linking the measured quantities to the actual state of the system.
As an example, the states of the digital thread stem from the analysis of the strains measured on the system in operation, the fail stresses of the materials acquired during coupon tests, and the duration of the different stages of the manufacturing process.

Finally, reliability-based design optimisation is used with gradient-based structural optimisation methods to compute the optimal design decisions based on the knowledge gathered and inferred by the digital thread, therefore providing decision aids to minimise the production cost of new generations of this VTP, based on the mass of this system and the time required for its manufacture. On the other hand, these decisions must ensure that the suggested design does not lead to structural failure of the system during normal operation. This is achieved by introducing failure criteria in the constraint function, based on static deformation and buckling considerations applied to the loads and material properties determined by the digital thread.

## 1.4. *Structure*

This paper starts by introducing the context related to data-driven applications dedicated to the design of structural components. Then, Section 2 provides the mathematical definition and the concepts required to implement a digital thread. Section 3 goes on to present the overall methodology and the architecture implemented in this study, while the parametric modelling tool used by the digital thread is described in Section 4, which outlines the geometric, structural, aerodynamic and manufacturing cost models used in this framework. A specific highlight is also made in Section 5 to specify the Bayesian filtering and the sensor measurement techniques meant for the uncertainty reduction which is at the core of this strategy, and to exhibit the resulting uncertainty reduction achieved on the structure investigated in this study. Thereafter, Section 6 summarises the mathematical formulation and the modelling of the optimisation problem, as well as the predominant phenomena and main assumptions affecting the underlying implementation choices, before presenting optimisation results and benefits enabled by the whole framework for the VTP use case. Finally, the conclusion summarises the key findings of the study and proposes directions for future work.

# 2. Towards a Digital Thread Implementation: Mathematical Definition and Requirements

## 2.1. *Concept*

The digital thread is a data-driven framework that seamlessly integrates information across the entire lifecycle of a product (design, manufacturing, operation, and maintenance), enabling real-time data exchange, enhanced decision-making, and continuous improvement through a connected digital environment.

This study applies the digital thread as a data-driven framework to optimise the design of complex aeronautical structures, in this case the vertical tail plane (VTP) of an A320-like aircraft, by integrating multidisciplinary models, Bayesian inference, and optimisation. It reduces uncertainty about the system's real state by merging lifecycle measurements, digital models, and design decisions, using a Bayesian filter to infer the real state of the system from essential data such as operational loads, material properties, and manufacturing geometries obtained from indirect measurements. The framework interacts seamlessly with structural, aerodynamic, and cost models to link measured quantities to the inferred states. Additionally, a reliability-based design optimisation (RBDO) approach [52] combines inferred knowledge with gradient-based methods to guide design decisions that minimise production costs while ensuring structural reliability by introducing failure criteria into the constraint function, based on static deformation and buckling considerations applied to the loads and material properties determined by the digital thread.
This integrated approach emphasises the digital thread’s capability to enhance decision-making across a product’s lifecycle by fusing data and multidisciplinary insights.

## 2.2. Definitions and Requirements

Following [22], the digital thread can be defined as a data structure containing simultaneously all the resources related to the given system at stage ($R_t$, see details in Section 4) and the probabilistic knowledge of the system at stage $t$ ($Q_t$, see further details in Section 5). Mathematically, the digital thread can therefore be expressed as a combination of these two sets of data:

$$ \mathcal{D}_t = (Q_t, R_t) \qquad (1) $$

Furthermore, this digital thread should be able to update these sets in order to take the time variation of the characteristics of the system of interest into account. In this respect, it should be able to update the probabilistic knowledge of the system, based on design decisions ($u_t$) or measurements ($z_t$) performed during the different phases of the product lifecycle. More specifically, this requires providing the digital thread with the following update capabilities, including design decisions and measurements:

$$ \mathcal{D}_{t+1} = \Phi_t(\mathcal{D}_t, u_t, z_t) \qquad (2) $$

Considering the relationship defined in Equation (1), this update capability can be defined as a combination of the following two independent update features, dedicated, namely, to the probabilistic knowledge and the resources related to the system:

$$ Q_{t+1} = \tau_t(\mathcal{D}_t, u_t, z_t) \qquad (3) $$

$$ R_{t+1} = \Psi_t(R_t, u_t) \qquad (4) $$

## 2.3. Probabilistic Knowledge

Among these two functions, on one hand, the uncertain state transition function ($\tau_t$) should update the probabilistic knowledge of the system ($Q_{t+1}$, see Section 5) by merging the content of the digital thread ($\mathcal{D}_t$) at the previous time step (and more specifically, the probabilistic knowledge available at the previous time step) with the additional information brought by observations (measurements) performed during the different phases of the lifecycle at stage $t$ ($z_t$). This first data fusion is achieved using Bayesian filtering.

## 2.4. Resources

On the other hand, the resources transition function ($\Psi_t$, see Section 4) should generate the resources available and model the system at the next time step ($R_{t+1}$), by updating the state of these resources at the current time step ($R_t$) in order to account for the new design “deployed” (e.g., update of the geometry of the system of interest, of the dimensions of its structural components, or of the architecture of its structure, based on the design decisions taken at this step).

The global transition function of the digital thread is then defined as a function ($\Phi_t$) combining these two features, in order to update correspondingly the internal state of the digital thread ($\mathcal{D}_t$) each time a new measurement ($z_t$) is observed or a new design decision ($u_t$) is made. This behaviour is summarised in Figure 3, which depicts how these transition functions interact with the digital thread, the data, and the measurements in order to update the current state of the digital thread.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Update process of the digital thread.</figcaption></figure>

## 2.5. Overall Data-Driven Process

In this respect, the overall data-driven process can therefore be viewed as a repeated iteration of a sequence that takes some design decisions, based on the current state of the digital thread, performs observations on the new design resulting from these design decisions, and finally updates the state of the digital thread based on the last design decisions and measurements. The corresponding behaviour is outlined in Figure 4.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Overall data-driven process sequence.</figcaption></figure>

# 3. Overall Methodology

The objective is to develop an approach allowing retrieval of information from different stages of the lifecycle of a system of interest (design, manufacturing, and operation...). This seeks to reduce the uncertainties related to the state and the operational conditions of the system, with the end goal of providing design decision aids, based on the extended knowledge acquired.

To meet this objective, a dedicated framework, presented in Figure 5, is derived from the architectural definition of the digital thread provided by [22]. It is made up of two core elements interacting with a virtual representation of the studied design, further referred to as the *geometry* in Figure 5:

1. **A data filtering application**, which aims at gathering all available information regarding the system. This encompasses observable data ($z_t$) collected during the different phases of the lifecycle of the design of interest and all the resources ($R_t$) exploited in relation with this system. Within this framework, observable data include measurements of the behaviour of the system during its operational lifespan, such as the mechanical strains measured on the structure during its operations ($\epsilon$). It also integrates observations of the fail stress ($\sigma^{\text{lim}}$) of the core material of the structure, which are obtained thanks to coupon tests performed during experiments prior to the design of a new generation of the system of interest. Finally, it includes measurements of the duration of each step required to manufacture the system ($\Delta t_{\text{man}}$). In addition, the filtering application aims at combining these measurements with information derived from models used along the system lifecycle. This includes simulations that predict the mechanical behaviour of the system and the aerodynamic loads resulting from the interactions between the system and its external environment in given conditions. Additionally, it incorporates models detailing the geometry and the structural characteristics of the studied design.
These models are then used to derive the relations that link the measurements defined above with the properties of interest of the system, such as the operational loads actually experienced by its structure or the real failure properties of the materials in use. This link is finally used by a Bayesian filter in order to infer information related to the internal state of the system using the measurements available, and then to improve the knowledge of the uncertain properties of the system (bias reduction and precision improvement) by merging this additional information with the prior knowledge of these uncertain data.

2. **A reliability-based design optimiser**, which relies on the updated probabilistic knowledge of the system (stored in $D_t$, the digital thread itself), in order to compute the optimal decision for the new generation of design to produce. The underlying design decisions are obtained using optimisation methods such as stochastic programming in order to identify the design that allows minimising a given cost function ($r_t$), such as the weight of the structure or the time required to manufacture the final product, while ensuring that the design suggested by the optimiser complies with constraints ($g_t \ge 1$), for example to ensure structural integrity of the optimised design during normal operation or for the sizing cases considered.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Detailed flow diagram of the digital thread implemented for the VTP.</figcaption></figure>

## Architectural Definition

More specifically, the implementation of the two core features defined above implies the use of the following functional blocks interacting together in the global framework depicted in Figure 5 (for each functional block enumerated below, see the bullet displaying the corresponding number on this flow chart):

1. The *virtual representation* of the system, used as the reference source, provides the different resources and models of the framework predicting the behaviour of the system in a given situation (experimental test, response to an external stimulus...). In this framework, the representation comprises, exclusively, the geometry of the system.
2. The real inputs ($y_t$), containing the actual values of all the uncertain properties characterising the system (operational loads, true properties of the material used to manufacture the final product, properties related to the manufacturing process...). These data are used by the resource space in order to model and predict the behaviour of the system in a specific context (such as off-line coupon tests performed during the early design phase of the product or measurements performed during representative periods of its operational life). However, they cannot be accessed by the rest of the framework, which considers the properties related to these data as uncertain inputs ($Y_t$, see further details in item 5).

3. The set of resources ($R_t$) gathering all the elements used in conjunction with the system of interest during at least one phase of its lifecycle. This includes numerical solvers (such as Finite Element Analysis software, allowing prediction of the structural behaviour of the system or the aerodynamic loads that it experiences), empirical models (for example, the models predicting the time required to manufacture a given composite structure), geometric engines (for example, the environments used to generate the meshes used by the structural and aerodynamic solvers or to compute the mass and geometric properties of the given design) and more generally software used to interface with the model or to post-process and display the content of the digital thread.

4. A measurement ($z_t$) feature, which converts the real properties of the system and its ideal behaviour into measurements by adding all the sensing non-idealities, such as measurement noise and potential bias, to the ideal data obtained thanks to the computations and predictions performed by the Resources ($R_t$) defined above. This block might be replaced by a feature importing data from real measurements carried out on the system during key periods of its lifecycle (for instance, measurements of strain gauges on the structure in operation).

5. A Bayesian filter, which relies on these measurements in order to update the probabilistic knowledge of the internal state of the system.

6. The digital thread itself ($D_t$), which merges all the information available to optimise the product knowledge, as well as a link to the resources ($R_t$), in order to be able to interact directly with these entities.

7. The stochastic optimiser, which computes the design decisions to be implemented, based on the current state of the digital thread. This optimal decision is obtained by starting from an initial design point, defined arbitrarily, and converging iteratively to the optimal decision, by computing the design decisions that maximise the expectancy of a given value function ($\psi_t(D_t)$) for different configurations of the digital thread. For that purpose, it relies on a value function that combines some cost considerations ($r_t$), such as the cost and complexity required to manufacture the real system and constraint functions ($g_t$), ensuring that the design suggested by the optimal decision will not lead to a failure of the system during its lifetime. It is worth noting that in Figure 5, $\psi_t^*(D_t)$ denotes the optimal value function, hence the minimal cost due to the transition from the current design of a system ($t$) to its next "generation" ($t+1$).

8. The features computing the design decisions ($u_t$) to be adopted based on the optimal solution computed previously, such as, for example, the ply thickness to manufacture for the different components of the structure, depending on the current value of $D_t$. These decisions are then shared with the reference geometry of the product (and in turn, indirectly with $R_t$), in order to update their respective representation of the current design according to these decisions and therefore to the new product to be generated (modify the composition of the structure of the design of interest, the dimensions or thicknesses of its structural components...).
The Real inputs ($y_t$), the Resource ($R_t$), and the Measurements block ($z_t$) have the purpose of providing virtual measurements of the behaviour of the system under its normal operation or during design experiments in order to evaluate the sensitivity and the potential benefits of implementing a digital thread with a chosen set of measurements, performed during the different phases of the system life. In addition, these virtual measurements are used to compute the optimal decision ($u_t^*$), based on the actual state of the digital thread. Furthermore, the modularity of this framework allows direct transposition of the same approach on real systems by replacing the generation of "virtualisation measurements" ($y_t$, $R_t$ and $z_t$) with measurements performed directly on a real system.

# 4. Parametric Wingbox Generation and Modelling

This section details the different models that are used by the digital thread to refine the knowledge of the design of interest. First, it introduces the environment used to generate and analyse the geometry and the internal structure of the system. Secondly, it describes the computation of the aerodynamic loads resulting from the interactions of the system with its external environment and the prediction of the mechanical behaviour of the structure under these conditions. Finally, it details the model and assumptions that are used to compute the time required by the different steps of the manufacturing process based on the properties computed in the models above.

## 4.1. *Need for a Common Data Structure*

In order to predict the behaviour and performance of the system of interest based on theoretical models, it is first and foremost necessary to create a data structure that contains all the key properties of this system. This structure aims at gathering all the information required by the models of the framework in a single database, to ensure that the models generated automatically rely on the same data, and therefore that the predictions of each model are in accordance with the definition of the system that is used in the other models.

In this respect, a dedicated data structure is defined. It is used to gather and store all the properties required by the different models, such as the external geometry of the system (required to perform mechanical and aerodynamic simulations), its internal structure (needed to accurately predict the mechanical behaviour of the system), or the properties of the external environment in which the product operates (which are used to estimate the external loads induced by this environment), constituting therefore both the geometric properties and the real inputs ($y_t$) of the digital thread (see Figure 5).

The information saved in this data structure is then shared throughout the digital thread as a unique data source characterising the system under investigation to ensure that all models generated by the resources ($R_t$) of the thread rely on the same definition of the system.

## 4.2. *Geometry Generation*

As most of the models of the framework require a detailed knowledge of the main geometric features of the system of interest, the first objective is to implement an application able to generate and analyse automatically the geometry of the system studied by the digital thread. More specifically, this requires the development of a framework able to link high-level parameters of a system to its detailed geometry. The objective of this framework is therefore to generate the reference geometry that will be used by the different models embedded in the digital thread.

To this end, the external geometry of the product is defined using the main shape parameters of the sequence of wing sections composing it. The wing geometry is defined through the different cross-sections of the wing geometry (airfoil profile and chord length)
and the relative positioning of these cross-sections (spanwise distance between these
sections, dihedral angle (ν), twist angle (defined at 50% of the local chord, in this framework)
($a_{tw}$), and sweep angle (defined at 25% of the local chord, for this implementation) ($\varphi_{25}\%$).
The detailed geometry is then generated by performing a linear lofting between each of
these cross-sections. In this respect, Figure 6a shows an example of geometry and structure
generated automatically by the geometric modeller for the VTP investigated in this paper.

For that purpose, the generation tool produces two different instances of this geometry.
The first one is modelled in Abaqus 2020, a Finite Element Analysis (FEA) program with
modelling and analysis capabilities for mechanical components and assemblies. This
software is used to run structural analyses on the geometry in order to determine the
mechanical behaviour of the system depending on the loads experienced by its structure.
In addition, it relies on the geometric modelling capabilities of Abaqus to compute the
detailed features of the system components (mass and dimensions of the structural parts,
surface of intersection between these components...). The second instance is generated
in OpenVSP 3.26.1, where it is combined with the properties of the flow surrounding the
system, in order to compute the aerodynamic loads experienced by the external structure
in these conditions.

The corresponding geometry is then extended with the definition of the internal
structure of the system (number of spars, ribs or stringers, position of these structural
components...), based on the data-set defining the architecture of the internal structure of
the system. More specifically, the setup creates the geometric entities corresponding to each
element of the structural system and adds them to the automatically generated reference
geometry. For that purpose, the spars are defined using multiple points, characterised by
their coordinates along the span of the geometry of interest, and along the local chord at
this spanwise position.

The spars are then instantiated by sequentially extruding shell elements along the
path defined by these different points in order to define a shell filling the interval between
the upper and the lower skins of the wing box. The stringers are then defined by extruding
an assembly of shell elements modelling the cross-section selected by the user (I, L, Z,
or T-shapes) for these components along the paths defined using the same method as
for stringers. Finally, the ribs are added to the structure, depending on their position
and the rotation of these ribs with respect to the axis normal to the plane formed by the
spanwise and chordwise direction of each wing section. In addition, lightening holes can
be integrated in the definition of the ribs. These holes are currently shaped as ellipses with
a major axis aligned with the local chord of the rib and defined by means of the position of
their centre and the dimensions of their (major and minor) axes.

For further details, Figure 6a,b show, namely, the geometry, the internal structure,
and the position where pressure fields are “probed” by the aerodynamic solver (the Cp
slices) of the reference system that serves as a use case, the VTP of an Airbus A320, for the
implementation defined in this work.

Once generated, the geometry and the internal structure of the system, resulting
from the automated framework defined in Section 4.2, are then fed to numerical solvers,
allowing the prediction of their mechanical behaviour and the aerodynamic loads that they
experience, depending on the operational conditions.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Reference model generated for the vertical tail plane use case.</figcaption></figure>

### 4.3. Aerodynamic Loads

A second geometry is generated in OpenVSP [53], a parametric aircraft geometry tool that includes a tool performing aerodynamic analyses (VSPAERO 6.2.0), using intermediate fidelity Computational Fluid Dynamics (CFD) methods, such as the Vortex Lattice Method (VLM) and the Panel Method. This instance is used to compute the aerodynamic pressure field experienced by the external envelope of the system of interest, according to the flow conditions in which the system operates. These computations are performed using a linear potential flow method (assuming inviscid and incompressible aerodynamic flow and neglecting therefore turbulence effects, separation, or complex flow phenomena [54]). Additionally, the compressibility effects are approximated through the Prandtl–Glauert correction [55], making it mainly suitable for low subsonic cases, such as those encountered in the majority of sizing conditions of the VTP, since those conditions are mainly related to flight phases at low Mach numbers, such as take-off and approach. Then, the software uses strip theory with viscous friction to estimate the viscous drag contribution by adding a skin friction component estimate computed using a simple plate drag model for each panel of the model [55]) and including lift-dependent viscous corrections that are derived from an empirical fit of NACA 0012 airfoil data [56] for all lifting surfaces of the model. Finally, a modified version of the Carlson Pressure Correlation [57] limits the local pressure and lift coefficients computed, using a semi-empirical 3D method depending on the local Mach and Reynolds of the flow.

The objective of the simulations performed by this environment is to compute the pressure field resulting from the interactions between the external geometry of the system and the free stream. This field is computed after defining the flow conditions (the angle of attack (α), the sideslip angle (β), and the Mach number (M)), configuring the simulation to run (number of wake iterations...) and defining the cross-sections of the external geometry where the distribution is to be computed. These sections are further referred to as the “Cp Slices” and are currently distributed uniformly along the span of the VTP (see the slices that were defined in Figure 6b). In this regard, Figure 7a depicts the pressure distribution that is computed by OpenVSP for one of these “Cp slices”, while Figure 7b shows the pressure field computed on the whole skin of the VTP.
<figure><img src="image_9.png" /><figcaption>(a) Cp distribution of the 19th "Cp slice" (z=2.6m).</figcaption></figure>

<caption>(b) 3D pressure distribution.</caption>

<caption>**Figure 7.** Pressure distributions computed by OpenVSP for the reference VTP.</caption>

It is worth noting that for higher Mach numbers, this tool also includes a second-order Karman-Tsien correction feature [56], which might prove more relevant than the Prandtl-Glauert correction for computations in transonic conditions. This may be necessary to use this environment for wing analysis, which would require, among others, the performance of analyses as cruising conditions. Those modelling assumptions might also produce slightly higher $C_{LMax}$ values than the ones of the actual VTP, which would lead to overestimation of the maximal aerodynamic loads at high angles of attack and might therefore result in slightly conservative optimisation results regarding the designed structure. On the other hand, this CFD solver is not able to model post-stall lift coefficients either, which could lead to less reliable results from the controllability point of view. This, however, is currently out of the scope of this study.

Regarding the fidelity level of the software selected for this study, its main objective was to provide a first proof of concept dedicated to preliminary design phases, where computational performance generally takes precedence over the accuracy of the results computed. However, this environment could benefit from the integration of higher-fidelity CFD methods, such as a viscous solver (e.g., using steady RANS) or a combined version of high- and low-fidelity CFD tools in order to validate and correct on-line the predictions of this first CFD tool, which would in turn increase the accuracy of the simulation results computed, while maintaining a moderate impact on computational performance.

Furthermore, the computation of the pressure field obtained from the aerodynamic model is performed without displacement transfer between the geometry and the aerodynamic solver (i.e., assuming a rigid link transfer from the geometry to the aerodynamic solver). Therefore, it does not take into account the impact of elastic phenomena on the pressure field computed. This assumption should not, however, significantly impact the computations performed by the current framework, which is intended to investigate the mechanical behaviour of stiff structures, such as a VTP made of carbon fibre reinforced polymers, under static loading conditions.

### 4.4. *Structural Simulations*

In this first model, the structure (Figure 6a) is assumed to be made of carbon fibre reinforced polymer (CFRP) materials. More specifically, the mechanical system is modelled using transverse isotropy and linear elasticity assumptions, using the MTM45-1/12K AS4 [58] material as the constitutive material of the structure generated by this environment (see Appendix A.1 for further details about the mechanical properties considered to model
this material). In addition, the structure is meshed using general purpose quadrilateral shell elements with finite membrane strains and reduced integration (S4R elements).

In order to reduce the complexity and the number of degrees of freedom generated by the structural modeller, the composite layups generated for the different shells of the systems are made of a symmetric stacking of four plies of variable thicknesses ($e_{0°}$, $e_{90°}$ and $e_{±45°}$), with the following layup: $[45°, -45°, 90°, 0°]_{sym}$ (see the corresponding layout in Figure 8), leading therefore to a zero membrane-bending coupling matrix ($B$) and thus to a decoupling of the bending and the membrane response of the corresponding section.

This layup is defined in such a way that it optimises the membrane stiffness and/or the critical buckling load of the laminate (for a given set of plies thicknesses). More specifically, the central plies of the layup (0° and 90°) aim at providing the membrane stiffness of the layup, while the plies with a fibre angle of 45° and -45° are located on the external plies of the layup (as far as possible from its symmetry plane), in order to maximise their impact on the buckling resistance of the layup. In addition, the thicknesses of the 45° and -45° plies are set equally, forcing the stacking sequence to be balanced, which results in a decoupling of the shear-stretch response of the layup ($A_{16} = A_{26} = 0$) and ensures (almost) symmetrical buckling behaviour of the shells made of this layup.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Layup composition.</figcaption></figure>

A different thickness is then set for each component of the structure. By default, one thickness is set independently for each of these sets: the upper skin, the lower skin, all the spars, the root rib, and the tip rib. Finally, the thickness of each of the intermediate ribs is linearly interpolated, based on the spanwise position of these ribs along the wing.

The boundary conditions of the model are then set by clamping the root rib of the wing geometry generated, while two independent pressure fields are used to apply load boundary conditions to the upper and the lower skins of the structure in order to model the aerodynamic loads that are experienced by the wing model. These pressure fields are computed depending on the flow conditions in which the vertical tail plane operates. The way these computations are performed is detailed further in Section 4.3. Finally, the whole structure is meshed using quad elements.

Static deformations (Figure 9a) and buckling (Figure 9b) simulations are then run on these models, in order to predict the corresponding behaviour of the structure for a given set of aerodynamic loads, derived from the current flight conditions of the system. The local stress ($\sigma$) fields resulting from the static analysis are then used to assess whether the structure might fail under these static conditions. This verification is achieved by means of the Tsai–Wu [59] failure criteria (see further details in Section 6.1.1), using the failure strengths defined in Appendix A.1. In addition, the eigenvalues of the first buckling modes of the structure ($\lambda$) are imported in order to check that its loading does not lead to buckling phenomena. These failure criteria are finally used to calculate the Reserve Factors (RF) defining the safety margin between the current simulation and the loads that would lead to a system failure (i.e., a failure occurs if $RF \le 1$). These factors are then aggregated in a
single vector (RF), containing therefore two different kinds of Reserve Factors, computed as follows:

$$RF_{\text{buckling}}^{i} = \lambda_{i} \qquad (5)$$

$$RF_{\text{static}}^{j} = \frac{1}{f(\sigma)} \qquad (6)$$

where $RF_{\text{buckling}}^{i}$ and $\lambda_{i}$ are the Reserve Factor and the eigenvalue associated to the ith buckling mode of the structure, $RF_{\text{static}}^{j}$ is the factor resulting from the static loading of the structure, evaluated at the jth integration point of the model, and $f(\sigma)$ is the value of the stress-based failure criterion computed using the Tsai-Wu failure criterion at this integration point. The synthesis of the interactions between the mechanical simulations and the post-processing operations described above is shown in Figure 10.

<figure><img src="image_11.png" /><figcaption>(a) Static deformation of the structure (b) Buckling shape (mode 4) predicted. Figure 9. Structural behaviour simulated by Abaqus for the reference VTP.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10. Flow diagram of the structural computations.</figcaption></figure>
## 4.5. *Manufacturing Properties*

One of the objectives of this digital thread is to predict the duration of the differ-
ent steps required to manufacture the system of interest, and therefore the underlying
manufacturing cost. This cost is currently assumed to depend only on the mass of mate-
rial and on the time required to manufacture the corresponding part. More specifically,
it is computed in the following way: cost = mass × raw cost + manufacturing dura-
tion (Σ ACCEM) × hourly rate, considering a raw cost of 50 €/kg and an hourly rate of
20 €/h. These durations are predicted thanks to the Advanced Composite Cost Estimating
Manual [51,60] (ACCEM), which identifies and summarises empirical models predicting
the duration of every step that can be involved in the manufacturing process of a composite
structure, depending on the geometric properties of this composite element (dimensions,
number of flanges, holes or curved bends...).

To this end, this model defines and pools all the manufacturing steps to be performed
according to the type of operation associated to these steps, such as layup, core preparation,
curing, consolidation, and finishing operations. It begins with the layup of the material to
form the components of the final assembly (ribs, stringers, skin...), which can be performed
manually, using an automatic tape laying machine, or via hand-assisted methods. Then, the
model defines the steps required to prepare the core material, in case of parts containing a
honeycomb core (these steps are currently not included in the implementation described
in this paper). The next manufacturing steps computed by the model are related to the
curing process of the composite material, which can be achieved using vacuum bagging in
an autoclave (curing over 1.034 bar) or an oven (curing at or below 1.034 bar), or by the use
of heating elements in thermal expansion moulds. Afterwards, the model considers the
joining and splicing operations required to build the assembly (for example, for applying
adhesive or handling the part). Finally, finishing operations can be executed on the cured
part. There are two main types of these operations. The first category is related to the
drilling of holes, required to match the specifications of the assemblies, in order to be able
to integrate some external systems (actuators, flight control surfaces...) into the structure
manufactured. The second type of finishing operations is used to trim a part to bring it to
its final dimension specifications.

The duration of all the manufacturing steps defined by the ACCEM is then computed
based on the geometric properties of the final assembly to build, along with parameters
related to the manufacturing process in use. To this end, the geometric engine of Abaqus
is used to calculate the main geometric properties of each component of the structure
investigated (total external perimeter, surface, volume... of the part). An example of the
corresponding properties of the reference VTP is shown in Figure 11a. Additionally, it
computes the geometric interfaces between these different parts (i.e., the areas of inter-
sections between each sub-component of the structure). This information is then sent to
this ACCEM implementation, and used along with user parameters defining the materials
and tools used to manufacture the structure (see description below) in order to predict the
time required to perform each step of the manufacturing process, the duration needed to
manufacture each sub-component of the structure, and finally the time to generate and
assemble the whole structure.

For that purpose, the user parameters defining the specific features of the manufactur-
ing process in use, as implemented in the digital thread, are the following:

* Ply deposition method: either manual, hand-assist, or automatic;

* Size of the tape used to manufacture the different layups;

* Speed at which the automatic ply deposition is performed (if applicable);

* Type of bags used for the curing processes: reusable or disposable; and
* Curing mean: autoclave or oven (using vacuum bag curing), heating elements, or oven (if thermal expansion moulding is used).

The model predicting the duration of each manufacturing step of the design is then computed automatically, assuming that the generation of each sub-component of the final system involves a layup generation [51], pp. 28–30 (the material is considered to be made of a single layup, without any core material embedded in the middle of that layup), an autoclave curing [51], pp. 47–78, and finally the trimming [51], pp. 28–30 of the part manufactured.

Each assembly operation performed to generate the final structure is then assumed to imply a splicing [51], p. 50 task, and then thermal expansion moulding [51], p. 63 of the junction area of the parts to be assembled. These operations are then used for each step of the final assembly, which is assumed to start with the conglomeration of all spar flanges (on both skin halves), the addition of every rib and spar web, and finally the joining of both skin halves.

The models and parameters resulting from the predictions of these different manufacturing steps are then merged in a single model ($H_t^m$) and input vector ($y_t^m$) in order to obtain the global manufacturing model predicting the duration of every step of the manufacturing process of the structure modelled for a given set of parameters. To illustrate these calculations, Figure 11a shows the geometric properties that are computed on the reference VTP. These dimensions are used by the ACCEM to compute the linearised manufacturing parameters of the model ($y_t^m$) and the transfer function ($H_t^m$), linking these parameters with the duration of the different steps required to manufacture the whole system. These data are then used to assess the duration of every step of the manufacturing process, before combining all the steps related to a given component in order to compute the time required to manufacture each component of the reference structure (Figure 11b).

<figure><img src="image_13.png" /><figcaption>Figure 11. Input parameters and resulting outputs of the ACCEM for the reference VTP.</figcaption></figure>

## 4.6. Global Framework

In a nutshell, the overall architecture of the multi-physical modelling and simulation environment comprises two geometric engines, generating automatically the geometry and structure of the system of interest, depending on its major characteristics. The resulting geometries are then fed simultaneously to the aerodynamic solver, in order to compute the pressure loads exerted by the external environment of the system, and to the mechanical solver, so as to compute the mechanical behaviour of the system in these aerodynamic conditions. The corresponding framework is summarised in Figure 12, which delineates the
interactions between the structural and aerodynamic solvers, as well as with the geometric and manufacturing models. Furthermore, it depicts the outputs resulting from all those analyses, along with the interfaces and file formats used to exchange seamlessly different types of information (such as geometry and loads) among these models.

In this respect, the models used by this environment are defined and analysed using primary parameters, related to the geometry and structure definition, as well as to the free stream conditions of the environment where the system operates and to the properties of the core material used to manufacture the system. These parameters are summarised in Table 1.

<figure><img src="image_14.png" /><figcaption>Figure 12. Global framework of the parametric modelling tool.</figcaption></figure>

<table><caption>Table 1. Summary of the parameters available for the geometric and structural generation.</caption><thead><tr><th>Geometric Parameters</th><th>Structural Components</th><th>Aerodynamic Configuration</th></tr></thead><tbody><tr><td><ul><li>Airfoil profile</li><li>Chord length</li><li>Span</li><li>Sweep angle (φ<sub>25%</sub>)</li><li>Twist angle (α<sub>tw</sub>)</li><li>Dihedral (ν)</li></ul></td><td><ul><li>Ribs<ul><li>Position</li><li>Holes</li></ul></li><li>Spars<ul><li>Position</li></ul></li><li>Stringers<ul><li>Position</li><li>Cross-section</li></ul></li></ul></td><td><ul><li>Angle of attack (α)</li><li>Sideslip angle (β)</li><li>Mach number (M)</li><li>Air density (ρ)</li><li>Free stream velocity (V<sub>∞</sub>)</li></ul><p><strong>Layup definitions</strong></p><ul><li>Ply thicknesses</li></ul></td></tr></tbody></table>
## 5. Uncertainty Reduction Using Bayesian Filtering

This section illustrates how the simulation results from the above models are used to infer the features of the measurements that would be performed on the real system during the different phases of its lifecycle (experiments performed during the design phase, structural deformation measured throughout the operational life...). Finally, it presents the way these measurements are processed in order to update the initial knowledge of the system state.

The main objective of this filter is to provide a feature able to update and improve the knowledge of the system state, based on data measured during the lifecycle of the system. For that purpose, the filter developed takes different kinds of measurements into account in order to update the current knowledge of the system state:

1. Direct measurements, which are directly representing the internal state of the system (e.g., fail-stresses of the constitutive material of the structure, measured during coupon tests) and can therefore be used directly in the Kalman filter with model-free methods.

2. Indirect measurements, using analytic models to link these measurements with the internal state of the system (as for example the ACCEM, which relies on empirical relations in order to link the duration of the different steps of the manufacturing process of a given structure to the manufacturing properties of this structure).

3. Indirect measurements, using numerical models trained by results from digital simulations (e.g., loads withstand by the structure and derived from strain measurements).

In order to achieve this objective, this application relies on Bayesian filtering using a variant of the discrete-time Kalman filter where the measurement step is performed before the prediction step, in order to compute a minimum variance unbiased estimator of the state variables of interest ($\hat{y}_t$). This section details the implementation of this filter, along with the system modelling assumptions used for that purpose. The notations and modelling assumptions of this section are partially derived from the ones detailed in [22].

### 5.1. Dynamic System Model

To implement this filter, the first step is to define the linear models characterising the state transition of the system ($y_t$, see Equation (7)), as well as the relation linking this state with the observations performed on the system ($z_t$, see Equation (8)). To this end, the state transition of the system is modelled using the following equation:

$$Y_{t+1} = F_t(R_t, u_t)Y_t + B_t(R_t)u_t + \zeta_t \quad (7)$$

where $F_t$ is the state transition function, $B_t$ is the input transition matrix, and $\zeta_t$ represents the process noise, which is assumed to be a white Gaussian noise of covariance $\Theta_t$. Therefore, $\zeta_t \sim N(0, \Theta_t)$.

In addition to that, the observation model is assumed to link the system state ($y_t$) to the measurements performed ($z_t$) through the following linear system, using the observation model functions ($H_t$ and $h_t$) that can be obtained from analytic models (such as the ACCEM) or by computing the numerical model linking these properties using system identification or surrogate modelling techniques. Furthermore, this model considers the impact of the observation noise ($\epsilon_t$), which is assumed to be a white Gaussian noise of covariance $\Lambda_t(\epsilon_t \sim N(0, \Lambda_t))$. The global observation model is therefore expressed as:

$$Z_t = H_t(R_t, u_t)Y_t + h_t(R_t, u_t) + \epsilon_t \quad (8)$$
At the present time, the assumption is made that the information does not vary between the different design stages ($F_t = I_{N_y}$) except for some errors coming from $\zeta_t$, and that the control inputs does not influence the evolution of the state of the system ($B_t = 0_{N_y}$).

## 5.2. Bayesian Filter [22], pp. 75–76

The objective of this filter is to use Bayes' theorem in order to define how the probabilistic knowledge of the system state should evolve, given observations performed (directly or indirectly) on some variables characterising a system. This is achieved by combining the prior belief of the system state before any observation is made ($p(y_t)$), the likelihood defining the probability of making a specific observation for a given state of the system ($p(z_t|y_t)$) derived from the observation models, and the knowledge of the stochastic properties of the sensors ($p(z_t)$, the evidence) in order to compute the posterior distribution of the variables of interest ($p(y_t|z_t)$), using the following relation:

$$p(y_t|z_t) = \frac{p(z_t|y_t) \cdot p(y_t)}{p(z_t)} \qquad (9)$$

For the Kalman filter, this is achieved using a sequence of update and prediction steps, where the filter merges information about the dynamics of the system and the observations made on the system in order to predict or update the estimates of its state variables.

Each time a measurement is performed, the updated uncertain input resulting from data assimilation is computed using Kalman theory, $\bar{Y}_t \sim N(\bar{\mu}_t, \bar{\Sigma}_t)$, starting from the prior distribution of the uncertain inputs, $Y_t \sim N(\mu_t, \Sigma_t)$, where the distribution of these two variables is assumed to be Gaussian. For the system defined in Section 5.1, the updated mean and covariance properties of these data are computed as:

$$\bar{\mu}_t = \mu_t + K_{Kalman} \cdot e_t \\ \bar{\Sigma}_t = (I - K_{Kalman} \cdot H_t) \cdot \Sigma_t \qquad (10)$$

where $K_{Kalman}$ is the Kalman gain, the weight given to the measurements and current state estimates, expressed as:

$$K_{Kalman} = \Sigma_t H_t^T (H_t \Sigma_t H_t^T + \Lambda_t)^{-1} \qquad (11)$$

$e_t$ is the pre-fit residual, i.e., the difference between the current predictions of the sensor model and the real measurements performed on the system.

$$e_t = (z_t - H_t \mu_t - h_t) \qquad (12)$$

The update between different design stages is then performed using the following relation:

$$\mu_{t+1} = F_{t+1}\bar{\mu}_t + B_t u_t \\ \Sigma_{t+1} = F_{t+1}\bar{\Sigma}_t F_{t+1}^T + \Theta_{t+1} \qquad (13)$$

This sequential behaviour of the Kalman filter can be modelled using Figure 13.
<figure><img src="image_15.png" /><figcaption>Figure 13. Block diagram of the Kalman filter.</figcaption></figure>

## 5.3. Initial Prior Distributions

The internal properties of the system (further referred to as the inputs, $y_t$) are assumed to be made of three different types of data: the loads applied on the structured (outlined using the *l* superscript), the properties related to the manufacturing process (indicated by the *m* superscript), and the properties of the core material of the system (labelled using the *a* superscript). The same notation is used for the other variables, properties, and model functions related to these type of data (e.g., $y_t^l$, $z_t^l$, and $H_t^l$, $h_t^l$ denote, namely, the inputs, measurements, and observation model related to the load properties of the system (internal loads, strain measurements and functions linking the strain measurements and the structural loads), $y_t^m$, $z_t^m$, and $H_t^m$, $h_t^m$ are used for the variables related to the manufacturing properties of the system, and $y_t^a$, $z_t^a$ and $H_t^a$, $h_t^a$ represent the variables modelling the material properties of the system). Moreover, these three kinds of data ($y_t^l$, $y_t^a$, and $y_t^m$) are assumed to be independent of each other.

In addition, the initial prior distribution of these properties $p(y_0 | I_0, u_0)$ is set to be Gaussian, i.e., $Y_0 \sim N(\mu_0, \Sigma_0)$ where $\mu_0 \in \mathbb{R}^{N_y}$ and $\Sigma_0 \in \mathbb{R}^{N_y \times N_y}$. For the material properties $(\mu_0^a, \Sigma_0^a)$ and the manufacturing parameters $(\mu_0^m, \Sigma_0^m)$, these initial properties are derived from the ones suggested by [22]. Moreover, the combination of the mean ($y_0$) and covariance ($\Sigma_0$) of these properties can be expressed as:

$$ \mu_0 = \begin{bmatrix} \mu_0^l \\ \mu_0^a \\ \mu_0^m \end{bmatrix} \qquad \Sigma_0 = \begin{bmatrix} \Sigma_0^l & 0 & 0 \\ 0 & \Sigma_0^a & 0 \\ 0 & 0 & \Sigma_0^m \end{bmatrix} $$

### 5.3.1. Material Properties

The prior distribution of the material strength properties at stage 0 is defined such that the mean of the initial distribution equals the 90th percentile of the actual material strength properties. The other material properties (elastic moduli and Poisson's ratios) are assumed to be certain variables perfectly known (unbiased). This assumption aims at approaching the B-value (value of the distribution such that at least 90% of the population equals or exceeds the material properties selected, with a 95% confidence interval) of these properties, in order to comply with the structural certification requirements from CS 25.613. The corresponding data are obtained from the “modified” specimen (to account for an additional variation in the material properties when it is produced in larger amounts over a lengthy period of time (to manufacture the real system), compared to the variability that could be observed by performing coupon tests on the initial specimens (not perfect independence of
the ingredients used to manufacture the composite, specimens manufactured within a shot period of time, which are generally not representative of the actual production material variability). The B-value properties are enumerated in [58], p. 36 for Room Temperature Dry (RTD) environmental conditions.

Therefore, the real properties of the core material selected in Section 4.4 are assumed to be ideal (constant) values of the failure stress of the core material of the structure under a tensile ($\sigma^t$) or compressive ($\sigma^c$) deformation, and are set to be:

$$y_t^a = \begin{bmatrix} \sigma_t^c \\ \sigma_t^t \\ \sigma_t^c \\ \tau_{lt}^{\lim} \end{bmatrix} = \begin{bmatrix} 1867 \\ -1398 \\ 48 \\ -185 \\ 65 \end{bmatrix} \quad [MPa]$$

While the initial mean ($\mu_0^a$) and covariance ($\Sigma_0^a$) of the prior distribution of these proper-

$$\text{ties is set as: } \mu_0^a = \begin{bmatrix} 1655 \\ -1234 \\ 31 \\ -167 \\ 58 \end{bmatrix} \quad [\text{MPa}] \quad \Sigma_0^a = \begin{bmatrix} 3.3 \times 10^4 & 0 & 0 & 0 & 0 \\ 0 & 3.3 \times 10^4 & 0 & 0 & 0 \\ 0 & 0 & 17 & 0 & 0 \\ 0 & 0 & 0 & 500 & 0 \\ 0 & 0 & 0 & 0 & 17 \end{bmatrix} \quad [\text{MPa}]^2$$

### 5.3.2. Manufacturing Parameters

The mean of the initial distribution assumed for the manufacturing parameters is computed assuming an initial overestimate of 20% of the material properties with respect to the actual manufacturing parameters of the real system ($\mu_0^a = 1.2 \cdot y_t^a$). Moreover, the covariance matrix of the manufacturing parameters is defined assuming independence between each manufacturing parameter, while the standard deviation of each parameter is assumed to equal 10% of the corresponding parameter. This can be re-written as:

$$(\Sigma_0^m)_{ij} = \begin{cases} [0.1 \cdot (y_0^m)_i]^2 & \text{if } i = j \\ 0 & \text{if } i \neq j \end{cases} \quad (14)$$

where $(y_0^m)_i$ is the ith manufacturing property of the structure to manufacture, calculated using the ACCEM, based on the reference geometry (see Section 4.5), and therefore

$$Y_0^a \sim \mathcal{N}(1.2 \cdot y_t^a, \Sigma_0^a) \quad (15)$$

### 5.3.3. Loads

The prior knowledge of these loads is computed by assuming an initial bias of 0.5 with respect to the actual pressure field experienced by the external skin of the structure ($\mu_0^l = 1.5 \cdot y_t^l$), in order to investigate the impact of the 1.5 safety factor between limit and ultimate loads (CS 25.303) on the convergence properties of the digital thread.

In addition, the loads applied to the external structure are assumed to be spatially correlated. Their initial covariance is therefore computed as:

$$\Sigma_0^l = k_{SE}(x, y) = \sigma^2 \cdot \exp \left[ -\frac{1}{2} \cdot \left( \frac{\|x - y\|_G}{\ell} \right)^2 \right] \quad (16)$$

The hyper-parameters of this kernel are $\ell$ and $\sigma$, where $\ell$ represents the length scale of the kernel considered, which defines how close two points have to be to influence each other
significantly, and $\sigma$ is the scale factor of the kernel, which defines the standard deviation of the load at a given position. Moreover, $\|\cdot\|_G$ denotes the geodesic distance between two points located on the surface of the reference geometry (i.e., the length of the shortest path connecting these two points along the external surface of the system).

For the use case presented in this paper, the following parameters are considered: $\ell = 0.1$ m and $\sigma^2 = 10^{-2}$. Therefore,

$$Y_0^l \sim \mathcal{N}(1.5 \cdot y_t^l, \Sigma_0^l) \qquad (17)$$

## 5.4. Observation Models

The measurements are assumed to be linear functions of the internal state of the system. In addition, these measurements are considered to be subject to an additive white Gaussian noise. The link between the internal state of the system and the measurements can therefore be re-expressed using the following relation:

$$Z_t = H_t(R_t, u_t)Y_t + h_t(R_t, u_t) + \epsilon_t \qquad (18)$$

where $H_t$ and $h_t$ are the observation model functions, derived from analytic relations or computed by identifying the linear model that links a set of inputs ($y_t$) with the corresponding measurements ($z_t$), and $\epsilon_t$ represents the additive measurement noise. It is assumed to be a white Gaussian noise with covariance $\Lambda_t^l$ ($\epsilon_t^l \sim \mathcal{N}(0, \Lambda_t^l)$). This measurement covariance is defined based on the assumptions set in [22], depending on the kind of measurements performed (see Table 2).

<table><caption>**Table 2.** Measurement variance for the example design problem.</caption><thead><tr><th>Component</th><th>Variance (σ<sup>2</sup>)</th></tr></thead><tbody><tr><td>Z<sub>t</sub><sup>l</sup></td><td>100 (μ-strain)<sup>2</sup></td></tr><tr><td>Z<sub>t</sub><sup>a</sup></td><td>[2, 10<sup>3</sup>] (MPa)<sup>2</sup></td></tr><tr><td>Z<sub>t</sub><sup>m</sup></td><td>4 (s)<sup>2</sup></td></tr></tbody></table>

### 5.4.1. Material Properties

The observation model of the material properties is defined as:

$$Z_t^a = H_t^a(R_t, u_t)Y_t^a + \epsilon_t^a \qquad (19)$$

where the observation model function ($H_t^a$) is the identity matrix ($H_t^a = I_5$), since material properties are direct measurements of the corresponding state of interest of the system ($y_t^a = z_t^a$) and $\epsilon_t^a \sim \mathcal{N}(0, \Lambda_t^a)$ is an additive noise term with covariance $\Lambda_t^a$. This covariance matrix is defined as:

$$\Lambda_t^a = \begin{bmatrix} 1000 & 0 & 0 & 0 & 0 \\ 0 & 1000 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 30 & 0 \\ 0 & 0 & 0 & 0 & 2 \end{bmatrix} [\text{MPa}]^2 \qquad (20)$$

### 5.4.2. Manufacturing Timestamps

This observation model is linking the geometry to the manufacturing duration required to create the component of assembly, measured during the manufacturing of the part using the following relation:

$$Z_t^m = H_t^m(R_t, u_t)Y_t^m + \epsilon_t^m \qquad (21)$$
where the observation model function ($H_t^m$) is the “transfer function” computed by linearising the ACCEM [60] based on the current properties of the geometry modelled (see Section 4.5) and $\epsilon_t^m \sim N(0, \Lambda_t^m)$ is an additive noise term with covariance $\Lambda_t^m$. An example of linearised transfer function computed by the ACCEM for a reduced set of the reference model (the first Rib of the structure modelled) is shown in Figure 14, while this covariance matrix is defined as

$$ \Lambda_t^m = (2s)^2 = \left(\frac{2}{3600}\right)^2 \cdot I_{N_{\text{timesteps}}} \quad (22) $$

<figure><img src="image_16.png" /><figcaption>Figure 14. Transfer function of the sensor model predicting the manufacturing properties of a Rib of the VTP depending on its geometry.</figcaption></figure>

### 5.4.3. Strain Sensors

This last observation model aims at linking the geometry, the aerostructural model, and the environment to the structural strains measured on the system during its operational life. This is achieved using the following relation:

$$ Z_t^l = H_t^l(R_t, u_t)Y_t^l + h_t^l(R_t, u_t) + \epsilon_t^l \quad (23) $$

where $H_t^l$ is the observation model function that links an external loading applied to the structure with the in-plane strains arising from this loading, $h_t^l$ is the observation model function linking the displacement boundary conditions imposed to the structure with the in-plane strains due to these conditions, and $\epsilon_t^l \sim N(0, \Lambda_t^l)$ is an additive noise term with covariance $\Lambda_t^l$ used to model the loads measurement noise. Similarly to [22], the noise covariance related to term is currently defined as:

$$ \Lambda_t^l = (10\mu\epsilon)^2 = (10^{-5})^2 \cdot I_{N_{\text{loads}}} \quad (24) $$

Furthermore, the loads observations are defined such that

$$ \epsilon_t = H_t^l \cdot C_p + h_t^l \quad (25) $$

Since $h_t$ represents the impact of the displacement boundary conditions on the in-plane strains measured by the different sensors, it can be recovered from simulations of the system constrained by the displacement boundary conditions, without any external loading.

The respective uncertain inputs ($y_t^l$) are the pressure coefficients ($C_p$) resulting from the aerodynamic loads applied to the structure, and the measurements ($z_t^l$) characterising these inputs are the in-plane strains ($\epsilon_t$) measured on the deformed structure. These data are ex-
$$ \text{pressed as: } \quad y_t^l = C_p = \begin{bmatrix} Cp^{\text{Slice1}} \\ Cp^{\text{Slice2}} \\ \vdots \\ Cp^{\text{SliceN}} \end{bmatrix} \in \mathbb{R}^{3N_{\text{Slices}} \cdot N_{C_p} \times 1} \qquad z_t^l = \varepsilon_t = \begin{bmatrix} \varepsilon_l \\ \varepsilon_t \\ \gamma_{lt} \end{bmatrix} \in \mathbb{R}^{3N_s \times 1} $$

The pressure distribution (Cp) to apply to the model is obtained from aerodynamic simulations, computing the pressure field applied to the structure in the specified free stream conditions. On the other hand, the ideal (noiseless) strain measurements are predicted using Finite Element Method (FEM) to simulate the structural deformation of the geometry of interest due to the pressure loads computed above.

The load sensor model function ($H_t^l$) is then approximated as the least-square solution of Equation (26), using the pressure (Cp) and in-plane strain ($\varepsilon_t$) fields computed above.

$$ (\varepsilon_t - h_t^l)^T = Cp^T \cdot (H_t^l)^T \qquad (26) $$

where $\varepsilon_l$, $\varepsilon_t$, and $\gamma_{lt}$ are the three in-plane strain component measurements expressed in the material frame (the frame aligned with the three principal material directions, including the fibre of the material, i.e., the (l, t, z) frame). These data are measured at the location of all the strain sensors, which are distributed linearly along the different components of the structure (or along the planar projection of the surface of these components, for the skins of the structure). An example of distribution of these sensors on the skin, on a spar, and on a rib of the VTP is shown in Figure 15a-c. These strain gauges are assumed to be located on the outer surface of the composite layups. The “virtual measurements” at these locations are then computed using a distance weighting algorithm to interpolate the strain field at these positions, based on the values of the strain fields at the integration points of the mesh.

<figure><img src="image_17.png" /><figcaption>Figure 15. Distribution of the sensors (“virtual” strain gauges) on the main components of the VTP.</figcaption></figure>

### 5.4.4. Overall Model

The global sensor model is then obtained by condensing the expressions calculated in Sections 5.4.1–5.4.3 into one single measurement model, as follows:

$$ Z_t = H_t(R_t, u_t)Y_t + h_t(R_t, u_t) + \varepsilon_t \qquad (27) $$

where

$$ H_t = \begin{bmatrix} H_t^l & 0 & 0 \\ 0 & H_t^a & 0 \\ 0 & 0 & H_t^m \end{bmatrix}, \quad h_t = \begin{bmatrix} h_t^l \\ 0 \\ 0 \end{bmatrix}, \quad \varepsilon_t = \begin{bmatrix} \varepsilon_t^l \\ \varepsilon_t^a \\ \varepsilon_t^m \end{bmatrix}, \quad \Lambda_t = \begin{bmatrix} \Lambda_t^l & 0 & 0 \\ 0 & \Lambda_t^a & 0 \\ 0 & 0 & \Lambda_t^m \end{bmatrix} \qquad (28) $$
### 5.5. Filtering Convergence

This section highlights the improvements in the knowledge of the state of the system enabled by the filtering feature detailed in Section 5. More specifically, Figure 16 depicts the evolution of the prior knowledge of the properties of the core material of the system at the different stages after performing two sequential measurements of these properties using coupon tests.

<figure><img src="image_18.png" /><figcaption>Figure 16. Convergence of the material properties along the different steps (t).</figcaption></figure>

On the other hand, the convergence of the manufacturing time predicted by the model after filtering of additional measurements performed during the manufacturing stage of the new design can be analysed by comparing Figure 17a, which depicts the initial prior distribution of the manufacturing timestamps assumed for the studied vertical tail plane, while Figure 17b shows the posterior knowledge resulting from the filtering of the a priori distribution with time measurements performed during the manufacturing process.

Equivalently, Figure 18a,b depict, respectively, the initial and updated pressure distribution for the 19th slice, and Figure 19a,b depict, respectively, the initial and updated covariance known for the uncertain loads.

<figure><img src="image_19.png" /><figcaption>Figure 17. Statistics ($\mu \pm 2\sigma$ bounds) of the manufacturing timestamps for the different components.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Statistics ($\mu \pm 2\sigma$ bounds) of the pressure distribution predicted for the 19th slice (z = 2.6 m).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Covariance known for the uncertain loads (pressure distribution).</figcaption></figure>

# 6. Optimisation of Composite Wing Structures

In this section, the models developed in Section 4 are used to compute the optimal wing structure that could be operated safely under these most critical conditions. The best design choice is made based on the level of uncertainty related to the knowledge of the real system. This section details the sequences of the original approach developed for the optimisation of the VTP structure thanks to the combined use of two powerful and complementary tools in such a way that the respective weaknesses of each tool are avoided in the global solution.

In this respect, a parametric optimisation environment has been developed for the VTP structure (see Figure 20). It relies on the parametric modelling feature described in Section 4 in order to generate the FEM model (that includes the geometry and the structure) of the system of interest, depending on its high-level structural parameters (see Table 1). For that purpose, the framework implemented relies on two different loops

1. **Outer loop:** This structural optimisation process (see Figure 20) is used to generate the models required to perform the function and constraint evaluations used by the thickness optimiser. A design of experiment method is used to manage this feature in order to feed the thickness optimisation loop with all the geometric and structural configurations to be investigated (e.g., by varying the number of ribs, the spanwise distribution of these ribs, the cross-section of the stringers, or the number of layups of independent thicknesses on the VTP skins).

2. **Inner loop:** This component of the framework performs the thickness optimisation (see Figure 20) that computes the optimal thicknesses of all the layups of the wing structure generated, ensuring that it does not fail under the aerodynamic loading resulting from knowledge of the uncertain inputs of the system. To achieve this, each
function evaluation of the optimiser relies on mechanical simulations using FEA, followed by a post-processing constraints check.

<figure><img src="image_22.png" /><figcaption>Figure 20. Parametric optimisation environment for the VTP structure.</figcaption></figure>

This optimisation framework is implemented using iSIGHT 2020 [61], a Low-code/No-code Integrated Modelling Environment that helps to streamline the engineering processes (digital thread) in a single engineering platform by allowing the automation of processes such as the simulation and optimisation of designs. This tool allows a seamless integration of interfaced applications, such as FEA, CAD, PLM, or scripting tools, as well as interfaces with programming languages. Moreover, it natively embeds tools, allowing assistance with decisions (design exploration and optimisation tools), uncertainty quantification (Monte Carlo, reliability techniques...), approximation features, and design of experiments (DoE) methods. This software supports a broad variety of tools and models, allowing the performance of multidisciplinary analyses and optimisation, making it therefore suitable for a seamless and scalable integration of the processes involved in the implementation of a digital thread. It also implements a broad variety of optimisation techniques (gradient-based, genetic, exploratory, direct search...).

In this respect, the structural optimisation performed for this study relies on multistage optimisation, where the first optimisation step relies on non-linear quadratic programming using the NLPQLP [62] algorithm, and the second step uses the Modified Method of Feasible directions (MMFD, [63]).

## 6.1. Constraints

### 6.1.1. Stress Constraint Under Static Deformation: Tsai-Wu Failure Criterion

The Tsai-Wu criterion is a quadratic failure criterion used for anisotropic composite materials with different compressive and tensile strengths. To ensure that the material remains within the allowable limits and therefore does not fail under the given loading conditions, it requires that:

$$ \left(\frac{1}{\sigma_l^t} + \frac{1}{\sigma_l^c}\right) \sigma_l + \left(\frac{1}{\sigma_t^t} + \frac{1}{\sigma_t^c}\right) \sigma_t - \frac{\sigma_l^2}{\sigma_l^t \sigma_l^c} - \frac{\sigma_t^2}{\sigma_t^t \sigma_t^c} + \left(\frac{\tau_{lt}}{\tau_{lt}^{\text{lim}}}\right)^2 - \frac{\sigma_l \sigma_t}{\sqrt{\sigma_l^t \sigma_l^c \sigma_t^t \sigma_t^c}} \le 1 \quad (29) $$

where $\sigma_k^t$ and $\sigma_k^c$ are, namely, the tensile and compressive failure strength in the three main directions of anisotropy of the material ($k \in \{l, t, z\}$, see Table A1), $F_i$ and $F_{ij}$ are
the components of the strength tensors of the second and fourth ranks ($i, j \in \{1, 2, \dots, 6\}$), and $\sigma_i$ are the stress components, expressed in contracted notation (i.e., $i \in \{1, 2, \dots, 6\}$) or using tensor notation (i.e., $i \in \{l, t, z, tz, lz, lt\}$).

### 6.1.2. Buckling Constraint

This constraint on the eigen values of the 40th first buckling modes aims at ensuring that no part of the optimal structure computed will buckle. In addition, following the recommendation of [64], it aims at applying a different threshold value for the eigenvalues of local and global buckling modes. The following constraint is therefore introduced: $RF \le 1$ for local buckling modes (limit load since these deformations are not detrimental for the structure). It is also recommended to include a constraint of $RF \le 1.5$ for the global buckling modes [64].

Since a post-buckling state of the different panels making up the wing structure would allow it to carry some additional loads without causing structural failure, the local buckling eigenvalues are considered as limit loads, indicating that no part of the structure can buckle. However, this buckling is an ultimate failure mode, and a panel can undergo a post-buckled state, where it can still support a load beyond its elastic buckling limit. These behaviours are highly non-linear and are not considered for the optimisation. In contrast, global buckling of a column (Euler buckling) is an ultimate failure mode, unlike buckling in a panel, which is constrained at its edge by other panels.

## 6.2. Thickness Constraints

The thickness constraints are distributed as follows:

*   **Ply fraction constraint:** a constraint is applied on the composition of each layup in order to ensure that the cumulative thickness of the plies of a given orientation equals at least 5% of the total thickness of the corresponding layup.

*   **Total thickness constraint:** a variation interval between a minimum and a maximum value for the thickness is considered to guarantee the manufacturability.

*   **Adjacency constraint:** a constraint is applied on the allowable thickness variation between two neighbouring sections of the skin.

*   **Symmetry constraint:** due to its nature, the VTP must react symmetrically to a load applied from mirror directions.

## 6.3. Sensitivity to the Number of Buckling Modes

This being stated, and keeping in mind the need to have a tool that is both precise and efficient from a calculation time point of view, a trade-off search must be carried out between the number of buckling modes and the computational time. Figure 21a below shows the impact of the buckling modes on the computational time, given various peak memory requirements, while Figure 21b gives the memory consumption related to the number of buckling modes. Impact on the computational time appears to be the most influential parameter and, considering the results, a trade-off value around 40 buckling modes seems a reasonable optimum value to consider.
<figure><img src="image_23.png" /><figcaption>Figure 21. Influence of the number of buckling modes requested and the maximal memory allocation on the computational performances and peak memory requirements.</figcaption></figure>

## 6.4. Optimisation Objectives

Optimisation has the goal of minimising operational costs through reducing the mass, and therefore the fuel consumption and environmental impact of the component, taking into account as a working hypothesis that the fabrication cost is directly related to the manufacturing duration of the component, which is in itself directly related to the component mass. Other hypotheses are not considered; for example, complex composite draping techniques that could impact the fabrication cost.

## 6.5. Decision Problem Components and Features

Based on the above definition of the optimisation problem, the main features and components of the digital thread involved in the overall optimisation framework can therefore be summarised as follows:

### 1. Inputs

$y_t = [y_t^l \ y_t^a \ y_t^m]^T \in \mathcal{Y}$, where:

(a) **Loads**: $y_t^l : \partial \mathcal{B}_t \to \mathbb{R}^5$

(b) **Material Properties**: $y_t^a \in \mathbb{R}^{N_p}$ with $N_p = 12$, based on MTM45-1/12K AS4 material [58]

(c) **Manufacturing Parameters**: $y_t^m \in \mathbb{R}^{N_v}$ where $N_v$ and $N_m$ are determined by the manufacturing process and associated parameters detailed in the Advanced Composite Cost Estimating Manual (ACCEM) cost model [51,60]

### 2. Measurements

$z_t = [z_t^l \ z_t^a \ z_t^m]^T \in \mathcal{Z}$, where:

(a) **Strains**: $z_t^l \in \mathbb{R}^{3N_s}$, for $N_s = 120$

(b) **Coupon Properties**: $z_t^a \in \mathbb{R}^{N_p}$

(c) **Manufacturing Timestamps**: $z_t^m \in \mathbb{R}^{N_m}$ for each of the $N_m$ steps of the manufacturing process.

### 3. Decision variables

Decision space: $\mathcal{U}$

$u_t = [u_t^z \ u_t^d \ u_t^r]^T \in \mathcal{U}$, where:

(a) **Thickness**: $u_t^z : \mathcal{B}_t^0 \to [e_{Min}, e_{Max}]$ ($\mathcal{B}_t^0 \subset \mathcal{B}_t$)

(b) **High-Level decisions (geometry)**: $u_t^d \in \{0, 1\}$

* Number of independent skin layups
* Number of ribs

* Stringer's cross-section

(c) Rib distribution: $u_t^r: B_t \mapsto \mathbb{R} \times \{0, 1\}$

* Distribution coefficient

* Distribution type

### Problem dimensions

$N_p$ Number of material properties (12)

$N_m$ Number of steps of the manufacturing process ($\propto$ ACCEM)

$N_s$ Number of strain sensors (120)

$N_v$ Number of manufacturing parameters ($\propto$ ACCEM)

## 4. Cost Functions

(a) Component Thicknesses: $\mathbf{x}^0 \in \mathcal{B}_t^0$

$$f_t^{\text{thick}}(\mathbf{x}^0, u_t) = e_t^{\text{tot}}(\mathbf{x}) \qquad (30)$$

where $f_t^{\text{thick}}: \mathcal{B}_t^0 \times \mathcal{U}_t \mapsto [e_{\text{Min}}, e_{\text{Max}}]$

Global cost function: The total cost associated to this metric is given by an integration of that metric over the component mid-plane $\mathcal{B}_t^0 \subset \mathcal{B}_t$

$$I_t^{\text{mass}}(u_t) = m(u_t) = \rho \int_{\mathcal{B}_t^0} f_t^{\text{thick}}(\mathbf{x}^0, u_t) \mathrm{d}a_t \qquad (31)$$

(b) Manufacturing timestamps

$I_t^{\text{man}}: \mathcal{U} \times \mathcal{Y} \mapsto [0, \infty[$ is the sum of the manufacturing timestamps predicted by the manufacturing model (ACCEM, see Section 4.5) for the system characterised by the digital thread.

Total cost function:

$$r_t(\mathcal{D}_t, u_t, y_t) = I_t^{\text{mass}}(u_t) + I_t^{\text{man}}(u_t, y_t) \qquad (32)$$

## 5. Constraint Function

(a) Structural failure: An arbitrary failure criterion for structural analysis can be put in the form $\mathcal{F}_t^m(\sigma_t, \epsilon_t, y_t^a)$, where $\sigma_t$ denotes a stress tensor with components that map to $\mathbb{R}^{3 \times 3}$ at stage $t$, $\epsilon_t$ denotes a strain tensor with components that map to $\mathbb{R}^{3 \times 3}$ at stage $t$, and $m$ denotes the $m$th mode of M total failure modes. For this study, the failure criteria that were investigated are the Tsai-Wu criteria [59] and the non-buckling failure criteria.

This expression is local and denotes the failure index at some spatial coordinate in $\mathcal{B}_t$. Failure occurs whenever $\mathcal{F}_t^m > 1$. The constraint function for failure criteria is then given by a maximisation over spatial coordinates on the component body and failure mode:

$$g_t(\mathcal{D}_t, u_t, y_t) = \max_{m \in [1, \dots, M]} \mathcal{F}_t^m(\sigma_t(\mathbf{x}, u_t, y_t), \epsilon_t(\mathbf{x}, u_t, y_t), y_t^a) \qquad (33)$$
Therefore:

$$ 
\begin{aligned}
\mathcal{F}_t^i = & \left( \frac{1}{\sigma_l^t} + \frac{1}{\sigma_l^c} \right) \sigma_l + \left( \frac{1}{\sigma_t^t} + \frac{1}{\sigma_t^c} \right) \sigma_t - \frac{\sigma_l^2}{\sigma_l^t \sigma_l^c} - \frac{\sigma_t^2}{\sigma_t^t \sigma_t^c} \\
& + \left( \frac{\tau_{lt}}{\tau_{lt}^{\lim}} \right)^2 - \frac{\sigma_l \sigma_t}{\sqrt{\sigma_l^t \sigma_l^c \sigma_t^t \sigma_t^c}} \quad \forall i \in \{1, \dots, n_{\text{layups}}\} \quad (34)
\end{aligned}
$$

$$ \mathcal{F}_t^{i+n_{\text{layups}}} = \frac{1.5}{\lambda_i} \qquad \forall i \in \{1, \dots, n_{\text{modes}}\} $$

## Problem Formulation

The mathematical expression of the approach described before is given below:

$$ \text{minimize} \quad \text{Weight}(e_{0^{\circ}}, e_{\pm 45^{\circ}}, e_{90^{\circ}}) + \text{ACCEM}(e_{0^{\circ}}, e_{\pm 45^{\circ}}, e_{90^{\circ}}) $$

$$ \text{subject to} \quad Ku = f $$

$$ (K + \lambda K_d) \cdot v = 0 $$

$$ \begin{array}{ll} \text{RF}_i^{\text{static}} \ge 1.5 & i = 1, \dots, n_{\text{sections}} \\ \lambda_i \ge 1 & i = 1, \dots, n_{\text{modes}} \end{array} \qquad (35) $$

$$ \text{pf}_i \ge 5\% \qquad i \in \{0, \pm 45, 90\} $$

$$ e_{\text{Min}} \le e_i^{\text{tot}} \le e_{\text{Max}} \quad i \in \{0, \pm 45, 90\} $$

$$ |e_{i+1}^{\text{tot}} - e_i^{\text{tot}}| \le \Delta t_{\text{adj}} \quad i = 1, \dots, n_{\text{sections}} - 1 $$

## XZ Symmetry

where pfᵢ is the ply fraction of a given ply of the layup (0°, ±45°, or 90°) defined in Figure 8, eMin = 2 mm, and eMax = 300 mm.

## 6.6. Vertical Tail Plane Optimisation Results

As an illustration of the main results, Figure 22 depicts some of the first local buckling modes of the optimised structure and Figures 23 and 24 depict the convergence of the optimisation step relying on NLPQLP constrained with 100 buckling modes and on the MMFD algorithm constrained with 40 buckling modes.

<figure><img src="image_24.png" /><figcaption>Figure 22. Buckling modes of the vertical tail plane resulting from the deterministic optimisation.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 23.** Convergence of the first optimisation step, relying on NLPQLP constrained with 100 buckling modes.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 24.** Convergence of the second optimisation step, relying on MMFD constrained with 40 buckling modes.</figcaption></figure>

The impact of the overall framework on the performances of the optimised structure, and more specifically, the benefit of the uncertainty reduction brought by the Bayesian filtering on the final optimisation of the VTP structure, is shown in Figures 25 and 26, which show, namely, the optimal VTP structure resulting from the optimisation process taking the initial and the updated knowledge of the uncertainties related to the state and inputs of the system into account. Based on these two figures, it can be observed that the impact of the uncertainties reduction of the structural inputs leads to an improvement in the distribution of the system response, allowing, in turn, a reduction of 10.5% of the total mass of the structure, with similar buckling performances and slightly lower reserve factors regarding the failure under static conditions. However, it can be observed that for both solutions, the constraint having the greatest impact on the solutions computed by the optimiser appears to be the non-buckling constraint.
<figure><img src="image_27.png" /><figcaption>**Figure 25.** Optimised structure with respect to the initial knowledge of the system state.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure 26.** Optimised structure with respect to the updated knowledge of the system state.</figcaption></figure>

# 7. Conclusions

## 7.1. Key Findings

In this paper, the digital thread architecture has been extended from the single component level, as previously considered in the literature, to the aero-structural level, exemplified with an application dedicated to an aircraft structure, the vertical tail plane of an Airbus A320 made of carbon fibre reinforced polymer (CFRP). This implementation shows promising results with respect to the uncertainty reduction capability provided by a digital thread integrating multi-disciplinary considerations (e.g., aerodynamics, structure mechanics, and manufacturing aspects) for such a structure, with a potential for cost reduction, time savings, and decision-making support for the system characterised by this digital thread, and with all of these parameters being of crucial importance for the future of complex structure design, particularly in aeronautics.

The key findings of this study are the following:

* Development of an automated framework for the parametric generation of wing geometric, structural, and aerodynamic models, taking into account the aerodynamic loads experienced by the system due to its surrounding environment instead of the arbitrary loads considered in other implementations from the literature.

* Prediction of the behaviour of mechanical systems in a given environment using mechanical simulations, buckling analyses, and modelling to calculate high-level properties, afterwards integrated with empirical and analytical models in order to automate lifecycle measurement predictions and manufacturing process durations for a complex aircraft structure.
* Recursive Bayesian estimation to combine predictions, models, and lifecycle measurements, in order to update the system state knowledge and reduce uncertainty in system core material properties, manufacturing parameters, and aerodynamic loads through data fusion.

* Optimisation of composite VTP structure by updating aerodynamic load knowledge through strain measurements and surrogate models, in order to allow efficient handling of complex systems and seamless digital thread integration via a modular, functional block interface.

## 7.2. *Ways Forward*

Based on these results, the next steps will logically address the following points:

* Implementation of the optimisation feature assisting the decision-making process based on improved system knowledge, as well as the unique and integrated design environment provided by the digital thread.

* Modelling of flight dynamics and flight control architecture interacting with the vertical tail olane, including multiple load cases, derived from the sizing cases enforced by the Certification Specifications for Large Aeroplanes (CS 25) for the VTP (gust loads, yaw manoeuvre, rudder control reversal...), as performed in [65].

* Multi-physics consideration to include the different subsystems interacting with the structure of the aircraft, studied with the final objective of extending the digital thread to the product level (complete aircraft).

* Consideration for non-linear data filtering.

* Accuracy of the simulated pressure field ($C_p$) to improve the results for the virtual data used to simulate the real data set. Consideration for multi-fidelity simulations or surrogate models instead of VSPAERO to generate the virtual measurement and the corresponding transfer functions.

* Multi-parameter measurements for characterisation of the system state behaviour to improve the prediction capabilities of the Bayesian filter.

* Finally, adapting the current methodology to more complex transitions of the system state between different stages (i.e., $F_t \neq I$ and $B_t \neq 0$) to investigate the benefit of integrating time-varying design properties or design decisions taken by the user according to this given methodology, considering that the information varies between one design generation and the next.

**Author Contributions:** Conceptualisation, S.d.L., C.B., E.B., J.J. and Y.G.; methodology, S.d.L.; software, S.d.L.; validation, S.d.L.; investigation, S.d.L.; resources, S.d.L. and C.B.; data curation, S.d.L.; writing—original draft preparation, S.d.L.; writing—review and editing, S.d.L., C.B. and Y.G.; visualisation, S.d.L.; supervision, C.B., E.B., J.J. and Y.G.; funding acquisition, E.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by Airbus through the CEDAR research chair ("Chair for Eco Design of AircRaft").

**Data Availability Statement:** Samples of the data presented in this study or supporting the findings of this article are available from the corresponding author upon reasonable request.

**Acknowledgments:** The authors would like to thank the CALMIP High Performance Computing centre for providing the computational resources required to conduct this study.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
## Abbreviations

The following abbreviations are used in this manuscript:

<table><tr><td>ACCEM</td><td>Advanced Composite Cost Estimating Manual</td></tr><tr><td>CS</td><td>Certification Specification</td></tr><tr><td>CFD</td><td>Computational Fluid Dynamics</td></tr><tr><td>CFRP</td><td>Carbon Fibre Reinforced Polymer</td></tr><tr><td>CSM</td><td>Computational Structural Mechanics</td></tr><tr><td>DoE</td><td>Design of Experiment</td></tr><tr><td>FEA</td><td>Finite Element Analysis</td></tr><tr><td>FEM</td><td>Finite Element Methods</td></tr><tr><td>MMFD</td><td>Modified Method of Feasible Directions</td></tr><tr><td>RBDO</td><td>Reliability-Based Design Optimisation</td></tr><tr><td>RF</td><td>Reserve Factor</td></tr><tr><td>VLM</td><td>Vortex Lattice Method</td></tr><tr><td>VTP</td><td>Vertical Tail Plane</td></tr></table>

## Nomenclature

The following nomenclature is used in this paper. It is, amongst other, derived from the one suggested in [22]:

<table><tr><td>C<sub>p</sub></td><td>Pressure coefficient</td></tr><tr><td>D<sub>t</sub></td><td>Digital Thread at stage t</td></tr><tr><td>e<sub>i</sub></td><td>Thickness of the ply with an orientation of i° (i ∈ {0, ±45, 90})</td></tr><tr><td>e<sub>t</sub></td><td>Pre-fit residual (innovation) at stage t</td></tr><tr><td>f</td><td>Force vector</td></tr><tr><td>g<sub>t</sub></td><td>Stage constraint function at stage t</td></tr><tr><td>H<sub>t</sub>, h<sub>t</sub></td><td>Observation (sensor) model functions</td></tr><tr><td>H<sup>a</sup><sub>t</sub>, h<sup>a</sup><sub>t</sub></td><td>Material properties model function at stage t</td></tr><tr><td>H<sup>l</sup><sub>t</sub>, h<sup>l</sup><sub>t</sub></td><td>Strain sensor model function at stage t</td></tr><tr><td>H<sup>m</sup><sub>t</sub>, h<sup>m</sup><sub>t</sub></td><td>Manufacturing timestamps model function at stage t</td></tr><tr><td>K</td><td>Stiffness matrix of the structure</td></tr><tr><td>K<sub>d</sub></td><td>Geometric stiffness matrix of the structure</td></tr><tr><td>l</td><td>Correlation length scale</td></tr><tr><td>I<sub>Ny</sub></td><td>Identity matrix of dimension N<sub>y</sub> × N<sub>y</sub></td></tr><tr><td>N</td><td>Normal distribution</td></tr><tr><td>p f<sub>i</sub></td><td>Ply fraction of the ply with an orientation of i° (i ∈ {0, ±45, 90})</td></tr><tr><td>Q<sub>t</sub></td><td>Probabilistic knowledge of the system at stage t</td></tr><tr><td>R<sub>t</sub></td><td>Resource available at stage t</td></tr><tr><td>r<sub>t</sub></td><td>Stage cost function at stage t</td></tr><tr><td>R<sub>t</sub></td><td>State of the resources available at stage t</td></tr><tr><td>u</td><td>Nodal displacement vector</td></tr><tr><td>u<sub>t</sub></td><td>Decision variable at stage t</td></tr><tr><td>Y<sub>t</sub></td><td>Uncertain inputs of the system at stage t</td></tr><tr><td>y<sub>t</sub></td><td>Real inputs of the system at stage t</td></tr><tr><td>Z<sub>t</sub></td><td>Random measurement variable at stage t</td></tr><tr><td>z<sub>t</sub></td><td>Measurements at stage t</td></tr><tr><td>α<sub>tw</sub></td><td>Twist angle of a wing section</td></tr><tr><td>ε<sub>t</sub></td><td>Measurement noise at stage t</td></tr><tr><td>ε<sub>i</sub></td><td>Strain components (contracted notation: ∀i ∈ {1, 2, ..., 6}, tensor notation: ∀i ∈ {l, t, z, tz, lz, lt})</td></tr><tr><td>λ<sub>i</sub></td><td>Eigenvalue of the ith buckling mode</td></tr><tr><td>Λ<sub>t</sub></td><td>Measurement noise covariance at stage t</td></tr><tr><td>μ<sub>t</sub></td><td>Multivariate mean vector</td></tr><tr><td>μ̃<sub>t</sub></td><td>Data assimilated input mean at the end of stage t</td></tr><tr><td>ν</td><td>Dihedral angle of a wing section</td></tr><tr><td>φ<sub>i</sub></td><td>Eigenvector of the ith buckling mode</td></tr></table>
$\varphi_{25\%}$ Sweep angle at 25% of the chord of a wing section

$\psi_t$ Value function ("cost-to-go") at stage $t$

$\Sigma_t$ Multivariate covariance matrix

$\Sigma_t$ Data assimilated input covariance at the end of stage $t$

$\sigma_i$ Stress components (contracted notation: $\forall i \in \{1, 2, \dots, 6\}$, tensor notation: $\forall i \in \{l, t, z, tz, lz, lt\}$)

$\sigma_i^c$ Compressive failure strength in the three main directions of anisotropy ($\le 0$), $\forall i \in \{l, t, z\}$

$\sigma_i^t$ Tensile failure strength in the three main directions of anisotropy, $\forall i \in \{l, t, z\}$

$\tau_{ij}^{\lim}$ Material shear strength in the three planes of symmetry, $\forall ij \in \{tz, lz, lt\}$

The following convention of notations is adopted:

$\bar{X}$ Variable resulting from the data assimilation of X

$X_l$ Material property in the fibre direction

$X_t$ Material property in the transverse direction

$X_t$ Random variable X at stage $t$

$x_t$ Realisation of the random variable X at stage $t$

$X_z$ Material property in the out-of-plane direction

$X^a$ Subset of X dedicated to material properties

$X^c$ Compressive material property

$X^l$ Subset of X dedicated to loads and strain measurements

$X^m$ Subset of X dedicated to manufacturing parameters and properties

$X^t$ Tensile material property

$X^*$ Optimal value of X

$(y_t^m)_i$ ith variable of the multivariate sample $y_t^m$

$(\Sigma_t^m)_{ij}$ ij entry of the multivariate covariance matrix $\Sigma_t^m$

$\|\cdot\|_G$ Geodesic distance operator

# Appendix A. Modelling Assumptions

## Appendix A.1. Material Properties

This section details the main properties that were considered for the CFRP material that was used for the structural model of the VTP (MTM45-1/12K AS4). Those values are summarised from the characterisation performed in [58].

<table><caption>Table A1. Properties of the MTM45-1/12K AS4 material from [58].</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td><strong>Density</strong></td><td></td><td>1600 kg/m<sup>3</sup></td></tr><tr><td colspan="3"><strong>Transversely isotropic elastic properties</strong></td></tr><tr><td>E<sub>l</sub></td><td>Young's modulus in the fibre direction</td><td>117 GPa</td></tr><tr><td>E<sub>t</sub> = E<sub>z</sub></td><td>Transverse and out-of-plane Young's Modulus</td><td>8.6 GPa</td></tr><tr><td>G<sub>lt</sub> = G<sub>lz</sub></td><td>In-plane and out-of-plane shear modulus</td><td>3.6 GPa</td></tr><tr><td>G<sub>tz</sub></td><td>Out-of-plane shear modulus</td><td>4.3 GPa</td></tr><tr><td>ν<sub>lt</sub> = ν<sub>lz</sub></td><td>Poisson's ratio</td><td>0.31</td></tr><tr><td colspan="3"><strong>Fibre failure stresses</strong></td></tr><tr><td>σ<sub>l</sub><sup>t</sup></td><td>Tensile stress limit in the fibre direction</td><td>1867 MPa</td></tr><tr><td>σ<sub>l</sub><sup>c</sup></td><td>Compressive stress limit in the fibre direction</td><td>-1398 MPa</td></tr><tr><td>σ<sub>l</sub><sup>t</sup></td><td>Tensile stress limit in the transverse direction</td><td>48 MPa</td></tr><tr><td>σ<sub>l</sub><sup>c</sup></td><td>Compressive stress limit in the transverse direction</td><td>-185 MPa</td></tr><tr><td>τ<sub>lt</sub><sup>lim</sup></td><td>In-plane shear strength at 5% strain (τ<sub>lt</sub><sup>5%</sup>)</td><td>65 MPa</td></tr><tr><td colspan="3"><strong>Fibre failure strains</strong></td></tr><tr><td>ε<sub>l</sub><sup>t</sup></td><td>Tensile strain limit in the fibre direction</td><td>14,250 με</td></tr><tr><td>ε<sub>l</sub><sup>c</sup></td><td>Compressive strain limit in the fibre direction</td><td>-11,100 με</td></tr><tr><td>ε<sub>l</sub><sup>t</sup></td><td>Tensile strain limit in the transverse direction</td><td>14,150 με</td></tr><tr><td>ε<sub>l</sub><sup>c</sup></td><td>Compressive strain limit in the transverse direction</td><td>-10,900 με</td></tr><tr><td>γ<sub>lt</sub><sup>lim</sup></td><td>In-plane shear strain limit</td><td>18,700 με</td></tr></tbody></table>
## Appendix A.2. Modelling Parameters

<table><caption>**Table A2.** Aerodynamic conditions under which the VTP structure has been optimised.</caption><thead><tr><th colspan="4">Aerodynamic Configuration</th></tr></thead><tbody><tr><td>α</td><td>Angle of attack</td><td>0°</td><td></td></tr><tr><td>β</td><td>Sideslip Angle</td><td>5°</td><td></td></tr><tr><td>M</td><td>Mach Number</td><td>0.3</td><td></td></tr><tr><td>ρ</td><td>Air density</td><td>1.225</td><td>kg/m<sup>3</sup></td></tr><tr><td>V<sub>∞</sub></td><td>Free stream velocity</td><td>102</td><td>m/s</td></tr><tr><td>Re</td><td>Reynold's Number</td><td>4 × 10<sup>7</sup></td><td></td></tr></tbody></table>

<table><caption>**Table A3.** Wing parameters used to model the VTP geometry, based on [66].</caption><thead><tr><th colspan="4">Geometric Parameters</th></tr></thead><tbody><tr><td>b</td><td>Span</td><td>6.26</td><td>m</td></tr><tr><td>c<sub>r</sub></td><td>Root Chord</td><td>6.421</td><td>m</td></tr><tr><td>c<sub>t</sub></td><td>Tip Chord</td><td>1.900</td><td>m</td></tr><tr><td>S</td><td>Area</td><td>21.5</td><td>m<sup>2</sup></td></tr><tr><td>φ<sub>25%</sub></td><td>Quarter-chord sweep angle</td><td>34°</td><td></td></tr><tr><td>α<sub>tw</sub></td><td>Twist angle</td><td>0°</td><td></td></tr><tr><td>ν</td><td>Dihedral</td><td>0°</td><td></td></tr><tr><td>ε</td><td>Taper Ratio</td><td>0.296</td><td></td></tr><tr><td>λ</td><td>Aspect Ratio</td><td>1.8</td><td></td></tr><tr><td>af</td><td>Airfoil sections</td><td>NACA 0012</td><td></td></tr></tbody></table>

<table><caption>**Table A4.** Overview of the parameters dedicated to the structural model of the VTP and initial conditions used by the structural optimiser.</caption><tbody><tr><td colspan="2"><strong>Structural parameters</strong></td></tr><tr><td>Number of ribs</td><td>9</td></tr><tr><td>Rib distribution</td><td>Linear</td></tr><tr><td>Rib holes</td><td>Ø</td></tr><tr><td>Skin sections of independent thicknesses</td><td>1</td></tr><tr><td>Number of spars</td><td>2</td></tr><tr><td>Spar positions</td><td>20% chord<br>70% chord</td></tr><tr><td colspan="2"><strong>Numerical parameters</strong></td></tr><tr><td>Number of buckling modes</td><td>100</td></tr><tr><td>Structural Mesh seed size</td><td>40 mm</td></tr><tr><td colspan="2"><strong>Initial layup thicknesses</strong></td></tr><tr><td>Spars Thickness</td><td>15 mm</td></tr><tr><td>Skin Thickness</td><td>2 mm</td></tr><tr><td>Root rib Thickness</td><td>10 mm</td></tr><tr><td>Tip rib Thickness</td><td>6 mm</td></tr><tr><td>Internal ribs thicknesses</td><td>6 + 4 · <sup>Ribpos</sup>/<sub>6.26</sub> mm</td></tr></tbody></table>

## References

1. Chaudhuri, A.; Jasa, J.; Martins, J.R.R.A.; Willcox, K.E. Multifidelity Optimization Under Uncertainty for a Tailless Aircraft. In Proceedings of the 2018 AIAA Non-Deterministic Approaches Conference, Kissimmee, FL, USA, 8–12 January 2018. [CrossRef]

2. Ledermann, C.; Hanske, C.; Wenzel, J.; Ermanni, P.; Kelm, R. Associative parametric CAE methods in the aircraft pre-design. *Aerosp. Sci. Technol.* **2005**, *9*, 641–651. [CrossRef]

3. Chell, B.W.; Hoffenson, S.; Blackburn, M.R. A Comparison of Multidisciplinary Design Optimization Architectures with an Aircraft Case Study. In Proceedings of the AIAA Scitech 2019 Forum, San Diego, CA, USA, 7–11 January 2019. [CrossRef]

4. Grose, D. Reengineering the aircraft design process. In Proceedings of the 5th Symposium on Multidisciplinary Analysis and Optimization, Panama City Beach, FL, USA, 7–9 September 1994. [CrossRef]

5. Bertoni, M.; Bertoni, A. Designing solutions with the product-service systems digital twin: What is now and what is next? *Comput. Ind.* **2022**, *138*, 103629. [CrossRef]

6. Singh, M.; Fuenmayor, E.; Hinchy, E.; Qiao, Y.; Murray, N.; Devine, D. Digital Twin: Origin to Future. *Appl. Syst. Innov.* **2021**, *4*, 36. [CrossRef]
7. Grieves, M.W. Completing the Cycle: Using PLM Information in the Sales and Service Functions. In Proceedings of the SME Management Forum, Troy, MI, USA, 31 October 2002.

8. Glaessgen, E.; Stargel, D. The Digital Twin Paradigm for Future NASA and U.S. Air Force Vehicles. In Proceedings of the 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, 20th AIAA/ASME/AHS Adaptive Structures Conference, 14th AIAA, Honolulu, HI, USA, 23–26 April 2012. [CrossRef]

9. Hazbon, O.; Gutierrez, L.; Bil, C.; Napolitano, M.; Fravolini, M. Digital Twin Concept for Aircraft System Failure Detection and Correction. In Proceedings of the AIAA Aviation 2019 Forum, Dallas, TX, USA, 17–21 June 2019. [CrossRef]

10. Kallou, E.; Gharbi, A.; Briceno, S.I.; Mavris, D.N. Analysis and Visualization of Failure Modes on the Wing Rib. In Proceedings of the 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Kissimmee, FL, USA, 8–12 January 2018. [CrossRef]

11. Korobenko, A.; Pigazzini, M.; Singh, V.; Kim, H.; Allaire, D.L.; Willcox, K.E.; Marsden, A.; Bazilevs, Y. Dynamic-Data-Driven Damage Prediction in Aerospace Composite Structures. In Proceedings of the 17th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Washington, DC, USA, 13–17 June 2016. [CrossRef]

12. Li, C.; Mahadevan, S.; Ling, Y.; Choze, S.; Wang, L. Dynamic Bayesian Network for Aircraft Wing Health Monitoring Digital Twin. *AIAA J.* **2017**, *55*, 930–941. [CrossRef]

13. West, T.D.; Pyster, A. Untangling the Digital Thread: The Challenge and Promise of Model-Based Engineering in Defense Acquisition. *Insight* **2015**, *18*, 45–55. [CrossRef]

14. Pang, T.Y.; Pelaez Restrepo, J.D.; Cheng, C.T.; Yasin, A.; Lim, H.; Miletic, M. Developing a Digital Twin and Digital Thread Framework for an ‘Industry 4.0’ Shipyard. *Appl. Sci.* **2021**, *11*, 1097. [CrossRef]

15. McGibney, A. Digital Thread in Smart Manufacturing. In Proceedings of the Leveraging Applications of Formal Methods, Verification and Validation, Rhodes, Greece, 22–30 October 2022; Margaria, T., Pesch, D., McGibney, A., Eds.; Lecture Notes in Computer Science; Elsevier: Cham, Switzerland, 2022; Volume 13704, pp. 179–183.

16. Abdel-Aty, T.A.; Negri, E. Conceptualizing the digital thread for smart manufacturing: A systematic literature review. *J. Intell. Manuf.* **2024**, *35*, 3629–3653. [CrossRef]

17. Kapteyn, M.G.; Pretorius, J.V.R.; Willcox, K.E. A probabilistic graphical model foundation for enabling predictive digital twins at scale. *Nat. Comput. Sci.* **2021**, *1*, 337–347. [CrossRef]

18. Kinard, D.A. F-35 Advanced Manufacturing and the Digital Thread. In Proceedings of the 2018 Aviation Technology, Integration, and Operations Conference, Atlanta, GA, USA, 25–29 June 2018. [CrossRef]

19. Kraft, E. HPCMP CREATE-AV and the Air Force Digital Thread. In Proceedings of the 53rd AIAA Aerospace Sciences Meeting, Kissimmee, FL, USA, 5–9 January 2015. [CrossRef]

20. Kraft, E.M. The Air Force Digital Thread/Digital Twin—Life Cycle Integration and Use of Computational and Experimental Knowledge. In Proceedings of the 54th AIAA Aerospace Sciences Meeting, San Diego, CA, USA, 4–8 January 2016. [CrossRef]

21. Kraft, E.M. Integrating Model-Based Systems Engineering and Uncertainty Quantification to Design Optimum Aircraft Test Campaigns. In Proceedings of the 2018 Applied Aerodynamics Conference, Atlanta, GA, USA, 25–29 June 2018. [CrossRef]

22. Singh, V. Towards a Feedback Design Process Using Digital Thread. Ph.D. Thesis, Massachusetts Institute of Technology, Boston, MA, USA, 2019.

23. Gharbi, A.; Briceno, S.I.; Mavris, D.N. Enabling the Digital Thread in Commercial Aircraft Companies. In Proceedings of the 2018 Aviation Technology, Integration, and Operations Conference, Atlanta, GA, USA, 25–29 June 2018. [CrossRef]

24. Oroz, J.; Roohi, Z.A.; Abelezele, S.; Fronk, G.; Al Fawares, R.; Pinon-Fischer, O.J.; Gharbi, A.; Mavris, D.N.; Petersen, M.; Karl, A.; et al. Implementing the Digital Thread—A Proof-of-Concept. In Proceedings of the AIAA SCITECH 2023 Forum, National Harbor, MD, USA & Online, 23–27 January 2023; p. 1405. [CrossRef]

25. Eskue, N. Digital Thread Roadmap for Manufacturing and Health Monitoring the Life Cycle of Composite Aerospace Components. *Aerospace* **2023**, *10*, 146. [CrossRef]

26. Fioriti, M.; Boggero, L.; Corpino, S.; Isyanov, A.; Mirzoyan, A.; Lombardi, R.; D'Ippolito, R. Automated Selection of the Optimal On-board Systems Architecture within MDO Collaborative Environment. In Proceedings of the 18th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Denver, CO, USA, 5–9 June 2017. [CrossRef]

27. Fioriti, M.; Boggero, L.; Corpino, S.; Prakasha, P.S.; Ciampa, P.D.; Nagel, B. The Effect of Sub-systems Design Parameters on Preliminary Aircraft Design in a Multidisciplinary Design Environment. *Transp. Res. Procedia* **2018**, *29*, 135–145. [CrossRef]

28. Gallard, F. Optimisation de Forme d'un Avion Pour sa Performance sur une Mission. PhD Thesis, École Doctorale Aéronautique-Astronautique, Toulouse, France, 2014.

29. Goertz, S.; Ilic, C.; Jepsen, J.; Leitner, M.; Schulze, M.; Schuster, A.; Scherer, J.; Becker, R.; Zur, S.; Petsch, M. Multi-Level MDO of a Long-Range Transport Aircraft Using a Distributed Analysis Framework. In Proceedings of the 18th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Denver, CO, USA, 5–9 June 2017. [CrossRef]

30. Hwang, J.T.; Jain, A.V.; Ha, T.H. Large-scale multidisciplinary design optimization—Review and recommendations. In Proceedings of the AIAA Aviation 2019 Forum, Dallas, TX, USA, 17–21 June 2019. [CrossRef]
31. Singh, V.; Willcox, K.E. Engineering Design with Digital Thread. In Proceedings of the 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Kissimmee, FL, USA, 8–12 January 2018. [CrossRef]

32. Curran, R.; Price, M.; Raghunathan, S.; Benard, E.; Crosby, S.; Castagne, S.; Mawhinney, P. Integrating Aircraft Cost Modeling into Conceptual Design. *Concurr. Eng.* **2005**, *13*, 321–330. [CrossRef]

33. Gharbi, A.; Sarojini, D.; Kallou, E.; Harper, D.J.; Petitgenet, V.; Rancourt, D.; Briceno, S.I.; Mavris, D.N. STAnDD: A Single Digital Thread Approach to Aircraft Detailed Design. In Proceedings of the 55th AIAA Aerospace Sciences Meeting, Grapevine, TX, USA, 9–13 January 2017. [CrossRef]

34. AIAA Digital Engineering Integration Committee. *Digital Twin: Reference Model, Realizations & Recommendations*; AIAA, AIA and NAFEMS Implementation Paper; Aerospace Industries Association: Arlington, VA, USA, 2023.

35. Singh, V.; Willcox, K.E. Decision-Making Under Uncertainty for a Digital Thread-Enabled Design Process. *J. Mech. Des.* **2021**, *143*, 091707. [CrossRef]

36. Renganathan, S.A.; Harada, K.; Mavris, D.N. Aerodynamic Data Fusion Toward the Digital Twin Paradigm. *AIAA J.* **2020**, *58*, 3902–3918. [CrossRef]

37. Pagan, J.E.; Huynh, D.; Schafer, S.R.; Pinon-Fischer, O.J.; Mavris, D.N. Revolutionary Affordable Architecture Generation & Evaluation—Application to a System of Attributable UAVs. In Proceedings of the 2018 AIAA Aerospace Sciences Meeting, Kissimmee, FL, USA, 8–12 January 2018. [CrossRef]

38. Price, M.; Raghunathan, S.; Curran, R. An integrated systems engineering approach to aircraft design. *Prog. Aerosp. Sci.* **2006**, *42*, 331–376. [CrossRef]

39. Tuegel, E.J.; Kobryn, P.; Zweber, J.V.; Kolonay, R.M. Digital Thread and Twin for Systems Engineering: Design to Retirement. In Proceedings of the 55th AIAA Aerospace Sciences Meeting, Grapevine, TX, USA, 9–13 January 2017. [CrossRef]

40. Zhang, Q.; Liu, J.; Chen, X. A Literature Review of the Digital Thread: Definition, Key Technologies, and Applications. *Systems* **2024**, *12*, 70. [CrossRef]

41. Stark, J. *Product Lifecycle Management*, 3rd ed.; Decision Engineering; Springer International Publishing: Cham, Switzerland, 2015; Volume 1. [CrossRef]

42. Azamatov, A.; Lee, J.W.; Byun, Y.H. Comprehensive aircraft configuration design tool for Integrated Product and Process Development. *Adv. Eng. Softw.* **2011**, *42*, 35–49. [CrossRef]

43. AIAA Digital Engineering Integration Committee. *Digital Twin: Definition & Value*; AIAA and AIA Position Paper; Aerospace Industries Association: Arlington, VA, USA, 2020.

44. Tuegel, E.J.; Ingraffea, A.R.; Eason, T.G.; Spottswood, S.M. Reengineering Aircraft Structural Life Prediction Using a Digital Twin. *Int. J. Aerosp. Eng.* **2011**, *2011*, 154798. [CrossRef]

45. Singh, V.; Willcox, K.E. Methodology for Path Planning with Dynamic Data-Driven Flight Capability Estimation. *AIAA J.* **2017**, *55*, 2727–2738. [CrossRef]

46. Pinon Fischer, O.J.; Matlik, J.F.; Schindel, W.D.; French, M.O.; Kabir, M.H.; Ganguli, J.S.; Hardwick, M.; Arnold, S.M.; Byar, A.D.; Lewe, J.; et al. Digital Twin: Reference Model, Realizations, and Recommendations. *Insight 2022*, 25, 50–55. [CrossRef]

47. Toriello, M. Creating the Digital Thread. *Insight 2022*, 25, 34–37. [CrossRef]

48. Abollado, J.R.; Shehab, E.; Bamforth, P. Challenges and Benefits of Digital Workflow Implementation in Aerospace Manufacturing Engineering. *Procedia CIRP 2017*, *60*, 80–85. [CrossRef]

49. AIAA Digital Engineering Integration Committee. *Digital Thread: Definition, Value and Reference Model*; AIAA, AIA and NAFEMS Implementation Paper; Aerospace Industries Association: Reston, VA, USA, 2023.

50. Tuegel, E. The Airframe Digital Twin: Some Challenges to Realization. In Proceedings of the 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, 20th AIAA/ASME/AHS Adaptive Structures Conference, 14th AIAA, Honolulu, HI, USA, 23–26 April 2012. [CrossRef]

51. LeBlanc, D.J.; Lorenzana, J.; Kokawa, A.; Bettner, T.; Timson, F. *Advanced Composite Cost Estimating Manual, Volume 1*; Technical Report AFFDL-TR-76-87/ADA041495; Northrop Corporation Aircraft Division: Hawthorne, CA, USA, 1976.

52. Okamoto, K.; Tsiotras, P. Optimal Stochastic Vehicle Path Planning Using Covariance Steering. *IEEE Robot. Autom. Lett.* **2019**, *4*, 2276–2281. [CrossRef]

53. McDonald, R.A.; Gloudemans, J.R. Open Vehicle Sketch Pad: An Open Source Parametric Geometry and Analysis Tool for Conceptual Aircraft Design. In Proceedings of the AIAA SCITECH 2022 Forum, San Diego, CA, USA & Virtual, 3–7 January 2022. [CrossRef]

54. Erickson, L.L. *Panel Methods: An Introduction*; Technical report NASA Technical Paper 2995; National Aeronautics and Space Administration: Moffett Field, CA, USA, 1990.

55. Ingraham, D.; Kinney, D. Towards Low-Noise Design of a Proprotor with the Vortex Lattice Method and Gradient-Based Optimization. In Proceedings of the AIAA SCITECH 2024 Forum, Orlando, FL, USA, 8–12 January 2024. [CrossRef]

56. Sheridan, C.N.; Pham, D.D.; Whiteside, S. Evaluation of VSPAERO Analysis Capabilities for Conceptual Design of Aircraft with Propeller-Blown Wings. In Proceedings of the AIAA Aviation 2021 Forum, Virtual Event, 2–6 August 2021. [CrossRef]
57. Valarezo, W.O.; Chin, V.D. Method for the prediction of wing maximum lift. *J. Aircr.* **1994**, *31*, 103–109. [CrossRef]

58. Clarkson, E. NCP-RP-2008-004 N/C: MTM45-1/12K AS4 145gsm 32\% RW Unidirectional Qualification Statistical Analysis Report; Technical Report; Wichita State University: Wichita, KS, USA, 2016.

59. Tsai, S.W.; Wu, E.M. A General Theory of Strength for Anisotropic Materials. *J. Compos. Mater.* **1971**, *5*, 58–80. [CrossRef]

60. Gutowski, T.; Hoult, D.; Dillon, G.; Neoh, E.T.; Muter, S.; Kim, E.; Tse, M. Development of a theoretical cost model for advanced composite fabrication. *Compos. Manuf.* **1994**, *5*, 231–239. [CrossRef]

61. Van Der Velden, A.; Wujek, B.; Koch, P. iSIGHT-FD: A Tool for Multi-Objective Data Analysis. In Proceedings of the 12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Victoria, BC, Canada, 10–12 September 2008. [CrossRef]

62. Schittkowski, K. NLPQLP: A Fortran Implementation of a Sequential Quadratic Programming Algorithm with Distributed and Non-Monotone Line Search; Report; Department of Computer Science, University of Bayreuth: Bayreuth, Germany, 2006. Available online: https://www.schittkowski.de/downloads/numerical/NLPQLP.pdf (accessed on 16 December 2024).

63. Vanderplaats, G.N. An efficient feasible directions algorithm for design synthesis. *AIAA J.* **1984**, *22*, 1633–1640. [CrossRef]

64. Qian, J.; Alonso, J.J. Design of Unconventional Aircraft Structures with Aeroelastic and Buckling Considerations. In Proceedings of the AIAA SCITECH 2022 Forum, San Diego, CA, USA & Virtual, 3–7 January 2022. [CrossRef]

65. Scarano, M.; de Longueville, S.; Jézégou, J.; Gourinat, Y.; Bénard, E. Automated sizing chain of aircraft structures: A Vertical Tail Plane use case. *J. Phys. Conf. Ser.* **2023**, 2526, 012031. [CrossRef]

66. Jenkinson, L.; Simpkin, P.; Rhodes, D. Civil Jet Aircraft Design—Aircraft Data File—Airbus Aircraft. 2001. Available online: https://booksite.elsevier.com/9780340741528/appendices/data-a/table-1/table.htm (accessed on 26 June 2024).

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.