Sara Croxford

Annemarie Rinder Stengaard
email: annemarie.rinder.stengaard@regionh.dk

Johanna Brännström

Lauren Combs

Nikos Dedes

Enrico Girardi

Sophie Grabar

Ole Kirk

Giorgi Kuchukhidze

Jeffrey V Lazarus

Teymur Noori

Anastasia Pharris

Dorthe Raben

Jurgen K Rockstroh

Daniel Simões

Ann K Sullivan

Dominique Van Beckhoven

Valerie C Delpech

Late diagnosis of HIV : An updated consensus definition

Keywords: definition, diagnosis, monitoring, HIV, late presentation

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

MANUSCRIPT INTRODUCTION

Late HIV diagnosis is associated with poor outcomes, an increased risk of ongoing HIV transmission and high healthcare costs. [START_REF]Late Presentation Working Groups in EuroSIDA and COHERE. Estimating the burden of HIV late presentation and its attributable morbidity and mortality across Europe 2010-2016[END_REF][START_REF]Recording and investigation of late HIV diagnoses: good practice position statement[END_REF] As such, late diagnosis remains a key public health metric in assessing the success of HIV testing programmes. In 2010, a consensus statement was published in which late presentation of HIV was defined as presenting for HIV care having a CD4 count <350 cells/µL or with an acquired immunodeficiency syndrome (AIDS)-defining event. [START_REF] Antinori | Late presentation of HIV infection: a consensus definition[END_REF] This definition was endorsed by the European Centre for Disease Prevention and Control (ECDC) and the World Health Organization (WHO) Regional Office for Europe and has been used across Europe for clinical research and public health monitoring for more than 10 years. Completeness of HIV surveillance data on CD4 count at diagnosis is now high for most countries reporting to ECDC and WHO. [START_REF]HIV/AIDS surveillance in Europe 2021: 2020 data[END_REF] In recent years, testing for HIV has expanded and frequency has increased across some populations and regions, particularly in relation to the roll out of pre-exposure prophylaxis (PrEP) programmes. This has resulted in an increasing number of people, particularly men who have sex with men (MSM), being diagnosed with HIV during seroconversion, when their CD4 count may be temporary low, (known as the 'seroconversion effect' [START_REF] Lodi | Symptomatic illness and low CD4 cell count at HIV seroconversion as markers of severe primary HIV infection[END_REF][START_REF] Lodi | Time from human immunodeficiency virus seroconversion to reaching CD4+ cell count thresholds <200, <350, and <500 cells/mm³: assessment of need following changes in treatment guidelines[END_REF]). Using the current definition of late HIV presentation, these individuals are incorrectly assigned as being diagnosed late. This issue of overestimation has already been raised by research groups in Belgium, [START_REF] Sasse | Late presentation to HIV testing is overestimated when based on the consensus definition[END_REF] Sweden [START_REF] Brännström | Deficiencies in the health care system contribute to a high rate of late HIV diagnosis in Sweden[END_REF] and the United Kingdom (UK). [START_REF] Kirwan | Reassessing the late HIV diagnosis surveillance definition in the era of increased and frequent testing[END_REF] This has led to correction factors being applied to the late diagnosis rate of specific subgroups. The magnitude of these correction factors depends on the reclassification criteria, population, country, and study period, but was estimated to be as high as 9% in Belgium. [START_REF] Sasse | Late presentation to HIV testing is overestimated when based on the consensus definition[END_REF] Therefore, a working group established under the EuroTEST Initiative, with the support of the ECDC, WHO Regional Office for Europe and European AIDS Clinical Society (EACS), decided to revisit this definition, reviewing the feasibility of incorporating data on markers of recent infection to enable better distinction between people diagnosed with HIV late and people recently acquiring HIV.

MATERIALS AND METHODS

The EuroTEST Initiative convened a working group of experts in HIV from Europe, including clinicians, epidemiologists, public health professionals and civil society, to review the existing late HIV diagnosis definition in January 2022. Multiple meetings were held for stakeholders to discuss possible updates to the definition based on previous research. [START_REF] Sasse | Late presentation to HIV testing is overestimated when based on the consensus definition[END_REF][START_REF] Brännström | Deficiencies in the health care system contribute to a high rate of late HIV diagnosis in Sweden[END_REF][START_REF] Kirwan | Reassessing the late HIV diagnosis surveillance definition in the era of increased and frequent testing[END_REF] The full list of experts consulted can be found in Appendix 1.

In May 2022, a short survey was developed by the EuroTEST HIV Late Diagnosis Definition Working Group and sent to all 53 national surveillance contact points for HIV in the WHO European Region by the ECDC and WHO (one per country). Responses were entered onto the Research Electronic Data Capture (REDCap) online tool, hosted at the Centre of Excellence for Health, Immunity and Infections (CHIP). [START_REF] Harris | The REDCap consortium: Building an international community of software platform partners[END_REF][START_REF] Harris | Research electronic data capture (REDCap)--a metadata-driven methodology and workflow process for providing translational research informatics support[END_REF] The main aim of this survey was to better understand the availability and flow of data on recent infection needed to be able to re-classify late HIV diagnoses and to monitor the modified late diagnosis indicator at a national and European level. Respondents were asked about baseline assessments carried out when people are diagnosed with HIV, data collection, data caveats and whether there are currently any adjustments made to national late HIV diagnosis figures to account for recent infection in their country. Submitted data was validated by the ECDC, where applicable. The full survey can be found in Appendix 2.

RESULTS

Definition

Late HIV diagnosis: A person first diagnosed with HIV with a CD4 count <350 cells/µl or with an AIDS-defining event, regardless of the CD4 cell count.

People with evidence of recent infection (i.e. being diagnosed during seroconversion) should be reclassified as 'not late'. Evidence of recent infection should be considered hierarchically; the individual must have either: i) laboratory evidence of recent infection (recent infection testing algorithm (RITA), p24 antigen), ii) a last negative HIV test within 12 months of HIV diagnosis or iii) clinical evidence of acute infection (e.g. seroconversion illness). People with evidence of being previously diagnosed, either abroad or elsewhere, should be excluded from the calculation of the proportion of the population diagnosed late. This definition has been adapted from the previously published definition (3) to describe late HIV diagnosis rather than late presentation, to focus on people newly diagnosed with HIV.

Twelve months was chosen as a pragmatic cut-off for timing of last negative HIV test to broadly align with current HIV testing recommendations. [START_REF]Public health guidance on HIV, hepatitis B and C testing in the EU/EEA -An integrated approach[END_REF] It is also known that most people who are seroconverting recover their CD4 count within a year,(13) typically within six months, and those who do not could rather be considered as 'fast progressors'. [START_REF] Lodi | Time from human immunodeficiency virus seroconversion to reaching CD4+ cell count thresholds <200, <350, and <500 cells/mm³: assessment of need following changes in treatment guidelines[END_REF][START_REF] Robb | Prospective study of acute HIV-1 Infection in adults in East Africa and Thailand[END_REF][START_REF] Kaufmann | Impact of early HIV-1 RNA and T-lymphocyte dynamics during primary HIV-1 infection on the subsequent course of HIV-1 RNA levels and CD4+ T-lymphocyte counts in the first year of HIV-1 infection. Sydney Primary HIV Infection Study Group[END_REF] Furthermore, research from the UK and Sweden shows that among people newly diagnosed with HIV with a negative test recorded, the vast majority had their last negative test within the most recent year, so lengthening the cut-off would have minimal effect. [START_REF] Kirwan | Reassessing the late HIV diagnosis surveillance definition in the era of increased and frequent testing[END_REF][START_REF] Brännström | EuroTEST HIV Late Diagnosis Definition Working Group Meeting[END_REF] The 12month cut-off was endorsed by members of the EACS Governing Board and antiretroviral therapy (ART) panel group.

We recommend people with evidence of being previously diagnosed with HIV should be excluded from the late diagnosis calculation as their current positive test is not their first and thus, they are not being newly diagnosed.

Survey to European countries

Overall, the response rate for the survey was relatively low, with respondents from only 18 countries participating: Albania, Belgium, Denmark, France, Georgia, Germany, Greece, Ireland, Italy, Liechtenstein, Luxembourg, Malta, Netherlands, Serbia, Slovakia, Spain, Sweden and the United Kingdom (UK). All questions were answered by respondents from all countries.

Baseline assessments at presentation: In all 18 countries, the standard of care for people newly diagnosed with HIV is to have a CD4 count taken and be asked about previous HIV diagnosis elsewhere (e.g. abroad). Seventeen (94%) respondents reported that people are assessed for clinical symptoms of AIDS-defining illnesses (clinical judgement based on symptoms and medical history), 16 (89%) that HIV testing history information is collected, 16 (89%) that the individual is assessed for clinical symptoms of seroconversion illness and 14 (78%) that laboratory testing for evidence for seroconversion (HIV polymerase chain reaction (PCR) or antigen positive but HIV antibody negative) is carried out. Only six (38%) country respondents reported people newly diagnosed with HIV are tested for recent infection (e.g.

avidity testing).

Clinical data flows: Not all country respondents reported being able to capture data needed to re-classify late HIV diagnoses at a national level (Table 1); 15 (83%) reported AIDS at diagnosis is able to be captured through national surveillance mechanisms, 15 (83%) CD4 count at diagnosis, 12 (67%) last negative HIV test, 12 (67%) previous HIV diagnosis, 11 (61%) information on seroconversion illness at diagnosis and only five (28%) incident HIV antibody test results. While the respondent for Belgium reported that laboratory findings of seroconversion are collected as part of national HIV surveillance with a coverage of almost 100%, the respondent for Ireland reported that all avidity testing has been paused due to resourcing issues. In Georgia, RITA surveys are conducted periodically, dependent on funding, but are not part of routine surveillance. In Sweden, incidence HIV antibody results are only available from HIV clinical cohort study data from Stockholm and Gothenburg. Only respondents from three countries, France, Ireland and the UK, reported that all surveyed indicators needed to apply the revised late HIV diagnosis definition are captured nationally.

Most respondents who reported at least one data item not being available at a national level, indicated that these data were available either locally (nine countries) or through cohort studies (three countries), highlighting potential for expanding data collection nationally and opportunities for collaboration. Data caveats: Respondents for each country were asked what data caveats would need to be considered if HIV surveillance data indicating recent infection were used to adjust late diagnosis figures at a national level. In terms of incident antibody testing, seven respondents reported testing was not carried out in their country, while five indicated the testing results were not collected centrally, one reported there was a significant reporting delay in receiving the results, one reported incomplete linkage between datasets and four that the data source did not cover all cases. The most common reasons for difficulty in providing last HIV test information were significant missing data (10 respondents), incomplete coverage (five respondents) and incomplete linkage between datasets (five respondents). Significant missing data was also a barrier reported by 10 respondents in collecting data on seroconversion illness; a further five respondents reported seroconversion illness data are not currently collected centrally.

Respondents from Spain, the UK, Denmark, France and Belgium reported that late HIV diagnosis figures are adjusted for recent infection in their country. The respondent from the Netherlands reported that figures had been adjusted previously and the respondent from Serbia reported that adjustments had been attempted previously. Fourteen respondents expressed an interest in being involved in future work to attempt to quantify the 'seroconversion effect' correction factor for MSM diagnosed late in their countries.

DISCUSSION

We present a revised consensus definition of late HIV diagnosis, in which people known to be diagnosed during seroconversion are reclassified as 'not late'. This represents a pragmatic approach to take into account observed increases in HIV testing frequency in Europe and this revision will ensure continued relevance of this longstanding key HIV metric for public health monitoring.

The application of this definition by national public health agencies and institutions is recommended but it is acknowledged this is dependent on reporting of surveillance data relating to late HIV diagnosis (CD4 count and AIDS at diagnosis) and evidence of recent infection (presence of seroconversion illness, last negative HIV test date or incident antibody testing). Our survey findings show this is currently only possible in some countries. However, it is important to note that respondents from only 18 countries responded to the survey, likely due to competing priorities, including coronavirus and monkeypox virus, with limited representation from Eastern European countries. Expansion of HIV surveillance mechanisms to facilitate data collection of these markers is essential to ensure uptake is comprehensive across Europe. In some circumstances, collaboration with HIV clinical cohort studies of people with HIV may be able to fill data gaps. Analyses from Spain and France show that in the absence of surveillance data, the use of cohort data is feasible for identifying the prevalence of recent infection at national level. [START_REF] Hernando | Recent infections among newly diagnosed cases of HIV infection in Spain, 2015-2016. National estimates using cohort data[END_REF][START_REF] Guillou | Reaching the second and third Joint United Nations Programme on HIV/AIDS 90-90-90 targets is accompanied by a dramatic reduction in primary Human Immunodeficiency Virus (HIV) infection and in recent HIV infections in a large French nationwide HIV cohort[END_REF] This consensus definition of late HIV diagnosis has been endorsed by the EACS and adopted by the ECDC and WHO Regional Office for Europe. The European Surveillance System (TESSy) for HIV currently includes a variable on acute HIV infection, [START_REF]TESSy -The European Surveillance System -HIV/AIDS Reporting Protocol and Analysis Plan 2021[END_REF] which will be adapted based on this work. This will facilitate comparisons between countries and assessment of trends over time across Europe where possible. It will also allow for a more accurate assessment of the effectiveness of national HIV testing programmes in reaching people who are underserved.

Respondents from some Western European countries reported already adjusting their national late HIV diagnosis figures for recent infection. [START_REF] Sasse | Late presentation to HIV testing is overestimated when based on the consensus definition[END_REF][START_REF] Brännström | Deficiencies in the health care system contribute to a high rate of late HIV diagnosis in Sweden[END_REF][START_REF] Kirwan | Reassessing the late HIV diagnosis surveillance definition in the era of increased and frequent testing[END_REF] With regard to the extent that these adjustments have been found to affect late HIV diagnosis rates, there was variation by country and sub-population. A study from Belgium, showed that, in 2012, late HIV diagnosis dropped by 9%, from 42%, to 33%, after reclassification based on reported recent infection by clinicians. [START_REF] Sasse | Late presentation to HIV testing is overestimated when based on the consensus definition[END_REF] In the UK, late HIV diagnosis dropped by 7%, from 49% to 42% in 2019, following reclassification based on RITA testing and/or a negative test within the last 24 months. [START_REF] Kirwan | Reassessing the late HIV diagnosis surveillance definition in the era of increased and frequent testing[END_REF] In both Belgium and the UK, reclassification was more frequent among MSM compared to people who acquired HIV through heterosexual contact. [START_REF] Sasse | Late presentation to HIV testing is overestimated when based on the consensus definition[END_REF] Data from the national Swedish InfCareHIV registry between 2017 to 2021, showed that late HIV diagnosis dropped from 55% to 52% after reclassification of people with primary HIV infection and evidence of a negative HIV test within one year of HIV diagnosis to 'not late'. [START_REF] Brännström | EuroTEST HIV Late Diagnosis Definition Working Group Meeting[END_REF][START_REF] Brännström | Email to: S. Croxford (UKHSA) regarding InfCareHIV study[END_REF] In conclusion, adoption of this revised consensus definition of late HIV diagnosis by national health agencies, institutions and researchers is needed to ensure consistent monitoring of access to HIV testing. Reclassification of individuals with recent infection will help to reduce overestimation of late HIV diagnosis estimates and better identify populations most at risk of poor HIV outcomes and areas for intervention, to further expand and target HIV testing in the era of elimination of HIV transmission. International public health bodies, such as the ECDC and WHO, should continue to work with countries to improve reporting of data needed to reclassify late HIV diagnoses. Furthermore, collaboration between agencies responsible for national HIV surveillance and HIV clinical cohorts should be strengthened to try to address gaps in data availability.

ACKNOWLEDGEMENTS

We would like to thank the ECDC and the WHO Regional Office for Europe for circulating the late HIV diagnosis survey to their European HIV Surveillance Network and the ECDC/WHO national HIV surveillance focal points for taking the time to complete it. We would also like to acknowledge the support of the EACS governing board and ART panel, the ECDC and the WHO Regional Office for Europe, as well as the members of the EuroTEST steering committee, not already listed as authors: Ben Collins (European Testing Week working group chair, ReShape/International HIV Partnerships, UK), Brian Gazzard (Imperial College, School of Medicine, HIV Research Director, Chelsea and Westminster Hospital, UK), Cary James, (World Hepatitis Alliance, UK), Daniela Rojas Castro (Coalition PLUS, France), Ferenc Bagyinszky (AIDS Action Europe, Germany), Francesco Negro (European Association for the Study of the Liver, Switzerland), Igor Karpov (Department of Infectious Disease, Belarus State Medical University, Belarus), Irith De Baetselier (Antwerp Institute of Tropical Medicine, Belgium), Jack S. Lambert (University College Dublin, Ireland), Jeffrey V. Lazarus (ISGlobal, Barcelona Institute for Global Health, Hospital Clínic -University of Barcelona, Spain), Jens D. Lundgren (Department of Infectious Diseases, University of Copenhagen and Centre of Excellence for Health, Immunity and Infections, Denmark), Jordi Casabona (Center for HIV/STI Epidemiological Studies of Catalonia, Spain), Jürgen Rockstroh (Department of Medicine I, University Hospital Bonn, Germany), Kira Grazava (TB Coalition Europe, Ukraine), Lella Cosmaro (Fondazione LILA Milan, Italy), Liudmyla Maistat (Medicines Patent Pool, Switzerland), Magnus Unemo, (WHO Collaborating Centre for Gonorrhoea and Other STIs, National Reference Laboratory for STIs, Örebro University Hospital, Sweden), Mark Vermeulen (Aidsfonds -Soa Aids Nederland, Netherlands), Mojca Matičič (University Medical Centre Ljubljana, Slovenia), Nino Tsereteli (Center for Information and Counseling on Reproductive Health -Tanadgoma, Georgia), Rajul Patel (Department of Genitourinary Medicine, Southampton University Hospitals and International Union against Sexually Transmitted Infections, UK), Tom Platteau (HIV/STI clinic of the Institute of Tropical Medicine in Antwerp, Belgium) and Yazdan Yazdanpanah (INSERM, France).

FUNDING

The EuroTEST Initiative receives funding from Gilead Sciences, Merck MSD, ViiV Healthcare and the ECDC through a framework contract unrelated to the content of this work (ECDC/2021/020). The funders had no input into the content of the manuscript.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available on reasonable request from the corresponding author. The data are not publicly available due to privacy restrictions.

AUTHOR CONTRIBUTIONS

All authors were involved in the development of the survey on late HIV diagnosis and contributed important intellectual content to this manuscript. All authors commented on the manuscript and approved the final draft. SC was responsible for creating the survey and carried out data analysis with LC and ARS. SC also drafted the manuscript, incorporated author comments, and was responsible for the final draft to be published. LC created the survey in REDCap. ARS and DR were responsible for convening the expert group on behalf of EuroTEST. VCD and DR contributed to the study conception. AS, DVB, GK, VD, AP, SC and TN provided public health and surveillance expertise; AS, SG, OK, EG and JB provided clinical input. DS and ND provided a community perspective.

DECLARATIONS OF INTEREST

SC, ARS, LC, JVL, JB, ND, SG, GK, AP, TN, DR, DS, AKS, DVB and VCD report no conflicts of interest to disclose. JKR reports honoraria for consulting or speaking at educational events from Abivax, Boehringer, Galapagos, Gilead, Janssen, Merck, Theratechnologies and ViiV, outside the submitted work. EG received grants support from Gilead Sciences and Mylan, and speaker's honoraria from Gilead Sciences, outside the submitted work. OK reports a research grant from Gilead, and non-financial support from Merck, outside the submitted work.


Charles Mareau
email: charles.mareau@ensam.eu

Thermodynamic framework for variance-based non-local constitutive models

Keywords: Thermodynamics, Constitutive laws, Non-locality, Damage, Plasticity 1 Introduction

The present work deals with the development of a framework dedicated to the construction of constitutive models with non-local internal variables. Such internal variables allow considering the impact of microstructural changes on the current state of a material point. Non-locality is introduced by considering, not only the spatial average, but also the spatial variance of an internal variable in constitutive relations. The proposed framework relies on continuum thermodynamics to construct the set of constitutive equations. Such a framework allows including some information regarding the spatial distribution of internal variables when constructing non-local models for thermomechanical applications. In contrast with gradient-type models, this strategy does not require additional equilibrium equations and boundary conditions. For the purpose of illustration, some numerical examples are presented. According to the numerical results, the proposed framework can be used to circumvent the difficulties associated with excessive spatial localization or to consider size effects.

materials. For the majority of constitutive models, the evolution equations associated with the different internal variables are purely local. Specifically, the evolution of the internal state of a material point solely depends on the current state of that material point. However, there are some situations for which it is needed to include some information regarding the spatial distribution of state variables. As discussed by Bazant and Jirásek [START_REF] Bazant | Nonlocal Integral Formulations of Plasticity and Damage: Survey of Progress[END_REF], there are different motivations for such constitutive models, which are often referred to as non-local. Indeed, non-local models allow (i) considering deviations from locality due to material heterogeneity at small scales, (ii) limiting strain localization resulting from softening, which is largely used in the context of damage mechanics, and (iii) capturing size effects associated with metal plasticity [START_REF] Nix | Mechanical properties of thin films[END_REF][START_REF] Fleck | Strain gradient plasticity: Theory and experiment[END_REF] or fracture of ceramic materials [START_REF] Walsh | Fracture of plain concrete[END_REF][START_REF] Jin | Experimental and numerical investigations on the size effect of moderate high-strength reinforced concrete columns under small-eccentric compression[END_REF].

Non-local constitutive models can be classified as either gradient-type or integral-type models 1 . Gradient-type models treat the spatial gradients of internal variables as additional state variables. The consequence is that the evolution of internal variables is governed by a set of equilibrium equations and some boundary conditions. Such equations and boundary conditions are generally derived from an extended principle of virtual power [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF][START_REF] Cermelli | Geometrically necessary dislocations in viscoplastic single crystals and bicrystals undergoing small deformations[END_REF][START_REF] Nguyen | Some remarks on standard gradient models and gradient plasticity[END_REF]. Gradienttype models have been largely used in the context of conventional plasticity, either to capture size effects [START_REF] Aifantis | The physics of plastic deformation[END_REF][START_REF] Fleck | A reformulation of strain gradient plasticity[END_REF] or to limit strain localization associated with softening [START_REF] Anand | A large-deformation gradient theory for elasticplastic materials: Strain softening and regularization of shear bands[END_REF]. Also, the role of Geometrically Necessary Dislocations (GND) on strain hardening can be incorporated in the framework of crystal plasticity by considering the spatial gradients of plastic shear strains as state variables [START_REF] Peng | Modeling dislocation absorption by surfaces within the framework of strain gradient crystal plasticity[END_REF][START_REF] Niordson | Computational strain gradient crystal plasticity[END_REF][START_REF] Wulfinghoff | A gradient plasticity grain boundary yield theory[END_REF]. The phase-field method, which is widely used for solving interfacial problems, also relies on the spatial gradients of some order parameters to evaluate surface energy. Many applications of the phase-field method to either brittle [START_REF] Miehe | Thermodynamically-consistent phase field models of fracture: variational principles and multi-field FE implementations[END_REF][START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF][START_REF] Borden | A higher-order phasefield model for brittle fracture: formulation and analysis within the isogeometric analysis framework[END_REF][START_REF] Gmati | A phasefield model for brittle fracture of anisotropic materials[END_REF] or ductile fracture [START_REF] Hernandez Padilla | A coupled phase-field model for ductile fracture in crystal plasticity[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF] have been proposed. When applied to fracture, the phase-field method allows circumventing the difficulties associated with excessive damage localization and considering the increase of surface energy resulting from crack nucleation and crack growth.

In contrast with gradient-type models, integral-type models use some nonlocal variables obtained after spatial integration of their local counterparts. Generally speaking, such non-local variables can be interpreted as the spatial averages of the corresponding local variables over the neighborhood of the material point of interest. Integral-type models have been largely used in the context of brittle [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Benvenuti | A thermodynamically consistent nonlocal formulation for damaging materials[END_REF][START_REF] Borino | A symmetric nonlocal damage theory[END_REF] or ductile fracture [START_REF] Leblond | Bifurcation effects in ductile metals with nonlocal damage[END_REF][START_REF] Andrade | A Ductile Damage Nonlocal Model of Integral-type at Finite Strains: Formulation and Numerical Issues[END_REF], mostly to limit damage localization. Some alternative formulations, using either an evolving length scale [START_REF] Pijaudier-Cabot | Nonlocal damage model with evolving internal length[END_REF] or an internal time [START_REF] Desmorat | Nonlocal models with damage-dependent interactions motivated by internal time[END_REF], have been proposed to consider damage-dependent interactions. As discussed by [START_REF] Jirásek | Comparison of integral-type nonlocal plasticity models for strain-softening materials[END_REF], integral-type models for plasticity-induced softening have also been developed. For instance, Bazant and Lin [START_REF] Bazant | Nonlocal yield-limit degradation[END_REF] proposed a non-local Mohr-Coulomb model that replaces the local plastic strain tensor by its non-local counterpart. The approach of Polizzotto et al. [START_REF] Polizzotto | A thermodynamically consistent formulation of nonlocal and gradient plasticity[END_REF], which is developed in a consistent thermodynamic framework, uses a non-local definition of the isotropic hardening variable. In the context of crystal plasticity, Gao and Huang [START_REF] Gao | Taylor-based nonlocal theory of plasticity[END_REF] and Counts et al. [START_REF] Counts | Predicting the HallPetch effect in fcc metals using non-local crystal plasticity[END_REF] developed a non-local formulation that accounts for the role of GND. For this purpose, they have used an integraltype of approach to evaluate the spatial gradient of the plastic deformation gradient tensor.

The present work aims at enriching integral-type constitutive models by incorporating not only the spatial average, but also the spatial variance of an internal variable in constitutive relations. The objective is to include additional information regarding the spatial distribution of an internal variable when constructing integral-type constitutive models. The proposed framework is developed in the context of continuum thermodynamics and can be applied to different types of problems, including plasticity and damage. In the first section, the definitions of the average and the variance of an internal variable are introduced. The general form of constitutive equations is detailed in the second section. Some illustrative examples are presented in the final section in either the context of damage mechanics, or crystal-plasticity.

Non-local internal variables

The development of constitutive relations, which is detailed in the next section, uses some internal variables to consider the impact of microstructural transformations on the state of a material point. The main idea of integral-type non-local models consists of assuming that the state of a material point (e.g. stress, internal energy) depends on the spatial distribution of state variables in its neighborhood. For this purpose, two different types of non-local variables are considered: the spatial average and the spatial variance. The definitions of these variables are briefly presented here.

Spatial average

Most integral-type models use some non-local variables that are obtained from the application of an averaging procedure. Specifically, for any local internal variable β, the spatial average µ β is defined according to [START_REF] Borino | A Symmetric Formulation for Nonlocal Damage Models[END_REF]:

µ β [X] = 1 W Vc w[X -X] β[X ] -β[X] dX + β[X] (1)
As indicated by the above equation, the non-local variable µ β is actually the first raw moment associated with the spatial distribution of the local variable β. The averaging procedure uses a weight function w that describes the interactions between material points with initial positions X and X . The region of space over which such interactions are considered is denoted by V c . The factor W is obtained by integrating the weight function over an infinite region of space V ∞ , that is:

W = V∞ w[X ]dX (2)
For the development of constitutive relations, it may be needed to perform a double-averaging operation. For any internal variable β, the corresponding double-averaged variable µ ββ is given by:

µ ββ [X] = 1 W Vc w[X -X] µ β [X ] -µ β [X] dX + µ β [X] (3)
In some situations, it may be useful to construct a gradient approximation of the spatial average. For this purpose, the internal variable at position X can be obtained from a Taylor series expansion at position X according to:

β[X ] = β[X] + k 1 k! X -X ⊗k • ∇ ⊗k ⊗ β[X] (4)
The above equation uses the following notations 2 :

X -X ⊗k def = X -X ⊗ • • • ⊗ X -X k (5)
∇ ⊗k ⊗ β def = ∇ ⊗ • • • ⊗ ∇ k ⊗β (6)
Provided that the material point of interest is away from any internal or external boundary, the application of the averaging operator of Equation ( 1) is not impacted by spatial gradients of odd order. Also, when the interaction volume is small, the spatial average µ β can be approximated from the second order derivative of the internal variable β with:

µ β [X] ≈ β[X] + 1 2 L : ∇ ⊗2 ⊗ β[X] (7)
with:

L = 1 W Vc w[X -X] X -X ⊗2 dX (8)
The above approximation has been used by [START_REF] Voyiadjis | Thermodynamic framework for coupling of non-local viscoplasticity and non-local anisotropic viscodamage for dynamic localization problems using gradient theory[END_REF] and [START_REF] De Borst | Gradient damage vs phase-field approaches for fracture: Similarities and differences[END_REF] to transform integraltype models into gradient-type models.

Spatial variance

When constructing constitutive relations, it might be necessary to provide additional information regarding the spatial distribution of an internal variable in the vicinity of a material point. For this purpose, higher order central moments of the internal variable β can be included in the list of state variables. In the following, our attention is restricted to the variance of the internal variable 2 The notation ∇ is used to denote the derivative with respect to the initial position X.

β (i.e. the second central moment). The spatial variance, which is denoted by ξ β , is obtained from:

ξ β [X] = 1 W Vc w[X -X] µ β [X] -β[X ] ⊗2 -µ β [X] -β[X] ⊗2 dX + µ β [X] -β[X] ⊗2 (9)
The spatial variance is a tensor whose rank is twice that of β. It measures the scatter associated with the spatial distribution of the internal variable around the material point of interest.

The gradient approximation of the variance is obtained by combining Equations (4) and [START_REF] Berbenni | A fast Fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals[END_REF]. Upon the condition that the point of interest is away from any internal or external boundary, the variance ξ β can be approximated as a quadratic function of the first gradient of the local internal variable β with:

ξ β [X] ≈ (β[X] ⊗ ∇) • L • (∇ ⊗ β[X]) (10)
This quadratic form is largely used in many gradient-type models (e.g. [START_REF] Forest | Strain gradient crystal plasticity : thermomechanical formulations and applications[END_REF][START_REF] Wulfinghoff | Equivalent plastic strain gradient enhancement of single crystal plasticity: theory and numerics[END_REF]).

The above result indicates that such gradient-type models can be perceived as approximations of variance-based integral-type models.

Weight function

The weight function w, which is needed for the evaluation of non-local internal variables, allows controlling the mutual interactions between material points with initial positions X and X . The weight function is conveniently expressed as the composition of two functions k and r such that:

w[X -X] = k • r[X -X] (11)
with:

r = X -X • R -2 • X -X (12)
The symmetric and positive semi-definite second order tensor R is a material parameter that provides information regarding the distances and directions over which interactions are considered. Many theories use an isotropic weight function. For this specific case, the interaction tensor R is defined from an interaction radius R (with R ≥ 0) according to:

R = R 1 (13)
A more general definition consists of using an anisotropic weight function, in which case the interaction tensor R takes the following form: where the radii R i and mutually orthogonal unit vectors r i allow defining an interaction ellipsoid. It is worth mentioning that the choice of the interaction tensor R should comply with the restrictions due to material symmetry. When the weight function w has bounded support, interactions are effective only in the region corresponding to the intersection of the interaction ellipsoid with the region V c over which the internal variable is defined (see Figure 1). If the weight function w has unbounded support, interactions are effective within the entire definition domain of the internal variable V c . The underlying idea is that different constitutive models, with different sets of internal variables, may be used when modeling the thermomechanical behavior of a body. In such a case, only the material points for which the definition of the internal variable makes sense are considered when evaluating the non-local variable.

R = i R i r i ⊗ r i (14) Vc[X ] X X Vc[X ] X Vc[X]

Field equations

In a thermomechanical context, the history of a material point with initial position X is described by the functions χ and θ defining the motion and temperature histories. Such functions allow determining the current position x and absolute temperature T with:

x = χ[X, t] and T = θ[X, t] (15)
The evolution of the current position and temperature fields is governed by some partial differential equations, which are commonly classified as either conservation or constitutive equations. In the following, conservation equations are first briefly recalled. The general framework used for the construction of constitutive relations is then detailed.

Conservation equations

For a closed system, mass conservation requires that the initial mass density , which corresponds to the mass per unit volume in the initial configuration, does not change with time [START_REF] Šilhavý | The Mechanics and Thermodynamics of Continuous Media[END_REF]:

˙ = 0 (16)
For linear and angular momenta to be conserved, the first Piola-Kirchoff stress tensor P must be such that:

P • ∇ + B = ẍ ( 17
)
P • F T = F • P T (18)
where F = ∂x/∂X is the deformation gradient tensor and B is the body force density. Also, according to the first law of thermodynamics, the evolution of the specific internal energy e is given by:

ė = 1 P : Ḟ - 1 ∇ • Q + r q + r nl ( 19
)
where Q is the heat flux density vector and r q is the specific heat supply. While the mass and momentum conservation equations display their classical form, it is important to mention that the local form of the first law of thermodynamics includes a non-conventional contribution that takes the form of a specific energy source r nl . This additional energy source, which is sometimes referred to as the non-locality residual [START_REF] Polizzotto | A thermodynamically consistent formulation of nonlocal and gradient plasticity[END_REF], results from the assumption of non-locality. Specifically, it represents the transfer of internal energy between interacting material points. For any material point, the specific source r nl must satisfy the following condition:

rnl [X] = 1 V c Vc r nl [X ]dX = 0 (20)
In the above equation, rnl denotes the average value of the non-locality residual over the region V c , which corresponds to the domain over which the constitutive model is adopted (see Figure 1). In general, this region is not the volume occupied by the whole body. The condition [START_REF] Fleck | A reformulation of strain gradient plasticity[END_REF] reflects the fact that a material point may transfer some internal energy to neighborhing material points but, when averaged over the volume V c , non-locality does not result in any net production or loss of energy. It is worth mentioning that the region over which the constitutive model is defined may coincide with the region occupied by the body of interest. For this specific but rather common situation, the average value of the non-locality residual rnl is identical for all material points, hence does not depend on the initial position. Finally, the second law of thermodynamics requires the specific entropy production source σ to be non-negative at any time and for any material point, which means that:

σ = ṡ + 1 ∇ • Q T - r q T ≥ 0 ( 21
)
where s is the specific entropy.

Constitutive equations

In a thermomechanical context, the list of state variables3 used for the development of constitutive relations includes the deformation gradient tensor F and its rate Ḟ , the absolute temperature T , its spatial gradient ∇T , a local internal variable β and its non-local counterparts µ β and ξ β . The constitutive model consists of some constitutive relations for the first Piola-Kirchoff stress tensor P , the heat flux density vector Q, the specific free energy a = e -sT and the specific entropy s. Such constitutive relations are formally written in the following form:

a = â F , Ḟ , T, ∇T, β, µ β , ξ β (22)
s = ŝ F , Ḟ , T, ∇T, β, µ β , ξ β (23)
P = P F , Ḟ , T, ∇T, β, µ β , ξ β (24)
Q = Q F , Ḟ , T, ∇T, β, µ β , ξ β (25)
Also, when some internal variables are introduced, the corresponding evolution equations must also be included in the constitutive model. These equations allow describing how the microstructure of a material point is affected by a thermomechanical process. The evolution equations take the following form:

β = β F , Ḟ , T, ∇T, β, µ β , ξ β (26)
It is worth mentioning that no evolution equation is needed for non-local internal variables. Indeed, non-local internal variables are directly connected to their local counterparts with their definition given by either Equation [START_REF] Aifantis | The physics of plastic deformation[END_REF] or Equation [START_REF] Berbenni | A fast Fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals[END_REF]. Also, though the finite strain framework is used here for the purpose of generality, the proposed framework can be used in the infinitesimal strain context by considering the infinitesimal strain tensor (and its rate) as external state variables instead of the deformation gradient tensor (and its rate) in constitutive relations. The specific dissipation source d = σT , which must be non-negative, is obtained by combining the above equation with the energy conservation equation [START_REF] Fleck | Strain gradient plasticity[END_REF] and the definition of the specific free energy:

d = 1 P : Ḟ + r nl -ȧ -s Ṫ - 1 ∇ ln[T ] • Q ≥ 0 ( 27
)
The specific free energy rate ȧ, which appears in the above inequality, is expressed as follows 4 :

ȧ = ∂a ∂F : Ḟ + ∂a ∂ Ḟ : F + ∂a ∂T Ṫ + ∂a ∂∇T •∇ Ṫ + ∂a ∂β • β+ ∂a ∂µ β • μβ + ∂a ∂ξ β ξβ (28)
Using the above relation, the expression of the specific dissipation source becomes:

d = 1 P - ∂a ∂F : Ḟ - ∂a ∂ Ḟ : F -s + ∂a ∂T Ṫ - ∂a ∂∇T • ∇ Ṫ - ∂a ∂β • β - ∂a ∂µ β • μβ - ∂a ∂ξ β ξβ + r nl - 1 ∇ ln[T ] • Q ≥ 0 ( 29
)
The above inequality should hold for any temperature and motion history.

Since the dissipation source exhibits a linear dependence with respect to F , Ṫ and ∇ Ṫ , the specific free energy must satisfy the following conditions:

∂a ∂ Ḟ = O, ∂a ∂T = -s and ∂a ∂∇T = 0 (30)
These conditions indicate that both the specific free energy and specific entropy should not depend on the deformation gradient rate and the temperature gradient, i.e.:

a = â F , T, β, µ β , ξ β (31)
s = ŝ F , T, β, µ β , ξ β (32)
The dissipation inequality therefore takes the following form:

d = 1 P ir : Ḟ + 1 b • β + 1 m • μβ + 1 z ξβ + r nl - 1 ∇ ln[T ] • Q ≥ 0 (33)
For the purpose of conciseness, the following quantities have been introduced:

P ir = P - ∂a ∂F (34) b = - ∂a ∂β (35) m = - ∂a ∂µ β (36) z = - ∂a ∂ξ β (37)
The dissipation inequality [START_REF] Mandel | Plasticit classique et viscoplasticit[END_REF] does not allow constructing the evolution equation associated with the internal variable β. Indeed, the contributions associated with the evolution of local and non-local internal variables are separated though these variables are not independent from each other. To circumvent this difficulty, it is convenient to introduce the quantities µ m , µ z and µ x such that (see Appendix A):

m • μβ = µ m • β ( 38
)
z ξβ = 2 (µ z • β -µ x ) • β (39)
To satisfy the above averaging conditions, the quantities µ m , µ z and µ x are defined according to:

µ m [X] = 1 W Vc w[X -X] m[X ] -m[X] dX + m[X] (40)
µ z [X] = 1 W Vc w[X -X] z[X ] -z[X] dX + z[X] (41)
µ x [X] = 1 W Vc w[X -X] z[X ] • µ β [X ] -z[X] • µ β [X] dX + z[X] • µ β [X] (42)
Using the properties ( 38) and ( 39), the dissipation inequality is rewritten as:

d = 1 P ir : Ḟ + 1 B • β - 1 ∇ ln[T ] • Q ≥ 0 (43)
with:

B = b + µ m + 2µ z • β -2µ x (44)
The above equation indicates that the dissipative force B driving the evolution of the internal variable β includes a local contribution as well as some non-local contributions resulting from the introduction of the spatial average and spatial variance in the list of state variables. Also, to obtain the above dissipation inequality, the following definition of the specific energy source r nl has been adopted:

r nl = m • μβ + z ξβ -µ m • β -2 (µ z • β -µ x ) • β (45)
Because of the equalities [START_REF] Niordson | Computational strain gradient crystal plasticity[END_REF] and [START_REF] Nix | Mechanical properties of thin films[END_REF], the condition ( 20) is automatically fulfilled for any material point and at any time for the body of interest. The dissipation inequality [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] is the starting point for the construction of evolution equations. Indeed, it indicates that the evolution of the state of a material point is controlled by the dissipative forces P ir , B and ∇ ln[T ]. In contrast with conventional (i.e. local) theories, the evolution of the internal variable β is driven by the dissipative force B that includes both local and non-local contributions when the proposed framework is adopted.

Numerical examples

For illustration purposes, some numerical examples using the proposed framework are presented in this section. Such examples aim at showing the possible advantages offered by the introduction of higher order non-local internal variables when constructing constitutive models for thermomechanical applications.

The following examples deal with crystalline materials. In the following, an isotropic interaction tensor (see Equation ( 13)) is used whatever the crystal structure is, mostly for simplicity reasons. The underlying assumption is that the interactions between different material points only depend on their relative distance. The function r therefore takes the following form:

r[X -X] = ||X -X|| R (46)
Also, a biquadratic function is used for k, that is5 :

k[r] = 15 16 1 -r 2 2 + ( 47
)
To solve the differential equations resulting from static equilibrium and compatibility, the spectral method [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] is used. This method allows evaluating the behavior of a volume element with periodic boundary conditions. The spectral method has been implemented in the context of finite strains according to the iterative procedure presented by Eisenlohr et al. [START_REF] Eisenlohr | A spectral method solution to crystal elasto-viscoplasticity at finite strains[END_REF]. The sole modification concerns the evaluation of non-local internal variables. Specifically, non-local variables are obtained from their local counterparts once constitutive relations have been integrated from a discrete convolution operation. Non-local variables are then used as an input for the integration of constitutive relations during the next iteration. This procedure is repeated until convergence is achieved. Also, for simplicity, the temperature field is assumed to be uniform and constant, which corresponds to the specific case of isothermal conditions.

Example 1: Damage localization

Constitutive relations For the present application, the construction of constitutive relations relies on the general framework of continuum damage mechanics [START_REF] Lemaitre | A course on damage mechanics[END_REF]. Specifically, a scalar internal variable β is introduced to represent the degradation of mechanical resistance of a material point. The damage variable is equal to one for a fully damaged material point. At the opposite, a zero value corresponds to the absence of damage. Assuming that fracture is controlled by brittle damage, the specific free energy a is decomposed into three contributions according to:

a = a e + a d + a th (48)
The elastic strain energy contribution a e is given by:

a e = 1 2 E : C : E ( 49
)
where C is the stiffness tensor and E is the Green strain tensor:

E = 1 2 F T • F -1 (50)
To consider closure effects, the impact of damage on stiffness properties is described according to:

C = (1 -β) C 0 if tr[E] > 0 (1 -β) C 0 + β P s : C 0 : P s if tr[E] ≤ 0 ( 51
)
where C 0 is the initial stiffness tensor. The above equation assumes that the spherical contribution to the stiffness tensor is recovered when the spherical strain is negative. For this purpose, Equation ( 51) uses the spherical projection tensor P s given by:

P s = 1 3 1 ⊗ 1 (52)
For the contribution of damage to free energy a d , different options (a-d) are explored:

a d = 1 A 2β -β 2 (53a)
a d = 1 A 2β -β 2 + 1 G 2 (β -µ β ) 2 (53b)
a d = 1 A 2β -β 2 + 1 G 2 (β -µ β ) 2 + ξ β (53c) a d = 1 A 2β -β 2 + 1 G 2 ξ β (53d)
In the above equations, A and G are material parameters controlling respectively the local and non-local contributions of damage to free energy. While the local contribution takes the form suggested by Wu [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF], different propositions regarding the non-local contribution are considered. The first option (a) corresponds to a purely local damage model. For the other options (b-d), the non-local contribution, which vanishes for a uniform damage field, is estimated from the spatial average µ β and/or the spatial variance ξ β of the damage variable. The case where the spatial average is the sole non-local internal variable (b) corresponds to the standard non-local approach. The options involving the spatial variance (c and d), alone or together with the spatial average, are presented to demonstrate the possibilities offered by the proposed framework. From a thermodynamic point of view, it is worth mentioning that, in comparison with the local model (a), an affine distribution of the damage variable leads to an increase of free energy for variance-based models (c or d). At the opposite, an affine distribution of the damage variable does not provide any additional contribution to free energy when the standard average-based model (b) is adopted. The thermal contribution to free energy, which is denoted by a th , does not need to be specified because isothermal conditions are assumed.

The state equations are established by differentiating the specific free energy with respect to the different state variables. In the absence of viscous contribution to the stress state (i.e. P ir = O), the first Piola-Kirchoff stress tensor is:

P = ∂a ∂F = F • (C : E) (54)
The differentiation of free energy with respect to the damage variable leads to the expression of the energy restitution rate b:

b = - ∂a ∂β = - 1 2 E : ∂C ∂β : E + 2A (β -1) (55a or 55d) b = - ∂a ∂β = - 1 2 E : ∂C ∂β : E + 2A (β -1) + G (µ β -β) (55b or 55c)
The driving force associated with the spatial average of the damage variable is given by:

m = - ∂a ∂µ β = 0 (56a or 56d) m = - ∂a ∂µ β = G (β -µ β ) (56b or 56c)
In a similar fashion, the driving force associated with the spatial variance is:

z = - ∂a ∂ξ β = 0 (57a or 57b) z = - ∂a ∂ξ β = - G 2 (57c or 57d)
Using the above state equations, the dissipation equality reduces to:

d = 1 B β - 1 ∇ ln[T ] • Q ≥ 0 (58)
with:

B = - 1 2 E : ∂C ∂β : E + 2A (β -1) (59a) B = - 1 2 E : ∂C ∂β : E + 2A (β -1) + G (2µ β -β -µ ββ ) (59b) B = - 1 2 E : ∂C ∂β : E + 2A (β -1) + 2G (µ β -β) (59c) B = - 1 2 E : ∂C ∂β : E + 2A (β -1) + G (µ ββ -β) (59d)
When the formulation is purely local (a), the dissipative force driving the development of damage does not contain any contribution limiting damage localization. For the models that use the spatial variance (c and d), the nonlocal contribution to the dissipative force has a clear meaning. Whatever the damage field is, the non-local contribution participates in damage diffusion by reducing the difference between the damage variable β and its spatial (single or double) average (µ β or µ ββ ). This non-local contribution can therefore be interpreted as a localization limiter. For the standard non-local approach (b), which considers the damage variable and its spatial average as internal variables, the non-local contribution also promotes damage diffusion. This nonlocal contribution tends to reduce the difference between the local damage variable β and a non-local variable (2µ βµ ββ ).

From a computational perspective, the mixed model (c), which uses both the variance and the average, requires less computational efforts than the standard non-local approach (b). Indeed, for the evaluation of the non-local contribution to the dissipative force, only one averaging operation needs to be performed. At the opposite, when the standard non-local approach (b) is adopted, both the single and the double-averaging operations must be applied to the damage variable to evaluate the non-local contribution to the dissipative force driving damage development.

In the context of standard materials, the evolution equations associated with the damage variable and the heat flux density vector are derived from a dissipation potential. Because of the restriction to isothermal conditions, the constitutive equation for the heat flux vector does not need to be established here. The dependence of the dissipation potential to the temperature logarithm gradient is therefore ignored. In the present case, the dissipation potential φ displays the following form:

φ = K 2 B + K 2 (1 -β) ( 60
)
where K is a viscosity parameter. The differentiation of the dissipation potential with respect to the dissipative force B leads to the following evolution equation for the damage variable:

β = ∂φ ∂B = B + K (1 -β) (61)
It is worth mentioning that the above evolution equation only allows damage growth (i.e. healing is not allowed).

Parameters and loading conditions

The volume element used for this application is shown in Figure 2. It is composed of 64 equiaxed randomly oriented grains, with an average grain size of 25 µm. To evaluate the impact of spatial discretization on numerical results, resolutions ranging from 32 3 to 128 3 voxels have been used for the application of the spectral method. Since the size of the volume element is 100 µm, the voxel size ranges from 1.6 µm to 0.5 µm. Also, to investigate the development of brittle damage, periodic boundary conditions corresponding to uniaxial tension with a constant axial strain rate of 10 -3 s -1 are prescribed to the volume element.

The material parameters used for the numerical simulations are listed in Table 1. An interaction radius R of 2 µm is used for the computation of nonlocal damage variables.

Results

The evolution of the macroscopic axial stress as a function of the macroscopic axial strain obtained for a resolution of 128 3 is plotted in Figure 3. The macroscopic behavior is qualitatively similar for the different models. Specifically, the axial stress first increases and drops rapidly when damage develops significantly. When the formulation is purely local, the maximum axial stress is lower than that obtained with non-local models. Indeed, when the non-local contribution is adopted, the development of damage is accompanied by a diffusion process that reduces the impact of the singularities resulting from damage development.

The total dissipated energy and peak stress are plotted as a function of resolution in Figure 4. According to the numerical results, some differences exist regarding mesh dependency. Specifically, in contrast with non-local formulations (b-d), the results obtained with the local formulation (a) strongly depend on mesh resolution. This difference is attributed to the localization of damage within a few voxels. Strong spatial damage localization, hence mesh dependency, is largely reduced when non-locality is introduced. This is particularly true for variance-based models (c or d) that exhibit the lowest dependency with respect to mesh resolution. The difference with the average-based model (b) is attributed to the fact that the non-local contribution to the dissipative force, which drives damage diffusion, is higher when the spatial variance is introduced.

Example 2: Polycrystalline plasticity

Constitutive relations To describe the elastic-viscoplastic behavior of a crystalline material point, the deformation gradient tensor is decomposed into elastic (superscript e) and plastic (superscript p) contributions according to:

F = F e • F p (62)
In addition to the plastic contribution to the deformation gradient tensor F p , the list of internal state variables used for the present constitutive model includes some isotropic hardening variables. In the following, the hardening variable attached to the sth slip system is denoted by β s and the corresponding spatial average and variance are respectively denoted by µ s β and ξ s β . The proposed model is similar in spirit to the non-local softening plasticity models of Vermeer and Brinkgreve [START_REF] Vermeer | A new effective non-local strain measure for softening plasticity. Localization and bifurcation theory for soils and rocks[END_REF] and Strömberg and Ristinmaa [START_REF] Strömberg | FE formulation of a nonlocal plasticity theory[END_REF] in the sense that non-locality is introduced via isotropic hardening variables. This approach contrasts with common strain gradient plasticity models (e.g. that consider the spatial gradients of plastic strains, rather than hardening variables.

The proposed set of state variables allows decomposing the specific free energy a into three contributions corresponding to elastic strain energy a e , defect energy a d and thermal energy a th such that:

a = a e + a d + a th (63)
The elastic contribution to free energy is given by:

a e = 1 2 E e : C : E e ( 64
)
where C is the stiffness tensor and E e is the elastic Green strain tensor:

E e = 1 2 F T e • F e -1 (65)
As for the damage model, different options (a-d) are explored for the contribution of defects to free energy:

a d = 1 2 s β s t H st β t (66a)
a d = 1 2 s β s t H st β t + 1 2 G s (β s -µ s β ) 2 (66b)
a d = 1 2 s β s t H st β t + 1 2 G s (β s -µ s β ) 2 + ξ s β (66c) a d = 1 2 s β s t H st β t + 1 2 G s ξ s β (66d)
In the above equation, H is the isotropic hardening matrix and G is a material parameter controlling the non-local contribution to free energy. The first option (a) corresponds to a purely local formulation while the other options (b-d) are non-local. For the latter formulations, the non-local contribution is either evaluated from the spatial averages (b and c) and/or the spatial variance (c and d) of the isotropic hardening variables. Since isothermal conditions are considered, the thermal contribution to free energy, which solely depends on temperature, does not need to be specified.

In the absence of viscous stress, the first Piola-Kirchoff stress tensor is given by:

P = ∂a ∂F (67) = F e • (C : E e ) • F -T p (68)
Depending on the retained option, one obtains the following state equations for the driving force associated with the hardening variable for the sth system:

b s = - ∂a ∂β s = - t H st β t (69a or 69d) b s = - ∂a ∂β s = - t H st + G(µ s β -β s ) (69b or 69c)
For each slip system s, the driving force associated with the spatial average of the hardening variable is:

m s = - ∂a ∂µ s β = 0 (70a or 70d) m s = - ∂a ∂µ s β = G(β s -µ s β ) (70b or 70c)
For the spatial variance of the hardening variable, the corresponding driving force is given by:

z s = - ∂a ∂ξ s β = 0 (71a or 71b) z s = - ∂a ∂ξ s β = - G 2 (71c or 71d)
Using the proposed framework, the expression of the specific dissipation source d reduces to:

d = 1 Σ : L p + 1 s B s βs - 1 ∇ ln T • Q (72)
where L p = Ḟ p • F -1 p is the plastic contribution to the velocity gradient. The above equation indicates that plastic flow is controlled by the Mandel stress tensor Σ [START_REF] Mandel | Plasticit classique et viscoplasticit[END_REF], whose expression is:

Σ = F T e • P • F T p (73)
The stress-like quantity B s is the conjugate force to the hardening variable β s . The application of relation [START_REF] Pijaudier-Cabot | Nonlocal damage model with evolving internal length[END_REF] leads to:

B s = - t H st β t (74a) B s = G 2µ s β -µ s ββ -β s - t H st β t (74b) B s = 2G µ s β -β s - t H st β t (74c) B s = G µ s ββ -β s - t H st β t (74d)
In the context of crystal plasticity, each slip system s is defined from two unit vectors: the slip plane normal n s and the slip direction m s . When plastic flow is the sole result of crystallographic slip, the plastic contribution to the velocity gradient L p is obtained from:

L p = s m s ⊗ n s ν s (75)
where ν s is the plastic shear strain rate for each slip system. To consider the size-dependence of the hardening rate, the critical resolved shear stress r s , which represents the resistance to plastic flow for the sth slip system, is defined according to:

r s = -B s (76a) r s = -B s 1 + L (β s -µ s β ) 2 (76b) r s = -B s 1 + L (β s -µ s β ) 2 + ξ s β (76c) r s = -B s 1 + L ξ s β ( 76d
)
where L is a material parameter. In contrast with local crystal-plasticity based models (a), the resistance to plastic flow includes a non-local contribution (bd) that results from inhomogeneous hardening. Such a contribution vanishes when the hardening variable is (at least locally) uniformly distributed, in which case the linear isotropic hardening rule is retrieved. To determine whether the conditions for plastic flow are met or not, it is convenient to introduce a yield function f s for each slip system s such that:

f s = |τ s | -r s (77)
with:

τ s = m s • Σ • n s (78)
The dissipation potential φ used to obtain the evolution equations takes the following form:

φ = 1 s K M + 1 f s K M +1 (79)
where K and M are some viscosity parameters. The plastic shear strain rate is thus given for each slip system s by:

ν s = ∂φ ∂τ s = f s K M sign [τ s ] (80)
Also, the evolution equations for the hardening variables obtained for the different options are: Parameters and loading conditions The volume element used to investigate the elasto-plastic behavior is the same as for the previous example (see Figure 2). To study size effects, three different grain sizes have been considered: 12.5 µm, 25 µm and 50 µm. The results obtained for an infinite grain size, which correspond to the local formulation, are also provided. The boundary conditions prescribed to the volume element correspond to a simple shear test with a constant shear strain rate of 10 -2 s -1 . The material parameters used for performing the numerical simulations are listed in Table 2. For the present application, the hardening variables have an initial value of 1.25%, which corresponds to an initial critical resolved shear stress of 30 MPa.

βs = ∂φ ∂B s = - ∂φ ∂r s ∂r s ∂B s = |ν s | (81a) βs = ∂φ ∂B s = - ∂φ ∂r s ∂r s ∂B s = |ν s | 1 + L (β s -µ s β ) 2 (81b) βs = ∂φ ∂B s = - ∂φ ∂r s ∂r s ∂B s = |ν s | 1 + L (β s -µ s β ) 2 + ξ s β (81c) βs = ∂φ ∂B s = - ∂φ ∂r s ∂r s ∂B s = |ν s | 1 + L ξ s β (81d)

Results

The evolution of the macroscopic tangential stress as a function of the macroscopic shear strain is plotted in Figure 5 for the different grain sizes.

According to the numerical results, the reduction of grain size leads to an increase of both the yield stress and the hardening rate. In the present context, this size effect is attributed to the inhomogeneous distribution of the hardening variable, which affects the resistance to plastic flow. It is worth noticing that, in comparison with the average-based model (b), the impact of grain size on the elasto-plastic behavior is more significant when the spatial variance (c or d) is introduced. Since the spatial average does not much contribute to size effects, the mixed (c) and variance-based (d) models provide similar results. The yield stress obtained from the non-local models (b, c or d) is plotted as a function of the inverse of the square root of the average grain size in Figure 6. According to the results, the yield stress tends toward a minimum value when the grain size approaches infinity, which is consistent with the Hall-Petch equation. This effect is due to the fact that the non-local contribution to hardening is negligible for large grain sizes. When the grain size is reduced, the dependence of the yield stress with respect to the average grain size is close to that given by the Hall-Petch equation when either the mixed (c) or variance-based (d) models are adopted. While the results obtained with the average-based model (b) are size-dependent, the role of grain size on the development of plasticity is not fully consistent with the Hall-Petch equation. Such results highlight the possible benefits of considering the spatial variance of hardening variables in the description of size effects.

It should be mentioned that, in contrast with gradient-type models (e.g. [START_REF] Haouala | Simulation of the Hall-Petch effect in FCC polycrystals by means of strain gradient crystal plasticity and FFT homogenization[END_REF], [START_REF] Berbenni | A fast Fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals[END_REF]), integral-type models are not well suited for very small grain sizes, i.e. inferior to the interaction radius R. Indeed, when the interaction volume is large in comparison with the grain size, the spatial average and spatial variance reach asymptotic values in the sense that they are not much affected by further grain size reduction. As a result, the impact of grain size on the yield stress is not correctly depicted for very small grain sizes. Additional developments are therefore needed to improve the description of size effects over a wide range of grain sizes with (average-or variance-based) integral-type models.

The spatial distribution of the cumulated plastic strain β eq = s β s obtained at the end of the shear test with the variance-based model (d) is shown in Figure 7 for different grain sizes. Plastic strain localization is less pronounced for the finest microstructure due to non-local effects. At the opposite, for the coarsest microstructure, significant plastic strains have been accumulated near grain boundaries as a result of internal stresses. In this latter case, non-local effects have a negligible impact on the spatial distribution of plastic deformation. The norm of ∇ × F p , which provides a quantitative measure of plastic deformation incompatibility, is also plotted in Figure 8 for the different grain sizes. This quantity is related to the density of GND [START_REF] Arsenlis | Crystallographic aspects of geometrically-necessary and statistically-stored dislocation density[END_REF]. According to the results, though the proposed model does not explicitly use GND densities, the reduction of grain size leads to higher plastic deformation incompatibilities, hence higher GND densities. Also, GND preferably accumulate near grain boundaries, which are the regions where significant plastic deformation incompatibilities are observed.

To illustrate the energetic implications of the present model, the evolution of the average internal energy density as a function of the macroscopic shear strain is plotted in Figure 9 for the variance-based model (d). According to the results, the grain size reduction leads to an increase of internal energy that is attributed to heterogeneous hardening.

Conclusions

In the present work, a general framework has been proposed to construct variance-based non-local constitutive models. In contrast with most integraltype constitutive models, which only use the spatial averages, the spatial variances of internal variables are treated as additional state variables. The proposed framework relies on continuum thermodynamics to construct the state equations and the evolution equations, which incorporate some non-local contributions. Such contributions allow considering the impact of the spatial distribution of internal variables on the thermomechanical behavior of solid materials without introducing additional equilibrium equations and boundary conditions. For the purpose of illustration, some numerical examples have been presented. According to the numerical results, the proposed framework can be used to circumvent the difficulties associated with excessive spatial localization or to consider size effects. Further work should focus on the construction of more advanced constitutive models within the proposed framework and their validation from experimental observations.

A Averaging properties of non-local variables

The construction of constitutive relations uses some specific averaging properties of nonlocal variables. These properties are briefly demonstrated here. Specifically, for any local internal variable β, the time derivative of the spatial average µ β is given by:

μβ [X] = 1 W Vc w[X -X] β[X ] -β[X] dX + β[X] (82)
For any tensor field m, whose rank is that of β, the above expression leads to:

m • μβ = 1 Vc 1 W Vc Vc w[X -X]m[X] • β[X ] -β[X] dX dX + 1 Vc Vc m[X] • β[X]dX (

Fig. 1

 1 Fig. 1 Schematic diagram showing the interaction volume for different material points.

Fig. 2

 2 Fig. 2 Polycrystalline volume element used for the application of the spectral method.

Fig. 3

 3 Fig.3Stress-strain diagrams obtained for a uniaxial tension test with the different damage models.

Fig. 4

 4 Fig.4Evolution of the total dissipated energy (top) and peak stress (bottom) as a function of the mesh resolution.

Fig. 5 Fig. 6

 56 Fig. 5 Stress-strain diagrams obtained for a simple shear test for different grain sizes with the average-based (top), mixed (middle) and variance-based (bottom) non-local models.

Fig. 7

 7 Fig. 7 Spatial distributions of the cumulated plastic strain obtained at the end of a simple shear test for different grain sizes with the variance-based non-local model.

Fig. 8

 8 Fig. 8 Spatial distributions of ||∇ × F p|| obtained at the end of a simple shear test for different grain sizes with the variance-based non-local model.

Table 1

 1 Material parameters for the brittle damage model. The Laue class corresponding to the crystalline solid is 6/mmm (hexagonal system). Only independent elastic constants are given.

			Stiffness
	C 1111	C 1122	C 3333	C 1133	C 1212
	163.7 GPa	36.4 GPa	63.5 GPa	53.0 GPa 38.8 GPa
		Damage		Radius
	A	G	K	R
	10 -1 MPa	1 MPa	5 × 10 -3 MPa.s	2 µm

Table 2

 2 Material parameters for the crystal plasticity model. The Laue class corresponding to the crystalline solid is m 3m (cubic system). Only independent elastic constants are given.

		Stiffness		Viscoplasticity
	C 1111	C 1122	C 1212	K	M
	168.3 GPa 122.1 GPa 75.7 GPa 10 MPa	20
		Hardening		Radius
	Hss	Hst	G	L	R
	200 MPa	200 MPa	10 GPa	1000	2 µm

For the purpose of conciseness, a single internal variable is considered here. The extension of the proposed framework for multiple internal variables is straightforward.

The inner product of two tensors whose rank is that of β (e.g. µ β ) is denoted by •. In a similar fashion, the inner product of two tensors whose rank is twice that of β (e.g. ξ β ) is denoted by .

The positive part of a variable (say x) is given by x + = (x + |x|)/2.

Since the weight function w is even, the above expression becomes:

Also, the time derivative of the spatial variance ξ β is:

For any symmetric tensor z, whose rank is twice that of β, one obtains that:

The evenness of the weight function allows re-writing the above expression as follows:

For the purpose of conciseness, it is convenient to introduce the tensor field x such that:

Using the above definition, expression (88) becomes: