PhD Anna Calissano
email: anna.calissano@inria.fr

| Theodore Papadopoulo

| Xavier Pennec

Samuel Deslauriers-Gauthier

Graph Alignment Exploiting the Spatial Organisation Improves the Similarity of Brain Networks

Every brain is unique, having its structural and functional organisation shaped by both genetic and environmental factors over the course of its development. Brain image studies tend to produce results by averaging across a group of subjects, under the common assumption that it is possible to subdivide the cortex into homogeneous areas while maintaining a correspondence across subjects. We investigate this assumption: can the structural properties of a specific region of an atlas be assumed to be the same across subjects? This question is addressed by looking at the network representation of the brain, with nodes corresponding to brain regions and edges to their structural relationships. Using an unsupervised graph matching strategy, we align the structural connectomes of a set of healthy subjects, considering parcellations of different granularity, to understand the connectivity misalignment between regions.

Firstly, we compare the obtained permutations with four different algorithm initializations: Spatial Adjacency, Identity, Barycenter and Random. Our results suggest that applying an alignment strategy improves the similarity across subjects when the number of parcels is above 100 and when using Spatial Adjacency and Identity initialization (the most 1 plausible priors). Secondly, we characterize the obtained permutations, revealing that the majority of permutations happens between neighbours parcels. Lastly, we study the spatial distribution of the permutations. By visualizing the results on the cortex, we observe no clear spatial patterns on the permutations and all the regions across the context are mostly permuted with first and second order neighbours. Every brain is unique, having its structural and functional organisation shaped by both genetic and environmental factors over the course of its development [START_REF] Thompson | Genetic influences on brain structure[END_REF]. This uniqueness directly translates into inter-subject variability in the location and extent of functional brain areas, the network organisation of structural connectivity, and the structure-function relationship [START_REF] Vigneau | Metaanalyzing left hemisphere language areas: phonology, semantics, and sentence processing[END_REF][START_REF] Sporns | The human connectome: a structural description of the human brain[END_REF][START_REF] Mueller | Individual variability in functional connectivity architecture of the human brain[END_REF][START_REF] Tate | Probabilistic map of critical functional regions of the human cerebral cortex: Broca's area revisited[END_REF].

In imaging studies which tend to produce results by averaging across a group of subjects, this variability introduces temporal and spatial variance in the localisation of specialized brain areas and their network connectivity. In particular, this variability undermines the common assumption that it is possible to subdivide the cortex into homogeneous areas while maintaining a correspondence across subjects. In other words, can we assume the structural and functional properties of a specific region of an atlas to be the same across subjects? This question is of particular interest to studies that represent the brain as a network, where nodes correspond to brain regions and edges to their structural or functional relationship [START_REF] Sporns | The human connectome: a structural description of the human brain[END_REF][START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF]. If there is a spatial variability across subject, then the nodes of the network may not be matched, and the networks of distinct subjects should not be compared directly.

A possible solution to this problem is to make no assumption on the correspondence between the nodes of the network and to try to derive it using a data driven approach. This solution acts at the connectivity matrix level and jointly addresses: a possible spatial misalignment between the parcels that occurred during the registration phase; and a possible structural connectivity difference across subjects. Aligning the nodes of two networks based on their connectivity properties is referred as graph matching or network alignment [START_REF] Conte | Thirty years of graph matching in pattern recognition[END_REF]. Given two networks, network alignment corresponds to finding a permutation of the nodes of one graph which minimizes a cost function -usually a dissimilarity measure between the graphs. As the problem is NP-hard, many different techniques have been proposed in the literature, from exact matching to inexact matching (see [START_REF] Conte | Thirty years of graph matching in pattern recognition[END_REF][START_REF] Emmert-Streib | Fifty years of graph matching, network alignment and network comparison[END_REF] for exhaustive reviews on the topic). Notice that network alignment differs from the fingerprint matching and the parcel matching problem [START_REF] Mars | Comparing brains by matching connectivity profiles[END_REF]Gallardo et al., 2018a) as it only relies on the intrinsic properties of the network, for example edge weights, instead of node specific features. The network alignment problem gained attention in the context of neuroscience, aiming at understanding the heterogeneity inter-subjects brain connectivity networks. However, current approaches were designed for the alignment of general networks, for example social networks, and do not typically impose spatial constraints on the node permutations. In the context of brain networks, arbitrary permutation may result in the large spatial displacement of brain areas which is highly improbable. Therefore, while local permutation explaining inter-subject variability is expected, permutation resulting in large displacement should be discouraged (e.g., permuting regions in the occipital and frontal lobes). Therefore, a graph alignment strategy promoting or enforcing local permutation may be more appropriate to brain network alignment. In Figure 1, we illustrate the difference between the identity matching, where the regions across subjects are considered as the same nodes, and the graph matching, where we search for a better matching of regions across subjects. Both different similarity measures as well as different matching strategies can be defined, creating a wide range of approaches [START_REF] Frigo | Network alignment and similarity reveal atlas-based topological differences in structural connectomes[END_REF][START_REF] Deslauriers-Gauthier | A unified framework for multimodal structure-function mapping based on eigenmodes[END_REF][START_REF] Rasero | Consensus clustering approach to group brain connectivity matrices[END_REF][START_REF] Frost | Improved correspondence of resting-state networks after macroanatomical alignment[END_REF][START_REF] Takerkart | Graph-based inter-subject pattern analysis of fmri data[END_REF][START_REF] Becker | Spectral mapping of brain functional connectivity from diffusion imaging[END_REF][START_REF] Weinstein | A simple permutation-based test of intermodal correspondence[END_REF][START_REF] Ravindra | Rigid graph alignment[END_REF]. In this work, we aimed at exploring the stability of group-wise parcellations when performing inter-subject analysis by looking at network alignment performances and characteristics on a set of control subjects. Given a set of structural connectivity matrices derived from a extrinsic connectivity parcellation (Gallardo et al., 2018b), we study the network alignment results. Among the different alignment strategies, we choose the Fast Approximate Quadratic (FAQ) algorithm [START_REF] Vogelstein | Fast approximate quadratic programming for graph matching[END_REF], as it represents a state of the art in the network alignment literature. It is also known to be well tailored for the analysis of brain connectivity networks [START_REF] Glaser | The development and analysis of integrated neuroscience data[END_REF][START_REF] Vogelstein | Fast approximate quadratic programming for graph matching[END_REF]. The analysis is conducted on the different priors of the matching procedure, by imposing different initial conditions to the graph matching algorithm, and at different parcels granularities. The proposed approach allows for the definition of the spatial initialization as a suitable initialization for structural connectivity data and the identification of the most stable and most unstable brain area in brain matching across subjects.

Our main contribution consists in showing with a purely data driven approach how inter-subject variability of structural connectivity graphs can be decreased using graph matching and with the correct initialization. We also characterized the area of the brain that are more often permuted, showing where the misalignment often resides.

| MATERIALS AND METHODS

| Brain Graph Matching

A weighted graph is a mathematical object G i = (N , E i , w i ) with a set of nodes N = {1, . . . , n }, a set of edges E = { (j , k ) ∈ {0, 1}, j , k ∈ N }, and a weight function w i : E → . In matrix term, a graph can be represented as a weighted adjacency matrix G i ∈ n×n . In the context of brain networks, the number of nodes n corresponds to the number of parcels per hemisphere and w often encodes the connectivity, either structural or functional, between parcels. Given the matrix representation, different metrics can be defined to compare graphs [START_REF] Donnat | Tracking network dynamics: A survey using graph distances[END_REF][START_REF] Frigo | Network alignment and similarity reveal atlas-based topological differences in structural connectomes[END_REF][START_REF] Bassett | Small-world brain networks revisited[END_REF][START_REF] Pandit | Traumatic brain injury impairs small-world topology[END_REF]. In this work, we are using the Frobenious norm:

d F (G 1 , G 2 ) = | |G 1 -G 2 | | 2 = n i =1 n j =1 (G 1 (i , j ) -G 2 (i , j )) 2
Note that the Frobenious norm is not invariant under node permutations, it thus allows for an analysis of the suggested nodes' permutations. In addition, it allows for a quick graph matching procedure [START_REF] Vogelstein | Fast approximate quadratic programming for graph matching[END_REF] and it has been broadly used in the context of structural connectivity comparison. Graph matching corresponds in finding an injective map from the nodes of the graph G i to the nodes of the graph G j . In geometrical terms, the nodes (parcels in this context) permutation is expressed using the matrix representation of the permutation group P = {P ∈ {0, 1} n×n : P T 1 = P 1 = 1}. Given an adjacency matrix G i , it can be permuted into potentially n! different adjacency matrices by applying the permutation action. The results is an equivalence class of graphs [G i ] = {P T G i P : P ∈ P }, obtained by permuting the rows and the columns of the adjacency matrices by matrix multiplication. If we are using the metric as the alignment cost function and the permutation matrix representation, aligning two networks G i and G j corresponds in finding the permutation matrix P i j minimizing:

d ( [G i ], [G j ]) = min P ∈P d F (G i , P T G j P ) (1)
This minimization problem corresponds to the graph matching problem. In addition, d

( [G i ], [G j ]
) is an induced metric on on the quotient space n×n /P obtained by applying the permutation action to the space of adjacency matrices.

This space is often called graph space (see for details [START_REF] Jain | Structure spaces[END_REF]Calissano et al. 2023). Note that the graph matching is only the simplest operation which can be performed in graph space, as the theoretical framework for the analysis of graphs is richer, including for example Fréchet mean (Jain and Obermayer, 2010), principal components (Calissano et al., 2023), regression [START_REF] Calissano | Graph-valued regression: Prediction of unlabelled networks in a non-euclidean graph space[END_REF][START_REF] Severn | Non-parametric regression for networks[END_REF], and testing [START_REF] Ginestet | Hypothesis testing for network data in functional neuroimaging[END_REF].

In the context of brain network analysis, graph matching represents a way to understand if there is a misalignment between subjects, and provides a way to correct it. A misalignment is detected when the optimal permutation matrix obtained by equation ( 1) is not the identity. In this situation, the similarity between subjects can be improved by permuting the labels of one subject, implying that the initial labels do not match. The optimal permutation matrix P to apply to the labels of subject j is also produced by Equation ( 1) and can thus be applied as preprocessing step before further analysis or comparison of the networks.

| FAQ algorithm for graph matching

Graph Matching is an NP-hard problem which has been largely explored in the literature (see [START_REF] Conte | Thirty years of graph matching in pattern recognition[END_REF] for a review). Among the different inexact graph matching algorithms, we selected the Fast Quadratic Assignment method proposed by [START_REF] Vogelstein | Fast approximate quadratic programming for graph matching[END_REF] as it represents the state of the art method and it has been used in the context of human and insect brain matching [START_REF] Frigo | Network alignment and similarity reveal atlas-based topological differences in structural connectomes[END_REF][START_REF] Winding | The connectome of an insect brain[END_REF]. The FAQ matching is a fast method for inexact graph matching based on the Frobenious norm between adjacency matrices and it allows to select an initialization (a prior) to the matching procedure, which is important in this context where not all the regions permutations are meaningful.

The idea behind the methodology is to relax the domain of permutation matrices from the discrete set P to the continuous convex hull of doubly stochastic matrices D = {P ∈ n×n : P T 1 = P 1 = 1, P ≥ 0}. Starting from the minimization of the Frobenious norm between two adjacency matrices G i , G j ∈ n×n , the problem is manipulated into the following minimization:

minimizing f (P ) = -trace(G i P G T j P T ) (2)
subject to P ∈ D

The local solution in found by gradient descend in the space D. To obtain a permutation matrix as output, the solution is then projected onto P. As many optimization problems, the FAQ algorithm relies on an initialization P 0 , which could be interpreted as a prior knowledge on the parcels' matching. Given the complexity of the problem, this prior knowledge has a significant impact on the quality of the final solution. In FAQ, such initialization matrices has to be a doubly stochastic matrix [START_REF] Gagniuc | Markov chains: from theory to implementation and experimentation[END_REF]. Given a permutation matrix, we obtain a doubly stochastic matrix using the Sinkhorn Knopp Algorithm [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF]. We are proposing four different initializations:

1. Barycenter: starting from a non informative doubly stochastic matrix where all entries have the same value. This corresponds to picking all the regions' permutations as equally probable;

2. Random: starting from a random doubly stochastic matrix. This encode guessing some permutation at random; 3. Identity: starting from the identity matrix. This initialization suggests that the current labels are correct;

4. Spatial Adjacency: starting from the spatial structure of the regions. This matrix is derived from the parcellation of the cortical surface where adjacent regions are considered as neighbors. Such initialization suggests a permutation with the spatial neighbours.

| Structural brain networks

Our dataset consists of 100 healthy young adults from the Human Connectome project, (U100 subject group) (Van Essen et al., 2012; [START_REF] Glasser | The minimal preprocessing pipelines for the Human Connectome Project[END_REF]. For each subject, a five-tissue-type image [START_REF] Smith | Anatomically-constrained tractography: improved diffusion MRI streamlines tractography through effective use of anatomical information[END_REF] was obtained using the Freesurfer pipeline (Fischl, 2012) invoked through Mrtrix3 [START_REF] Tournier | MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation[END_REF]. A response function was estimated for the white matter, gray matter, and cerebrospinal fluid using a maximal spherical harmonic order of 8 for all tissues [START_REF] Jeurissen | Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data[END_REF]. The fiber orientation distribution functions (fODFs) were then computed using the multi-shell multi-tissue constrained spherical deconvolution algorithm [START_REF] Jeurissen | Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data[END_REF]. Finally, the fODFs were used as input for probabilistic anatomically constrained tractography performed with the iFOD2 algorithm [START_REF] Smith | Anatomically-constrained tractography: improved diffusion MRI streamlines tractography through effective use of anatomical information[END_REF] seeding from the gray matter -white matter interface and obtaining a total of five million streamlines per subject. The structural brain connectivity network are built from the parcellation introduced in Gallardo et al.

(2018b), as we are interested in group-wise parcellation. This parcellation is built based on 138 tractograms built from co-registered cortical surfaces from Human Connectome Project [START_REF] Van Essen | The human connectome project: a data acquisition perspective[END_REF]. The resulting parcellation is defined on the average FreeSurfer surface (MSMAll) with 32k vertices. A clustering procedure offers different regions granularity by selecting a threshold and sequentially aggregating regions with similar connectivity structure. Thus, we study the parcellation at different granularity from 50 up to 1000 nodes, to assess the threshold where the original alignment across subjects is inexact.

| Experiments Set Up

The aim is to study how graph matching behaves when comparing healthy subjects' brain structural graph. By analysing the graph matching results at different parcellations and with different initialization, we can infer which areas of the brain most suffer from misalignment errors in the parcellation strategy. To evaluate the quality of the alignment, we use the ratio between F opt , the objective function obtained FAQ graph matching, and F i d the objective function obtained with the identity matching (i.e., leaving the graphs with the current nodes labels).

F i j opt (P ) = trace(G i P G T j P T ) F i j i d = trace(G i I T G j I ) (3)
Given P ∈ P, the permutation matrix estimated via FAQ and I ∈ n×n the identity matrix of the original alignment,

F i j opt ( P ) F i j i d
If this ratio is below 1, the observations are better aligned with the identity matching. If the ration is above 1, the observations are better aligned if a permutation is applied. F i j opt and F i j i d refers to the alignment between subject i and subject j . The choice of pairwise matching rather than template matching is driven by two main factors: The computation of an average is non trivial in the case of node unlabelled graphs (see [START_REF] Kolaczyk | Averages of unlabeled networks: Geometric characterization and asymptotic behavior[END_REF]Calissano et al. (2023) for further details); The pairwise matching also offers an higher number of matching and permutation matrices, for a more detailed and robust understanding of the misalignment between subjects.

| Comparing different initialization

For all the parcellations and hemispheres, we are aligning all the graphs G i with all the other graphs G j , i , j = 1, . . . , 100 computing via FAQ the optimal permutation matrices Pij . The aim of this first experiment is to see how the initialization parameter changes the performance of the graph matching with respect to the identity matching.

| Comparing different permutations

Once the graph matching is performed, we are interested in understanding how the regions are matched by looking at the optimal permutation matrices Pij . We divide the type of permutations into three categories: Self Matching -when the parcel is permuted with itself; Neighbours Matching -when the parcel is permuted with a neighbouring parcel; and Others -when the permuted regions are not direct neighbors. In Figure 2, we provide a conceptual visualization of the type of permutation between regions.

F I G U R E 2 Types of regions matching: (Top) Self Matching representing regions which are permuted with themselves; (Middle) Neighbours Matching representing regions permuted with their spatial neighbours on the brain; (Bottom) Others representing permutations with other regions rather then themselves or neighbours.

| Spatial distribution of permutations

The matching results can be also illustrated on cortical surface to visualize the areas which are more often matched.

Specifically, for each region we computed its permutation rate, differentiating over the spatial distance of the proposed permutation. We denote r n p the rate of permutation with an n th neighbor. The symbol r 0 p therefore corresponds to the rate of "no permutation" (or self permutation), r 1 p the rate of permutation with a direct neighbor, etc. This visualization provides insights on the spatial misalignment across subjects and if certain regions are more permuted than other. It is also possible to visualise individual matching results by computing the rate of permutation of one region with all others across all subjects. This provides a region by region view of the stability of regions and the most likely permutations.

| RESULTS

| Comparing different initialization strategies

In Figure 3, we illustrate the distribution of the ratio for every parcellation and every initialization. The barycentric and random initialization are systematically outperformed by the original regions matching, indicating that it is preferable to perform no matching. Contrarily, the identity and the spatial initialization outperform the original regions matching for every atlas with a number of regions at or above 100, suggesting inter-subjects alignment. Furthermore, the improvement increases with the number of regions. For example, at 100 parcels with the barycenter initialization the percentage of permuted parcels is 21% and it doubles at 1000 parcels, reaching 58.5%. The results are in line with the expectations, as the most significant priors on the brain regions matching are the identity and the spatial rather then the random and the barycenter, as it is more plausible to have either none or only local misalignment between regions in the brain. The plot shows how the alignment with identity or spatial initialization can decrease the misalignment error inter-subjects. F I G U R E 3 For each parcellation {P 50, ..., P 1000}, we compute all the subject-to-subject matching. We plot the distribution of the ratio between the graph matching objective function and the identity matching objective function. When the ratio is above 1, the graph matching is preferable over the identity matching. Different algorithm initializations are considered (Spatial Adjacency, Identity, Barycenter and Random).

| Comparing different permutations

The percentage of regions belonging to Self Matching, Neighbours Matching, and Others is illustrated in Figure 4.

When the graph matching was outperformed on average by the identity matching, indicating a poor matching, we report the results with a dashed line. For all four initialization strategies the percentage of permuted nodes increases as the number of regions increases. For the identity and the spatial initializations, the percentage of Self Matched nodes is higher then the percentage of matched nodes. Among the matched nodes, the percentage of nodes matched with Neighbours is predominant with respect to regions permuted with others. These are plausible results, as we expect the regions to be permuted with neighbours due to local misalignment during the parcellation procedure. In addition, the identity initialization show an higher percentage of self permuted regions while the spatial initialization an higher percentage of neighbours matched regions.

| Spatial organisation of the permutations

To understand the spatial distribution of the different permutations (self, other, neighbours), we represent the permutations obtained with the Adjacency Initialization on the cortex. We chose this initialization as it is the best performing one according to the previous results. Figure 5 illustrates the permutation rates r 0 p , r 1 p , and r 2 p , for atlases with a 100, 300, and 1000 regions. In agreement with the results of Figure 4, the majority of permutation are local, with their rate quickly decreasing with increased pair distance. No spatial organisation is noticeable, indicating that all regions have the same permutation distribution. Examples of permutation for specific regions are illustrated in Figure 6 for similarly located regions across the atlases with 100, 300, and 1000 regions. We observe a high locality of permutation, with the rate of permutation rapidly decreasing with increased distance from the region of interest. As the number P50 100 0 P200 P400 P600 P800 P1000 F I G U R E 4 Type of permutations. For each parcellation, the proportion of self permuted parcels, parcels permuted with neighbours and parcels permuted with other parcels. The lines are dotted when the graph matching is under-performing the identity matching. By looking at both the Identity Initialization and the Spatial Adjacency Initialization we notice an increase in the number of permuted parcels when the number of parcels increase, maintaining an higher proportion of parcels permuted with neighbours.

F I G U R E 5 Surface visualization of the permutation type. For each regions, we considered the rate at which a regions was permuted with any other region (i.e. not a self permutation), with a region that is not a 1st degree neighbor, and with a region that is not a 1st or 2nd degree neighbor. Warmer and colder colors indicate a higher and lower permutation rate, respectively. Gray indicates a permutation rate below 1%. The three rows correspond to a brain parcellation with 100, 300, 1000 regions. of regions of the atlas is increased, the permutations become more local and stable around the region of interest.

The increased locality for atlases with more parcels is expected as the number of parcels increases the parcel size is reduced, leading to permutations with short distances. This is essentially a consequence of the neighborhood being defined via parcel adjacency instead of cortical distance.

| DISCUSSION

Studies representing the brain as a network often assume a correspondence between the nodes across the subjects.

However, due to inter-subject variability affecting the structural and functional organisation of the cortex, the labels of the nodes may not be matched across subjects, leading to incorrect comparisons. In this work, we proposed a network alignment strategy that exploits the spatial organisation of the cortical surface and thus promotes local permutation of node labels. When applied to the structural connectivity matrix of healthy young adults, our proposed alignment reduced inter-subject variance.

F I G U R E 6 Given the brain region marked with the black dot, we plot the frequency of permutation with other regions for three atlas granularity (100, 300, 1000 regions). The frequencies are computed considering the whole set of subject-to-subject permutations obtained with the spatial initialization. The considered region is permuted with gray regions less than 1% of the times and with the other regions as indicated by the legend: warmer colors indicate high frequency of permutation, colder colors indicate low frequency of permutation.

In more detail, we considered the structural brain networks of 100 healthy subjects from the HCP, obtained from a state-of-the-art diffusion MRI tractography pipeline. The nodes of the network corresponded to the regions of a groupwise structural parcellation of the cortex with the number of regions varied from 50 to 1000. Among the different available network alignment techniques, we selected the FAQ matching as it allowed us to include priors on the spatial organisation of the cortex. We choose a spatial adjacency matrix as initialization of the algorithm, encoding the spatial relations of the parcels. Compared to other initialization strategies, using the adjacency matrix as a starting point systematically increased the similarity between subjects for atlases with more than 50 regions, as illustrated in Figure 3. On the other hand, using a random or barycentric initialisation lead to a decrease of inter-subject similarity.

Note that our baseline for similarity is the identity matrix, which corresponds to the assumption that the atlas is matched across subject made in current studies. Our results suggest that, for a healthy cohort, network alignment should be performed before any comparison of the networks across subjects. The whole framework is implemented in a notebook in the geomstats python package1 .

The importance of exploiting the spatial organisation of the cortex when performing network alignment is further explored by looking at the individual permutations. Specifically, we considered the distance between permuted regions, in terms of the neighborhood. As illustrated in Figure 4, most permutation are either with the same label (no permutation r 0 p ) or with a direct neighbor (r 1 p ). As the number of regions increase, there is an increase in the rate of long range permutation (r n p with n > 1), reaching 21% for 1000 regions. This phenomenon can be explained by the size of the regions of the atlas which decreases as the number of parcel increases. In other words, r 1 p and r 2 p may correspond to the same cortical distance on two different atlases. These results suggest that a better notion of neighborhood could be estimated with the distance between cortical regions rather than their adjacency in terms of cortical mesh. Such approach would require the selection of an arbitrary distance threshold.

To understand the spatial distribution, the permutation frequencies for each region of the brain were illustrated in Figure 5. Overall, no specific patterns where identified and all regions presented a similar permutation distribution.

From the inter-subject variability point of view, this supports the idea that the variance is not region specific and is instead uniform over the cortex. As illustrated in Figure 6, the permutation patterns were also consistent over different resolutions of the atlas. Regions were in general permuted with close neighbors, again highlighting the importance of the neighborhood information in the graph matching problem. While our proposed strategy favors local permutations, it does not forbid long range ones and they were observed in some cases. Because we considered structural connectivity, the permutation of these nodes can indicate that they fulfil similar roles in terms of structural network organisation, or at least that they have similar connectivity to the rest of the brain. This could occur, for example, in tightly interconnected subnetworks with few connections to other subnetworks.

Previous studies [START_REF] Tian | A deep learning-based multisite neuroimage harmonization framework established with a traveling-subject dataset[END_REF]; [START_REF] Chamberland | On the origin of individual functional connectivity variability: the role of white matter architecture[END_REF] have reported that interindividual variability is heterogeneous across the cortex, contrasting with the relatively homogeneous permutation maps of our results. Indeed, because our study relies on very similar data and processing, it is reasonable to expect our connectivity matrices to have the same variability profile. However, our goal in performing graph alignment is to improve the global similarity or the matrices and our approach does not seek to specifically correct for regions with high variability. In addition, the permutations identified by the graph alignment do not all have the same impact on the variability. Some permutations may greatly reduce it (e.g. permutations in regions of high variability) while others will have a small impact. Our results in Figure 3 indicate that the variability is indeed systematically reduced when using the spatial initialization (a ratio above 1 corresponds to reduced variability). However, due to noise, approximation in the atlas, and general model mismatch, not all variability can be compensated by the graph alignment and some interindividual subject will remain and may still be heterogeneous. In brief, the maps representing the permutation rate in Figure 5 should not be interpreted as interindividual variability as the relationship between the two is not trivial. Finally, our approach could be optimized to promote not only neighborhood permutations, but also favor regions of high variability susceptible to correction by permutation.

In summary, our results show that brain network alignment increases the inter-subject similarity of structural connectivity, considered through the structural connectome. Notice that the underlying uncertainty might be due to different factors, including the group-wise parcellation misalignment, and the subjects specific differences. A comparison of the registration error in the groupwise parcellation and the graph matching results might decouple the source of uncertainty. In the context of healthy cohorts, the reduction in variability is advantageous as we expect all subjects to have normal, and therefore similar, connectivity. Extending these results to pathology is challenging as the intersubject variance may be a feature of the disease and should not be corrected explicitly. A possible solution would be to apply network alignment using a different alignment targets per class, as is common in registration. Finally, our results considered only structural connectivity and may not generalize to functional connectivity or scenarios where both are considered simultaneously, e.g. structure-function mapping. In such cases, the alignment strategy should be modified to consider all available data to ensure a consistent permutation across the different networks.

| CONCLUSION

Considering a set of structural connectivity matrices of healthy subjects, this study shows how the cross-subjects variability can be decreased by applying graph matching techniques. We compare cross-subjects similarity using different initialization of the FAQ graph matching algorithm. The improvement is significant when the algorithm is initialized using the identity or spatial structure of the regions and when the number of parcels increases. By looking at the results, the suggested matching are mostly between neighbours, showing local misalignments due to groupwise parcellation strategy. These preliminary results highlight the potential of using graph matching prior to group connectivity studies. Applying graph matching between subjects using graph matching technique improve the crosssubjects similarity. F I G U R E 8 Right Hemisphere: Visualization of the ration between objective functions for each parcellations and each subjects.

F

  I G U R E 1 Left: Conceptual visualization of identity matching and graph matching in the context of structural brain connectivity. Right: Conceptual visualization of graph embedding in the context of identity matching and graph matching.

  Left Hemisphere: Visualization of the ration between objective functions for each parcellations and each subjects.

https://github.com/geomstats/geomstats

Funding information European Research Council (ERC) (grant agreement G-Statistics No. 786854)

Conflict of Interest

There are no relevant financial or non-financial competing interests to report.

The data that support the findings of this study are openly available from the Human Commectome Project at https://www.humanconnectome.org/.


Hugo Nicolas Irenav

Paolo Perali

Matthieu Sacher

Patrick Bot
email: patrick.bot@ecole-navale.fr.

Boundary Element Method Analysis of 3D Effects and Free-Surface Proximity on Hydrofoil Lift and Drag Coefficients in Varied Operating Conditions

Keywords: hydrofoil, free surface, aspect ratio, potential flow, boundary element method (BEM)

The use of hydrofoils to enhance ship performance raises the scientific issue of free-surface proximity, which is important to consider during the design stage, to feed velocity prediction programs, for instance. Typically, the flow over a shallowly submerged hydrofoil is characterized by the Froude number, the submergence depth-to-chord ratio, the angle of attack, and geometric parameters of the lifting surface. Among these parameters, the present paper investigates the influence of the wing aspect ratio on the lift and drag coefficients of hydrofoils operating near a free surface. For this purpose, rectangular wings with an H105 profile at 2 • angle of attack and aspect ratios ranging from 4 to 20 are systematically analyzed using a 3D boundary element method. The free surface is modeled using a linearized Neumann-Kelvin boundary condition. Chord-based Froude numbers of 0.5, 1.1, and 6.3 are studied. The submersion depth is swept between 0.1 and 30 times the foil chord length. The evolution of the normalized lift and drag coefficients with respect to the foil submersion and the aspect ratio is discussed in detail. Flow velocity is shown to play a significant role in the evolution of the lift and drag coefficients with submersion depth, close to the free surface, for all the aspect ratios. Its influence gets reduced by moving away from the free surface. The critical submersion depth, where the free-surface effects cease, is found to increase with higher flow velocity and aspect ratio. Furthermore, both positive and negative correlations between the force coefficients and the aspect ratio are identified, depending on the operating conditions. It is found that when the proximity to the free surface either enhances or impairs a force coefficient relative to its value in unbounded flow, increasing the aspect ratio amplifies this effect. Overall, this study confirms the effectiveness of steady boundary element methods for simulating the flow around hydrofoil wings in the vicinity of a free surface and contributes to further understanding the influence of geometric parameters on hydrofoil performance.

NOMENCLATURE c

Hydrofoil chord length

[m] C D Drag coefficient [-] C D∞ Drag coefficient in unbounded flow [-] C L Lift coefficient [-] C L∞ Lift coefficient in unbounded flow [-] F r
Froude number [-] F r c

Chord-based Froude number (F r c = U 0 / √ gc) [-] g Gravitational acceleration [m s -2

INTRODUCTION

Hydrofoils enhance the performance of marine vehicles, including sailboats and commercial ships, by reducing hydrodynamic resistance and improving seakeeping. Whether partially submerged or fully submerged, these appendages often operate close to the water's surface. Their proximity to the free surface is a critical factor to consider in the design stage [START_REF] Bal | High-Speed Submerged and Surface Piercing Cavitating Hydrofoils, Including Tandem Case[END_REF], as recently demonstrated through the implementation of a velocity prediction program (VPP) for the America's Cup AC50 and AC75 yachts [START_REF] Patterson | Development of a Six Degree of Freedom Velocity Prediction Program for the Foiling America's Cup Vessels[END_REF]). Yet, the complex physics involved is challenging to capture. The tools available for addressing the free-surface problem consist of simulations of varying fidelity levels, including methods such as lifting-line theory (LLT), vortex lattice method (VLM), panel method, and Reynolds-averaged Navier-Stokes (RANS) solvers, as well as experimental approaches. When it comes to simulations, the choice of fidelity level often results from a trade-off between computational efficiency and model accuracy.

Previous studies have observed that the proximity to the free surface alters the flow field compared to unbounded flow, which can significantly impact hydrofoil performance [START_REF] Zhu | Dynamics of a Three-Dimensional Oscillating Foil Near the Free Surface[END_REF]. Specifically, the low pressure on the suction side of a lifting surface acts to deform the free surface, resulting in the generation of transverse waves whose characteristics are dictated by the flow regime [START_REF] Parkin | Pressure Distribution on a Hydrofoil Running near the Water Surface[END_REF]. [START_REF] Parkin | Pressure Distribution on a Hydrofoil Running near the Water Surface[END_REF] further illustrated that the presence of a free surface leads to two distinct flow regimes, corresponding to low and high chord-based Froude numbers (F r c ). Therefore, the Froude number (F r) emerges as a dominant parameter in the problem, alongside the submergence depth-to-chord ratio (h/c), and other parameters such as the angle of attack (α) [START_REF] Kwag | Numerical Simulation of Free-Surface Flows 3-D Submerged Hydrofoil by N-S Solver[END_REF].

The aforementioned transition in the wave patterns between low and high Froude numbers has been extensively documented in the scientific literature, primarily on 2D profiles [START_REF] Bai | A Localized Finite-Element Method for Two-Dimensional Steady Potential Flows with a Free Surface[END_REF][START_REF] Yasko | Boundary Element Method for a Hydrofoil near the Free Surface[END_REF][START_REF] Chen | A Vortex Based Panel Method for Potential Flow Simulation around a Hydrofoil[END_REF][START_REF] Carabineanu | The Study of the Potential Flow Past a Submerged Hydrofoil by the Complex Boundary Element Method[END_REF]. At high Froude numbers, a smooth transverse wave forms above the foil, significantly diminishing lift as the body nears the free surface [START_REF] Daskovsky | The Hydrofoil in Surface Proximity, Theory and Experiment[END_REF]. In contrast, low Froude numbers disrupt the free surface by gravity waves and hydraulic jumps [START_REF] Parkin | Pressure Distribution on a Hydrofoil Running near the Water Surface[END_REF]. In this flow regime, the flow behavior depends on both the submergence depth-to-chord ratio and the Froude number. In close proximity to the free surface, a standing wave crest can manifest slightly ahead of the foil, with another crest further aft along the chord, resulting in a decrease in lift [START_REF] Laitone | Limiting Pressure on Hydrofoils at Small Submergence Depths[END_REF][START_REF] Pernod | Free-Surface Effects on Two-Dimensional Hydrofoils by RANS-VOF Simulations[END_REF]. Increasing the submersion depth increases the wavelength [START_REF] Karim | Numerical Simulation of Free Surface Water Wave for the Flow around NACA 0015 Hydrofoil Using the Volume of Fluid (VOF) Method[END_REF]. [START_REF] Uddin | Application of Volume of Fluid (VOF) Method for Prediction of Wave Generated by Flow around Cambered Hydrofoil[END_REF] have also demonstrated that the wave amplitude is influenced by the submergence depth, with increased steepness as the free surface is approached, potentially leading to wave-breaking [START_REF] Duncan | The Breaking and Non-Breaking Wave Resistance of a Two-Dimensional Hydrofoil[END_REF]. In cases of wave-breaking, unsteady breaking waves introduce fluctuations in both lift and drag [START_REF] Prasad | Numerical Simulation of Free Surface Flows around Shallowly Submerged Hydrofoil by OpenFOAM[END_REF]. By further reducing the Froude number, the free surface behaves akin to a rigid wall [START_REF] Faltinsen | Hydrodynamics of High-Speed Marine Vehicles[END_REF]. In such instances, flow acceleration occurs above the suction side of the hydrofoil, analogous to a ground effect [START_REF] Giesing | Potential Flow about Two-Dimensional Hydrofoils[END_REF]. Consequently, hydrofoil lift increases until a critical value of h/c is reached, above which it rapidly diminishes when approaching too close to the free surface [START_REF] Filippas | Hydrodynamic Analysis of Flapping-Foil Thrusters Operating beneath the Free Surface and in Waves[END_REF].

The variation of the lift coefficient (C L ) with the Froude number is characteristic of the challenge posed by the free-surface proximity. The impact of the alterations in the wave patterns between low and high Froude numbers on the lift coefficient was demonstrated by [START_REF] Yeung | A Hybrid Integral-Equation Method for Steady Two-Dimensional Ship Waves[END_REF] through a hybrid integral-equation method. The analyses conducted by [START_REF] Yeung | A Hybrid Integral-Equation Method for Steady Two-Dimensional Ship Waves[END_REF] have since served as benchmarks for validating numerical models for 2D and 3D hydrofoils operating beneath a free surface [START_REF] Bal | A Potential Based Panel Method for 2-D Hydrofoils[END_REF][START_REF] Bal | Numerical Analysis of 2-D and 3-D Cavitating Hydrofoils Under a Free Surface[END_REF]. [START_REF] Kennell | A Second-Order Theory for the Potential Flow About Thin Hydrofoils[END_REF] additionally noted that the lift coefficient attains its maximum value when operating at a low Froude number in proximity to h/c = 1, as observed in their case. Furthermore, [START_REF] Hough | Froude Number Effects on Two-Dimensional Hydrofoils[END_REF] reported that the lift coefficient remains nearly constant above a chord-based Froude number of 10. Consequently, [START_REF] Hough | Froude Number Effects on Two-Dimensional Hydrofoils[END_REF] designated the range of interest for high-speed operating conditions as

F r c ∈ [3 -10].
Concerning drag, a peak in the drag coefficient is observed in close proximity to the free surface, with respect to h/c [START_REF] Bal | Numerical Analysis of 2-D and 3-D Cavitating Hydrofoils Under a Free Surface[END_REF]. This peak's occurrence is independent of the flow regime, but its magnitude decreases with increasing flow velocity [START_REF] Xie | Performance Analysis of 3D Hydrofoil under Free Surface[END_REF]. The significance of drag mitigation due to flow velocity is underscored by Plotkin's observation that the drag coefficient (C D ) experiences significant damping at low Froude numbers [START_REF] Plotkin | The Thin-Hydrofoil Thickness Problem Including Leading-Edge Corrections[END_REF]. Ultimately, the influence of the free-surface proximity on the lift and drag coefficients diminishes above a specific submergence threshold [START_REF] Ali | Numerical Study of Free Surface Effect on the Flow Around Shallowly Submerged Hydrofoil[END_REF]. Above this threshold, the submersion depth can be regarded as effectively infinite [START_REF] Hoque | Simulation of Water Wave Generated by Shallowly Submerged Asymmetric Hydrofoil[END_REF].

Regarding the impact of the angle of attack on the hydrodynamic characteristics of shallowly submerged hydrofoils, the lift coefficient follows a linear progression with the angle of attack (until reaching stall), similar to the behavior observed in unbounded flow [START_REF] Wadlin | A Theoretical and Experimental Investigation of the Lift and Drag Characteristics of Hydrofoils at Subcritical and Supercritical Speeds[END_REF][START_REF] Binns | The Effect of Heel Angle and Free-Surface Proximity on the Performance and Strut Wake of a Moth Sailing Dinghy Rudder T-foil[END_REF]. [START_REF] Ni | Performance of a Hydrofoil Operating Close to a Free Surface over a Range of Angles of Attack[END_REF] demonstrated that motion-induced surface waves can enhance hydrofoil performance by delaying stall, particularly at high angles of attack. Moreover, the evolution of the drag coefficient exhibits a quadratic relationship with respect to the lift coefficient, as discussed by [START_REF] Binns | The Effect of Heel Angle and Free-Surface Proximity on the Performance and Strut Wake of a Moth Sailing Dinghy Rudder T-foil[END_REF]. [START_REF] Binns | The Effect of Heel Angle and Free-Surface Proximity on the Performance and Strut Wake of a Moth Sailing Dinghy Rudder T-foil[END_REF] also revealed that increasing the heel angle does not impair the lift-to-drag ratio of a full T-foil and highlighted ventilation inception under extreme sailing conditions, characterized by α, h/c, and F r c .

As for hydrofoil geometry, early studies on its impact on the hydrodynamic characteristics employed a lifting-line approach [START_REF] Michel | Hydrofoil Handbook Volume II -Hydrodynamic Characteristics of Components[END_REF]. More recently, an iterative BEM developed by [START_REF] Bal | A BEM for the Prediction of Free Surface Effects on Cavitating Hydrofoils[END_REF] was used for investigating the effects of curved wing tips and taper ratio on the performance of 3D shallowly submerged cavitating hydrofoils under various operating conditions at chordbased Froude numbers up to 1.6 and 0.8, respectively (Bal, 2023a;[START_REF] Bal | The Taper Ratio Influence on the Performance of 3-D Cavitating Hydrofoils Moving under a Free Surface[END_REF]. These studies have highlighted the considerable influence of hydrofoil geometry in the context of free-surface proximity.

Among the typical parameters describing the geometry of lifting surfaces, there remains limited knowledge concerning the influence of the wing aspect ratio (Λ) on the performance of hydrofoils near a free surface, as no thorough investigation into its influence over a broad range of varied operating conditions was conducted. [START_REF] Thiart | Vortex Lattice Method for a Straight Hydrofoil Near a Free Surface[END_REF] developed a VLM for finite-span hydrofoil wings operating beneath a free surface to explore the effects of aspect ratio, taper ratio, sweep angle, and dihedral angle on hydrofoil lift and drag. [START_REF] Thiart | Generalized Vortex Lattice Method for Prediction of Hydrofoil Characteristics[END_REF] focused on rectangular wings with aspect ratios up to 10 under fixed operating conditions, with the exception of the chord-based Froude number, which ranged from 0 to 20. In a separate study, a panel method was employed to compute the lift coefficient for foils with aspect ratios of 4, 5, and 6 at low Froude numbers and a submergence depth-to-chord ratio of 1 [START_REF] Xie | Performance Analysis of 3D Hydrofoil under Free Surface[END_REF]. Using a similar level of fidelity, [START_REF] Ozdemir | Wave Field Generated by Finite-Span Hydrofoils Operating beneath a Free Surface[END_REF] assessed the wave field and hydrodynamic performance of shallowly submerged hydrofoils with aspect ratios ranging from 4 to 8, operating at low and moderate Froude numbers and with h/c equal to 0.75, 1, and 1.5. The comparisons between 2D and 3D cases revealed that, although performance tends to improve with increasing aspect ratio, the behavior of shallowly submerged 3D hydrofoils can differ significantly from that of their 2D counterparts, depending on the Froude number.

Therefore, a gap remains in the understanding of the effects of the aspect ratio on the performance of submerged hydrofoils, particularly at high Froude numbers (F r c > 3) commonly encountered in practical foiling scenarios. Thus, this paper introduces a method to explore how the wing aspect ratio and the submersion depth affect the lift and drag coefficients of hydrofoils operating at low, moderate, and high Froude numbers. Section 2 outlines the methodology and establishes the numerical framework, which consists of modeling and meshing wing geometries with different aspect ratios, and solving the fluid problem for different Froude numbers and submersion depths, using a 3D BEM approach. Section 3 details the numerical setup. In Section 4, the numerical framework is applied to the parametric case studies. The lift and drag coefficients computed for each wing, under each operating condition, are analyzed with respect to submersion depth and aspect ratio.

NUMERICAL FRAMEWORK

Problem Formulation

Figure 1 provides a descriptive representation of the flow around a shallowly submerged hydrofoil. In Figure 1, U 0 is the flow velocity in the x-direction and η is the free-surface elevation. The submersion depth is defined at the mid-chord of the centerline, coinciding with the center of rotation of the angle of attack. The submergence depth-to-chord ratio is defined between this point and the free surface at rest, z = 0.

Geometry Modeling and Meshing

The appendages were generated using LEModeler (https://github.com/hugo-nicolas/LEModeler), an open-source code for modeling parametric hydrofoil geometries developed for the purpose of this study, using the Rhinoceros® CAD engine. The geometries were meshed using structured grids refined at both the leading and trailing edges. To capture the 3D effects, the spanwise panel distribution was refined to leave smaller square cells at the foil tips.

Hydrodynamic Model

The fluid problem was addressed using PUFFIn (https://www.ensta-bretagne.fr/fr/optifoil), a 3D potential flow solver based on a BEM approach. The solver PUFFIn is detailed in [START_REF] Perali | Comparaison de Conditions de Surface Libre Linéarisées Pour l'étude d'un Hydrofoil Submergé à l'aide d'une Approche Potentielle[END_REF] and the theoretical framework for steady simulations used in the present work is summarized here.

The fluid flow is considered incompressible and irrotational, and viscous effects are neglected. The velocity field (u) is the superposition of the undisturbed velocity (U 0 ) with the gradient of the perturbation potential (φ). It satisfies Laplace's equation,

∇ 2 φ = 0. (1)
Applying the second Green identity to this equation yields the boundary integral equation [START_REF] Katz | Low-Speed Aerodynamics[END_REF],

φ = - 1 4π S ï σ 1 r -µ n • ∇ Å 1 r ãò dS, ( 2
)
where σ is the source strength and µ is the doublet strength such that

σ = - ∂φ ∂n , µ = -φ. (3)
S is the domain boundary, which typically consists of the body surface (S B ), the wake surface (S W ), and the free-surface boundary (S F S ). n is the outward pointing vector normal to S. Under the assumption that the wake should not support hydrodynamic loads, the source distribution on the wake surface is null. Consequently, Equation 2 can be rewritten as

φ = - 1 4π S B +S F S ï σ 1 r -µ n • ∇ Å 1 r ãò dS + 1 4π S W ï µ n • ∇ Å 1 r ãò dS, (4)
with the source strength given by the non-penetration condition on the hydrofoil,

σ = -U 0 • n. ( 5
)
PUFFIn models the free surface using either a linearized Neumann-Kelvin, a symmetry condition, or an antisymmetry boundary condition. The symmetry and antisymmetry boundary conditions result from the asymptotic behaviors of the Froude number F r -→ 0 and F r -→ ∞, respectively [START_REF] Newman | Marine Hydrodynamics[END_REF]. In the scope of this paper, the simulations were run using the Neumann-Kelvin formulation linearized at z = 0,

(U 0 • ∇) 2 φ + g ∂φ ∂z = 0, ( 6
)
where g is the gravitational acceleration.

To solve the problem numerically, the boundaries are discretized using quadrilateral elements. The distribution of sources and doublets is assumed constant on each element. Writing Equation 4 at the geometric center of each element yields the general matrix form [START_REF] Filippas | Hydrodynamic Analysis of Flapping-Foil Thrusters Operating beneath the Free Surface and in Waves[END_REF])

A(η)   µ B σ F S   = B(η)   σ B µ F S   + W (η) µ W , ( 7
)
where the free-surface elevation is given by

η = - 1 g U 0 • ∇φ. (8)
σ B and µ B denote the strengths of sources and doublets on the body, respectively. Similarly, σ F S and µ F S represent the strengths of sources and doublets on the free surface. µ W corresponds to the doublet strengths on the wake panels. σ B is determined by Equation 5 while the other strength vectors are unknown. A(η), B(η), and W (η) can be computed for constant-strength singularity panels. Considering the linearized free-surface boundary condition 6, the problem is solved by inverting A(η) once at the beginning of the computation, using a lower-upper decomposition.

A Lagrangian approach is employed to construct the wake surface during the computation. This process is iterative in time, meaning that the position of a wake panel at the current time step is determined based on the velocity values from the previous time step. Since the flow is considered inviscid, the potential on a wake panel convected by the flow remains constant. Thus, µ W only depends on the strengths of the doublets on the panels located downstream of the trailing edge, which are computed by imposing pressure equality between the suction side and the pressure side, known as the Kutta condition, at the trailing edge of the hydrofoil.

The pressure distribution on the body is yielded by the Bernoulli relation. From there, the hydrodynamic forces and moments acting on the hydrofoil are computed by integrating the pressure over the body surface. The simulations end when a residual criterion on the last N iterations is met.

Remarks on the Hydrodynamic Model

The numerical model applies to steady simulations. As noted by [START_REF] Prasad | Numerical Simulation of Free Surface Flows around Shallowly Submerged Hydrofoil by OpenFOAM[END_REF], wave-breaking may occur at low submergence, leading to unsteadiness and fluctuations in both lift and drag. However, since PUFFIn considers inviscid and irrotational flow, it does not capture this phenomenon. Additionally, unsteadiness in the wake of the hydrofoil can arise in close proximity to the free surface, particularly at low Froude numbers. Extremely close proximity to the free surface can result in significant changes in wave steepness, affecting hydrofoil immersion and making the convergence of the simulations more challenging. Most importantly, the model does not consider multiphase flow phenomena like cavitation or ventilation, which limits the ability to draw quantitative conclusions at very low immersion levels. Therefore, this paper does not provide conclusions for these extreme operating conditions. Nevertheless, despite these limitations, the results presented in Section 4 offer qualitative insights into the evolution of the lift and drag coefficients in the vicinity of the free surface, aligning with observations in the existing literature.

NUMERICAL SETUP

Numerical Parameters

The numerical framework was applied to hydrofoil wings with an aspect ratio of 4, 6, 10, and 20, operating at different Froude numbers and submersion depths. In particular, three chord-based Froude numbers equal to 0.5, 1.1, and 6.3 were analyzed. They respectively correspond to a low, moderate, and high Froude number. The lifting surfaces were modeled rectangular to avoid accumulating effects induced by more complex geometries and rather focus on the influence of the aspect ratio. Their chord length was fixed to 0.125 m. Consequently, each Froude number corresponds to a given flow velocity, and the changes in the aspect ratio were achieved by adjusting the wing span. The angle of attack remained constant and equal to +2 • . The cross-section profile was an H105 profile [START_REF] Speer | The BASILICUS Project -Return of the Cruising Hydrofoil Sailboat[END_REF], which is close to most sections used in foil design. Thus, the zero-lift angle of attack was -2.35 • . Concerning submersion depth, a wide range of submergence depth-to-chord ratios between 0.1 and 30 were swept. Water density was assumed constant and equal to 1025 kg/m 3 . Finally, an infinite depth of water was considered to avoid the effects of shallow water. The operating conditions investigated in this paper are summarized in Table 1.

Table 1. Operating conditions investigated.

F r c

[-] U 0 [m s -1 ] α [ • ] h/c [-] 0.5 0.6 2 0.1-30 1.1 1.2 2 0.1-30 6.3 7.0 2 0.1-30
The size of the free surface domain was determined by the foil's chord length. Upstream of the foil was left a distance of 5 chords while downstream of the foil was left 40 chord lengths. Port and starboard, the distance to the foil tips was 20 times the chord length such that the hydrofoil was centered in the transverse direction of the flow. The residual criterion for stopping the hydrodynamic simulations was defined such that the maximum deviation on the forces is 0.005 N for the last N = 20 iterations.

Mesh Convergence Study

The mesh topology was first established on the hydrofoil with an aspect ratio of 4 and further adapted to the lifting surfaces of higher aspect ratios while keeping the same mesh topology as Λ = 4. For this reason, the mesh convergence analysis was conducted on the wing with an aspect ratio of 4. The convergence study was conducted at a Froude number of 6.3 and a submergence depth-to-chord ratio of 20. The results of the mesh convergence study on the lift coefficient are shown in Figure 2.

To balance cost and accuracy, the mesh of 1280 panels (16 panels spanwise and 40 along the chord on each side) was selected as the reference. Figure 3 shows a 3D view of the reference geometry.

RESULTS

Influence of the Free-Surface Proximity

The lift and drag coefficients computed for each wing, under each operating condition in Table 1, were normalized by their respective value in unbounded flow (C L∞ and C D∞ ). The results of these simulations are presented in Figure 4, with Figure 4 In Figure 4 (a), (b), and (c), the lift coefficient plummets when approaching too close to the free surface, no matter the flow velocity. As discussed in Section 2.4, the values of the force coefficients in close proximity to the free surface (typically h/c < 0.5 here) must be taken with caution since unsteadiness may occur. Nevertheless, the significant reduction in the lift coefficient observed in the present work is in agreement with the literature, giving confidence in the use of steady BEM to provide an overview of hydrofoil performance near the free surface [START_REF] Chen | A Vortex Based Panel Method for Potential Flow Simulation around a Hydrofoil[END_REF].

Figure 4 (a), (b), and (c) illustrates that the characteristic evolution of the lift coefficient between low and high Froude numbers is well captured by the simulations. Notably, a portion of the normalized lift coefficient exceeds 1 for the chord-based Froude numbers 0.5 and 1.1. The occurrence of such a gain in the lift coefficient was also identified in other studies at low Froude numbers [START_REF] Pernod | Free-Surface Effects on Two-Dimensional Hydrofoils by RANS-VOF Simulations[END_REF]. Figure 4 (a) exhibits that, in the present work, this gain in the lift coefficient occurs near h/c = 1 at the chord-based Froude number of 0.5. It enhances the lift coefficient by up to 15.8 % compared to its value in unbounded flow for the wing with an aspect ratio of 20, which suggests a potential benefit of operating in the vicinity of the free surface for maximizing the lift coefficient, under specific operating conditions. Figure 4 (b) illustrates that the amplitude of the gain is smaller at F r c = 1.1 than at F r c = 0.5. The gain also occurs at a greater h/c compared to F r c = 0.5, close to h/c = 3 for the larger wings. At the chord-based Froude number of 6.3, it is made clear from Figure 4 (c) that the gain in the lift coefficient has disappeared completely. Therefore, the proximity to the free surface systematically reduces the lift coefficient at the higher Froude number, which is consistent with other studies when inertial forces largely dominate [START_REF] Daskovsky | The Hydrofoil in Surface Proximity, Theory and Experiment[END_REF]. for h/c < 0.5. The occurrence of the peak preceded by the reduction in the drag coefficient is delayed to greater submergence by increasing the flow velocity. Consequently, Figure 4 (f) shows that the reduction in the drag coefficient common to F r c = 0.5 and F r c = 1.1 has disappeared at F r c = 6.3. Furthermore, the amplitude of the peak is reduced with increasing speed, as shown by [START_REF] Xie | Performance Analysis of 3D Hydrofoil under Free Surface[END_REF].

Figure 4 (a), (b), and (c) illustrates that no matter the aspect ratio, the lift coefficient systematically converges to its value in unbounded flow by increasing h/c. Similar conclusions apply to the drag coefficient based on Figure 4 (d), (e), and (f). In particular, the disappearance of the free-surface effects with increasing h/c is achieved from a threshold value of h/c, whose existence is largely observed in the literature [START_REF] Hoque | Simulation of Water Wave Generated by Shallowly Submerged Asymmetric Hydrofoil[END_REF]. Above this threshold, increasing the aspect ratio increases the lift coefficient and decreases the drag coefficient according to the known effect of the aspect ratio in unbounded flow. Comparing the convergence of the normalized lift and drag coefficients with respect to h/c at each flow velocity reveals that the threshold h/c for achieving behavior akin to infinite submersion depth increases as the chord-based Froude number increases. It is also greater with a larger aspect ratio.

Influence of the Aspect Ratio

The aspect ratio demonstrates variable effects on the normalized force coefficients in Figure 4, depending on the influence of the free-surface effects on the lift and drag coefficients. These variations include both positive and negative correlations between the aspect ratio and the normalized force coefficients. The characterization of these correlations consists in the main finding of the present work and is discussed below.

Figure 4 (a), (b), and (c) illustrates that increasing the aspect ratio reduces the normalized lift coefficient when C L /C L∞ 1 (negative correlation) and increases it when C L /C L∞ 1 (positive correlation). The same conclusion may be drawn from Figure 4 (d), (e), and (f) for the drag coefficient. Therefore, when the proximity to the free surface either enhances or impairs a force coefficient relative to its value in unbounded flow, increasing the aspect ratio amplifies the effect. This, however, does not appear systematic based on the numerical simulations when the normalized lift and drag coefficients reach a value of 1 under the influence of the free surface.

To provide a clearer representation of the effects of the aspect ratio when the normalized force coefficients approach a value of 1 near the free surface, Figure 5 presents the data from Figure 4 focusing on submergence depth-to-chord ratios up to 4. Figure 5 illustrates the existence of ranges of operating conditions where the correlation between the normalized force coefficients and the aspect ratio is unclear, such as for the lift coefficient at h/c between 0.25 and 2 and F r c = 0.5. In such instances, the correlation is considered transient. The existence of ranges with transient correlation may result from numerical uncertainties, leading to variations in the results, notably at low Froude numbers.

In light of the changes in the evolution of the normalized lift and drag coefficients versus h/c discussed in Section 4.1, the submersion depth at which the normalized lift and drag coefficients reach a value of 1 increases with the chord-based Froude number. Consequently, the influence of the aspect ratio on the normalized force coefficients is characterized by both the Froude number and the submergence depth-to-chord ratio. 2 (b) respectively summarize the ranges of operating conditions where the correlation between the aspect ratio and the normalized lift and drag coefficients is positive ("+"), negative ("-"), or transient ("-→"). Table 2 helped select the test cases shown in Figure 6 to analyze quantitatively the effects of the aspect ratio on the lift and drag coefficients in varied operating conditions. F r c h/c 0.1 0.25 0.0 1 0.0 1.5 0.0 2 0.0 0.0 0.0 3 . . . 30

0.5 - | -→ | + 1.1 - | -→ | + 6.3 - (b)
F r c h/c 0.1 0.25 0.0 1 0.0 1.5 0.0 2 0.0 0.0 0.0 3 . . . 30

0.5 + | -→ | - 1.1 + | -→ | - 6.3 +
Correlation between a given normalized force coefficient and the aspect ratio:

negative + positive -→ transient submersion depths shown in Figure 6 were selected based on Table 2. Figure 6 (a) and (d) corresponds to h/c = 1, where a gain in the normalized lift coefficient was observed in Figure 4 (a). In Figure 6 (a), a negative correlation between the normalized lift coefficient and the aspect ratio can be observed for the chord-based Froude numbers of 1.1 and 6.3. Conversely, a positive correlation between the normalized drag coefficient and the aspect ratio can be observed in Figure 6 (d) for the three Froude numbers. As indicated in Table 2 (a), the correlation between the aspect ratio and the lift coefficient is transient at h/c = 1 and F r c = 0.5, which justifies the non-monotonic evolution of C L /C L∞ versus Λ for F r c = 0.5 in Figure 6 (a). The submergence depth-to-chord ratio of 3 was selected as the intermediate case, where the effects of the free surface remain significant yet reduced. At h/c = 3, a wing with a larger aspect ratio provides a greater lift coefficient and a smaller drag coefficient, for the chord-based Froude numbers of 0.5 and 1.1. This is shown in Figure 6 (b) and (e), respectively through the positive correlation of the normalized lift coefficient with the aspect ratio and the negative correlation of the drag coefficient with the aspect ratio, at F r c = 0.5 and F r c = 1.1. The correlations with the aspect ratio are respectively negative and positive for the normalized lift and drag coefficients, at F r c = 6.3. Finally, Figure 6 (c) and (f) shows the evolution of the normalized lift and drag coefficients when only a small proportion of the free-surface effects persists. The correlations between the aspect ratio and the force coefficients at h/c = 5 are consistent with those observed at h/c = 3.

The significant disparity in the lift and drag coefficients at h/c = 1 illustrated in Figure 6 (a) and (d) aligns with the observations made in Section 4.1 regarding the significant impact of the Froude number near the free surface. Comparing the evolution of the lift and drag coefficient at each submergence depth-to-chord ratio in Figure 6 reveals that the curves C L /C L∞ and C D /C D∞ for F r c = 0.5 and F r c = 1.1 tend to superimpose by increasing h/c, for all the aspect ratios analyzed in the present work. Thus, it may be concluded that, for all the aspect ratios, the effects of flow velocity on the lift and drag coefficients decrease and disappear, eventually, by moving away from the free surface.

CONCLUSIONS

The present paper investigates the effects of both the free surface and the aspect ratio on the hydrofoil lift and drag coefficients. To do so, simulations on fully submerged rectangular wings of variable aspect ratios under different flow velocities and submersion depths are performed using a 3D boundary

Figure 1 .

 1 Figure 1. Schematic description of the fluid problem.

Figure 2 .Figure 3 .

 23 Figure 2. Mesh convergence on the lift coefficient for the wing with an aspect ratio of 4. The chord-based Froude number is 6.3 and the submergence depth-to-chord ratio equals 20.

  (a), (b), and (c) showing the normalized lift coefficient, and Figure 4 (d), (e), and (f) the normalized drag coefficient with respect to the submergence depth-to-chord ratio. The trends observed in Figure 4 are contextualized within the existing literature.

Figure 4 Figure 4 .

 44 Figure 4 (d) shows that, at the chord-based Froude number of 0.5, the normalized drag coefficient first decreases before peaking by approaching the free surface. A similar behavior may be observed at the intermediate Froude number in Figure4(e), with a noticeable reduction in the drag coefficient

Figure 6 Figure 6 .

 66 Figure 6 illustrates how the lift and drag coefficients, normalized by their values in unbounded flow, evolve with the aspect ratio for different submersion depths and chord-based Froude numbers. The

Table 2 (

 2 a) and Table

Table 2 . Summary of the correlations between the wing aspect ratio and the normalized (a) lift coefficient and (b) drag coefficient in varied operating conditions.

 2

(a)

Downloaded from http://onepetro.org/JST/article-pdf/8/01/183/3331105/sname-jst-2023-10.pdf/1 by guest on 04 December 2023

ACKNOWLEDGEMENTS

The authors are especially grateful to Galen W. Ng for providing valuable feedback on the manuscript draft.

element method with a linearized Neumann-Kelvin boundary condition for the free surface. The difference between low and high Froude numbers is well captured by the flow solver, which demonstrates the applicability of steady panel methods to the problem of fully submerged hydrofoils of variable aspect ratios operating near a free surface.

Several conclusions may be drawn from the present work:

• The threshold submersion depth above which the free-surface effects cease is delayed with both flow velocity and aspect ratio. • Increasing the aspect ratio can enhance or impair the lift and drag coefficients, depending on the Froude number and the submergence depth-to-chord ratio. When the proximity to the free surface either enhances or impairs a force coefficient relative to its value in unbounded flow, increasing the aspect ratio amplifies this effect. • At low submergence levels, flow velocity significantly influences the lift and drag coefficients, for all the aspect ratios. Its influence diminishes as immersion depth increases, eventually reaching a point where it becomes negligible.

To provide a comprehensive understanding of the variations in the lift and drag coefficients discussed in this paper, additional analyses of local quantities (such as the pressure coefficients, the free-surface elevation, and the effective angle of attack) are essential. Moreover, investigating the influence of downwash and the transition phenomenon identified in this study should contribute to offering more meaningful and physical interpretations of the results.

Overall, this study contributes to enhancing the understanding of how geometric parameters impact the performance of closely submerged hydrofoils, in view of performing multi-fidelity optimizations, for example.