Louis Toulouse
email: louis.toulouse@etu.univ-littoral.fr

Patrick Mucci

Gautier Thierry Pezé

Influence of grade of

Keywords: children, adolescents, obesity, treadmill, VO 2max

come

Influence of grade of obesity on the achievement of VO2max using an incremental treadmill test in youths

The purpose of this study was to analyze the influence of grade of obesity on the probability of achieving a VO 2 plateau and threshold secondary criteria for verifying VO 2max during a treadmill walk test in youths with obesity. Therefore, 72 youths with obesity (aged 8-16) performed an incremental treadmill walk test to exhaustion during which oxygen uptake (VO 2 ), minute ventilation (VE), heart rate (HR) and rating of perceived exertion were continuously measured. HR corresponding to a "hard" level of perceived exertion was reported and expressed as a percentage of the predicted HR max . The rate of achievement of criteria for validation VO 2max (VO 2 plateau; HR>95% theoretical HR max ; RER>1.0; rating of perceived exertion ≥ "hard") was compared between participants with grade I and grade II obesity. 37% of the participants achieved a VO 2 plateau and 23% achieved both an HR>95% and RER >1.0. Youths with grade II obesity had lower minute ventilation (p<0.01) tended to be more likely to reach an HR>95% (OR = 0.33 ; P=0.06) and a "hard" rating of perceived exertion than grade I (OR = 4.5 ; P=0.07). However, there was no influence of grade of obesity on the achievement of VO 2 plateau, and RER>1.0. A higher grade of obesity was related to a lower peak in minute ventilation normalized by fat-free mass, lower peak heart rates and premature reaching of a "hard" rating of perceived exertion.

Therefore, youths with higher grades of obesity might be less likely to display maximal heart rate during an incremental treadmill test, probably due to higher ventilatory constraints and higher ratings of perceived exertion.

INTRODUCTION

Cardiorespiratory fitness is one of the best predictors of all-cause mortality in adults [START_REF] Barry | The Joint Association of Fitness and Fatness on Cardiovascular Disease Mortality: A Meta-Analysis[END_REF] , especially when it is normalized by fat-free mass [START_REF] Imboden | Cardiorespiratory Fitness Normalized to Fat-Free Mass and Mortality Risk[END_REF] . Cardiorespiratory fitness is assessed by the maximal oxygen uptake (VO 2max ), which is the highest value of VO 2 clinically attainable during a progressive cycling or treadmill test conducted to exhaustion [START_REF] Poole | Determinants of oxygen uptake. Implications for exercise testing[END_REF] . Since walking is more representative of free-living physical activities, treadmill walk tests are mostly used to assess the VO 2max in youth with obesity [START_REF] Gutin | Effects of exercise intensity on cardiovascular fitness, total body composition, and visceral adiposity of obese adolescents[END_REF][START_REF] Marinov | Ventilatory efficiency and rate of perceived exertion in obese and non-obese children performing standardized exercise[END_REF][START_REF] Breithaupt | Using the oxygen uptake efficiency slope as an indicator of cardiorespiratory fitness in the obese pediatric population[END_REF] . During these tests, the speed is constant and usually between 3.5 and 5 km.h -1 , which is close to the preferred walking speed in this population [START_REF] Peyrot | Why does walking economy improve after weight loss in obese adolescents?[END_REF] . The intensity is increased by increasing the slope, amounting to 1 or 2% per minute [START_REF] Gutin | Effects of exercise intensity on cardiovascular fitness, total body composition, and visceral adiposity of obese adolescents[END_REF][START_REF] Marinov | Ventilatory efficiency and rate of perceived exertion in obese and non-obese children performing standardized exercise[END_REF] . To verify the achievement of VO 2max (i.e., a maximal aerobic effort) during this type of test, it was proposed that youths with obesity should reach at least one of the three following physiological criteria: 1) the presence of a VO 2 plateau; 2) a heart rate > 200bpm or 95%HR max ; 3) a respiratory exchange ratio > 1.0. According to these criteria, Gutin et al. found that 66 out of 80 participants (~ 83%) reached VO 2max [START_REF] Gutin | Effects of exercise intensity on cardiovascular fitness, total body composition, and visceral adiposity of obese adolescents[END_REF] . Using the same protocol, a study 6 found that 29% of adolescents with grade II obesity achieved a VO 2 plateau. However, the rates of achievement of secondary threshold criteria were not specified. Finally, one study 5 did not specify the criteria [START_REF] Gutin | Effects of exercise intensity on cardiovascular fitness, total body composition, and visceral adiposity of obese adolescents[END_REF] used to verify the achievement of VO 2max , which does not allow to determine the extent to which their treadmill protocol elicited a maximal aerobic effort. Since VO 2max remains the gold standard for monitoring aerobic fitness in youths with obesity [START_REF] Thivel | High-intensity interval training in overweight and obese children and adolescents: systematic review and metaanalysis[END_REF] , it is necessary to ensure that a maximal aerobic effort has been made during the test.

However, there is still a lack of transparency regarding the criteria used to ensure that a maximal aerobic effort has been made during treadmill walk tests in youths with obesity. Furthermore, no study has investigated the influence of grade of obesity on the probability of achieving VO 2max during this type of protocol. Thus, further studies are needed to better understand the mechanisms associated with the non-achievement of a VO 2 plateau and secondary threshold criteria for verifying VO 2max in youth with obesity.

Therefore, this study aimed to assess the influence of grade of obesity on the probability of achieving a VO 2 plateau and secondary threshold criteria for verifying VO 2max during an incremental treadmill walk test in youth with obesity.

METHODS

Participants

72 youths with obesity were included in this study. The inclusion criteria were an age between 8 and 16 years old and a BMI superior to the 85 th percentile of the IOTF standards for BMI, defining grade I obesity in youths [START_REF] Cole | Extended international (IOTF) body mass index cut-offs for thinness, overweight and obesity[END_REF] . These youths are oriented by their practitioner towards obesity management networks (including the hospital center and associations) where they are offered to follow care adapted to their needs. During their inclusion appointment to the program, the youths participants and their parents received an oral and written explanation of the experimental procedure and the objectives of the study. For youths, the information notice was written in a language adapted to their level of understanding. After a reflection time of at least one week following the inclusion visit, the youths and the voluntary parents signed a written consent summarizing the information relating to the study.

General procedure

The experiments took place in the department of pulmonary medicine of the hospital center in the presence of a specialist doctor. Participants were called in the morning between 10 a.m. and 12 p.m. or in the afternoon between 3 p.m. and 7 p.m., at least two hours after their last food intake to avoid digestive discomfort during the test. Parents were asked to make sure that their youths did not exercise intensively during the 48 hours preceding the appointment. Before the experiments, the participants were measured in standing and seated positions using a stadiometer. Then, the body weight and composition were measured using bioimpedance-analysis (Tanita DC-360, Japan).

These measurements allowed us to obtain an estimate of their body composition (fatfree mass and fat mass) and to calculate their BMI z-score based on the IOTF [START_REF] Cole | Extended international (IOTF) body mass index cut-offs for thinness, overweight and obesity[END_REF] and WHO 10 standards. Participants who had a BMI z-score between 1 and 2 (not included) according to the IOTF standards were categorized into "grade I obesity", whereas those who had a BMI z-score > 2 were categorized into "grade II obesity". From the anthropometrical measurements, the Age of Peak Height Velocityage of PHV (also called "Maturity Offset") was estimated using the following sex-specific equations provided by Mirwald et al. [START_REF] Mirwald | An assessment of maturity from anthropometric measurements[END_REF] :

Then, age from the peak height velocity (age from PHV) was estimated by subtracting the estimated age of PHV by the chronological age. A negative value indicates that the PHV has not occurred yet [START_REF] Mirwald | An assessment of maturity from anthropometric measurements[END_REF] . The maturiy offset and age from PHV were used as an indicator of the maturity of the participants.

Incremental exercise

The incremental treadmill protocol consisted of performing a 3-minute warm-up at 0%

slope and a speed between 3.0 to 4.5 km.h -1 so that the participant reached a heart-rate (HR) steady-state at 60 ± 5% (110-130 bpm) of the predicted maximum HR estimated with the formula of Tanaka [START_REF] Tanaka | Age-predicted maximal heart rate revisited[END_REF] [HRmax = 208 -0.7 * age]. Normalizing the starting speed by the HR was intended to ensure an equitable starting intensity between individuals and optimal exercise duration (> 5 minutes excluding warm-up) for all participants [START_REF] Midgley | Criteria for determination of maximal oxygen uptake: a brief critique and recommendations for future research[END_REF] . During the first 5 minutes, the slope increased by 2% per minute and the speed was kept constant. From the sixth minute, the slope was maintained at 10% and the speed increased by 0.3 km.h -1 for a starting speed between 3.0 km.h -1 and 3.7 km.h -1 or by 0.4 km.h -1 for a starting speed between 3.7 and 4.5 km.h -1 . This protocol was selected because high treadmill grades would be associated with an increased failure to achieve VO 2max [START_REF] Midgley | Challenging a dogma of exercise physiology: does an incremental exercise test for valid VO 2 max determination really need to last between 8 and 12 minutes?[END_REF] . The test ended when the participant declared that he was no longer able to continue the test despite the encouragement of the experimenter.

Measurement and processing of data

During exercise, HR was measured using a recently validated optical heart rate monitor (Polar OH1, Finland) [START_REF] Hermand | Validation of a Photoplethysmographic Heart Rate Monitor: Polar OH1[END_REF] positioned at the left wrist. The peak heart rate (HR peak ) observed during the test was the highest value averaged on 3 points. The VO 2 was measured breath-by-breath using a portable respiratory gas analyzer (Cosmed K4b², Italy) calibrated before the exercises following the manufacturer's recommendations.

The VO 2peak and peak in minute ventilation (VE peakBM ) were reported as the highest values over a 30s moving-averaged interval, as recently recommended in youths [START_REF] Blanchard | Determination of peak cardiorespiratory fitness parameters in children: which averaging method should we use?[END_REF] .

Rating of perceived exertion (RPE) was continuously assessed using the CR-10 Borg scale previously used in youths with obesity [START_REF] Marinov | Ventilatory efficiency and rate of perceived exertion in obese and non-obese children performing standardized exercise[END_REF] . This scale describes the difficulties of breathing from level 1 ("not at all") to 10 ("maximal"). The HR measured when the participant declared a level of 5/10 on the scale -corresponding to the "hard" item (HR RPEhard ) -was reported and expressed as a percentage of the individual predicted HR max . We chose this indicator to get information about the level of physiological stress that each participant was able to tolerate before the effort was "hard", a term that has more clinically significance in youths [START_REF] Williams | CERT: a perceived exertion scale for young children[END_REF] .

Criteria for verifying VO 2max

Following previous studies in youths with obesity [START_REF] Gutin | Effects of exercise intensity on cardiovascular fitness, total body composition, and visceral adiposity of obese adolescents[END_REF][START_REF] Breithaupt | Using the oxygen uptake efficiency slope as an indicator of cardiorespiratory fitness in the obese pediatric population[END_REF] , the present cardiorespiratory-based criteria for verifying VO 2max were used: primary criterion = presence of a VO 2 plateau; secondary threshold criteria =HR>95% of theoretical maximal HR and RER > 1.0.

Participants were considered as VO 2max achievers when at least the VO 2 plateau or both the secondary threshold criteria were met. The presence of a VO 2 plateau was characterized when the increase in VO 2 during the last completed minute of effort was strictly less than 50% of the increase observable during the penultimate minute of effort.

This method based on the individual workload/VO 2 relationship was recently recommanded [START_REF] Niemeyer | Oxygen uptake plateau: calculation artifact or physiological reality?[END_REF] and used to detect VO 2 plateaus in adults with obesity [START_REF] Wood | Vo 2max in overweight and obese adults: do they meet the threshold criteria?[END_REF] . A fourth criterion based on the achievement of an RPE ≥ 5 (item "hard") was also used. This criterion was used as a measurement of the psychological tolerance to the treadmill test.

Conversely, the criteria based on volitional exhaustion was not used because the latter is determined rather by the experimenter than the participant, limiting the possibilities of inter-study comparisons.

Statistical analysis

The data were analyzed on the R software (version 3.6.2). The normality and homogeneity of the distributions were checked respectively with the Shapiro-Wilk test and the F test of Snedecor. Firstly, the clinical characteristics and peak exercise data of participants were compared between those who had grade I and those who had grade II obesity using Student t-tests and Cohen's d. Secondly, the number and rate of participants who met the criteria for verifying VO 2max were calculated and reported according to the grade of obesity. Chi-Square tests were performed to analyze the differences in the rate of achievement of each criterion according to the grade of obesity. If significant differences were found, logistic regression analysis was performed to examine the association between the grade of obesity and the achievement of each criterion. The significance threshold was set at P<0.05.

Results

Participants

Discussion

The main purpose of this study was to study the influence of the grade of obesity on the probability to meet the primary and secondary threshold criteria for verifying VO 2max during a treadmill test. The secondary purpose was to measure the cardiorespiratory and perceptual responses to this test in the youth with obesity, which is an under-evaluated population. The main findings of this study are that 1) youths with grade II obesity were less likely to display an HR>95%HR max and 2) youths with grade II obesity were more likely to indicate an RPE reflecting a "hard" level of effort than grade I. Furthermore, youth with grade II obesity displayed lower peaks in minute ventilation normalized by the fat-free mass than grade I.

Adults suffering from chronic diseases have a lower ability to achieve VO 2max during a progressive test conducted until exhaustion [START_REF] Poole | Measurement of the maximum oxygen upta e o 2max o 2peak is no longer acceptable[END_REF] . However, in youths with obesity, current studies provide little data regarding the achievement of criteria for verifying VO 2max [START_REF] Gutin | Effects of exercise intensity on cardiovascular fitness, total body composition, and visceral adiposity of obese adolescents[END_REF][START_REF] Marinov | Ventilatory efficiency and rate of perceived exertion in obese and non-obese children performing standardized exercise[END_REF][START_REF] Breithaupt | Using the oxygen uptake efficiency slope as an indicator of cardiorespiratory fitness in the obese pediatric population[END_REF] .

Therefore, knowledge about the individual's factors associated with the probability of achieving VO 2max remains limited. Although no statistically significant, the present results in a relatively large sample of youth with obesity (n=72) suggest that the grade of obesity may be associated with a lower probability of meet an HR>95%HR max and a higher probability of meet a "hard" level of RPE. Thus, youths with higher grades of obesity might perceive the test as more difficult. This result may be of significance since perceptual responses play a major role during a progressive test conducted until exhaustion [START_REF] Poole | Measurement of the maximum oxygen upta e o 2max o 2peak is no longer acceptable[END_REF] . The present data would be supporting those of our predecessors 5 who found that RPE was higher in grade II than grade I obesity at the end of a treadmill test.

They hypothesized that these differences could be due to an increased aerobic energy cost, a greater CO 2 turnover, and a lower anaerobic threshold in grade II obesity. In the present study, the peak in minute ventilation normalized by body mass was significantly lower in youths with grade II than grade I obesity, and this differences persisted even after removing differences in body composition. This could mean that ventilatory adaptations to exercise were lower in grade II vs. grade I. Interestingly, it has been previsouly suggested that lower cardiorespiratory fitness in youths with obesity could be partly explained by ventilatory factors [START_REF] Mendelson | Low Cardiorespiratory Fitness is Partially Linked to Ventilatory Factors in Obese Adolescents[END_REF] . The excess adipose tissue at the thoracic level would constrain them to breathe at low lung volumes, increasing the work of respiratory muscles to the detriment of locomotor muscles [START_REF] Mendelson | Low Cardiorespiratory Fitness is Partially Linked to Ventilatory Factors in Obese Adolescents[END_REF] . Nevertheless, ventilatory factors would have a little impact on cardiorespiratory fitness in youths [START_REF] Mendelson | Low Cardiorespiratory Fitness is Partially Linked to Ventilatory Factors in Obese Adolescents[END_REF] . Accordingly, no differences were found between grade I and grade II regarding VO 2peak after removing differences in body mass and composition. Therefore, the present data would be in accordance with the previous hypothesis 5 that ventilatory factors are involved in the higher feeling of difficulty of youths with grade II obesity, but has limited influence on VO 2peak [START_REF] Mendelson | Low Cardiorespiratory Fitness is Partially Linked to Ventilatory Factors in Obese Adolescents[END_REF] . Otherwise, breathing difficulties could also explain that youths with grade II obesity were less likely to met an HR>95%HR max . It should be noted however that the relative HR at which the RPE was 5/10 ("hard") was higher than what could be expected, supporting previous data showing that RPE is underestimated in youths with obesity [START_REF] Marinov | Ventilatory efficiency and rate of perceived exertion in obese and non-obese children performing standardized exercise[END_REF][START_REF] Thivel | Pediatric Obesity and Perceived Exertion: Difference Between Weight-Bearing and Non-Weight-Bearing Exercises Performed at Different Intensities[END_REF][START_REF] Thivel | Reproducibility of the intermittent Spartacus Run Test in obese adolescents[END_REF] .

Other factors may explain that youths with higher grades of obesity perceive their effort as more difficult. Firstly, obesity has many deleterious consequences on the musculoskeletal function of the lower limbs in youths [START_REF] O'malley | A pilot study to profile the lower limb musculoskeletal health in children with obesity[END_REF] . Thus, lower limb pain sensations may contribute to the higher RPE in those who have obesity. Considering the low RPE scores found at the end of an incremental running or walking fitness test [START_REF] Marinov | Ventilatory efficiency and rate of perceived exertion in obese and non-obese children performing standardized exercise[END_REF][START_REF] Mendelson | Low Cardiorespiratory Fitness is Partially Linked to Ventilatory Factors in Obese Adolescents[END_REF][START_REF] Thivel | Pediatric Obesity and Perceived Exertion: Difference Between Weight-Bearing and Non-Weight-Bearing Exercises Performed at Different Intensities[END_REF] , it could be relevant to use an RPE scale based on muscular sensations. Indeed, current RPE scales are based on cardiorespiratory sensations [START_REF] Quinart | Validation of a new scale to measure perceived exertion in obese adolescents: the Childhood Obesity Perceived Exertion (Cope-10) Scale[END_REF] , assuming that cardiorespiratory function is the limiting factor during an aerobic fitness test. Conversely, the combination of both a "local" and "overall" RPE scale would make it possible to dissociate the chest sensations from the lower limb sensations, which could help clinicians to better identify the reason for stopping the test. This approach has already been used in healthy children using an incremental cycling test and showed that "leg RPE" was significantly higher than "overall RPE" and "chest RPE". [START_REF] Robertson | Children's OMNI scale of perceived exertion: mixed gender and race validation[END_REF] Secondly, the lack of experience in vigorous-intensity exercises (due to poor engagement in physical activities) is known to negatively affect RPE measurements [START_REF] Belanger | Do obese children perceive submaximal and maximal exertion differently?[END_REF] . However, physical activity levels were not measured in this study, limiting further interpretation of current data.

The probability of achieving VO 2max using a treadmill protocol according to the grade of obesity has previously been studied in women [START_REF] Wood | Vo 2max in overweight and obese adults: do they meet the threshold criteria?[END_REF] . In accordance with the present results in youths, the authors found no differences between grade I and grade II obesity women regarding the achievement of a VO 2 plateau. However, they also reported no differences regarding the achievement of secondary threshold criteria (i.e., HR>95%HR max , RER>1.15, RPE>18/20, and blood lactate >8 mmol.L -1 ), contrary to the present study.

Therefore, the lower ability of grade II obesity individuals to meet an HR>95%HR max during a treadmill test could only concern youths. According to a recent study [START_REF] Berntsen | Do Obese Children Achieve Maximal Heart Rate during Treadmill Running? Sports[END_REF] , treadmill tests would elicit ~4% lower HR peak than actives games in children with obesity. This study reinforces the idea that a high proportion of youths with obesity is not achieving a maximal effort during a treadmill test.

The VO 2 plateau used to be the primary criteria for verifying VO 2max , but it is still debated whether the occurrence of VO 2 plateau is only influenced by the sampling intervals, or if physiological factors could be involved [START_REF] Astorino | Alterations in VOmax and the VO plateau with manipulation of sampling interval[END_REF][START_REF] Beltrami | High prevalence of false-positive plateau phenomena during VO2max testing in adolescents[END_REF][START_REF] Niemeyer | Oxygen uptake plateau: calculation artifact or physiological reality?[END_REF] . In adults with obesity [START_REF] Wood | Vo 2max in overweight and obese adults: do they meet the threshold criteria?[END_REF] , the VO 2 plateau was reached by 46% of the participants, which is higher than the 37% found in our study using the same definition of a VO 2 plateau, but close to previous studies in youths with [START_REF] Breithaupt | Using the oxygen uptake efficiency slope as an indicator of cardiorespiratory fitness in the obese pediatric population[END_REF] or without [START_REF] Sansum | Validity of the Supramaximal Test to Verify Maximal Oxygen Uptake in Children and Adolescents[END_REF] obesity. Since the range of VO 2 plateau can vary from 30 to 50% between the studies [START_REF] Breithaupt | Using the oxygen uptake efficiency slope as an indicator of cardiorespiratory fitness in the obese pediatric population[END_REF][START_REF] Sansum | Validity of the Supramaximal Test to Verify Maximal Oxygen Uptake in Children and Adolescents[END_REF][START_REF] Rowland | Oxygen uptake plateau during maximal treadmill exercise in children[END_REF] , the present differences are probably not meaningful. Furthermore, comparisons between studies might be not available since a different sampling interval [START_REF] Astorino | Alterations in VOmax and the VO plateau with manipulation of sampling interval[END_REF] can significantly affect the results regarding the VO 2 plateau.

Limitations and perspectives

This study has limitations. Although the methodology was based on the most recent recommandations [START_REF] Niemeyer | Oxygen uptake plateau: calculation artifact or physiological reality?[END_REF][START_REF] Marinus | Altered gas-exchange at peak exercise in obese adolescents: implications for verification of effort during cardiopulmonary exercise testing[END_REF] , the criteria that should be used to verify VO 2max are still debated.

However, the use of supramaximal tests to verify the achievement of VO 2max is not clinically appropriate in youths with obesity. Furthermore, the present results could have been strengthened by physical activity data and the use of an RPE scale based on the muscular feeling of pain. These data could have strengthened the interpretation of the results regarding the mechanisms associated with the achievement of a VO 2 plateau or secondary threshold criteria for verifying VO 2max . Despite these limitations, the present findings are of clinical importance since performing a maximal cardiorespiratory effort is the gold standard for measuring VO 2max [START_REF] Poole | Determinants of oxygen uptake. Implications for exercise testing[END_REF] . According to the criteria used in this study and those of predecessors [START_REF] Gutin | Effects of exercise intensity on cardiovascular fitness, total body composition, and visceral adiposity of obese adolescents[END_REF][START_REF] Marinov | Ventilatory efficiency and rate of perceived exertion in obese and non-obese children performing standardized exercise[END_REF][START_REF] Breithaupt | Using the oxygen uptake efficiency slope as an indicator of cardiorespiratory fitness in the obese pediatric population[END_REF] , a high proportion of youths with different grades of obesity (approximating 50% in our study) would not VO 2max using a treadmill walking test. Thus, clinicians should be aware that VO 2max may be underestimated in this population when it is measured using these protocols.

Conclusion

The rate of achievement of an HR>95% HR max and RPE>5/10 ("hard") was respectively independent of the grade of obesity in a sample of 72 youths. These criteria were achieved by a minority of participants (15 to 64%), suggesting that a maximal cardiorespiratory effort was not achieved by the participants during the test. Youths with grade II obesity displayed a lower peak in minute ventilation normalized to fat-free mass and a premature feeling of difficulty than grade I during progressive walk tests, which may be attributable to higher ventilatory constraints and premature feeling of leg

table 1

 1 ' characteristics, cardiorespiratory, and perceptual responses are shown in

	. Youths with grade II obesity were less mature (P<0.05, d=0.5), heavier
	(d=1.25) and had a higher BMI z-score (d=4.47) and percentage of body fat (d=2.38)
	than grade I (P<0.001 for each parameter). VE peak divided by body mass (P<0.001,
	d=1.34), fat-free mass (P<0.05, d=0.72) and VO 2peak divided by body mass (P<0.01,
	d=1.03) were lower in youths with grade II obesity than grade I. Conversely, no

*** Table

2

here ***

Table 1

 1 Participants' characteristics, cardiorespiratory and perceptual responses by grade of obesity subgroups.

	IOTF Grade I obesity	IOTF Grade II obesity
	(n=15)	(n=57)

Table 2 :

 2 Number and rate of participants having reached the criteria for verying VO 2max and achieved VO 2max according to these criteria, by grade of obesity and results of

	logistic regressions.
		VO 2max	VO 2 plateau	HR>95%	RER>1.0	RPE>5/10
		Yes	No	Yes	No	Yes	No	Yes	No	Yes	No
	Grade	n=7	n=8	n=4	n=11	n=9	n=6	n=8	n=7	n=12	n=3
	I	47%	53%	27%	73%	60%	40%	53%	47%	80%	20%
	Grade	n=33	n=24	n=23	n=34	n=19	n=38	n=39	n=18	n=54	n=3
	II	58%	42%	40%	60%	33%	67%	68%	32%	95%	5%
	X² tests	P=0.44	P=0.33	P=0.06	P=0.27	P=0.07
	OR	NS	NS	0.33	NS	4.5

Abbreviations : VO 2 = oxygen uptake ; HR = heart rate ; OR = odds-ratio; RER = respiratory exchange ratio ; RPE = rating of perceived exertion ; X² = chi-square tests

muscular pain. Thus, our study supports the idea that a dose-response relationship exists between the grade of obesity and the psychophysiological responses to an incremental treadmill test in children.


Bo Wang

Julien Réthoré

Katerina E Aifantis
email: kaifantis@ufl.edu

Capturing the stress evolution in electrode materials that

Keywords: Electrodes, Phase transformation, Chemo-mechanical coupling, Concentration-dependent elastic modulus

It is found that the DISs arise from the inhomogeneous volume expansions resulting from Li concentration gradients and the hydrostatic stress facilitates the diffusion of Li-ions under elastic deformation while it hinders diffusion in the plastic case. When the elastic modulus is reduced the magnitude of the diffusion-induced stress decreases but the strain increases under elastic deformation whereas the opposite occurs for the plastic case. Furthermore, if the electrode is assumed to undergo strain softening during plastic deformation, smaller stresses and higher plastic strains are predicted than when strain hardening is assumed. The models highlight the importance of chemo-mechanical coupling effects, concentration-dependent material properties and plastic deformation on diffusion-induced stresses. These findings render prospective insights

Introduction

Lithium ion batteries (LIBs) are currently widely used in portable electronics and will continue to play a significant role in next generation electric vehicles, utility grids and electric/hybrid-electric airplanes. To fulfill the practical requirements, LIBs should have an optimal electrochemical and mechanical performance, i.e, high energy density, high power density, and long cycle life. The operation of LIBs involves the continuous diffusion of Li-ions between the anode and cathode and therefore the performance of LIBs depends largely on the electrode active materials (materials that react with Li + ).

The cathode is the source of lithium ions (Li + ), and the anode should have the ability to host the Li + during charge and vice versa for the discharge process. Lithium insertion and de-insertion are accompanied by changes in lattice spacing. As a result, the lattice structure can be altered, and phase transformations can occur, resulting in volume changes. Consequently, diffusion-induced stresses (DISs) are generated. The energy density, power density and cycle life of LIBs are closely related to DISs. Excessive stresses result in pulverization of the electrode materials and their detachment from the current collector. As a result, fracture of the electrode and capacity decay occurs.

For instance, commercially used cathode materials (LiMn 2 O 4 and LiCoO 2 ) for LIBs have achieved only half of the theoretical capacity, due to the mechanical stresses that arise from phase transformations and changes in the crystal structures/lattice parameters during the charging/discharging processLiu et al. (2004); [START_REF] Venkatraman | Capacity of layered cathode materials for lithium-ion batteries-a theoretical study and experimental evaluation[END_REF]; [START_REF] Amatucci | Cobalt dissolution in licoo2based non-aqueous rechargeable batteries[END_REF]; [START_REF] Aifantis | An ideal elasticity problem for li-batteries[END_REF]. Similarly, high-capacity anode materials (Si and Sn) cannot retain their initial capacity ( 900-4000 mAh/g) due to the fracture that results from the lithiation-induced volume expansions (300%-400%) [START_REF] He | Shape evolution of patterned amorphous and polycrystalline silicon microarray thin film electrodes caused by lithium insertion and extraction[END_REF]; [START_REF] Wang | Two-phase electrochemical lithiation in amorphous silicon[END_REF]. Therefore, understanding the stress evolution process is crucial for revealing the degradation mechanisms of LIBs and thus, improving their performance. Initial models considered the stress evolution [START_REF] Aifantis | An ideal elasticity problem for li-batteries[END_REF] and fracture [START_REF] Aifantis | Cracking in si-based anodes for li-ion batteries[END_REF][START_REF] Aifantis | Fracture of nanostructured lithium batteries[END_REF] in active electrode particles by using a purely mechanical approach. Since then, multiple models have been developed to determine the DIS in material systems. The phenomenon of DIS was first studied by employing an analogy between thermal stress and DIS to analyze the transverse stresses introduced by solute lattice contraction of boron and phosphorus in a thin silicon plate during mass transport [START_REF] Prussin | Generation and distribution of dislocations by solute diffusion[END_REF]. A number of analytical solutions to DIS problems in elastic media of simple spherical and cylindrical geometries were provided in [START_REF] Li | Physical chemistry of some microstructural phenomena[END_REF], while the DISs arising from inhomogeneous concentration in materials was in investigated in [START_REF] Larcht'e | The effect of self-stress on diffusion in solids[END_REF]; [START_REF] Larche | The interactions of composition and stress in crystalline solids[END_REF]. The effect of DISs on diffusion in hollow cylinders [START_REF] Lee | Diffusion-induced stresses in a hollow cylinder[END_REF]; [START_REF] Wang | Effect of chemical stress on diffusion in a hollow cylinder[END_REF] and square sandwich composites [START_REF] Ko | Chemical stresses in a square sandwich composite[END_REF] has also been studied. Later on, numerous researchers followed this analogy approach to investigate the DISs in nanowire electrodes Deshpande et al. (2010a), spherical electrode particles Verbrugge & Cheng (2009); Cheng & Verbrugge (2010[START_REF] Cheng | The influence of surface mechanics on diffusion induced stresses within spherical nanoparticles[END_REF], cylindrical electrodes [START_REF] Deshpande | Effects of concentrationdependent elastic modulus on diffusion-induced stresses for battery applications[END_REF]; [START_REF] Song | Diffusion induced stresses in cylindrical lithium-ion batteries: analytical solutions and design insights[END_REF] and layered electrode plates [START_REF] Zhang | Diffusion induced stress in layered li-ion battery electrode plates[END_REF].

For intercalation-induced stresses, one-dimensional (1D) models have been developed by to investigate the stress generation due to Li + insertion into a spherical carbon anode [START_REF] Christensen | Stress generation and fracture in lithium insertion materials[END_REF], but instead of thermal analogy, the governing equations were derived based on lattice deformation analysis. Given that analytical solutions to diffusion equations under specific conditions are available [START_REF] Crank | The mathematics of diffusion[END_REF]; [START_REF] Jaeger | Analytical modeling and simulation of electrochemical charge/discharge behavior of si thin film negative electrodes in li-ion cells[END_REF], analytical expressions have been obtained to study the stress evolution in spherical electrode particles Cheng & Verbrugge (2009, 2010).

Most of these studies are centered on single-phase reactions, assuming that the electrode particle phase does not change during Li + insertion/de-insertion.

However, it's well known that lithium-ion insertion often leads to the forma-tion of different phases. Phase transformations have been experimentally observed during Li + insertion/de-insertion and could play an important role in the stress evolution during the discharging/charging process Shin & Pyun (1999a,b); [START_REF] Zhang | Moving boundary model for the discharge of a licoo2 electrode[END_REF]. Multiple models have been constructed to capture the effect of phase transformations on DISs. The "shrinking core" model [START_REF] Srinivasan | Discharge model for the lithium ironphosphate electrode[END_REF] considers that the core and shell would consist by lithium-rich and lithium-poor phases, respectively, during Li + insertion and vice-versa for de-insertion. This model leads to an abrupt change (jump) in the concentration at the phase boundary which also results in a stress discontinuity at the interface between the two phases. By using this model and extending the theoretical framework for a single-phase reaction [START_REF] Christensen | Stress generation and fracture in lithium insertion materials[END_REF], the stress generation in a spherical LiMn 2 O 4 particle was estimated Christensen & Newman (2006a). This "shrinking core" model has also been applied to study other spherical electrode particles and other configurations [START_REF] Deshpande | Diffusion induced stresses and strain energy in a phase-transforming spherical electrode particle[END_REF]; [START_REF] Liu | Stress fields in hollow core-shell spherical electrodes of lithium ion batteries[END_REF]; [START_REF] Esmizadeh | An investigation of the critical conditions leading to deintercalation induced fracture in twophase elastic electrode particles using a moving interphase core-shell model[END_REF]. Due to the sharp interface, concentration and stress discontinuities appear at the interface. It is noted that the moving phase boundary needs to be tracked and specific boundary conditions need to be applied during insertion and de-insertion. Therefore, this sharp interface model poses numerical difficulties for practical applications. A possible alternative approach would be phase field modeling of the diffusion. By employing a diffusive interface model, the ion distributions can be modeled with a continuous concentration field. As a result, the reaction front does not need to be explicitly tracked. Such phase field models have been successfully applied to study the stress generation of LiFePO 4 [START_REF] Zhang | Sodium ion batteries particles: Phase-field modeling with coupling of cahn-hilliard equation and finite deformation elasticity[END_REF] andLiMn 2 O 4 Huttin & Kamlah (2012) particles.

Phase field models, aside from the work in Hu et al. (2019a) focus on the DISs of cathode materials. Due to the generally small volume expansions that cathodes experience during lithiation, infinitesimal/small strain continuum mechanics formulations can be used. For example, LiFePO 4 shrinks approximately 7% after full delithiation [START_REF] Yamada | Phase change in li x fepo4[END_REF]; [START_REF] Wang | Cracking causing cyclic instability of lifepo4 cathode material[END_REF], and therefore cathodes are assumed to undergo linear elastic deformation during the interca-lation/deintercalation process [START_REF] Woodford | electrochemical shock" of intercalation electrodes: a fracture mechanics analysis[END_REF]; Hao & Fang (2013a,b).

However, all high capacity anodes, such as Si, are well known for large volume expansions during lithiation, which can reach 420% upon maximum lithiation (Si 4.4 Li). The deformation is so large that plastic deformation can take place.

Although linear elastic studies of the DISs Jagannathan & Chandran (2014); [START_REF] Tsagrakis | Gradient elasticity effects on the twophase lithiation of lib anodes[END_REF]; [START_REF] Golmon | Stress generation in silicon particles during lithium insertion[END_REF] 2014). An important finding in these aforementioned studies is that the tangential stress transits from compressive to tensile during lithiation when plasticity is taken into account, which is believed to be the cause of the experimentally observed fracture of Sn nanoparticles [START_REF] Aifantis | Capacity fade in sn-c nanopowder anodes due to fracture[END_REF], Si spherical particles [START_REF] Liu | Size-dependent fracture of silicon nanoparticles during lithiation[END_REF] and Si nanowires [START_REF] Ryu | Size-dependent fracture of si nanowire battery anodes[END_REF] during lithiation. Elastic models predict a compressive hoop stress on the outer surface at all times while plastic models can capture this transition (compression in the core to tension in the shell) since the hoop tension is the crack driving force of the observed surface cracking [START_REF] Huang | Stress generation during lithiation of high-capacity electrode particles in lithium ion batteries[END_REF]; [START_REF] Liu | Size-dependent fracture of silicon nanoparticles during lithiation[END_REF]; [START_REF] Chen | A phase-field model coupled with large elasto-plastic deformation: application to lithiated silicon electrodes[END_REF]; [START_REF] Chang | Phase field model for two-phase lithiation in an arbitrarily shaped elastoplastic electrode particle under galvanostatic and potentiostatic operations[END_REF].

Despite the fact that there have been a lot of models on DISs of electrode materials, a comprehensive understanding has yet to be achieved. One aspect is the chemo-mechanical (CM) coupling effect. Initial works, particularly those using the "shrinking core" model, did not consider CM coupling. Subsequent models considered the coupling effect between chemical diffusion and mechanical deformation [START_REF] Christensen | Stress generation and fracture in lithium insertion materials[END_REF]; [START_REF] Liu | Stress fields in hollow core-shell spherical electrodes of lithium ion batteries[END_REF]; [START_REF] Zhang | Sodium ion batteries particles: Phase-field modeling with coupling of cahn-hilliard equation and finite deformation elasticity[END_REF]; Hao & Fang (2013a,b); [START_REF] Tsagrakis | Gradient elasticity effects on the twophase lithiation of lib anodes[END_REF]; Zhao et al. (2011a); [START_REF] Chen | A phase-field model coupled with large elasto-plastic deformation: application to lithiated silicon electrodes[END_REF]. However, few works have taken into considera-tion the change that occurs in the elastic modulus of electrode materials during Li + insertion. A decrease or increase in the elastic modulus during Li-insertion (lithation), as a function of Li-concentration, have been observed experimentally [START_REF] Sethuraman | In situ measurement of biaxial modulus of si anode for li-ion batteries[END_REF]; [START_REF] Hertzberg | Ex-situ depth-sensing indentation measurements of electrochemically produced si-li alloy films[END_REF] and calculated through first principle calculations [START_REF] Maxisch | Elastic properties of olivine li x fepo 4 from first principles[END_REF]; [START_REF] Qi | Lithium concentration dependent elastic properties of battery electrode materials from first principles calculations[END_REF][START_REF] Qi | Threefold increase in the young's modulus of graphite negative electrode during lithium intercalation[END_REF]. Table 1 shows the Young's modulus of the lithiated and delithiated phases of some electrode materials. The effect of a concentration dependent elastic modulus on DIS was systematically examined using linear and gradient elasticity [START_REF] Natarajan | Size effects in nanostructured li-ion battery cathode particles[END_REF], showing that there was a negligible difference in DISs when the variation of the elastic modulus due to Li insertion was only 10% and a singlephase reaction during Li-insertion occurred. In such cases, constant moduli may be used for modeling DISs in a single-phase battery system. However, very few studies examined the effects of the concentration-dependent material properties on DISs in phase-transformation electrode materials, and there isn't much work comparing plastic and elastic models. Herein we present a new formulation that can account for chemomechanical coupling, concentration dependent elastic moduli for two phase transitions and large volume expansions, for both electrodes that undergo either pure elasticity or also plasticity during lithiation.

By performing such a comprehensive study of the DIS, more accurate predictions of the stress fields in the electrodes can be provided, allowing for a better selection of the materials and microstructure that will result in improved electrochemical performance. It should be noted that this work can be applied to other electrochemical systems that operate under the diffusion of ions, such as Na-ion batteries.

Phase field modeling of DISs

As mentioned earlier, a core-shell model is often used to model the diffusion process of electrode materials that exhibit phase segregation into Li-rich and Li-poor phases, respectively, during Li-ion insertion and de-insertion. It's often assumed that during insertion the phase boundary moves in the radial direction For the cathode with small volume variations, the deformation is assumed to remain in the elastic stage. The system (total) free energy of the domain Ω under consideration is

Ψ = Ω ψdΩ (1)
where Ψ is the free energy density, i.e., free energy per unit reference volume and is comprised of two parts:

ψ = ψ (ε, c, ∇c) = ψ e + ψ c (2)
where ψ e is the bulk elastic energy density solution, c is the actual species concentration, given in terms of moles per unit reference volume (the molar concentration per unit volume in the reference configuration), and ψ c is the energy contribution due to changes in the concentration solution. The strain and displacement are related as:

ε = 1 2 ∇u + (∇u) T ( 3
)
where u is the displacement field. The elastic strain energy density is

ψ e (ε e ) = 1 2 ε e : C : ε e (4)
where the elastic strain tensor is ε e = ε-ε c . ε c is the chemically induced strain tensor and is given by :

ε c = α (c * -c * 0 ) I ( 5
)
where I is a 2 nd -order identity tensor, C is a 4 th -order elasticity tensor, given by C = λδ ij δ kl + Λ (δ ik δ jl + δ il δ jk ), δ ij describes the Kronecker delta, λ and Λ are the Lamé constants, and c * is the relative concentration. The normalized concentration c * is used to characterize the relative saturation level of Li in an electrode. It's normalized by the maximum concentration c max in the Li/electrode system, ranging from 0 (no ions) to 1 (full ion insertion), given by:c * = c/c max . c * 0 is a reference concentration, usually set as 0. α is a dilatation coefficient (analogous to the thermal expansion coefficient), given by:

α = 3 1 + Ξc max -1 . = 1 3 Ξc max ( 6
)
where Ξ is the partial molar volume of Li-ions. Due to the presence of a sharp interface, the energy contribution due to a concentration gradient on diffuse interfaces cannot be neglected. It should be noted that for phase transformation materials, there is a sharp interface, however, phase field modeling uses diffusive interfaces with a finite thickness (but very thin) to approximate it, resulting in a continuous concentration. In traditional sharp interface models, the thickness is assumed to be zero, and as a result, the concentration is discontinuous across the interface. The chemical energy density can be split into two parts:

ψ c = ψ c (c, ∇c) = ψ ch + ψ cg (7)
where ψ ch is the homogeneous chemical free energy density (free energy density of a homogeneous system of uniform concentration) and ψ cg is the chemical gradient energy density due to the concentration gradient ∇c. The chemical stored free energy (chemical contribution to the stored energy) is assumed to have the form:

ψ ch = ψ ch (c) = RT c max [c * ln c * + (1 -c * ) ln (1 -c * ) + χc * (1 -c * )]+µ 0 c max c * . ( 8
)
Here R is the gas constant, T is the absolute temperature, and RT denotes the specific molar energy. µ 0 is a reference value of the chemical potential of the diffusing species, set as 0 since it will not affect the diffusion behavior. χ is the constant partial molar volume (the volume of a mole of the species , parameter χ indicates the convexity of the energy) which is dimensionless and controls the shape of the double-well energy function characterizing the Li-rich and Li-poor phases. This form of the homogenous energy density favors a separation of the system towards c * = 0 and c * = 1, i.e., minimization of the Gibbs free energy.

To ensure phase separation, χ > 2 . The chemical gradient energy density is given by:

ψ cg = κ 2 |∇c| 2 = κ 2 |∇c * | 2 (9)
depending on the concentration gradient and representing an interfacial energy. κ = κc max 2 and κ is the gradient energy coefficient, with units of energy per unit volume times a length squared. The parameter κ modulates the contribution of the large concentration gradient at the two-phase interface to the free energy.

When the elastic property of the material is strongly influenced by concentration, based on rule of mixtures from composites, it can be approximated by a simple linear relation:

E (c * ) = E 0 + kc * ( 10
)
where E is the concentration-dependent Young's modulus, and E 0 is the Young's modulus of the delithiated phase (i.e. initial elastic modulus), and k is a positive or negative parameter allowing for a reduction or increase in the modulus during lithiation. The chemical potential µ is defined as the variational derivative of the system free energy:

µ = δψ δc = ∂ψ ∂c -∇ • ∂ψ ∂∇c = µ 0 + RT ln c * 1 -c * + RT χ (1 -2c * ) -Ξσ h -c max ∇ • (κ∇c * ) - 1 2c max dS ijkl dc * σ ij σ kl (11)
where σ h = 1 3 tr(σ) is the hydrostatic stress and S ijkl is the concentrationdependent compliance tensor. Therefore, the diffusional potential is related to the Li concentration and the local stress state, as well as the concentration gradient and the changing material properties due to intercalation. For linear elastic and isotropic materials, the compliance tensor is given by:

S ijkl = 1 -2υ 3E δ ij δ kl + 1 + υ 2E δ ik δ jl + δ il δ jk - 2 3 δ ij δ kl ( 12
)
where υ is Poisson's ratio. As a result,

dS ijkl dc * σ ij σ kl = 1 E 2 dE dc * (υσ ii σ kk -(1 + υ) σ jl σ jl ) (13)
The Li transport in the host materials can be described by the diffusion equation,

as ∂c * ∂t + ∇• j = 0 in Ω ( 14
)
where j is the flux, and Ω is the bulk of the electrode. The gradient of the chemical potential is the chemical driving force for diffusion, i.e. the gradient of the diffusional potential of Li µ serves as the driving force for Li diffusion:

j = -M ∇µ ( 15
)
where M is the mobility tensor, which is generally a function of Li concentration c. For isotropic materials, it's assumed to have the form:

M = c * (1 -c * ) M I ( 16
)
where M is the solute mobility and is independent on time. The diffusion coefficient D is given by:

D = M RT (17)
To solve for the diffusion, suppose the electrode is under galvanostatic operation and the particle is charged galvanostatically, with a uniform lithium flux through its free surface according to

M c * (1 -c * ) ∂µ ∂r | r=ρ = i ( 18
)
Here i is the applied current density (per unit area), and ρ is the radius of the electrode particle. The variational boundary condition is given by:

κ∇c • n c = 0 on ∂Ω c (19)
with n c being the outward unit normal to the concentration boundary ∂Ω c .

The weak form of the diffusion equation can be given by: 20)

Ω ∂c * ∂t vdV = ∂Ω ivdS - Ω M c * (1 -c * )∇µ∇vdV (
Ω µwdV = Ω RT ln c * 1 -c * + 2χ (1 -c * ) wdV - Ω κc max ∇c * ∇wdV - Ω Ξσ h wdV - Ω 1 2c max dS ijkl dc * σ ij σ kl wdV (21)
where v and w are arbitrary test functions. Considering a spherical particle during the lithiation process, the mechanical equilibrium equation under spherically symmetrical deformation is given by:

dσ r dr + 2 (σ r -σ θ ) r = 0 (22)
where σ r is the radial stress, σ θ is the hoop stress and r is the radial distance.

The constitutive law (stress-strain relation) for an isotropic spherical electrode during lithiation is given by:

σ r = E (1 + υ) (1 -2υ) (1 -υ) du dr -αc * + 2υ u r -αc * σ θ = E (1 + υ) (1 -2υ) u r -αc * + υ du dr -αc *        (23)
Here ε r = du dr is the radial strain and ε θ = u r is the tangential strain and u is the radial displacement. The boundary conditions for this problem are: zero radial displacement at the center and zero radial stress on the outer surface, i.e., u | r=0 = 0 and σ r | r=ρ = 0 . As a result, when E is not dependent upon c, then the stress field can be given analytically by [START_REF] Hu | Capturing the differences between lithiation and sodiation of nanostructured tis2 electrodes[END_REF]:

σ r = 2αE 1 -υ 1 ρ 3 ρ 0 c * r 2 dr - 1 r 3 r 0 c * r 2 dr σ θ = 2αE 1 -υ 1 ρ 3 ρ 0 c * r 2 dr + 1 2r 3 r 0 c * r 2 dr - c * 2        (24)
For a general case, the Young's modulus may not be a constant. In such cases, the analytical solution cannot be obtained, and finite element or other numerical methods need to be employed. The weak form of the mechanical equilibrium equations can be given by:

ρ 0 q dσ r dr + 2 (σ r -σ θ ) r r 2 dr = 0 ( 25
)
where q is a test function. Integrating by parts and applying boundary conditions, one can finally obtain:

ρ 0 σ r dq dr + 2σ θ q r r 2 dr = 0 (26)
By substituting Equation 23 into Equation 26, we have the weak form in terms of the displacement as follows:

ρ 0 q (λ + 2Λ) r 2 d 2 u dr 2 + 2 (λ + Λ) r du dr -u -(3λ + 2Λ) αr 2 dc * dr dr = 0 (27)

Numerical implementation of the phase field model

In order to seek the solutions, the finite element method is used. Due to spherical symmetry, the problem can be reduced to a one-dimensional problem.

The weak form cannot be solved directly. Instead, we have to discretize in both the temporal and spacial space. Therefore, Equation 20 and Equation 21 can be rewritten as:

Ω c n+1 * -c n Ω µ n+1 wdV = Ω RT ln c n * 1 -c n * + 2χ (1 -c n * ) wdV - Ω κc max ∇c n * ∇wdV - Ω Ξ (σ h ) n wdV - Ω 1 2c max dS ijkl dc * (σ ij ) n (σ kl ) n wdV (29)
where the subscripts n + 1 and n are the corresponding values at the time step t n+1 and t n respectively, with ∆t n = t n+1 -t n being the n th time step. The hydrostatic stress at t n is given by:

(σ h ) n = 1 3 [(σ r ) n + 2 (σ θ ) n ] (30)
And Equation 13at t n can be given by:

dS ijkl dc * (σ ij ) n (σ kl ) n = 1 E 2 dE dc * [(σ r ) n + 2 (σ θ ) n ] 2 -(σ r ) 2 n + 2 (σ θ ) 2 n (31)
By using this kind of mixed formulation, the problem becomes to find c n+1 * and µ n+1 given that c n * and µ n are known. The built-in DOLFIN Newton solver in FEniCS Dol is used to solve the nonlinear equations. As for the mechanical equilibrium equation, the built-in NonlinearVariationalSolver in FEniCS is employed. The basic idea is to calculate the Gateaux derivative of the nonlinear form F (u; q) = 0 (corresponding to Equation 26). Thus, the Jacobian and the trial function will be automatically computed and solved. By solving the diffusion equation (Equation 28and Equation 29), the concentration and the chemical potential can be obtained. Then the radial displacement can be obtained by solving the weak form of the mechanical equilibrium equation (Equation 27). The unknown field variables (c * , µ, u) can be solved by using this staggered scheme at every time step. In the present simulations, the following dimensionless/normalized variables are used as shown in Table 3, where t 0 = ρ 2 D is the characteristic time, F is Faraday's constant, and z is the valence number of a Li-ion (z =1). As an example, the stress evolution of a spherical LiFePO 4 electrode is studied. The material properties of LiFePO 4 electrode are given in Table 2 Zhang & Kamlah (2018). The initial radius of the electrode is ρ = 100nm. Normalized radius is 1 and total number of the mesh elements of the unit interval is 100. (It follows that the model results can be applied to any electrode material since the variables are normalized. The general trend is of interest here and not the specific values.) The parameters are as follows: E 0 =124.5 GPa , υ=0.24,

Numerical results

Some

M * = 1, i * = 1, κ * = 0.0004, k = 0, Ξ=3.49×10 -6 m 3 /mol, χ = 2.6, c max =
2.29 × 10 4 mol/m 3 , T =300 K, ∆t=0.0005, α= 0.0221. These parameters are used to study the evolution of the concentration, radial stress, hoop stress, hydrostatic stress, von Mises stress, as well as the evolution of the hoop stress on the outer surface, as shown in Figure 1. Only the DISs during the lithium insertion process are considered and a constant Young's modulus is assumed for the two phases during lithiation. due to the fact that a higher stress is induced for materials with an increasing modulus as a function of Li-concentration when compared to the cases where the modulus is either not dependent on concentration or it decreases as a function of it. Also, the stresses predicted from the coupling case are smaller than those from the non-coupling case when the same material parameters are used. This is because the coupling model predicts a smaller miscibility gap (concentration differences between the Li-rich phase and Li-poor phase).

To gain a better understanding of the physics behind the stress evolution, let's consider the lithiation process. During Li-insertion, the outer shell contains more Li-ions than the inner core and, as a result, the core expands to a lesser degree. The expansion is constrained by the exterior region, leading to show similar trends as those of the von Mises stress distribution. An interesting finding is that the magnitudes of the strains are higher when the modulus de-creases during ion insertion, as opposed to when it increases with concentration.

Temporal distributions of the stresses in the center ( Figure 3 to the traction-free boundary conditions on the outer surface, the radial stress on the outer surface vanishes (i.e., is always zero). As a result, cracks are not likely to occur and propagate during lithiation.

DISs for elasto-plastic materials

Phase field model for elasto-plastic materials

In alloy-based electrode systems, such as Li/Si, a large volumetric change may be induced (over 400% for Li/Si), which can give rise to plastic deformation. If infinitesimal changes occur in the stress and chemical composition, a unique infinitesimal change of strain is produced, and the change is independent of the speed applied. In this case, we can say that the chemo-plasticity model is rate-independent. The model can be formulated from an energy perspective. Generally speaking, the formulation can be either in terms of stress space or strain space. The constitutive relation (flow rule) for the plastic strain rate is usually given in terms of the stress space. As a result, the stress-based formulation is employed here.

The Gibbs free energy ψ = ψ (σ, c * , ξ) is a function of the stress tensor (σ), chemical composition (c * ) and internal (intrinsic) variable (ξ). In this notation, the chemoelastic part of the strain tensor is given by:

ε ec = ∂ψ ∂σ (32)
The yield condition in stress space can be expressed as:

f (σ, c * , ξ) = 0 (33)
And the chemoplastic strain rate is given by:

εp = γ ∂g ∂σ ( 34
)
where g = g (σ, c * , ξ) is the plastic flow potential and γ is the rate of the plastic 300 multiplier. Since plastic deformation is an irreversible process, the multiplier must be non-negative. For the associative flow rule, g = f .

Kuhn-Tucker conditions must be satisfied during the deformation process, which are stated as:

f ≤ 0, γ ≥ 0, γf = 0 (35)
This consistency condition is also equivalent to:

ḟ ≤ 0, γ ≥ 0, γ ḟ = 0 (36)
The hardening behavior of the material is determined by h and it's a function of c * , and the history variable (ϑ). The hardening rule can be expressed as:

ξ = γh (c * , ϑ) (37)
and the hardening parameter H o is defined as:

H o = - ∂f ∂ξ h (c * , ϑ) (38)
The chemoelastic strain rate is given by:

εec = ∂ 2 ψ ∂σ ⊗ ∂σ : σ + ∂ 2 ψ ∂σ∂c * ċ * (39)
The general form of ψ is as follows:

ψ = ψ e (σ) + ψ ch (σ, c * ) + ψ p (c * , ϑ) (40)
The three terms represent the energy contribution from mechanical elastic part, chemical elastic part, and plastic part respectively. The mechanical elastic part and plastic part are given as:

ψ e (σ) = 1 4Λ tr σ 2 - λ 3λ + 2Λ tr 2 (σ) (41) ψ p (c * , ϑ) = 1 2 Hζ 2 + σ 0 y ζ (42)
Linear isotropic hardening is assumed with H being the hardening modulus, and

ζ is a strain-like internal variable (normally it's equivalent to plastic strain). As a result, the chemoelastic part of the strain rate could be given as:

εec = 1 2Λ I - λ 3λ + 2Λ I ⊗ I : σ + ∂ψ ch ∂c * ċ * (43)
For a homogenous and isotropic material, a linear expansion is assumed,

ψ ch (σ, c * ) = η (c * ) tr (σ) (44)
η is given by:

η = α∆c * = α (c * -c * 0 ) (45)
where c * 0 is a reference concentration, usually set as 0. The consistency condition is ḟ = 0 , i.e., ḟ = ∂f ∂σ

: σ + ∂f ∂ξ : ξ + ∂f ∂c * ċ * = 0 (46)
The consistency parameter is given by: γ = ċ *

H 0 + 2Λ = ∂f ∂σ : C : ε + ∂f ∂c * ċ * 2 3 H + 2Λ ( 48
)
When plasticity is taken into consideration, the constitutive relations for a spherical particle can be given by:

σ r = E (1 + υ) (1 -2υ) (1 -υ) du dr -ε p r -αc * + 2υ u r -ε p θ -αc * σ θ = E (1 + υ) (1 -2υ) u r -ε p θ -αc * + υ du dr -ε p r -αc *        (49)
And the deviatoric parts are:

σ r = 2 3 (σ r -σ θ ) , σ θ = σ φ = - 1 3 (σ r -σ θ ) (50)
As a result, the effective stress/von Mises stress is:

σ e = 3 2 s : s = |σ r -σ θ | (51)
When linear hardening is assumed, the yield function f and yield stress σ are given by:

f (σ, ξ) = √ s : s - 2 3 σ (ξ) (52) σ (ξ) = σ (ε p ) = σ 0 y + H 0 εp ( 53
)
where εp is the equivalent plastic strain, defined as below:

εp = 2 3 ε p : ε p (54)
The plastic strain components are given by:

εp r = γ σ r s , εp θ = γ σ θ s ( 55
)
For the associative flow rule,

∂g ∂σ = ∂f ∂σ = sign(σ r ) 1 √ 6 diag (2, -1, -1) (56) ∂f ∂c * = ∂f ∂σ : ∂σ ∂c * + ∂f ∂σ 0 y ∂σ 0 y ∂c * + ∂f ∂H 0 ∂H 0 ∂c * = 2 ∂Λ ∂c * ∂f ∂σ : ε -2 ∂Λ ∂c * ∂f ∂σ : ε p - 2 3 ∂σ 0 y ∂c * + ∂H 0 ∂c * εp (57)
When the yield function does not depend on the concentration field the consistency parameter can be determined by:

γ = √ 6Λ H + 3Λ sign(σ r ) ∂ u ∂r - u r Θ( σ e σ -1) (58)
With Θ being the Heaviside unit step function, whose value is 1 for any positive variable and 0 otherwise. And σ e = |σ r -σ θ | is the effective/ von Mises stress.

Numerical implementation of the phase field model

The problem becomes to seek solutions of c n+1 * , µ n+1 , u n+1 , and γ n+1 at the time step t n+1 given that c n * , µ n , u n , and γ n at the time step t n are known.

Since the weak form of the diffusion equation and the mechanical equilibrium equation are the same as in the elastic case, the same approach can be applied to the elasto-plastic case. The only difference comes from the plastic deformation.

The current model does not take kinematic hardening and viscous effects into account, i.e. only the quasi-static problem is solved. When the time step is small enough, the increment of the consistency parameter can be given by:

∆γ n = √ 6Λ H + 3Λ sign((σ r ) n ) ∂ (∆u n ) ∂r - ∆u n r Θ( (σ e ) n σn -1) ( 59
)
where

∆u n = u n+1 -u n , (σ r ) n = 2 3 [(σ r ) n -(σ θ ) n ], (σ θ ) n = -1 3 [(σ r ) n -(σ θ ) n ], (σ e ) n = |(σ r ) n -(σ θ ) n |, σn = σ 0 y + H (ε p ) n , (ε p ) n = 2 3 (ε p ) n : (ε p ) n , s n = 2 3 (σ e ) n , (∆ε r p ) n = ∆γ n (σr) n s n , (∆ε θ p ) n = ∆γ n (σ θ ) n s n .
As a result,

γ n+1 = ∆γ n + γ n , (ε r p ) n+1 = (∆ε r p ) n + (ε r p ) n , (ε θ p ) n+1 = (∆ε θ p ) n + (ε θ p ) n ,
and the other quantities at time step t n+1 can also be computed. The mechanical equilibrium equation is essentially the same as for the elastic case, and the same solver is used. Once the internal iterations for solving the weak form of the mechanical equilibrium equation (Equation 26) converge, one increment of the mechanical field can be obtained by forward Euler integration method. Thus, the mechanical equilibrium equation can be solved, and the stress field can be obtained incrementally.

Numerical results

Here the stress evolution in a spherical Si electrode with radius 100 nm is studied. It has been shown that both the elastic modulus and the yield stress of Si are reduced with continuous Li-insertion [START_REF] Sethuraman | In situ measurement of biaxial modulus of si anode for li-ion batteries[END_REF], however, this concentration-dependence has not been accounted for theoretically.

Furthermore, previous studies used either linear hardening or viscoplasticity to study the stress evolution. However, few looked at strain softening effects, which most likely take place based on the experimental observations Sethuraman et al.

(2010); [START_REF] Soni | Diffusion mediated lithiation stresses in si thin film electrodes[END_REF]. In the present study, both strain hardening and strain softening are considered. For the sake of simplicity, concentration-dependence of the modulus is not considered in this section (but in the next one), and therefore here the elastic modulus is kept the same for all cases. Normalized variables as in the elastic case are employed. The normalized radius is 1 and total number of the element mesh of the unit interval is 100. The other parameters are given by: E 0 =90 GPa , υ=0.24, M * = 1, i * = 1, κ * = 0.0004, χ = 2.6, k = 0 (same Young's modulus for both the core and the shell), σ 0 y = 0. gradually increase from the core towards the shell, and as soon as they exceed the yield stress, plastic deformation begins to take place.If there's no plastic deformation, the stress evolution will be the same as in the elastic deformation case, i.e., the maximum von Mises stress will occur near the core-shell interface.

355 However, once it exceeds the yield stress, there will be plastic deformation. As The von Mises stresses show abrupt changes near core/shell interfaces as they increase first, then decrease and increase again quickly. The von Mises stress towards the shell continues to increase for the strain hardening case while it decreases for the strain softening case.

As for the plastic deformation, the plastic strain components are zero in the core regions. The radial component of the plastic strain for the shell is always positive while the hoop component is always negative. The equivalent plastic strain continues to increase as lithiation continues, and is higher in the strain softening case. It's interesting to notice that there's a relative plateau in the equivalent plastic strain profiles (Figure 5(h), and Figure 6(h)). This is due to the stress drops near the core-shell interface. The von Mises stress near the phase boundary decreases quickly and then increases rapidly again.

However, according to the von Mises yield criteria, plastic deformation can continue to occur only if the von Mises stress exceeds the effective yield stress.

Thus, there won't be plastic strain increment over the stress drop regions. The radial and hoop stresses in the center (Figure 7(a)) are equal due to spherical symmetry. For the initial stages(up to 1000 ∆t), the stress profiles (Figure 7(a))

are almost the same and all in tension. However, they all undergo a transition from tension to compression and the magnitude is the highest when it's strain hardening without CM coupling and the magnitude is the smallest when it's strain softening with CM coupling. The hoop stress is the driving force for crack propagation and on the outer surface (Figure 7(b)) it is compressive during the early lithiation stages and transits from compression to tension at later stages. The hoop stresses on the outer surface are almost the same for the initial stages when strain hardening/softening is considered either with or without CM coupling. However, as soon as the electrode begins to yield, the hoop stress on the outer surface keeps increasing for the strain hardening case while it decreases for the strain softening case.

To have a better understanding of the stress evolution, the constitutive relation (Equation 49) may be revisited. The stress state depends on both lithiation-induced strain components and plastic strain components. In the elastic case, the hoop stress in the shell is always compressive. When plastic deformation is considered, the stress state depends on the competition between the strain components. The lithiation-induced strain is isotropic and positive during insertion for the elastic case. While for plasticity, the core remains elastic and yielding occurs in the shell region, with the radial plastic strain always being positive, (as can be seen from Figure 6

ε p θ = -1 2 ε p r .
This means that the plastic deformation leads to radial expansion and tangential shrinkage. As a result, the hoop stress can be given by: pression to tangential tension explains the surface cracking of electrodes that occurs experimentally upon maximum lithiation.

σ θ = E (1 + υ) (1 -2υ) u r -αc * + υ du dr -αc * + E 2 (1 + υ) ε p r (
It should be pointed out that this model is based on the small deformation framework, so the numerical results cannot be compared quantitatively with experimental measurements. For quantitative comparison purposes, models with finite strain deformation are more appropriate. Our conclusions are in fact consistent with those obtained using finite strain [START_REF] Chen | A phase-field model coupled with large elasto-plastic deformation: application to lithiated silicon electrodes[END_REF]. Even though the predictions of our model are not quantitative, the main effects are captured and the model results can provide a complete overview of the different scenaria that can occur depending on the model parameters, and provide a comparison between the different cases, something which is more difficult to do in finite strain and has not been done before.

Comparison between elastic and plastic deformation

To consider the effects of plasticity on DISs, the stress evolution under elastic deformation and elasto-plastic deformation in a spherical Si electrode during lithiation is calculated. To make the comparison less complicated, CM coupling is not considered here. As a result, the same concentration distribution (Figure 8(a)) is obtained when the same set of parameters are employed. At early stages (100 ∆t), phase equilibrium has not been reached. As a result, the Li concentration is higher near the outer surface and there are almost no Li-ions in the interior region. As lithiation continuous (900 ∆t, and 2000 ∆t), a core-shell-like structure forms with the Li concentration being close to 1 and 0, respectively. The core and shell are separated by a sharp interface. As lithiation progresses, the interface also moves towards the center of the particle. And the Li concentration in the shell increases slightly. A reduction in the elastic modulus during Li-insertion is also taken into account (based on Equation 57).

The parameters are the same as in the previous section except that k = -5 9 E 0 and H = 0.01E 0 . For plastic deformation, only strain hardening is considered.

When the material properties of the core and shell are taken as being the same, the stress evolutions are as shown in Figure 8 a reduction in the modulus with increasing concentration is considered, much smaller stresses are predicted in the elastic case, while for the plastic case, higher stresses are predicted. The stress evolution over time is also investigated, and it should be noted that there is a visible difference. The radial and hoop stresses (Figure 9(g)-(h)) in the center are always positive for elastic deformation while they are only positive for the initial lithiation stages and then become negative at later stages under elasto-plastic deformation. The hoop stress on the outer surface is always compressive in the elastic case while for the plastic case it's compressive only for the initial lithiation stages and transits from compression to tension at later stages. The stresses predicted from the elastic model are much higher than those the from plastic model, and the elastic model cannot explain the experimentally observed radial cracks during lithiation. This is due to the fact that the hoop stress near the outer surface is the crack driving force and compression cannot cause cracks to grow. As can be seen in Figure 8(h), the elastic model predicts a compressive hoop stress all the time during lithiation while the plastic model predicts compression only at early stages and it predicts tension at late stages, which allows for crack initiation and is consistent with experimental evidence showing crack formation right before Li de-insertion. As a result, the plastic model is more appropriate for high capacity electrodes coupled with large volume changes and phase transformations.

Conclusions

In this work, phase field models accounting for DISs in spherical active materials that undergo phase-transformations are developed. For electrodes with relatively small volume variations, elastic models can be employed while for electrodes with large volume changes, plastic models are preferred. The models accounts for the effect of phase. The sharp phase boundary is naturally captured by the phase field model, and the concentration field is obtained by a mixed formulation of the fourth-order Cahn-Hilliard equation. DISs are obtained by solving the variational form of the mechanical equilibrium equations.

It is found that the DISs arise from the inhomogeneous volume expansions resulting from Li concentration gradients and the hydrostatic stress facilitates the diffusion of Li-ions under elastic deformation while it hinders diffusion in the plastic case. For elastic electrodes, higher radial, hoop, hydrostatic, and von Mises stresses are predicted for the coupling cases when compared with the ones without CM coupling at the same time step. This is due to the fact that the coupling model predicts a smaller miscibility gap (concentration differences between the Li-rich phase and Li-poor phase), thus, less volume change and DISs are induced. For elasto-plastic solids, it's the opposite. The coupling model predicts a higher miscibility gap and the electrodes experience a slightly higher von Mises stress and equivalent plastic strain when it's compared with the ones without CM coupling at the same time step. When the lithiation results in a lower elastic modulus, as opposed to having the elastic properties independent of concentration, a decrease in DISs but an increase in strain occur under elastic deformation. The opposite occurs for the plastic case. Under elastic deformation, the radial stress is always tensile while the hoop stress is tensile in the core region and compressive in the shell. For the plastic case, the radial stress shows a transition from tension in initial stages to compression at later stages. The hoop stress in the core region also shows a similar trend while the hoop stress in the shell shows a transition from compression to tension, which is believed to be the cause of experimentally observed surface cracking during the Li insertion process. Furthermore, if strain softening due to plastic deformation is assumed, smaller stresses and higher plastic strains are predicted than the strain hardening case. To sum up, concentration-dependent material properties due to Li insertion and hardening behavior of the material due to plastic deformation plays a significant role on DISs in spherical phase transformation electrodes.

By taking these factors into consideration, more accurate predictions of the DISs can be achieved, thus providing an improved theoretical basis and insight for designing next-generation mechanically stable phase transforming electrode materials.

  in these anodes can provide some insights on the stress evolution, plastic models Hu et al. (2010); Zhao et al. (2011b) would be more reasonable since the volume change is large enough that yielding can occur. Both in-situ measurementsSethuraman et al. (2010) and first-principle studies Zhao et al. (2011a) of the stress evolution in thin-film amorphous Si anodes during lithiation/delithiation have revealed that yielding and plastic flow occur. Elastoplastic and viscoplastic models have been adopted to study the DISs in Si and Sn anodes during lithiation Chen et al. (2017); Huang et al. (2013); Liu et al. (2012); Chen et al. (

  Figure 1(a) shows the radial distributions of the normalized Li concentration at lithiation times 500∆t, 2000∆t and 4000∆t, respectively. A sharp interface is present between the Li-rich and Li-poor phases. At the initial stage (500∆t), the differences in the concentration profiles are not distinct, and the stresses and strains ( Figure 1(b)-(h)) have similar radial distributions. As lithiation continues 2000∆t and 4000∆t, the Li-poor phase tends to have a higher concentration when compared with the case without CM coupling, and the phase interface tends to move faster in the coupling case. The radial stress (Figure 1(b)), is almost constant in the core (Li-poor) region and gradually decreases to zero on the outer surface. The hoop stress (Figure 1(c)) is also constant in the core-region and transits from tensile in the core to compressive in the shell. The hydrostatic stress (Figure 1(d)) shows a similar trend. The Li-poor region is under constant tension while the Li-rich shell is in compression. From Equation (21), it can be seen that the compressive

Figure 1 :

 1 Figure 1: Radial distribution of (a) Concentration profiles, (b)Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress at time steps 500∆t, 2000∆t, and 4000∆t with/without chemo-mechanical coupling effects. The evolution of (f) Normalized hoop stress on the outer surface. All the stress components are normalized by the Young's modulus of the delithiated electrode material. The unit for the time step is ∆t. The solid and dashed lines are for with and without chemo-mechanical coupling effect respectively.

Figure 2 :

 2 Figure 2: Radial distribution of (a) Concentration profiles, (b)Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress, and (f) Radial strain, (g)Hoop strain and (h) Equivalent strain at time step 2000 ∆t with chemo-mechanical coupling effects under elastic deformation. Constant modulus, reduced modulus and increased modulus during ion insertion are considered. All the stress components are normalized by Young's modulus of the delithiated electrode material. The dotted and dashed lines are for a 10% reduction or increase in the elastic modulus during ion-insertion respectively. As for the solid lines, a constant elastic modulus is assumed.

Figure 3 :

 3 Figure 3: Radial distribution of (a) Normalized hoop stress in the particle center, and (c) Normalized hoop stress at the outer surface at time step 2000 ∆t with chemo-mechanical coupling effects under elastic deformation. Both of the modulus increasing or decreasing during ion insertion are considered. All the stress components are normalized by the Young's modulus of the delithiated electrode material. The dotted and dashed lines are for a 10% reduction or increase in the elastic modulus during ion-insertion respectively. As for the solid lines, a constant elastic modulus is assumed.

  (a)) and on the outer surface (Figure 1(f) and Figure 3(b)) are recorded. Radial stress in the particle center is equal to hoop stress in the particle center due to spherical symmetry. It can be seen that the radial/hoop stress in the center increase faster in the early lithiation stages than in the later stages. While for the hoop stress on the outer surface, the magnitude increases quickly during early lithiation stages and then gradually decreases. The radial and hoop stresses in the center are tensile while the hoop stress on the outer surface is always compressive. Due

  |x|) denoting the positive part of x. For J 2 plasticity and associative flow rule, ∂f ∂σ : I = s s : I = tr(s)

  05E 0 Chen et al. (2014), H = ±0.01E 0 (corresponds to strain hardening and softening respectively) Chen et al. (2014), Ξ=3.1×10 -5 m 3 /mol, c max =3.67 × 10 5 mol/m 3 , T =300 K, ∆t=0.0005, α= 3 √ 4.2-1=0.613 (corresponding to 420% volume expansion after full insertion). At the early stage of lithiation (Figure 4), the differences in concentration, stresses and strain distributions are small. The concentration profiles (Figure 4(a)) almost coincide with each other. The radial stresses (Figure 4(b)) are almost constant in the core region and gradually decrease to zero as they approach the shell; the hoop stresses (Figure 4(c)) transit from tension in the core to compression in the shell; and the hydrostatic stresses (Figure 4(d)) show similar trends as the hoop stresses. The resulting von Mises stresses (Figure 4(e))

Figure 4 :

 4 Figure 4: Radial distribution of (a) Concentration profiles, (b)Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress, (f) Radial plastic strain, (g) Hoop plastic strain and (h) Equivalent plastic strain at time step 100 ∆t with/without chemo-mechanical coupling ('w cp' / 'w/o cp') effects under elastoplastic deformation. Both strain hardening and strain softening cases are considered. All the stress components are normalized by Young's modulus of the delithiated electrode material.

  Figure 5(c), and Figure 6(c)) and hydrostatic stress (Figure 4(d), Figure 5(d), and Figure 6(d)) show similar trends for the three cases considered: strain hardening without CM coupling ('w/o cp'), strain hardening with CM coupling ('w cp') and strain softening without CM coupling. They all (stress/strain profiles) remain constant in the core regions with strain hardening without CM coupling being the largest and strain softening without coupling being the lowest in magnitude. The radial stress is compressive in both the core and shell, while the hoop and hydrostatic stresses are compressive in the core and transit from compression near the interface to tension in the shell exterior. This hydrostatic stress state hinders the diffusion of Li-ions, which is the reason for the slower migration of the reaction front. There are stress drops in the radial stress, hoop stress and hydrostatic stress profiles near the core-shell interface due to the occurrence of plastic deformation. The von Mises stress profiles (Figure 5(e), and Figure 6(e)) also show stress drops but the situation is a little bit different. The von Mises stresses (Figure 4(e), Figure 5(e), and Figure 6(e)) are almost zero in the core regions. As a result, the core remains elastic and no plastic deformation (Figure 4(f)-(h), Figure 5(f)-(h), and Figure 6(f)-(h)) occurs there.

Figure 5 :

 5 Figure 5: Radial distribution of (a) Concentration profiles, (b)Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress, (f) Radial plastic strain, (g) Hoop plastic strain and (h) Equivalent plastic strain at time step 1500 ∆t with/without chemo-mechanical coupling ('w cp' / 'w/o cp') effects under elastoplastic deformation. Both strain hardening and strain softening cases are considered. All the stress components are normalized by Young's modulus of the delithiated electrode material.

  (f), Figure 4(f), Figure 5(f)) the hoop plastic strain always being negative (as can be seen from Figure 6(g), Figure 4(g), Figure 5(g)), and

  60)The first part on the right-hand side of Equation 60 is the solution for the hoop stress under elastic deformation only, arising from mismatch strains due to the concentration gradient. The region close to the outer surface (Li-rich region/shell) contains significantly more Li-ions than the inner region (Li-poor region/core). As a result, the Li-rich regions expand to a larger degree. However, this expansion is constrained by the surrounding materials and the hoop stress in the shell is compressive initially. Due to the phase transformation, large stresses develop near the phase interface. When the effective stress (von Mises stress) is large enough, plastic deformation emerges near the phase interface. The second part of Equation 60 represents the contribution due to plastic deformation and it's always positive in the shell. As lithiation proceeds, the plastic region propagates and plastic deformation accumulates. When the stress component due to plastic deformation surpasses the component due to the concentration gradient, the hoop stress in the shell undergoes a transition from compression to tension. This is believed to be the cause of experimentally observed cracks during the lithiation process of spherical Si/Sn electrode particles. For example in paper Aifantis et al. (2012) it was shown that Sn and Sn nanoparticles severely fractured upon the first lithiation. Significant fracture was also observed in Si nanoparticles Liu et al. (2012); McDowell et al. (2012) during the first cycle upon Li-insertion. Particularly, in Liu et al. (2012); Mc-Dowell et al. (2012) the cracks were observed to propagate inward along the radial direction during the lithiation process. The transition of tangential com-

Figure 6 :

 6 Figure 6: Radial distribution of (a) Concentration profiles, (b)Normalized radial stress, (c) Normalized hoop stress, (d) Normalized hydrostatic stress, (e) Normalized von Mises stress, (f) Radial plastic strain, (g) Hoop plastic strain and (h) Equivalent plastic strain at time step 3000 ∆t with/without chemo-mechanical coupling ('w cp' / 'w/o cp') effects under elastoplastic deformation. Both strain hardening and strain softening cases are considered. All the stress components are normalized by Young's modulus of the delithiated electrode material.

Figure 7 :

 7 Figure 7: The evolution of (a) Normalized hoop stress in the particle center and (b) Normalized hoop stress at the outer surface under elasto-plastic deformation. All the stress components are normalized by Young's modulus of the delithiated electrode material. The unit for time step is ∆t.

  (b),Figure 8(d),and Figure 8(f). At the initial stage, the deformation of the electrode remains elastic. The radial and hoop stresses are the same in the core region. The radial stress gradually decreases to zero towards the outer surface, and there is an abrupt transition in the hoop stress from tension in the core to compression in the shell. The maximum von Mises stress occurs near the interface. Under elasto-plastic deformation, yielding begins to occur near the outer surface during the initial lithiation stages, and the stress evolutions are almost the same as in the previous section. The hoop and radial stresses in the core transits from tension to compression at later stages, while the opposite trend is seen for the shell. The von Mises stress is almost zero in the core and increases to a local maximum first, then decreases to a local minimum and finally keeps increasing towards the shell. When a reduction in the elastic modulus with increasing concentration is considered, similar distributions (Figure 8(c), Figure 8(e), and Figure 8(g)) are found for the stresses. Comparison between the DISs with and without concentration-dependent elastic moduli are also made under elastic (Figure 9(a)-(c)) and plastic (Figure 9(d)-(f)) deformation, respectively. It is found that when

Figure 8 :

 8 Figure 8: (a) Concentration profiles at time step 100 ∆t, 900 ∆t, and 2000 ∆t respectively. Radial distribution of normalized radial stress (σr), normalized hoop stress (σ θ ) and normalized von Mises stress (σv), with constant modulus, under elastic and elasto-plastic deformation at time steps of (b) 100 ∆t, (d) 900 ∆t, and (f) 2000 ∆t respectively. Radial distribution of normalized radial stress, normalized hoop stress and normalized von Mises stress, with reduction in modulus with increasing concentration, under elastic and elasto-plastic deformation at time steps of (c) 100 ∆t, (e) 900 ∆t, and (g) 2000 ∆t respectively. (h)The evolution of Normalized hoop stress on the outer surface under elastic ('Ela') and elasto-plastic ('Pla') deformation. All the stress components are normalized by the Young's modulus of the delithiated electrode material. The unit for time step is ∆t. For the stress profiles ((b)-(g)), the solid lines are for the elastic deformation case and the solid lines with asterisks are for the elasto-plastic deformation case.

Figure 9 :

 9 Figure 9: The effects of a reduced elastic modulus with lithiation on the radial distribution of (a)Normalized radial stress (σr), (b) Normalized hoop stress (σ θ ), and (c) Normalized von Mises stress (σv) under elastic deformation. The effects of a reduced elastic modulus with lithiation on the radial distribution of (d) Normalized radial stress, (e) Normalized hoop stress, and (f) Normalized von Mises stress under elasto-plastic deformation. The evolution of (g) Normalized radial stress and (h) Normalized hoop stress in the center. All the stress components are normalized by the Young's modulus of the delithiated electrode material. The unit for time step is ∆t. For the stress profiles ((a)-(f)), the solid lines are for the constant elastic modulus case and the dashed lines are for the case with reduction in modulus with increasing concentration.

Table 1 :

 1 Young's modulus of some active electrode materials in lithiated and

		delithiated phases
		Delithiated	Method	Young's	Lithiated	Young's
		phase		modu-	phase	modu-
				lus(GPa)		lus(GPa)
		CoO 2	Computation	59.8	LiCoO 2	264
			Qi et al.
			(2014)
		C(graphite)	Computation	32	LiC 6	109
			Qi et al.
			(2010)
		Si(amorphous) Computation	96	Li 3.75 Si	41
			Shenoy et al.
			(2010)
		Si(crystaline) Experiment	90	Li 3.75 Si	12
			Hertzberg
			et al. (2011)
		Sn(crystaline) Computation	51	Li 3.5 Sn	24.7
			Stournara
			et al. (2012)
		and the expansion is isotropic. To model the sharp phase boundary, the Cahn-
		Hilliard type diffusion equation can be used. The effect of chemo-mechanical
	150	coupling on DISs can be incorporated into the model via a variational formu-
		lation approach. As for the electric field, the conductivities of the electrode
		materials are high enough in most cases so that the electric potential inside the
		electrode particles are regarded as uniform, therefore, the influence of electric
		field on DIS can be neglected Zuev & Tsvetkov (2017).
		Both, the concentration and temperature fields are governed by the diffusion
		equation. Therefore, it's naturally to deduce that a change in concentration
		will cause DIS in a similar manner as a change in temperature would lead to

  typical phase-transformation cathode materials are LiFePO 4 , LiCoO 2 and LiMn 2 O 4 . Although all these materials undergo phase transformations, it should be noted that LiFePO 4 has a clear concentration jump during lithiation between the different phases and parameter χ in Equation 8 must be greater than 2 to ensure phase separation. However, LiCoO 2 and LiMn 2 O 4 undergo a structural transformation and there is no clear concentration jump. In such cases, χ must be less than 2 and in particular, it can be set as 0.

Table 2 :

 2 Material properties of LiFePO 4

	Name	Symbol and unit	Value
	Diffusion coefficient	D[m 2 /s]	10 -14
	Young's modulus	E [GPa]	124.5
	Poisson's ratio	ν	0.25
	Partial molar volume	Ξ [m 3 /mol]	2.9× 10 -6
	Gradient coefficient	κ [J/m]	5.02× 10 -10
	Expansion coefficient	α	2.21× 10 -2
	Maximum concentra-	c max [mol/m 3 ]	2.29× 10 4
	tion

Table 3 :

 3 Normalized/dimensionless variables

	Name		Symbol	Normalized variable
	Normalized gradient co-	κ *		κc max / RT ρ 2
	efficient
	Normalized mobility	M *	M RT t 0 /ρ 2
	Normalized time	t *		t/t 0
	Normalized radial dis-	r *		r/ρ
	tance
	Normalized stress com-	σ ij	*	σ ij /E 0
	ponents
	Normalized	chemical	µ *		µ/ (RT )
	potential
	Normalized ion flux	i *		iρ/ (zF Dc max )
	Normalized partial mo-	Ξ
	lar volume

* ΞE/ (RT )

It takes approximately 6000 ∆t to complete the lithium insertion process, i.e. reach 100% lithiation of the active particle.

Acknowledgments

The authors are grateful to the National Science Foundation for supporting this work through the CMMI grant (CMMI-1762602).