Varun Warrier

Claire S Leblond

Freddy Cliquet

Aims-2-Trials Leap

Thomas Bourgeron

Simon Baron-Cohen

Polygenic scores for intelligence, educational attainment and schizophrenia are differentially associated with core autism features, IQ, and adaptive behaviour in autistic individuals

à 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.

ADI-R (e.g. for ADI restricted and repetitive behaviour: Beta = 0.06, 95%CI = 0.09 -0.02). Most of these associations were also significant when restricting it to males only or to individuals with IQ > 70. We find limited evidence for heterogeneity between cohorts.

Conclusion and relevance

We identify specific and different patterns of association between PGS for the three phenotypes and core and associated autism features. This provides greater resolution into the shared genetics between autism and the three phenotypes, and suggests one method to investigate heterogeneity in autism and co-occurring conditions.

Introduction

Autism is a set of complex and heterogeneous neurodevelopmental conditions characterised by difficulties in social communication, unusually repetitive and restricted behaviour and interests, difficulties adjusting to unexpected change, and altered sensory processing (1,2). Autism also frequently co-occurs with other medical, neurodevelopmental, or mental health conditions (3,4). Autism is substantially heritable, with twin and familial heritabilities between 50-90% [START_REF] Colvert | Heritability of autism spectrum disorder in a UK population-based twin sample[END_REF][START_REF] Tick | Heritability of autism spectrum disorders: A meta-analysis of twin studies[END_REF][START_REF] Sandin | The familial risk of autism[END_REF][START_REF] Bai | Association of Genetic and Environmental Factors With Autism in a 5-Country Cohort[END_REF]. Significant progress has been made in identifying common genetic variants [START_REF] Grove | Identification of common genetic risk variants for autism spectrum disorder[END_REF] and genes enriched for rare and protein truncating variants associated with autism [START_REF] Sanders | Insights into Autism Spectrum Disorder genomic architecture and biology from 71 risk loci[END_REF][START_REF] Satterstrom | Novel genes for autism implicate both excitatory and inhibitory cell lineages in risk[END_REF][START_REF] Yuen | Whole-genome sequencing of quartet families with autism spectrum disorder[END_REF]. Additionally, genome-wide association studies (GWAS) have identified shared genetics between autism and a number of other phenotypes, including schizophrenia (r g = 0.24±0.04), educational attainment (r g = 0.20±0.02), and intelligence (r g = 0.22±0.03) [START_REF] Grove | Identification of common genetic risk variants for autism spectrum disorder[END_REF][START_REF] Weiner | Polygenic transmission disequilibrium confirms that common and rare variation act additively to create risk for autism spectrum disorders[END_REF].

Autism is unlikely to be a single condition but a composite of multiple, partly dissociable core features (e.g. social interaction and communication difficulties, and restricted interested and repetitive behaviour) and associated features like difficulties in adaptive behaviour and a wide IQ range [START_REF] Waterhouse | ASD Validity[END_REF][START_REF] Shuster | Review of Factor Analytic Studies Examining Symptoms of Autism Spectrum Disorders[END_REF][START_REF] Mandy | Research Review: What is the association between the social-communication element of autism and repetitive interests, behaviours and activities?[END_REF][START_REF] Frazier | Validation of Proposed DSM-5 Criteria for Autism Spectrum Disorder[END_REF][START_REF] Warrier | Social and non-social autism symptoms and trait domains are genetically dissociable[END_REF][START_REF] Greenberg | Testing the Empathizing-Systemizing theory of sex differences and the Extreme Male Brain theory of autism in half a million people[END_REF]. These core and associated features contribute to heterogeneity in autism and are associated with differential outcomes in autistic individuals [START_REF] Pickles | The adult outcome of children referred for autism: typology and prediction from childhood[END_REF][START_REF] Anderson | Predicting young adult outcome among more and less cognitively able individuals with autism spectrum disorders[END_REF]. Therefore, investigating the genetics of autism as a single phenotype may provide incomplete information about prognosis and outcomes in autistic individuals.

While schizophrenia, educational attainment, and intelligence are all genetically correlated with autism, we still do not know which core and associated features of autism these phenotypes are genetically correlated with. In other words, are polygenic scores (PGS) for these phenotypes associated with all core and associated autism features, or are they associated with only some of the core and associated features? Finally, how do other sources of heterogeneity such as participants' IQ or sex affect these associations? This can provide better resolution of the shared genetics between various core and associated features and other co-occurring conditions, and identify axes to potentially stratify autistic individuals along to provide better care and support.

To address these questions, we investigated the association of PGS for three wellpowered GWAS (intelligence, educational attainment and schizophrenia) with 12 different core and associated autism features in 2,512 -3,681 autistic individuals from three cohorts. These 12 features represent the core aspects of autism such as restricted and repetitive behaviour and interests, and social interaction and communication difficulties, and associated features such as measures of IQ and adaptive behaviour. We additionally conduct sensitivity analyses by stratifying for males and individuals with IQ > 70, and investigate differences in PGS associations between cohorts.

Methods Overview of the cohorts

We conducted PGS analyses in three cohorts of autistic individuals: The Simons Simplex Collection (SSC) [START_REF] Fischbach | The Simons Simplex Collection: A Resource for Identification of Autism Genetic Risk Factors[END_REF], the Autism Genetic Resource Exchange (AGRE, CHOP sample) (N max = 1,200) [START_REF] Geschwind | The autism genetic resource exchange: a resource for the study of autism and related neuropsychiatric conditions[END_REF], and the AIMS-2-TRIALS LEAP sample (N max = 262) [START_REF] Charman | The EU-AIMS Longitudinal European Autism Project (LEAP): clinical characterisation[END_REF]. Participants in all three cohorts were recruited across multiple sites. For each cohort, we restricted our analyses to autistic individuals who passed the genetic quality control and who had the relevant phenotype, resulting in a total sample size ranging from 2,512 -3,681. We compare the demographic characteristics of these cohorts in Table 1. Further details of these cohorts are provided in the eMethods 1.

Insert Table 1 here Genetic quality control, imputation, and polygenic score generation Full details of the genetic quality control are provided in eMethods2. Briefly, quality control was conducted for each cohort separately, by array. We excluded participants with genotyping rate < 95%, excessive heterozygosity (±3 standard deviations from the mean), who were not of European ancestry as identified using multidimensional scaling, had mismatched genetic and reported sex, and families with Mendelian errors > 10%. SNPs were excluded with genotyping rate < 10%, or deviated from Hardy-Weinberg equilibrium (pvalue < 1x10 -6 ). Imputation was conducted using Michigan Imputation Server (25) using the 1000 genomes Phase 3 v5 as the reference panel (26) (for AGRE and SSC) or using HRC r1.1 2016 reference panel (for AIMS-2-TRIALS). SNPs with minor allele frequency < 1% or an imputation r 2 < 0.6 were excluded.

GWAS used to generate PGS were identified a priori using statistical power analyses (27) using https://eagenetics.shinyapps.io/power_website/. Details of parameters used are provided in eTable 1. We assumed a SNP heritability of 0.15 for the core and associated autism features and a sample size of 3,000. We conducted power analyses for seven relatively well-powered GWAS that are genetically correlated with autism. The GWAS for autism [START_REF] Grove | Identification of common genetic risk variants for autism spectrum disorder[END_REF] was excluded as this overlaps with the SSC and AGRE datasets. We subsequently restricted our analyses to only GWAS that had around 80% statistical power at a genetic correlation of 0.25 or greater than 60% power at a genetic correlation of 0.2, which is similar to the genetic correlation between autism and educational attainment, intelligence and schizophrenia. Three phenotypes met our criteria: 1. Educational attainment (28) (N = 766,345, excluding the 23andMe dataset), 2. Intelligence (29) (N = 269,867), and 3. Schizophrenia (30) (56,418 cases and 78,818 controls). All three GWAS are genetically correlated with autism, and PGS for these GWAS are over-transmitted from parents to autistic children [START_REF] Weiner | Polygenic transmission disequilibrium confirms that common and rare variation act additively to create risk for autism spectrum disorders[END_REF], suggesting that the observed genetic correlations with autism are not due to other confounding factors such as ascertainment bias.

PGS were generated for three phenotypes using PRSice2 (31), using independent SNPs (clumping r 2 of 0.1 and 250 kb) present in both the GWAS and the testing datasets.

Given the limited testing sample size and the number of tests conducted, we chose a priori pvalue thresholds that explained maximum variance in the three phenotypes: p ≤ 1 for the educational attainment and intelligence GWASs [START_REF] Lee | Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals[END_REF][START_REF] Savage | Genome-wide association meta-analysis in 269,867 individuals identifies new genetic and functional links to intelligence[END_REF], and p ≤ 0.1 [START_REF] Lam | Comparative genetic architectures of schizophrenia in East Asian and European populations[END_REF] for the schizophrenia GWAS. The total number of SNPs included are provided in eTable 2.

Autism core and associated features

We identified 12 autism core and associated features that are: 1. Widely used in studies related to autism; 2. Are a combination of parent-, self-, other-report and performance-based measures to investigate if reporter status affects the PGS association, 3. Collected in all three cohorts; and 4. Cover a range of core and associated features in autism. Details of each feature is provided in the eMethods 3. The core features are:

1 4. Non-verbal IQ (NVIQ) * Verbal and non-verbal communication scores were not separately calculated in the AIMS2-TRIALS-LEAP dataset. We thus used a total communication score and meta-analysed it with the verbal and non-verbal communication scores in the other two cohorts.

Statistical analyses

For each cohort, PGS were regressed against the autism features with sex and the first 15 genetic principal components as covariates in all analyses, with all variables standardised (mean = 0, standard deviation = 1). Additionally, age of participants in the AIMS-2-TRIALS LEAP cohort, and genotype array in the SSC cohort were included as covariates. Age at completion of each measure was only available in the AIMS-2-TRIALS LEAP dataset. IQ type and ADOS module were included as additional covariates. For SSC and AIMS-2-TRIALS LEAP as only unrelated individuals were included, a linear model was used. AGRE includes a large number of multiplex families with multiple autistic individuals per family. Thus, to retain all participants, we used linear mixed effects regression with family ID included as a random intercept to account for relatedness between individuals.

Effect sizes were meta-analysed across the three cohorts using inverse variance weighted meta-analyses with the following formula:

ݓ ൌ ܧܵ/1 ଶ ܧܵ ௧ ൌ ඥ1/Σ ݓ ܽݐ݁ܤ ௧ ൌ Σ ߚ ݓ /Σ ݓ
Where ߚ is the standardized regression coefficient of the PGS, ܧܵ is the associated standard error, and ݓ is the weight. p-values were calculated from Z scores ܽݐ݁ܤ( ௧ ܧܵ/ ௧ ሻ.

Variance explained was estimated (R 2 ) by squaring the ܽݐ݁ܤ ௧ . Given the high correlation between the autism features and phenotypes, we used Benjamini-Hochberg False Discovery Rates to correct for multiple testing (q-value < 0.05).

Core and associated autism features are quantitatively and qualitatively different between sexes [START_REF] Lai | Sex/gender differences and autism: setting the scene for future research[END_REF]. Further, autistic individuals with co-occurring ID are enriched for rare and de novo protein-truncating variants (PTV) [START_REF] Satterstrom | Large-Scale Exome Sequencing Study Implicates Both Developmental and Functional Changes in the Neurobiology of Autism[END_REF][START_REF] Sanders | Insights into Autism Spectrum Disorder Genomic Architecture and Biology from 71 Risk Loci Article Insights into Autism Spectrum Disorder Genomic Architecture and Biology from 71 Risk Loci[END_REF], suggesting that some of the variation in these phenotypes may be attributed to this class of genetic variants in these individuals [START_REF] Bishop | Identification of Developmental and Behavioral Markers Associated With Genetic Abnormalities in Autism Spectrum Disorder[END_REF][START_REF] Lord | Damaging de novo mutations diminish motor skills in children on the autism spectrum[END_REF]. So, we conducted sensitivity analyses by restricting the samples to only males (N = 3,040 -2,125) or individuals with IQ > 70 (N = 2,319 -1,801). We did not test for association between FSIQ and VIQ and the three PGS in the AGRE in the two sensitivity analyses, as fewer than 20 individuals had scores on these measures. We followed the same statistical analysis pipeline for the two sensitivity analyses, and identified significant associations using Benjamini-Hochberg FDR correction (q-value < 0.05).

To understand how the features correlated with each other, we conducted Pearson's correlation analyses followed by hierarchical clustering for each cohort (eMethods 4).

All analyses were conducted in R. Data and software availability is provided in eMethods 5.

Ethics

We received ethical approval to access and analyse de-identified genetic and phenotypic data from the three cohorts from the University of Cambridge Human Biology Research Ethics Committee.

Results

Phenotypic correlation in the three cohorts

Pearson's correlation produced a similar correlation profile across all three cohorts with two clusters -a core autism cluster with all core autism features, and an associated feature cluster with all three IQ measures and VABS (Figure 1, eFigures 1 -3, eTables 3 -5). The two clusters were largely negatively correlated with each other. Within the autism cluster, hierarchical clustering suggested that instrument type had primacy over domain. For instance, the ADOS subscales clustered closely together, and as did the two parent-report measures (SRS and RBS).

Insert Figure 1 here Correlation analyses of the correlation coefficients across pairs of features identified high correlations for all three pairs (SSC-AGRE: r = 0.89, p-value < 2x10 -16 ; SSC-AIMS-2-TRIALS LEAP: r = 0.92, p-value < 2x10 -16 ; and AGRE -AIMS-2-TRIALS LEAP: r = 0.92, p-value < 2x10 -16 ) (Figure 2).

Insert Figure 2 here Polygenic score analyses PGS for intelligence were positively associated with all measures of IQ and VABS with concordant effect direction in all cohorts: FSIQ (Beta = 0.08±0.02, q-value = 6.69x10 -4 , R 2 = 0.63%), VIQ (Beta = 0.06±0.02, q-value = 5.21x10 -3 , R 2 = 0.37%), NVIQ (Beta = 0.09±0.02, q-value = 3.82x10 -5 , R 2 = 0.77%), and VABS (Beta = 0.06±0.02, q-value = 1.91x10 -2 ) (Figure 3

and eTable 6).

Similarly, PGS for educational attainment were positively associated with all measures of IQ and VABS with concordant effect direction in all cohorts: FSIQ (Beta = 0.06±0.02, q-value = 2.36x10 -3 , R 2 = 0.30%), VIQ (Beta = 0.06±0.01, q-value = 6.69x10 -4 , R 2 = 0.32%), NVIQ (Beta = 0.05±0.01, p-value = 7.57x10 -3 , R 2 = 0.21%), and VABS (Beta = 0.04±0.02, q-value = 4.53x10 -2 , R 2 = 0.20%) (Figure 3 and eTable 6). There were no statistically significant differences in effect sizes between the PGS for education and intelligence for all measures of IQ.

Interestingly, PGS for educational attainment were also positively associated with scores on the verbal and non-verbal communication subscales of the ADI, though the effect directions were not concordant in the three cohorts (ADI NVC: Beta = 0.05±0.02, q-value = 3.23x10 -2 , R 2 = 0.26% ADI VC: Beta = 0.05±0.02, q-value = 2.69x10 -2 , R 2 = 0.27%) (Figure 3 andeTable 6). Higher scores on the ADI are associated with more difficulties with communication.

PGS for schizophrenia were positively associated with two ADI subscales (ADI RRB: Beta = 0.06±0.02, q-value = 5.21x10 -3 , R 2 = 0.39%, ADI NVC: Beta = 0.05±0.02, q-value = 3.23x10 -2 , R 2 = 0.32%) (Figure 3 and eTable 6). The effect direction was concordant in all three cohorts for ADI RRB but not for ADI NVC (eTable 6).

Insert Figure 3 here

Assessing heterogeneity across cohorts and polygenic scores

We next conducted binomial sign test to investigate effect direction concordance between the three cohorts. Formal tests of heterogeneity are likely to be biased given the small number of studies included [START_REF] Hippel | The heterogeneity statistic I(2) can be biased in small metaanalyses[END_REF], so we assessed heterogeneity using sign tests assessing effect directions. 32 out of 36 tests had concordant effect direction between the SSC and the AGRE cohorts (p-value = 1.94x10 -6 , binomial sign test), 21 out of 30 between the SSC and AIMS-2-TRIALS (p-value = 0.04), and 22 out of 30 between the AGRE and the AIMS-2-TRIALS dataset (p-value = 0.01), suggesting largely similar effect directions across cohorts.

We also investigated the concordance of effect direction between the three sets of PGS for the 12 features tested using the meta-analysed effect direction. All 12 tests had concordant effect direction between the PGS for educational attainment and intelligence (pvalue = 4.8x10 -4 , binomial sign test), but only 6 of 12 has concordant effect direction between PGS for schizophrenia and both educational attainment and intelligence (p-value = 1). This is expected given the high genetic correlation between intelligence and educational attainment (r g = 0.72±0.011), but the low genetic correlation between schizophrenia and intelligence (r g = -0.21±0.025) and educational attainment (r g = 0.01±0.017).

Sensitivity analyses

As the AIMS-2-TRIAL LEAP dataset significantly differed from the AGRE and the SSC datasets in terms of fraction of female autistic or autistic individuals with IQ < 70, we conducted two sensitivity analyses (Methods). In the males-only analyses, PGS for intelligence were associated with FSIQ and NVIQ (eFigure 4 and eTable 7). PGS for educational attainment were associated with the two ADI communication subscales, RBS, FSIQ, and VIQ. Finally, PGS for SCZ were associated with ADI: RRB.

In the IQ > 70 analyses, PGS for intelligence and educational attainment were associated with FSIQ, VIQ, and NVIQ. Notably, the variance explained in this subset of individuals for these three phenotypes were 2 -3 times higher than the variance explained in the primary analyse. Additionally, PGS for educational attainment were associated with scores on the VABS, RBS, and ADI: NVC. Finally, PGS for schizophrenia were significantly associated with ADI: RRB and NVIQ (eFigure 5 and eTable 8). In both sensitivity analyses, the effect directions were concordant for all significant associations.

Discussion

Given demographic differences between the three cohorts, we first compared the patterns of phenotypic correlations for the core and associated autism features in the three cohorts, and identified a similar pattern in all three datasets. The core features clustered together, and the measures of intelligence and adaptive behaviour clustered together. This was confirmed by high Pearson's correlation between pairs of features across the three cohorts.

PGS for the three phenotypes were differentially associated with the features tested. Specifically, PGS for educational attainment and intelligence were positively associated with all three measures of IQ and VABS. The shared genetics between VABS and educational attainment and intelligence is unsurprising, given the modest to substantial phenotypic positive correlations between VABS and measures of IQ in the three cohorts (r = 0.36 to 0.68). Sensitivity analyses in autistic individuals with IQ > 70 confirmed the association between the measures of IQ and PGS for educational attainment and intelligence. Previous research has suggested that protein truncating rare and de novo variants are negatively associated with measures of IQ and scores on the VABS in autistic individuals [START_REF] Sanders | Insights into Autism Spectrum Disorder genomic architecture and biology from 71 risk loci[END_REF][START_REF] Satterstrom | Large-Scale Exome Sequencing Study Implicates Both Developmental and Functional Changes in the Neurobiology of Autism[END_REF][START_REF] Robinson | Genetic risk for autism spectrum disorders and neuropsychiatric variation in the general population[END_REF], though these sets of variants are enriched in individuals with IQ < 70. In contrast, in our analyses, PGS for educational attainment and intelligence explained a greater proportion of variance in measures of IQ in autistic individuals with IQ > 70 compared to all autistic individuals, suggesting a more prominent role for PGS in measures of IQ in autistic individuals without intellectual disability. Taken together, our analyses suggest that in autistic individuals, similar to the general population [START_REF] Lee | Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals[END_REF][START_REF] Savage | Genome-wide association meta-analysis in 269,867 individuals identifies new genetic and functional links to intelligence[END_REF][START_REF] Ganna | Ultra-rare disruptive and damaging mutations influence educational attainment in the general population[END_REF], genetic variants across the allelic frequency spectrum contribute to differences in IQ and related measures, but the relative contribution of different sets of genetic variants may differ based on intellectual disability.

PGS for educational attainment were also associated with higher scores on the two ADI communication subdomains (verbal and non-verbal communication), which is in the opposite direction to what's been observed between social and communication difficulties and educational attainment in the general population (r g = -0.30± 0.11, P = 0.007) [START_REF] Warrier | Social and non-social autism symptoms and trait domains are genetically dissociable[END_REF]. Sensitivity analyses restricting to either males or individuals with IQ > 70 also identified significant associations between PGS for educational attainment and ADI communication subdomains with concordant effect direction across all three datasets. This inversion in effect direction between social and communication difficulties and educational attainment in the typical population and autistic individuals must be further investigated. In the two sensitivity analyses PGS for educational attainment were negatively associated with scores on the RBS, though this was only nominally significant in the non-stratified analyses. The variance explained by the PGS for educational attainment for RBS-R in the subgroups was 2 -2.5x the variance explained in the primary analyses, highlighting, yet again, heterogeneity within the autism spectrum based on IQ and sex.

PGS for schizophrenia, on the other hand, were positively associated with scores on ADI: RRB and ADI: NVC. However, only the association with ADI: RRB was significant in the two sensitivity analyses. In the general population, there is some evidence to suggest that PGS for schizophrenia are associated with social and communication difficulties [START_REF] St Pourcain | ASD and schizophrenia show distinct developmental profiles in common genetic overlap with population-based social communication difficulties[END_REF][START_REF] Riglin | Schizophrenia risk alleles and neurodevelopmental outcomes in childhood: a population-based cohort study[END_REF]. However, to our knowledge, this is the first report linking PGS for schizophrenia with restricted and repetitive behaviour, though RRBs are elevated in individuals with schizophrenia [START_REF] Evans | Development of Two Dimensional Measures of Restricted and Repetitive Behavior in Parents and Children[END_REF].

The sample size for each individual cohort was limited, and the number of cohorts were small excluding the possibility of conducting formal tests of heterogeneity between cohorts. However, analyses of effect direction of the three sets of PGS revealed broad similarities between cohorts, with notably high concordance in effect direction between the AGRE and the SSC cohorts. This is notable given the substantial differences in demographic characteristics between the cohorts, and potential differences in the SNPs used to generate the polygenic scores for the three cohorts.

Comparing the physical correlations to the effect direction of the PGS further helps to elucidate the differences between the three sets of PGS. Whilst PGS for SCZ is positively associated with the autism core features, and negatively with the associated features, the PGS for educational attainment and IQ were positively associated with most features tested. Notably, PGS for schizophrenia were significantly associated with some autism core features, PGS for intelligence were significantly associated with all associated features, whereas PGS for educational attainment were significantly associated with all associated features and some autism core features suggesting differential association. However, we note that the statistical significance of these associations will be influenced by sample size. In other words, larger sample sizes may identify other associations between the PGS and the features.

Our study also provides a cautionary note about instrument bias which is reflected in the PGS associations. For example, restricted and repetitive behaviours are primarily measured here using RBS, ADI-RRB and ADOS-RRB. The phenotypic correlations between these are low to modest (r : 0.11 -0.48). Further, SCZ PGS are only associated with ADI-RRB, and PGS for educational attainment are only associated with RBS (in the sensitivity analyses). Together this suggests that idiosyncratic and systematic differences (e.g. reporter bias, number of questions) makes it challenging to compare different core and associated features even if they are meant to capture the same underlying latent trait. Phenotypically, this manifests as modest correlations, and genetically this suggests different underlying genetic architectures for these features, both of which have been demonstrated before [START_REF] Mallard | Not just one p: Multivariate GWAS of psychiatric disorders and their cardinal symptoms reveal two dimensions of cross-cutting genetic liabilities[END_REF][START_REF] Knott | Living with ASD: how do children and their parents assess their difficulties with social interaction and understanding?[END_REF][START_REF] Johnson | Discrepancies between self-and parentperceptions of autistic traits and empathy in high functioning children and adolescents on the autism spectrum[END_REF].

Limitations

Whilst the largest study to date, we were statistically well powered to investigate the shared genetics only with three phenotypes, and were not well-powered to explore the covariance explained at multiple different p-value thresholds in the training datasets. Larger datasets may identify additional significant associations between the current PGS and features tested. Further, whilst well powered, PGS still do not account for all the SNP heritability. We also do not examine the role of rare and low-frequency variants in this study, which are also associated with autism and the three phenotypes [START_REF] Satterstrom | Large-Scale Exome Sequencing Study Implicates Both Developmental and Functional Changes in the Neurobiology of Autism[END_REF][START_REF] Ganna | Ultra-rare disruptive and damaging mutations influence educational attainment in the general population[END_REF][START_REF] Singh | Rare loss-offunction variants in SETD1A are associated with schizophrenia and developmental disorders[END_REF]. Finally, each many phenotypes included in the study (e.g. SRS and RBS-R) can be further divided into additional subdomains based on factor analyses or evaluating the items included. A comprehensive investigation across all subdomains is warranted, yet is only possible with additional datasets that provide sufficient statistical power.

Conclusions

Core and associated autism features are differentially associated with PGS for intelligence, educational attainment, and schizophrenia. We find limited evidence for heterogeneity between cohorts despite differences in demographic characteristic, though stratifying for males and IQ > 70 identified additional associations with PGS. Our study provides insights into the underlying heterogeneity in autism, and provides greater resolution to the shared genetics between autism and intelligence, educational attainment, and schizophrenia.

Figure 1 :

 1 Figure 1: Heatmap of the phenotypic correlations in the three cohorts

Figure 2 :

 2 Figure 2: Correlations of correlations among pairs of autism core and associated features across the three cohorts

Figure 3 :

 3 Figure 3: Polygenic score association of the three phenotypes with the autism core features, measures of IQ, and adaptive behaviour


Tomi Kinnunen
email: tomi.kinnunen@uef.fi

Andreas Nautsch
email: andreas.nautsch@eurecom.fr

Nicholas Evans

Xin Wang

Massimiliano Todisco

Héctor Delgado

Junichi Yamagishi

Aik Kong

Lee

Visualizing Classifier Adjacency Relations: A Case Study in Speaker Verification and Voice Anti-Spoofing

Keywords:

Whether it be for results summarization, or the analysis of classifier fusion, some means to compare different classifiers can often provide illuminating insight into their behaviour, (dis)similarity or complementarity. We propose a simple method to derive 2D representation from detection scores produced by an arbitrary set of binary classifiers in response to a common dataset. Based upon rank correlations, our method facilitates a visual comparison of classifiers with arbitrary scores and with close relation to receiver operating characteristic (ROC) and detection error trade-off (DET) analyses. While the approach is fully versatile and can be applied to any detection task, we demonstrate the method using scores produced by automatic speaker verification and voice anti-spoofing systems. The former are produced by a Gaussian mixture model system trained with VoxCeleb data whereas the latter stem from submissions to the ASVspoof 2019 challenge.

Introduction

Whether it be for challenge results summarization [START_REF] Nautsch | Asvspoof 2019: spoofing countermeasures for the detection of synthesized, converted and replayed speech[END_REF][START_REF] Lee | Two decades into speaker recognition evaluation -are we there yet?[END_REF][START_REF] Greenberg | Two decades of speaker recognition evaluation at the National Institute of Standards and Technology[END_REF][START_REF] Nagrani | Voxceleb: Large-scale speaker verification in the wild[END_REF] or the analysis of classifier fusion, some means to compare alternative classifiers can often provide illuminating insight into robustness, coherence, and generalisation [START_REF] Ramos | Likelihood ratio data to report the validation of a forensic fingerprint evaluation method[END_REF]; (dis)similarity, and complementarity. Depending on prior knowledge (and data) available, one could compare • the underlying working principles, features, architectures, training data properties -classifier metadata,

• empirical performance (e.g. accuracy) on a common set of evaluation data, to name two possibilities. Often we want to understand how the two are related -what is the impact of specific modeling choices upon performance. A common approach is to fix some parameters while varying others to find out their impact on performance. Whenever one implements the classifier and has full control over the experiment, this white-box approach is relatively straightforward. There are, however, situations where one may not have access to exhaustive classifier details and configuration settings but would still like to learn about similarities between alternative solutions. One such black-box scenario are public machine learning challenges and technology benchmarks. Whether it be NIST speaker recognition evaluations [START_REF] Lee | Two decades into speaker recognition evaluation -are we there yet?[END_REF][START_REF] Greenberg | Two decades of speaker recognition evaluation at the National Institute of Standards and Technology[END_REF], VoxCeleb [START_REF] Nagrani | Voxceleb: Large-scale speaker verification in the wild[END_REF], ASVspoof [START_REF] Nautsch | Asvspoof 2019: spoofing countermeasures for the detection of synthesized, converted and replayed speech[END_REF], or indeed any other competitive campaigns, participants typically run their in-house systems on a common evaluation set and submit predictions (e.g. scores) for unlabeled data along with a system description. Since source codes or models are often not required, there can be uncertainty as to why specific challenge entries outperform others. The authors' personal motivation for the presented work stems partially from difficulties encountered in our efforts to link classifier properties to their performance in a recent challenge [START_REF] Nautsch | Asvspoof 2019: spoofing countermeasures for the detection of synthesized, converted and replayed speech[END_REF]. We wanted to find out what can be learnt about classifier differences based on detection scores.

To this end, we propose a method for computing the distance between binary classifiers based on the scores produced for common evaluation data (Fig. 1). We impose a stronger notion of classifier (dis)similarity than the established methods of detection error trade-off (DET) [START_REF] Martin | The DET curve in assessment of detection task performance[END_REF] and receiver operating characteristics (ROC) analyses: to be considered identical, a pair of classifiers must agree not only in their DET profiles but relative ordering of individual trials. Similar principles are the basis for a number of statistical significance tests but our perspective is in visualizing classifiers beyond DET/ROC plots.

As Fig. 1 suggests, we conjecture that distances derived from scores may provide information on (possibly unknown) classifier metadata differences. After describing the background and methodology, we report an example application of our methodology to automatic speaker verification (ASV) and voice anti-spoofing. The former includes classifiers constructed by the authors (known classifier metadata) while the latter represents the less controlled case of submissions to the ASVspoof 2019 challenge. Note, however, that the proposed methods are not in any way specific or limited to the speech domain. To this end, we provide an open-source reference implementation 1 .

Trade-offs in Characterizing Classifiers

There are various ways to describe a classifier. The coarsest description might include the general model class -such as 'a linear classifier' or 'a deep neural network.' A more refined description might detail the features and their dimensionality, the type and number of layers, the software toolkit, loss function, The behavior of classifiers can be controlled using a limited number of parameters meaningful to the experimenter. Not all the details, however, may be available to an evaluator. We propose to visualize classifier adjacencies based on their scores to complement other traditional analyses. The adjacencies can be informative of the adjacency of the underlying control parameters.

or training corpus to name a few possibilities.

Depending on the target group (e.g. general public, scientists, students, software engineers) one may prefer different levels of detail. What is in common is that the descriptions are designed for communicating key principles rather than enabling perfect reproducibility. The description is a human-readable tuple, h. For instance, our automatic speaker verification (ASV) systems described below are characterized by triplets h = (TrData, F, G) where TrData ∈ {VoxCeleb, LibriSpeech} indicates the training corpus while F and G are integers that indicate the number of mel filters and Gaussians, respectively. Each choice (attribute value in h) impacts results but there is no expectation for such triplets to be a complete description; it is a description intended for a particular audience (here, Interspeech attendee).

At the other extreme, the most precise description is a packaged/trained model -a list of all parameters, θ ∈ R D . For instance, θ could represent all the weights of a deep neural network (DNN), D being potentially in the order of millions. If scoring is deterministic (as is usually the case), θ and the functional form g(X; θ) completely specify how an arbitrary input X maps to an output (predicted class label or score), X ↦ Ŷ = g(X; θ). This leads to perfect reproducibilityone obtains the same scores, and hence also results, every time. Unlike h, however, θ is not intended to be comprehensible to the human: eyeballing the numerical values of DNN weights (an exercise which could take a while) is not particularly helpful in understanding how the system works or why it may work/fail in specific cases. Another problem is incomparability of the parameter vectors of different models. They could have different dimensionalities and there may be no meaningful correspondence between dimensions.

One more level of description are the scores that can be used to derive detection error trade-off (DET) profiles [START_REF] Martin | The DET curve in assessment of detection task performance[END_REF] and metrics such as the equal error rate (EER). Unfortunately, being based on score distributions, DETs or EERs are not informative of classifier responses to individual trials. A common compromise is to define a limited number of evaluation conditionssubsets of trials with characteristics useful for diagnostics (e.g. 'telephony', 'short-duration', 'attack S10', 'females'). Performance breakdown by condition is useful in analysing trends.

In summary, there are different abstraction levels to characterize classifiers. We aim at learning about their differences in terms of human-targeted descriptions without the precise knowledge of the attribute values themselves. We conjecture that classifiers that are similar in their h descriptions may produce similar outputs when executed on common data. In this case, one can learn about classifier design differences from their scores. To this end, we propose a method to compute pairwise classifier distances for common evaluation data. The distances can then be used for different purposes, such as agglomerative clustering or visualisation, to help in demonstrating trends or diagnosing classifier complementarities.

Proposed Classifier Adjacency Visualizer

Assumed input data

Let X and Y denote some input space and (binary) labels, respectively. We use

D = {(X i , Y i )} N i=1 ∼ i.i.d.
p(X, Y ) to denote a labeled (supervised) evaluation set of N trials. In an ASV set-up, each X i corresponds to an enrollment-test pair and Y i ∈ {0, 1} indicates whether the speaker identities in the two utterances are the same (target trial) or different (nontarget trial). In voice anti-spoofing, in turn, each X i corresponds to a test file and Y i indicates whether X i represents bonafide audio or a spoofing attack (such as computer-generated speech).

We assume a set G = {g 1 , . . . , g M } of M binary classifiers (detectors), g j ∶ X → R, each of which assigns a numeric membership value (detection score) to each trial, with the convention that high scores indicate support towards the positive class (e.g. target speaker). Unless we construct the classifiers ourselves, we do not have access to the scoring functions {g j } but only to their outputs for the shared evaluation data D (in turn, in a challenge setting, only the organizer knows the ground truth values {Y i } while the participants produce predictions on the test trials {X i } blindly). We use s i,j = g j (X i ) to denote the response of the jth classifier to the ith trial. The data used for our modeling consists of the responses of all classifiers to all the test trials -an N × M matrix S = [s 1 , . . . , s M ] where

s j = [s 1,j , . . . , s N,j ]
T is a column vector of size N × 1.

Optionally, one may have additional trial metadata. In this case, we assume that each trial can be uniquely assigned to one member of a mutually exclusive set of groups (conditions), c ∈ {1, 2, . . . , C}. In our experiments, for instance, groups in the ASV application correspond to the specific pair of speaker identities in the trial. We denote the number of trials in group c by N c . As the groups are mutually exclusive, N = ∑ C c=1 N c . We use either the original N × M score matrix or the smaller

C × M matrix S ′ = [µ 1 , . . . , µ M ]
where each group is represented using its average score. Here, µ j = [µ 1,j , . . . , µ C,j ] T is a column vector of size C × 1 for the jth classifier.

Classifier Adjacency

To identify similar (and dissimilar) systems, we define a distance function between classifiers, indicated by D(i, j), where i, j ∈ {1, . . . , M }. It is computed either from S or S ′ . We first note that, since the numerical range of detection scores can be arbitrary, a direct comparison of scores (e.g. using Euclidean distance) is usually not meaningful. In the ASV field, 'correction' for the range of scores is typically addressed through calibration (e.g. [START_REF] Brümmer | The BOSARIS toolkit user guide: Theory, algorithms and code for binary classifier score processing[END_REF]) that converts arbitrary scores into calibrated log-likelihood ratios (LLRs). Here, however, since score calibration is not our main concern, we opt for more straightforward scale-and translation-invariant distance computation. To this end, we define classifier distance through order statistics as opposed to raw scores. This choice stems from the knowledge that any two classifiers which yield the same order of trial scores on common data will share the same detection error trade-off (DET) curve and equal error rate (EER) [8, p.81] 2 . We use Kendall's τ [START_REF] Kendall | A new measure of rank correlation[END_REF][START_REF]The treatment of ties in ranking problems[END_REF][START_REF] Knight | A computer method for calculating Kendall's tau with ungrouped data[END_REF] as the rank correlation measure between classifiers i and j:

τ (i, j) = N con -N dis (N con + N dis + T i ) × (N con + N dis + T j ) , (2)
where N con and N dis are, respectively, the number of concordant and discordant pairs between i and j. Further, T i is the number of ties only in i, and T j is the number of ties only in j.

Concordance of a pair of trials means that the sort order agrees across the classifiers. For instance, if the a th and b th trial scores are ordered s a,i < s b,i in classifier i, the trial pair (a, b) is cordordant with classifier j, if similarly s a,j < s b,j . Likewise, if (s a,i > s b,i ) AND (s a,j > s b,j ), pair (a, b) is again corcondant. If concordance is not satisfied, the pair (a, b) is discordant. A tie within either classifier means that there identical scores produced for different trials. With real-valued scores stored in float or double precision this is generally a rare case.

An intuitive feeling of ( 2) can be developed by first assuming there are no ties (T i = T j = 0). Now, if the two classifiers place all trials into the same sort order, N dis = 0 and τ (i, j) reaches the maximum value of 1. The other extreme is obtained when the sort order of trials of one classifier is the reverse of the other. This corresponds to N con = 0 and yields the minimum value of -1. Finally, τ = 0 indicates lack of statistical association between the two. To sum up, τ (i, j) takes values in [-1, 1] (similar to Pearson correlation) and can be considered as a degree of agreement in the sort order of trials. Whenever τ (i, j) = 1, the classifiers have identical DET or ROC profiles.

For the purposes of visualization, we map τ into distances as D(i, j) = 1 2 (1 -τ (i, j)). Hence, identical systems (in the sense of their Kendall τ ) map to a distance of 0 while reverse ordering yields a maximum distance of 1. Finally, visualizations of classifier adjacency relations are obtained using classical (non-metric) multidimensional scaling (MDS). Each classifier is represented by a point in 2D space so that betweenclassifier distances approximate those in the given distance matrix (here, the M × M matrix containing all pairwise distances).

Experimental set-up

We demonstrate the proposed methods with two different tasks: automatic speaker verification (ASV) and voice anti-spoofing. The former consists of classifiers constructed by the authors with controlled parameters. The latter consists of countermeasures submitted to the ASVspoof 2019 challenge. Even if the system descriptions are known, we cannot interact with these systems nor do we know all their implementation nuances. 2 The converse does not necessarily hold: having the same DET curve or EER does not imply that the trial rankings must be the same.

Automatic speaker verification

The analysis of ASV methods uses the publicly available Vox-Celeb dataset. We use the standard trial list of 40 test speakers in VoxCeleb1 [START_REF] Nagrani | VoxCeleb: a largescale speaker identification dataset[END_REF]. It consists of 18860 target trials and the same number of nontarget trials. We form the groups (see Section 3) from the cross-product of speaker identities in the trials. With 40 speakers, this yields (40 ⋅ 41)/2 = 820 groups.

We use Gaussian mixture model-universal background model (GMM-UBM) based ASV system with an MFCCs frontend. While there are better ASV methods, we opt for the computationally light GMM-UBM as we want to generate a large number of classifiers. The focus of this study is in classifier adjacency rather than performance.

Following standard practice, extracted MFCCs are augmented with delta features and processed with RASTA and utterance-level cepstral mean and variance normalization (CMVN). For MFCCs, we include filterbanks comprising 12, 16, 20 and 24 filters. For UBMs, we vary the number of components in the powers of 2 between 2 and 1024. For UBM training, we select two different datasets: VoxCeleb2 [START_REF] Nagrani | Voxceleb: Large-scale speaker verification in the wild[END_REF] and LibriSpeech [START_REF] Panayotov | Lib-riSpeech: an ASR corpus based on public domain audio books[END_REF]. In total, we have 80 ASV systems.

Anti-spoofing

The analysis of voice anti-spoofing methods, or spoofing countermeasures (CMs), uses ASVspoof 2019 challenge submission entries. Details of the dataset [START_REF] Wang | Asvspoof 2019: A largescale public database of synthesized, converted and replayed speech[END_REF] and challenge results [START_REF] Nautsch | Asvspoof 2019: spoofing countermeasures for the detection of synthesized, converted and replayed speech[END_REF] are provided elsewhere. Here we focus on aspects relevant for the novel classifier adjacency analysis.

The challenge data consists of two different scenarios. The evaluation set of the logical attack (LA) scenario contains 13 different text-to-speech or voice conversion attacks (labeled A07. . . A19). The physical attack (PA) scenario, in turn, consists of simulated replay attacks from 27 different environments and 9 replay configurations. Both sets also contain additional bona fide (human speech) utterances. The participants processed the corresponding audio files blindly to obtain detection scores. In the LA scenario, the 3-class evaluation set consists of 5 370 bonafide/target, 1 985 bonafide/non-target, and 63 882 spoofed trials. For the PA scenario, the corresponding numbers are 12 960 bonafide/target, 5 130 bonafide/non-target, and 116 640 spoofed trials. For the LA scenario, the groups (see Section 3) are formed from the 13 attacks plus the bonafide class (a total of 14 groups). For the PA scenario, we used all crosscombinations of 27 environments × (9 replay configurations + 1 bona fide) to define 270 groups.

Results

The proposed method is demonstrated in a controlled-classifier setting (ASV) and in a setting of unconstrained classifiers (CMs).

Speaker Verification

Fig. 2 illustrates ASV results with varying number of filters, mixtures, and training data. The indicated C min llr performance [START_REF] Brümmer | The BOSARIS toolkit user guide: Theory, algorithms and code for binary classifier score processing[END_REF][START_REF] Brümmer | Application-independent evaluation of speaker detection[END_REF] is demonstrated for the lowest and highest number of mixture components to provide nominal guidance only. For the top-5 systems, C min llr varies between 0.466 and 0.481 (for visual clarity, we do not show all the C min llr values in Fig. 2). Differences in classifier parameters expose observable trends. For 12 filters, decision boundaries resulting from different training data become closer for as few as 4 components. Similar trajectories are followed as the number of components increases (as does the discrimination performance). Then, while the LibriSpeech trajectory continues, the VoxCeleb trajectory makes a sudden jump into the area of the top classifiers (at 512/1024 components and 20/24 filters). Adjacent classifiers (to the top ones) have less filters but more components; the top classifiers use more filters but less components. Our visuals supplement conventional DET plots: they provide insights into how decision boundaries are affected by data and parameters.

Anti-spoofing

Figs. 3 and4 show visualisations for the LA and PA scenarios, respectively. Note that the location of points in each panel for each dataset (left-to-right) are identical. In each panel, circles indicate a CM, whereas diamonds signify top-5 CMs. Colours indicate performance and frontend/backend meta-data.

The plots give an indication of the diversity among systems. Differences between the top-5 systems are reasonably representative of the full set of systems. The diversity may imply that the top-performing single systems are complementary. This hypothesis is supported by fusion results discussed in [START_REF] Nautsch | Asvspoof 2019: spoofing countermeasures for the detection of synthesized, converted and replayed speech[END_REF] which shows that the combination of top-performing single LA and PA systems leads to improved performance. One can presume that the reason for the ineffective fusion strategies for PA primary (a) Raw/FE (b) Grouped/FE (c) Raw/BE (d) Grouped/BE Figure 4: Visualisation of CMs on ASVspoof 2019 PA task; single systems. Frontends (FEs) and backends (BEs) are compared for raw and grouped scores. Legend as in Fig. 3. system submissions might be due to lack of complementary of the systems involved in the ensemble.

Adjacency from raw vs. grouped scores

Finally, we address classifier adjacency through score averaging by groups and subsequent use of the proposed method. Figs. 2 and4 illustrate a comparison of the proposed method for raw and grouped scores. Classifier adjacency is viewed through an alternative lens. In the ASV task, a clear separation between the two training is apparent. The VoxCeleb-trained classifiers (comparatively in-domain) become more alike. By using raw scores, trajectories by number of components (design parameters) are more alike. Feature resolution (the number of filters) corresponds to clear trajectories. In the uncontrolled CM setting, frontends as well as backends are more alike after grouping: the separation of non-GMM and GMM classifiers becomes oblivious.

Originally, we had two motivations for the grouped approach. First, the highly-compressed per-group average was hypothesized to give potentially a stable high-level classifier signature. Second, there is a computational advantage for very large trial sets. However, one needs additional metadata (and domain knowledge) to define the trial grouping. As we see, raw vs. grouped scores produce different findings. We prefer not to recommend either variant without a further study. A related approach, convex optimisation of error trade-offs 3 , resonates with the Bayesian decision framework; see use of the ROC convex hull instead of the corresponding 'steppy' profile [START_REF] Brümmer | The BOSARIS toolkit user guide: Theory, algorithms and code for binary classifier score processing[END_REF].

Conclusions

We proposed a simple approach to visualise classifier adjacency on common data. An example application under controlled conditions show that the tool renders classifiers adjacent, when they are similar in metadata. In the absence of unified metadata, analyses for uncontrolled, challenge conditions prove more challenging. The proposed tool nonetheless reveals intriguing, new visual insights into classifier adjacency not provided by any existing tools.

Acknowledgements

This work was supported by a number of projects and funding sources: VoicePersonae, supported by the French Agence Nationale de la Recherche (ANR) and the Japan Science and Technology Agency (JST) with grant No. JPMJCR18A6; Academy of Finland (proj. 309629); Region Grand Est, France.

Figure 1 :

 1 Figure1: The behavior of classifiers can be controlled using a limited number of parameters meaningful to the experimenter. Not all the details, however, may be available to an evaluator. We propose to visualize classifier adjacencies based on their scores to complement other traditional analyses. The adjacencies can be informative of the adjacency of the underlying control parameters.

Figure 2 :

 2 Figure 2: Visualisation of ASV similarity on VoxCeleb1 test set. The ASV systems are toy classifiers to show case how the proposed visualisation method behaves in a controlled setting. Shown C min llr values serve nominal orientation.

Figure 3 :

 3 Figure 3: CM similarity on ASVspoof 2019 LA test set. The CM systems are single systems (no fusions) submitted by participants; an uncontrolled setting is showcased (raw scores).The ASVspoof 2019 primary metric is shown, the so-called min t-DCF[START_REF] Kinnunen | Tandem assessment of spoofing countermeasures and automatic speaker verification: Fundamentals[END_REF]. Colours show different frontend transforms & feature types and backend classifiers.

In preliminary experiments, we used the proposed method in combination with the pool adjacent violators algorithm[START_REF] Brümmer | The PAV algorithm optimizes binary proper scoring rules[END_REF]. System dependently, 37720 ASV scores are reduced to 61 to 98 groups of same rank: visual differences are minimal; we spare comparisons.


Nattawut Chaomuang
email: nattawut.chaomuang@irstea.fr

Denis Flick

A simplified heat transfer model of a closed refrigerated display cabinet

Keywords: Closed refrigerated display cabinet, Simplified model, Heat transfer, Temperature

A simplified heat transfer model for a closed refrigerated display cabinet was developed based on a zonal approach. The developed model allows the prediction of time-averaged air and load temperatures at various positions in the display cabinet. These predicted values were compared with measured air and load temperatures. Good agreement was found between them with the maximum difference of 0.5°C for every studied position and the overall mean absolute error of 0.2°C. Air and load temperatures vary with positions. High temperature was observed at the front (highest value at the front top) and low at the back (lowest value at the back bottom). Based on the simulation, the cabinet thermal performance can be enhanced by minimizing the size of door gaps.

INTRODUCTION

It is well known that air infiltration across an air curtain is one of the major drawbacks of open refrigerated display cabinets. It accounts for about 66-77% of total refrigeration loads [START_REF] Gaspar | Experimental analysis of the thermal entrainment factor of air curtains in vertical open display cabinets for different ambient air conditions[END_REF] and also causes temperature heterogeneity in the cabinets (Laguerre et al., 2012b). High temperature is mostly at the front areas and low temperature at the rear [START_REF] Evans | Temperature and energy performance of refrigerated retail display and commercial catering cabinets under test conditions[END_REF]. As temperature is one of the most important factors affecting food quality and safety, minimizing the air infiltration would prolong product shelf life. Installation of doors on the open cabinets is one of the simple and effective solutions for the air infiltration reduction. Many achievements were demonstrated by various studies: energy saving potential [START_REF] Faramarzi | Performance and energy impact of installing glass doors on an open vertical deli/dairy display case[END_REF][START_REF] Fricke | Energy use of doored and open vertical refrigerated display cases[END_REF], better temperature homogeneity [START_REF] Chaomuang | Experimental analysis of heat transfer and airflow in a closed refrigerated display cabinet[END_REF][START_REF] Lindberg | Vertical display cabinets without and with doors-a comparison of measurements in a laboratory and in a supermarket[END_REF] and slower quality decay (Atilio de [START_REF] Atilio De Frias | Improving spinach quality and reducing energy costs by retrofitting retail open refrigerated cases with doors[END_REF]. As a result, the use of display cabinets with doors is increasingly used in retail stores worldwide. Closed display cabinets will account for about 75% of all operating display cabinets by 2020 in France [START_REF] Rpf | Fermeture des meubles réfrigérés: du positif pour les clients[END_REF]. Because door installation can modify airflow pattern and temperature distribution inside the cabinet [START_REF] Faramarzi | Performance and energy impact of installing glass doors on an open vertical deli/dairy display case[END_REF], an understanding of mechanism of heat transfer and airflow in the equipment becomes necessary to ensure food quality and safety. Although numerical simulation with CFD (Computational Fluid Dynamics) models is widely used to investigate the thermal phenomena in many refrigerating equipment [START_REF] Smale | A review of numerical models of airflow in refrigerated food applications[END_REF], this approach entails some limitations because of the requirement of calculation time and computational resources. A simplified heat transfer model based on a zonal approach becomes a complementary option as demonstrated by Laguerre et al. (2012a). They applied such a model to assess the influence of operating conditions and equipment design on temperature of products stored in an open display cabinet under steady state. This work on a closed display cabinet is a continuity of the previous studies of which the objective is to present a simplified heat transfer model based on the zonal approach. This developed model enables the understanding of the heat transfer phenomena in the equipment as well as the prediction of internal temperature distribution. This developed heat transfer model will be further integrated with the other models already developed by our teams for several refrigerating equipment. By this mean, it allows the prediction of time-temperature history of food products from a production plant to a household refrigerator. Furthermore, these models linked with predictive microbiological and/or quality models can be used as an evaluation tool for consumer risk [START_REF] Duret | Identification of the significant factors in food quality using global sensitivity analysis and the accept-andreject algorithm. Part II: Application to the cold chain of cooked ham[END_REF].

MATERIALS AND METHOD

Experimental setup

Temperature measurement was carried out in a closed refrigerated display cabinet loaded with test packages made of methylcellulose. This cabinet was installed in a test room in which the ambient temperature was controlled at 19.5°C. The lights in the test room and in the cabinet were switched on to represent the real-use condition in supermarket and the cabinet doors were always closed during the experiment. As shown in Fig. 1, calibrated thermocouples (T-type, accuracy ±0.2°C) were placed at various positions in the cabinet for the measurement of air and load (core and surface) temperatures. Perforation pattern of the back panel (PBP) of the studied cabinet is also depicted in this figure. According to the total plate area over each shelf, there are 2.3% of holes at the top shelf, 8.3% at the bottom shelf and 3.8% at all intermediate shelves.The temperatures were recorded every 10 s using a data logger for 24 hours. The time-averaged values were then calculated over 5 h of the quasi-steady state during which the defrosting period was excluded from the calculation. These results were then used for the model validation. heat exchanger at which the thermostat sensor is located (𝑇 𝑡ℎ ) and the external air temperature (𝑇 𝑒𝑥𝑡 ) are the input parameters which are assumed constant. To generalize the model, dimensionless temperature (𝜃 𝑖 ) based on these temperatures is defined as follows:

𝜃 𝑖 = 𝑇 𝑖 -𝑇 𝑡ℎ 𝑇 𝑒𝑥𝑡 -𝑇 𝑡ℎ Eq. (1)

Zone definition

Based on the zonal approach, heat exchange zones at different positions in the cabinet were established (Fig. 2). Two main domains were considered, the first one represents the heat transfer in the rear duct where cold air exchanges heat with the ambient air through the rear wall while the second one represents the heat transfer in the shelf space.

Heat balance equations

For every zone in both domains, algebraic heat balance equations were derived to describe the heat transfers in the display cabinet. Starting from Domain I and Position 1 (denoted by EZ [1,1], Fig. 2), the supply air (mass flow rate 𝑚̇𝑡 and temperature 𝜃 𝑡ℎ ) exchanges heat with the external air (𝜃 𝑒𝑥𝑡 ) through the rear wall by conduction and by convection (overall heat transfer coefficient 𝑈 𝑜𝑣𝑤 ). The air temperature then increases to 𝜃 𝑎𝑟𝑑1 . By integration along the rear duct located at this position, the heat balance equation becomes: Eq. ( 4)

𝜃 𝑎𝑟𝑑1 = 𝛼 𝑟𝑑1 𝜃 𝑡ℎ + (1 -𝛼 𝑟𝑑1 )
In addition to this convective heat exchange with air, the back load also exchanges heat with the front one by conduction (heat transfer coefficient 𝑈 𝑙𝑙1 ) because of the contact between these two loads.

𝜃 𝑎𝑠ℎ21 = 𝜃 𝑎𝑠ℎ11 + 𝛽 𝑠ℎ11 (𝜃 𝑙𝑓1 -𝜃 𝑙𝑏1 ) where 𝛽 𝑠ℎ11 = 𝑈 𝑙𝑙1 𝑆 𝑙𝑙1 𝑚̇𝑎 𝑠ℎ1 𝐶 𝑝𝑎 Eq. ( 5) A similar heat balance was developed for other shelves. The air then mixes with other two air streams at MZ [2,1,1] (Fig. 2) including air from DAG (mass flow rate 𝑚̇𝑎 𝑟𝑑5 and temperature 𝜃 𝑎𝑟𝑑7 ) and a part of the recirculating air (mass flow rate 𝛿 1 𝑚̇′ 𝑎𝑠ℎ1 and temperature 𝜃 𝑎𝑠ℎ51 ).

𝑚̇𝑎 𝑠ℎ1 ′ 𝜃 𝑎𝑠ℎ31 = 𝑚̇𝑎 𝑟𝑑5 𝜃 𝑎𝑟𝑑7 + 𝑚̇𝑎 𝑠ℎ1 𝜃 𝑎𝑠ℎ21 + 𝛿 1 𝑚̇𝑎 𝑠ℎ1 ′ 𝜃 𝑎𝑠ℎ51 where 𝑚̇𝑎 𝑠ℎ1 ′ = 𝑚̇𝑎 𝑠ℎ1 +𝑚̇𝑎 𝑟𝑑5 1-𝛿 1
Eq. ( 6)

This mixed air then exchanges heat with the front load, EZ [2,3,1], and its temperature becomes 𝜃 𝑎𝑠ℎ41 .

𝜃 𝑎𝑠ℎ41 = 𝛼 𝑠ℎ21 𝜃 𝑎𝑠ℎ31 + (1 -𝛼 𝑠ℎ21 )𝜃 𝑙𝑓1 where 𝛼 𝑠ℎ21 = exp (-

𝑈 𝑙𝑓1 𝑆 𝑙1 𝑚̇𝑎 𝑠ℎ1 ′ 𝐶 𝑝𝑎 )
Eq. ( 7)

In addition to this convective heat exchange with the air, the front load also exchanges heat with the back load by conduction (heat transfer coefficient, 𝑈 𝑙𝑙1 ) and receives heat from the conversion of the room and cabinet lighting powers (𝑞 𝑙𝑖 ). At steady state: Eq. ( 9)

𝜃 𝑎𝑠ℎ41 = 𝜃 𝑎𝑠ℎ31 + 𝛽 𝑠ℎ21
Then, the mixed air exchanges heat with the front load, EZ [2,3,2], and its temperature becomes 𝜃 𝑎𝑠ℎ42 .

𝜃 𝑎𝑠ℎ42 = 𝛼 𝑠ℎ22 𝜃 𝑎𝑠ℎ32 + (1 -𝛼 𝑠ℎ22 )𝜃 𝑙𝑓2 where 𝛼 𝑠ℎ22 = exp (-

𝑈 𝑙𝑓2 𝑆 𝑙2 𝑚̇𝑎 𝑠ℎ2 ′ 𝐶 𝑝𝑎 )
Eq. ( 10)

The load located at the front of Shelf 2 is subjected to the convection with the mixed air, the conduction with the adjacent load at the back (heat transfer coefficient 𝑈 𝑙𝑙2 ) and the heat conversion of visible radiation from the room lighting power (𝑞 𝑙𝑖2 ), thus

𝜃 𝑎𝑠ℎ42 = 𝜃 𝑎𝑠ℎ32 + 𝛽 𝑠ℎ22 (𝜃 𝑙𝑏2 -𝜃 𝑙𝑓2 ) + 𝑞 ′′ 𝑙𝑖2 where 𝛽 𝑠ℎ22 = 𝑈 𝑙𝑙2 𝑆 𝑙𝑙2 𝑚̇𝑎 𝑠ℎ2 ′ 𝐶 𝑝𝑎 ; 𝑞 ′′ 𝑙𝑖2 = 𝑞 𝑙𝑖2 𝑚̇𝑎 𝑠ℎ2 ′ 𝐶 𝑝𝑎 (𝑇 𝑒𝑥𝑡 -𝑇 ̅ 𝑡ℎ )
Eq. ( 11)

After exchanging heat with the front load, a part of the air at Shelf 1 (mass flow rate 𝛿 1 𝑚̇𝑎 𝑠ℎ1 and temperature 𝜃 𝑎𝑠ℎ41 ) flows toward the doors and is submitted to heat losses through the doors (overall heat transfer coefficient 𝑈 𝑜𝑣𝑑 ). This air also interacts with the air infiltration through the door gaps (mass flow rate 𝑚̇𝑔 𝑎𝑝 ). These exchanges at EZ [2,4,1] result in an air temperature increase to 𝜃 𝑎𝑠ℎ51 .

𝜃 𝑎𝑠ℎ51 = 𝛼 𝑠ℎ31 𝜃 𝑎𝑠ℎ41 + (1 -𝛼 𝑠ℎ31 )𝜃 𝑒𝑥𝑡 where 𝛼 𝑠ℎ31 = exp (-

𝐶 𝑒𝑞 𝑆 𝑑1 𝛿 1 𝑚̇𝑎 𝑠ℎ1 ′ 𝐶 𝑝𝑎 )
Eq. ( 12)

For simplification and as a first approach, it is considered that there is no air recirulation and air infiltration through the door gaps for the other shelves. A part of the air at these shelves, e.g. Shelf 2, is mixed with the air from the top shelf (mass flow rate (1 -𝛿 2 )𝑚̇𝑎 𝑟𝑑4 and temperature 𝜃 𝑎𝑠ℎ41 ) and exchanges heat with the external ambient air through the glass doors (overall heat transfer coefficient 𝑈 𝑜𝑣𝑑 ).

𝜃 𝑎𝑠ℎ52 = 𝛼 𝑠ℎ32 𝜃 𝑎𝑠ℎ41 + (1 -𝛼 𝑠ℎ32 )𝜃 𝑒𝑥𝑡 where 𝛼 𝑠ℎ32 = exp (- 13) This air then mixes with the air which exchanged heat with the load located at the front of the shelf at MZ [2,2,2].

𝑈 𝑜𝑣𝑑 𝑆 𝑑2 (1-𝛿 2 )𝑚̇𝑎 𝑟𝑑4 𝐶 𝑝𝑎 ) Eq. (
𝑚̇𝑎 𝑟𝑑3 𝜃 𝑎𝑠ℎ62 = 𝑚̇𝑎 𝑠ℎ2 ′ 𝜃 𝑎𝑠ℎ42 + (1 -𝛿 2 )𝑚̇𝑎 𝑟𝑑4 𝜃 𝑎𝑠ℎ52
Eq. ( 14)

A similar heat and mass balances at EZ [2,4,2] and MZ [2,2,2] were applied to Shelf 3 to 5. Input parameters used in the developed model were obtained from several sources: physical measurement, calculation based on literature and material properties and model fitting as summarized in Table 1. The simulation was run on MATLAB environment.

RESULTS AND DISCUSSION

Fig. 3a shows the predicted and measured air and load temperatures inside the closed display cabinet for the ambient temperature of 19.5°C. It was observed that the air temperature in the rear duct increased progressively from the duct upstream to the DAG (Fig. 3a) because of heat losses through the rear wall. As a part of this cold air also infiltrated into the cabinet through PBP, it explains higher temperature of back load at the top shelf compared to the lower shelves. Similarly, higher temperature of the load at the front top shelf was observed compared to that at the lower shelves. On the same shelf, the front load temperature was relatively higher (average difference of about 1.4°C) compared to the one at the back because of heat losses through the glass doors. The maximum difference (2.1°C) at the top shelf can be explained by additional heat loads from the warm air infiltration through the door gaps and the heat conversion of visible radiation from the room and cabinet lightings. The highest and lowest load temperature positions were located at the front of the top shelf and the back of the bottom shelf, respectively.

As shown in Fig. 3b, the calculated air and load temperatures agreed with the measured ones with an overall mean absolute error of about 0.2°C. A slight over-estimated air temperature at the front can be explained by the measurement position which may not be a good representation of the entire air across the front area. Nevertheless, the load temperature is well predicted. The maximum difference of the calculated and measured load temperatures was less than 0.5°C for every studied position. It is worth to remind that the developed model does not take into account the air infiltration due to door openings because the airflow during door openings is too complex to be predicted by the simplified model. The validated model was further used to investigate the effect of design parameters on the load temperature in the closed display cabinet including rear wall insulation thickness (𝑒) and air infiltration rate through the door gaps (𝑚̇𝑔 𝑎𝑝 ). As shown in Fig. 4a, doubling insulation wall thinkness resulted in a slight decrease of load temperatures (about 0.2°C). Load temperatures could be notibly decreased (up to 0.7°C) when the air infiltration through the door gap reduced by half (Fig. 4b). These findings suggested that the temperature performance of closed refrigerated display cabinets could be improved by minimizing the door gaps of the cabinets.

CONCLUSIONS

A steady state simplified heat transfer model was developed for a closed refrigerated display cabinet. Based on a zonal approach, the display cabinet was decomposed into various heat exchange zones allowing the prediction of air and load temperatures at different positions. In the shelves-space storage, time-averaged load temperatures at 10 positions (front and back of the five shelves) were predicted and compared with the experimental results. It was observed that the load temperatures were well predicted by the model whereas the air temperature in the air curtain was slightly over-estimated. Possible explanation of the discrepancy between the model and the experiment is that the air temperature was measured only at one position which may not be representative of the average air temperature inside the curtain . Nevertheless, the load temperature ranges and trends in the closed display cabinet are well represented by the proposed model with a short calculation time (less than 1 s). The overall mean absolute error was less than 0.2°C, and the maximum difference of the calculated and measured load temperatures never exceeds 0.5°C for every studied position.

Moreover, this validated model can be used to study the influence of display cabinet design on the thermal performance.

Figure

  Figure 1: Experimental setup for the temperature measurement

  𝜃 𝑒𝑥𝑡 where 𝛼 𝑟𝑑1 = exp (-𝑈 𝑜𝑣𝑤 𝑆 𝑤1 𝑚̇𝑡𝐶 𝑝𝑎 ) Eq. (2) A part of air flows into the shelve space through the perforated back panel (PBP) while the remaining part (mass flow rate 𝑚̇𝑎 𝑟𝑑 ) continues flowing towards the discharge air grille (DAG), along which the air temperature increases successively to 𝜃 𝑎𝑟𝑑7 . The heat balance of EZ [1,2…7] was therefore carried out in the same manner as that of EZ [1,1].

Figure 2 :

 2 Figure 2: Schematic view of the simplified heat transfer model of a closed refrigerated display cabinet.

Figure 3 :

 3 Figure 3: (a) Temperature distribution obtained from the numerical simulation (𝑻 𝒏𝒖𝒎 ) and measurement (𝑻 𝒆𝒙𝒑 ) in a closed display cabinet for ambient temperature of 19.5°C (b) comparison between the predicted and measured air and load temperatures.

Figure 4 :

 4 Figure 4: Numerical results showing the variation of air temperature at discharge air grille (DAG) and load temperatures with (a) wall insulation thickness and (b) air infiltration rate through the door gaps

  Domain II is composed of five shelves. Each shelf has four heat exchange zones (EZ [2,1…4,1…5]) and one or two air mixing zones (MZ [2,1…2,1…5]). At Position 1 Shelf 1, i.e. EZ [2,1,1] (see Fig.2), a part of the air (mass flow rate 𝑚 ̇𝑎𝑠ℎ1 = 𝛾 1 𝑚̇𝑡 and temperature 𝜃 𝑎𝑟𝑑5 ) flows horizontally through PBP at the mid height of Shelf 1. At the steady state, the air temperatures before and after flowing through PBP are the same: ], the air (temperature 𝜃 𝑎𝑠ℎ11 ) exchanges heat by convection (heat transfer coefficient 𝑈 𝑙𝑏1 ) with the back load (temperature 𝜃 𝑙𝑏1 ). The air temperature becomes 𝜃 𝑎𝑠ℎ21 .

	𝜃 𝑎𝑠ℎ11 = 𝜃 𝑎𝑟𝑑5		Eq. (3)
	At EZ [2,2,1𝜃 𝑎𝑠ℎ21 = 𝛼 𝑠ℎ11 𝜃 𝑎𝑠ℎ11 + (1 -𝛼 𝑠ℎ11 )𝜃 𝑙𝑏1 where 𝛼 𝑠ℎ11 = exp (-	𝑈 𝑙𝑏1 𝑆 𝑙1 𝑚̇𝑎 𝑠ℎ1 𝐶 𝑝𝑎	)

  Unlike Shelf 1, different heat balance is applied to the other shelves. The heat transfer mechanism on Shelf 2 was selected as a sample of the model development. The air flowing through PBP (mass flow rate 𝑚̇𝑎 𝑠ℎ2 and temperature 𝜃 𝑎𝑠ℎ22 ) mixes with a part of the air curtain from the upper shelf (mass flow rate 𝛿 2 𝑚̇𝑎 𝑟𝑑4 and temperature 𝜃 𝑎𝑠ℎ41 ).

		(𝜃 𝑙𝑏1 -𝜃 𝑙𝑓1 ) + 𝑞 ′′	𝑙𝑖1 where 𝛽 𝑠ℎ21 =	𝑈 𝑙𝑙1 𝑆 𝑙𝑙1 𝑚̇𝑎 𝑠ℎ1 ′ 𝐶 𝑝𝑎	; 𝑞 ′′	𝑙𝑖1 =	𝑚̇𝑎 𝑠ℎ1 ′	𝑞 𝑙𝑖1 𝐶 𝑝𝑎 (𝑇 𝑒𝑥𝑡 -𝑇 ̅ 𝑡ℎ )	Eq. (8)
	𝑚̇𝑎 𝑠ℎ2 ′	𝜃 𝑎𝑠ℎ32 = 𝑚̇𝑎 𝑠ℎ2 𝜃 𝑎𝑠ℎ22 + 𝛿 2 𝑚̇𝑎 𝑟𝑑4 𝜃 𝑎𝑠ℎ41 where 𝑚̇𝑎 𝑠ℎ2 ′	= 𝑚̇𝑎 𝑠ℎ2 + 𝛿 2 𝑚̇𝑎 𝑟𝑑4

Table 1 Values of input parameters used in the simplified heat transfer model

 1

	Parameter		Value Source		Value Source
	Mass flow rate (kg•s -1 )	𝑚̇𝑡	0.1120 Measurement	𝑚̇𝑔 𝑎𝑝 0.0051 Model fitting
	Air distribution	𝛾 1	0.02	Measurement	𝛿 1	0.90	Model fitting
	coefficient (-)	𝛾 2	0.02		𝛿 2	0.70
		𝛾 3	0.04		𝛿 3	0.70
		𝛾 4	0.05		𝛿 4	0.70
		𝛾 5	0.36		𝛿 5	0.70
	Heat transfer	𝑈 𝑙𝑏	15	Model fitting	𝑈 𝑜𝑣𝑤	0.77	Calculation
	coefficient (W•m -2 •K -1 )	𝑈 𝑙𝑓 𝑈 𝑙𝑙	24 3.4	Model fitting Calculation	𝑈 𝑜𝑣𝑑 𝐶 𝑒𝑞	1.02 17.6	Calculation Calculation
	Exchange area* (m 2 )	𝑆 𝑤1 𝑆 𝑤2…5 0.313 Measurement 0.216 Measurement	𝑆 𝑐 𝑆 𝑑1…4	0.500 Measurement 0.313 Measurement
		𝑆 𝑤6	0.156 Measurement	𝑆 𝑑5	0.431 Measurement
	Lighting power (W)	𝑞 𝑙𝑖1	10	Calculation	𝑞 𝑙𝑖2…5	0.25	Calculation
	*The cabinet is 1.25 m in length.


Mouna Merai

Steven Duret
email: steven.duret@irstea.fr

Onrawee Laguerre

Denis Flick
email: denis.flick@agroparistech.fr

Modelling the fluid flow inside a refrigerated trailer loaded with pork carcasses

Keywords: Aeraulics, Carcass, Computational Fluid Dynamics, Refrigerated vehicle

When pork carcasses are transported in a refrigerated trailer, all of them should be correctly ventilated. CFD simulations were performed to assess the homogeneity of air circulation around the carcasses. From an Xray scanner of a carcass, a simplified geometry was deduced. The air volume inside a semitrailer containing 430 pork half-carcasses was then meshed (21.10 6 cells). The Reynolds Averaged Navier Stokes equations were solved with the k-ε turbulence model. The flow pattern is very complex: air is blown at high velocity on the top of the container front, forming a wall jet. There is an air entrainment by this jet whose section increases whereas its velocity decreases. Part of the blown air flows throughout the 'porous' zone constituted by the carcasses while another part reaches the gap between the carcasses and the doors at the rear of the container. Finally air flows back toward the outlet. Zones of low velocity between half-carcasses were identified.

INTRODUCTION

When pork carcasses are transported in a refrigerated trailer, all of them should be correctly ventilated. But in fact, it has be shown both by experimentation and computational fluid dynamics (CFD) that air flow can be very heterogeneous in refrigerated trailers. Up to now, the studies concerned empty trailers [START_REF] Moureh | Airflow characteristics within slot ventilated enclosures[END_REF], trailers loaded with impermeable or porous pallets [START_REF] Tapsoba | Airflow patterns in a slot-ventilated enclosure partially loaded with empty slotted boxes[END_REF][START_REF] Moureh | Air velocity characteristics within vented pallets loaded in a refrigerated vehicle with and without air ducts[END_REF]) or a single carcass [START_REF] Pham | Finite element model for beef chilling using CFD-generated heat transfer coefficients[END_REF][START_REF] Kuffi | CFD modelling of industrial cooling of large beed carcasses[END_REF].

The present study concerns the case of a long trailer (13.3 m) loaded with 430 pork half-carcasses. The air inlet and outlet are both located at the front of the container, this generally implies a better ventilation in the front part than in the rear part. The air is blown at high velocity (11 m/s) near the ceiling of the container front, forming a turbulent wall jet. The particularity of this configuration in comparison with previous studies is the very complex geometry of the load. A numerical approach is here presented.

Many CFD studies of air flow in refrigeration equipment use the Reynolds Averaged Navier Stokes (RANS) approach. The objective is to predict the time-averaged mean velocity (v). Most often the Reynolds stresses, due to the velocity fluctuations (v') are expressed with the Boussinseq approximation by introducing the turbulent viscosity (ν t ). This way, the eddies are not directly solved but only characterized by their turbulent kinetic energy (k). Following Kolmogorov's point of view, mechanical energy of the mean flow is transformed into turbulent kinetic energy (because of the mean velocity gradients) and the energy of the smallest eddies is transformed into heat (molecular agitation) with a specific dissipation rate (ε). For complex geometries, an algebraic estimation of the turbulent viscosity cannot be used; the most popular approach relates ν t to k and ε (taking into account indirectly of the turbulence length scale l t ). Different versions are available which differ by some energy redistribution terms and near wall treatment (standard k-ε. realizable k-ε, k-ε RNG, k-ω SST). Near wall treatment is especially important if heat transfer has to be predicted. Sometimes the Boussinesq approximation cannot be applied because eddies are highly anisotropic and more sophisticated model (e.g. RSM) have to be used [START_REF] Defraeye | CFD Modeling of flow and scalar exchange of spherical food products; turbulence and boundary layer modelling[END_REF].

Near the walls, v varies very rapidly in a nonlinear way, so that a very fine mesh should be employed near the walls. But this becomes very expensive in term of memory and computing time. Therefore, assuming that the flow is almost parallel to the wall in the first cell near a wall, the boundary layer theory can be employed to estimate the wall shear stress. In the software the employed near wall relations are called 'wall functions'.

The present study aims to predict the homogeneity of the fluid flow around the carcasses. No precise velocity field or estimation of the heat transfer intensity is expected; therefore the standard k-ε model (with wall functions) was used. From an Xray scanner of a carcass, a simplified geometry was deduced as shown on Fig. 2.

GEOMETRY AND MESH

The simplified geometry was developed in order to have the same volume ensuring the same porosity of the zone occupied by the carcasses. The horizontal section at a given height is also about the same for the real carcass and the simplified geometry, ensuring almost the same convergent or divergent flow in the vertical direction. The air volume inside a semitrailer containing 430 pork half-carcasses was then meshed using 21 009 342 cells (Fig. 3). Near the top wall, the mesh is prismatic in order to better take into account the wall jet development.

MAIN HYPOTHESES AND BOUNDARY CONDITIONS

The fluid was assumed incompressible and isothermal. Properties of air at 4°C were used.

The steady state Reynolds Averaged Navier Stokes (RANS) equations with the standard k-ε model of turbulence and the standard (logarithmic) wall function (smooth walls) were solved using the finite volume method (Fluent).

The inlet conditions are: -section 1.100 m x 0.16 m (hydraulic diameter: 0.28 m) -normal mean velocity: 11 m/s (Reynolds number: 2.2 10 5 ) -turbulence intensity: 10 % -turbulence length : 7% of hydraulic diameter At outlet, uniform pressure was imposed (without viscous stress).

The simulation converged after 3500 iterations which needed 48 hours of computation on a 12 cores computer with 130 Go RAM.

RESULTS AND DISCUSION

Fig. 4 present some computed stream lines inside the loaded trailer (the carcasses are not represented in this figure but where taken into account in the simulation). This highlights the complexity of the 3D flow pattern. Globally air blown near the ceiling forms a wall jet. Its section increases whereas its velocity decreases (1). It adheres to the top wall (by Coanda effect) nearly up to the rear of the trailer. Part of the blown air flows throughout the 'porous' zone constituted by the carcasses while another part reaches the gap between the carcasses and the doors at the rear of the container (2). Finally air flows back toward the outlet near the bottom (3). There is an air entrainment by the jet which explains an upward flow on the front of the trailer (4). Fig. 5a presents the computed velocity magnitude in the symmetry plan. This velocity field can be compared outside the zone of carcasses with the experimental one (Fig. 4b) which was obtained on a reduced scale trailer loaded with reduced scale carcasses made of polyurethane foam [START_REF] Merai | Experimental characterization of airflow inside a refrigerated trailer loaded with carcasses[END_REF]) by 2D Laser Doppler Velocimetry (LDV). The main difference is that the predicted velocity decreases less rapidly along the jet axis than the measured one. This can be explained by several reasons. The simulated geometry does not take into account the presence of the hooks and the extremity of the feet. Also the heights at which the carcasses are hanged is not exactly the same in the experiment (and in real configuration) to the contrary of simulation. All these aspects could contribute to produce turbulence and increase the jet diffusion in the lateral and vertical directions. The k-ε model is also known to predict incorrectly flows with adverse pressure gradients. This explained that for of an empty trailer, the k-ε was unable to predict the experimentally observed separation of the jet from the wall [START_REF] Moureh | Airflow characteristics within slot ventilated enclosures[END_REF] whereas the RSM model, which better takes into account of turbulence anisotropy, predicted it well. More realistic geometry and more advanced turbulence models will be tested in future work.

Nevertheless, even taking in mind the limitations of the present numerical results, the interest of the CFD approach is that, contrarily to the experimental approach, it is possible to estimate the air flow around the carcasses. Fig. 7 presents the velocity magnitude at different position between the half carcasses. Some zones of low velocity between half carcasses can be identified. In plane (a) notably, a nearly stagnant zone at mid length of the trailer was predicted because the simulated jet wall continued flowing along the top wall until the rear of the trailer. In fact, according to the experimental velocity evolution on the jet axis, is seems that less air reaches the trailer rear and more air flows downward through the carcasses.

Fig. 8 present the vertical velocity component at mid height of the trailer. One could imagine that the carcasses form a porous medium of relative high pressure drop and that since the air is blown on the top, a relative homogenous downward airflow would occur throughout the carcass zone. But in fact, near the trailer rear the flow is indeed downward but near the front, the entrainment effect of the wall jet induces a suction effect and air flows upward throughout the carcasses in the front zone. In the intermediate region, there is sometimes neither downward and neither upward flow leading to low velocity zones. In the zones of low velocity, the refrigeration will be of low intensity. Thus, if the carcasses are loaded while there are still warm at core (e.g. 15°C at core and 7°C at surface), they could take a long time to reach everywhere the recommended temperature (<7°C). Even the surface temperature could slightly increase. This can be problematic for sanitary risk (pathogenic microorganism growth). In fact, the smallest velocities are observed in the loan region of the carcasses (at mid height of the carcasses). Since the loan thickness is small (compared to the jam for example) the loan is generally rapidly cooled down in the cold room of the slaughterhouse. So, microbial risk should be limited for the loan part. In opposite, the jam core is often the warmest part at the exit of the slaughterhouse but the velocity is rather high all along the trailer at the top of the carcasses (jam part). The shoulder regions are also quite well ventilated.

CONCLUSIONS

Numerical simulations were performed to determine the air flow distribution in a refrigerated trailer loaded with pork half carcasses using the standard k-ε model. In the studied configuration, some low velocity zones were located in the loan region of the carcasses but this should not be problematic for sanitary risk.

These numerical results are specific of the trailer geometry (inlet and outlet section and location, number of rails) and of the load (here 215 pork carcasses about 75 kg each). The numerical approach could be used to see the influence of the trailer configuration, including the effect of air ducts, and of load (lower or higher number of carcasses) on the airflow distribution.

Nevertheless, the numerical model should be improved, because it overestimates the flow rate which reaches the rear of the trailer. In further work, a more realistic geometry and more sophisticated turbulence model (e.g. RSM or k-ω SST) will be tested and the simulation will also integrate unsteady heat transfer.

Fig. 1

 1 Fig.1 presents the dimensions of the trailer and of the zone occupied by the carcasses.

Figure 1 :

 1 Figure 1: Dimensions of the trailer (real scale)

Figure 2 :

 2 Figure 2: Simplified geometry of the carcasses

Figure 3

 3 Figure 3 View of the mesh

Figure 4 :

 4 Figure 4 : Stream lines issued from the inlet section

Figure 5 :

 5 Figure 5: velocity magnitude in the symmetry plan (red color represents velocities ranging from 6m/s to 11 m/s which is the value at inlet)

Figure 6 :

 6 Figure 6 : Velocity along the jet axis

Figure 7 :

 7 Figure 7 : Velocity magnitude between half carcasses

Figure 8 :

 8 Figure 8: Vertical velocity (positive upward) at mid height of the trailer

ACKNOWLEDGEMENTS

This research work was carried out within the framework of a DIM ASTREA grant by the Regional Council of Ile-de-France (France).