Monica Del

Carmen Fong Lopez

M Rigal

L Rigal

G Vilarem

V Vandenbossche

Influence of temperature and soda concentration in a thermo-mechanochemical pretreatment for bioethanol production from sweet corn coproducts

Keywords:

Sweet corn co-product

Introduction

The rise in global population necessities and upsurge in food and energy demands lead to the current environmental crisis. For instance, an increment in Green House Gas emissions forecasts temperature rises and climate perturbations [START_REF] Schleussner | Science and policy characteristics of the Paris Agreement temperature goal[END_REF]. Therefore, to ensure sustainable food security, preservation of soil quality, the biodiversity and water reserves need to be addressed [START_REF] Godfray | Food security: the challenge of feeding 9 billion people[END_REF]. Additionally, with the increasing scarcity of fossil fuel, development of sustainable alternatives for transport energy are required [START_REF] Leggett | The implication for climate change and peak fossil fuel of the continuation of the current trend in wind and solar energy production[END_REF][START_REF] Solomon | Biofuels and sustainability[END_REF]. Transport energy demand solutions comprise bioethanol from renewable resources [START_REF] Zabed | Bioethanol production from renewable sources: current perspectives and technological progress[END_REF]. Bioethanol from lignocellulosic feedstocks (non-food raw material), such as agricultural residues, provide lower GHG emissions and fossil energy demand. However, its environmental impact significantly depends on land use conditions [START_REF] Cherubini | Crop residues as raw materials for biorefinery systemsa LCA case study[END_REF].

Lignocellulosic bioethanol is obtained by a sequential unit operations i.e.: pretreatment, saccharification, fermentation, and distillation [START_REF] Zabed | Fuel ethanol production from lignocellulosic biomass: an overview on feedstocks and technological approaches[END_REF]. At present, overcoming the technical barriers at each unit operation dictates the economic feasibility of the production process. Concerning the pretreatment phase, optimization to increase enzymatic hydrolysis and fermentation efficiency while maintaining low operational costs are the major challenges [START_REF] Chiaramonti | Review of Pretreatment Processes for Lignocellulosic Ethanol Production, and Development of an Innovative Method. Biomass and Bioenergy[END_REF]. The biomass complexity in its physiochemical structure limits the enzymatic performances to deconstruct the carbohydrate matrix (Himmel et al., 2007). Thus, various physical, chemical, physicochemical and biological technologies have been applied to lignocellulose to enhance the enzymatic accessibility; [START_REF] Alvira | Pretreatment technologies for an efficient bioethanol production process based on enzymatic hydrolysis: a review[END_REF][START_REF] Chiaramonti | Review of Pretreatment Processes for Lignocellulosic Ethanol Production, and Development of an Innovative Method. Biomass and Bioenergy[END_REF][START_REF] Putro | Pretreatment and conversion of lignocellulose biomass into valuable chemicals[END_REF][START_REF] Sun | Hydrolysis of lignocellulosic materials for ethanol production: a review[END_REF]. Currently, few technologies such as steam explosion, dilute acid, and alkali pretreatment have reached demonstration and commercialization phases, and their economic competitiveness are under study [START_REF] Irena | Innovation Outlook Advanced Biofuels 2016 Summary | Biofuel | Gasification.Innovation Outlook Advanced Biofuels 2016 Summary | Biofuel[END_REF].

Alkaline pretreatment technology has been extensively studied using different reagents, such as sodium hydroxide [START_REF] Chen | Understanding of Alkaline Pretreatment Parameters for Corn Stover Enzymatic Saccharification[END_REF], ammonia [START_REF] Prior | Hydrolysis of ammonia-pretreated sugar cane bagasse with cellulase, β-Glucosidase, and hemicellulase preparations[END_REF], potassium hydroxide ( [START_REF] Sharma | Potential of potassium hydroxide pretreatment of switchgrass for fermentable sugar production[END_REF]) and lime [START_REF] Kaar | Using lime pretreatment to facilitate the enzymic hydrolysis of corn stover[END_REF]. These alkali catalysts are known to disrupt the ester bonds between lignin and hemicelluloses, thereby providing partial or total lignin solubilization, resulting in the production of sugar oligomers and polymers. Furthermore, alkaline addition diminishes cellulose crystallinity by swelling which increases the internal surface area. In terms of process, temperature pretreatment https://doi.org/10.1016/j.indcrop.2019.03.044 hemispherical dishes with conical holes was used on module 6 to enable the filtrate to be collected. Feedstocks were fed into the extruder's first module using a feeder model number KCL-KT20 (K-tron, Clextral, France). Two piston pumps (DKM SUPER J PP 6.35, DKM SUPER K PP 16, Clextral, France) were used to inject, respectively, alkaline solutions of sodium hydroxide (NaOH) and acid solutions of sulfuric acid (H2SO4).

The conditions of extrusion were based on [START_REF] Vandenbossche | Suitability assessment of a continuous process combining thermomechano-chemical and bio-catalytic action in a single pilot-scale twin-screw extruder for six different biomass sources[END_REF] study at pilot scale. In this previous study, the pretreatment was conducted using sodium hydroxide at 10% (w/w) and barrel temperature of 100 °C, the neutralization was carried out with phosphoric acid at 3.5% (w/w)( [START_REF] Vandenbossche | Suitability assessment of a continuous process combining thermomechano-chemical and bio-catalytic action in a single pilot-scale twin-screw extruder for six different biomass sources[END_REF]). In this case, the screw rotation of the laboratory scale extruder was set at 300 rpm. The feeder screw discharged continuously the sweet corn coproduct at 1.8 Kg/h to obtain 1.6 Kg of dry sweet corn co-product per hour. Both pumps of reagents were fixed to maintain liquid solid ratios of 1.0 (Kg/Kg) in alkaline site and 4.2 (Kg/Kg) in acid introduction which correspond to a flow rate of 1.65 kg/h for alkali solution and 5.80 Kg/h for sulfuric acid solution. The pH of pretreated biomass was controlled to ensure that biomass was neutralized. The pH of extrudate and filtrate aimed were 5 and 4 respectively. Both extrudate and filtrate samples were stored in plastic containers at -20 °C.

Pretreatment and experimental design

To study the influence of NaOH and temperature in the pretreatment, the NaOH loading per g of dry matter of sweet corn co-product and the internal temperature monitored by the probe 1 (TP1) were evaluated. Pretreatment conditions were inspired by a 2 factors Dohelert experimental design. The experimental design consisted of a center at 120 °C and 6% variation steps of 60 °C for temperature and 2% for NaOH/SCC. The responses analyzed were the composition of extrudate (insoluble fiber, celluloses, hemicelluloses and lignin), hemicelluloses removal, the filtration efficiency parameters (extrudate dry matter and filtrate dry matter) and the enzymatic hydrolysis (enzymatic conversion and sugar yield).

The statistical analyses were performed with NemrodW®. Once the data was collected, the responses were normalized, then fitted to the empirical equations (second order polynomial regression equations) for identification of main variables:

= + + + + + Y a a X a X a X X a X a X O 1 1 2 2 12 1 2 11 1 2 22 2 2 (3.1)
Where the response is Y, the variables X 1 and X 2 represent the normalized parameters, alkaline loading and the internal temperature measured, respectively. The responses are correlated to the center of the domain a 0 , linear (a 1 and a 2 ), interactions (a 12 ) and quadratic (a 11 and a 22 ) coefficients. The coefficients a 1 and a 2 corresponding to the main effects of NaOH/SCC and internal temperature respectively.

Analytical methods

Extrudate analysis

= - × HR H in initial SCC H in extrudate H in initial SCC (( ))/( ) 100
(3.2) = × CY C inextrudate / C initial SCC (( ) ( )) 100 (3.3)
The French standard NF V 03-322 was used to determine moisture content and mineral content. An estimation of the three parietal constituents contained in the extrudates was performed using the ADF-NDF method from Van Soest and Vine (Soest et al., 1968). This method allows the quantification of insoluble fraction (IF), cellulose, hemicelluloses and lignin contained in the biomass. Insoluble fraction is obtained after neutral detergent wash and corresponds to the structural components in the fiber. All determinations were carried out in triplicate with a standard deviation < 1.5% for all measurements. Using presents advantages such as an affordable cost, high conservation of cellulosic material, low enzymatic production and fermentation inhibitors such as HMF and furfural [START_REF] Kim | A review on alkaline pretreatment technology for bioconversion of lignocellulosic biomass[END_REF][START_REF] Modenbach | Sodium Hydroxide Pretreatment of Corn Stover and Subsequent Enzymatic Hydrolysis: an Investigation of Yields, Kinetic Modeling and Glucose Recovery[END_REF][START_REF] Rodrigues | A molar basis comparison of calcium hydroxide, sodium hydroxide, and potassium hydroxide on the pretreatment of switchgrass and miscanthus under high solids conditions[END_REF].

Therefore, by coalescing alkaline pretreatment with twin screw extrusion, a continuous, adaptable and a high solid efficient process can be achieved [START_REF] Vandenbossche | Thermomechanical and thermo-mechano-chemical pretreatment of wheat straw using a twin-screw extruder[END_REF]. As a pretreatment, extrusion is characterized by shearing and pressure forces applied by the screw rotation, thus, leading to particle size reduction and fibrillation of the biomass [START_REF] Zhang | Pretreatment of corn stover for sugar production with combined alkaline reagents[END_REF]. The combined effects of extrusion with chemical pretreatment have been previously underlined: higher particle size reduction, solubilization of biomass components and sugar yields by enzymatic hydrolysis [START_REF] Duque | Extrusion as a pretreatment for lignocellulosic biomass: fundamentals and applications[END_REF]. In a recent study, Liu et al., 2018 compared alkaline pretreatment (AP) and alkaline twinscrew extrusion (ATSE) pretreatment on corn stover. According to the authors, mixing and chemical penetration in the extruder combined with the fibrillation enhance the lignin solubilization and therefore, the accessibility to cellulose and hemicelluloses of enzymes [START_REF] Liu | Comparative study of two different alkali-mechanical pretreatments of corn stover for bioethanol production[END_REF]. Comparison of alkali and alkali extrusion on corn residues was also analyzed by [START_REF] Khatri | Determination of optimal biomass pretreatment strategies for biofuel production: investigation of relationships between surface-exposed polysaccharides and their enzymatic conversion using carbohydrate-binding modules[END_REF]. Results indicate a loss of hemicelluloses at the biomass surface with both alkaline pretreatments, but lower crystalline cellulose was detected in alkali-extruded biomass [START_REF] Khatri | Determination of optimal biomass pretreatment strategies for biofuel production: investigation of relationships between surface-exposed polysaccharides and their enzymatic conversion using carbohydrate-binding modules[END_REF] Earlier studies have executed a continuous thermo-mechano-chemical-technology using twin screw extrusion [START_REF] Vandenbossche | Bio-catalytic action of twinscrew extruder enzymatic hydrolysis on the deconstruction of annual plant material: case of sweet corn co-products[END_REF] with successful pilot scale tests using six different biomass sources [START_REF] Vandenbossche | Suitability assessment of a continuous process combining thermomechano-chemical and bio-catalytic action in a single pilot-scale twin-screw extruder for six different biomass sources[END_REF]. This novel process includes different steps within the extruder, namely; alkaline pretreatment, neutralization, filtration, and enzymatic impregnation. Alkali pretreatment in the extruder prepares the biomass to have greater accessibility to the enzymes; whereas a neutralization phase acidifies the biomass and filtration step eliminates the yeast and enzymatic inhibitory compounds (salts and soluble molecules) [START_REF] Deparis | Engineering tolerance to industrially relevant stress factors in yeast cell factories[END_REF]. The results from the previous studies demonstrate the efficacy of alkaline pretreatment and the potential of bioextrusion to reach a viable techno-economic process to produce bioethanol. Thereby, in the present work, the influence of temperature and alkaline concentration in reactive extrusion is analyzed using the design of experiments (DOE) with an objective to produce an accessible and cellulose-rich substrate with a lower chemical load. For this matter, the parameters determined by DOE were applied to extrusion using sweet-corn co-product.

Materials and methods

Materials

Feedstock

Dehydrated sweet corn (Zea mays L. saccharata) co-products (SCC) were obtained from industrial corn canneries and were provided by SARL Soupro+ (Castelmoron sur Lot, France). It was milled using a hammer mill fitted with a 6 mm screen.

Twin-screw extrusion

The extrusion process was performed using a twin-screw extruder (Evolum 25, Clextral, France), composed of 10 modules with a length of 100 mm each. The extruder is configured to combine three different processing steps: alkaline pretreatment, neutralization, and filtration of the biomass (Fig. 1). The screw diameter is 25 mm. Four types of screws were used: a trapezoidal double-thread screw used for the feeding zone, conveying double-thread, screws bilobe paddle and reversed pitch double-thread screws used to produce transport, mixing, and shearing effects, respectively, along with the different zones of the process. A heating band thermoregulated modules and cooled by water circulation. Two internal temperature probes were introduced, in module 3 to monitor the temperature of alkali impregnation and module 8 to ensure matter temperature below 50 °C A filter section consisting of six these data and the extrudates flow rates, the mass balance of insoluble fiber, cellulose, hemicelluloses, and lignin were determined. The hemicellulose removal and the cellulose conservation rates were calculated by the following two equations (Eqs. (3.2) and (3.3)):

Where H stands for hemicelluloses and C to cellulose. The hemicelluloses and cellulose are expressed in mass.

Enzymatic hydrolysis

Enzymatic hydrolysis of extrudates were conducted at a consistency of 5% (W/W) on 0.5 g (extrudate dry matter) at 50 °C using a shaken hot-water bath and 25 mL glass Pyrex tubes for 48 h. Citrate-phosphate buffer (pH = 5) was added to maintain the pH at 5, and 0.02% sodium azide was used in the mixture to avoid microbial contamination. The commercial enzymatic cocktails used for the study were provided by Novozyme (Franklintong, USA), Cellic Ctec 2 and hemicellulases NS22002. The protein content of each cocktail was determined by the Lowry method (Biorad Protein assay kit), and cellulase activity was determined using the method reported by Ghose (Ghose, 2009). The protein content of Cellic Ctec 2 was 254.9 g/L and cellulase activity of 123.3 FPU/mL, and NS22002 protein content was 36.9 g/L. The Cellic Ctec 2 loading for saccharification assays was 15 FPU/g of the dry extrudate and was complemented NS22002 at 10% (protein basis).After 48 h of hydrolysis, the samples were diluted in 50 mL total volume. Afterwards, they were centrifuged and the supernatant was filtered through a 0.22 μm membrane. The hydrolysate was analyzed with the HPIC (High-performance ionic chromatography, Dionex, France) system using a Carbopac PA1 column. The efficiency of enzymes to release glucose and xylose was studied by glucose and xylose conversion rates, whereas the sugar yields evaluated the performance of the process to produce both monosaccharides using Eq. 3. Where m is the mass in grams. The mass of dry SCC was calculated as the ratio (Q ext /Q SCC )*0.5, where Q ext is the flow rate of oven-dry extrudate and Q SCC is the flow rate of dry matter of sweet corn coproduct. Data presented in this paper are the averages of the results from triplicate experiments with standard deviation lower than 5%.

Results and discussion

Raw biomass analysis

The corn husk and cob form the sweet corn co-product. After harvesting, it was dried by the supplier to a final dry matter of 89.9%. Characteristics in terms of composition and enzymatic hydrolysis are presented in Table 1. SCC has a large lignocellulosic fraction, 75.2% formed by a large portion of carbohydrates, hemicelluloses (39%) and cellulose (31%) and a low lignin amount (5%). Soluble sugars content was analyzed. Glucose concentration was 0.013 g/g of biomass and traces of xylose, arabinose and galactose were detected at minor concentrations. Due to its high carbohydrate content and low lignin content, SCC is a promising substrate to produce bioethanol by enzymatic hydrolysis and co-fermentation. In terms of availability, canned industry in south-west of France reported production capacity between 280,000-390,000 t/year of SCC. As it is currently used for animal feed, it is considered as a waste, hence, a high potential feedstock for upcycling with economic viability.

Effect of soda loading and temperature

Extrudate and filtrate composition

According to [START_REF] Alvira | Pretreatment technologies for an efficient bioethanol production process based on enzymatic hydrolysis: a review[END_REF][START_REF] Alvira | Pretreatment technologies for an efficient bioethanol production process based on enzymatic hydrolysis: a review[END_REF], as bioethanol is potentially produced by co-fermentation of C5 and C6 sugars, an economically viable pretreatment is based on the cellulose and hemicelluloses recovery. Simultaneously, minimizing the effect of hemicelluloses and lignin as a barrier for saccharification. Alkaline pretreatments are generally known for the disruption of lignin structure and solvation of acetic groups in hemicelluloses, thus, producing solubilization of part of the lignocellulose [START_REF] Kim | A review on alkaline pretreatment technology for bioconversion of lignocellulosic biomass[END_REF]. The analysis of extrudate composition at different conditions was conducted to monitor the effect of NaOH/SCC and temperature on the biomass carbohydrate content and lignin removal. The results are summarized in Table 2 and the statistical analysis in Table 3 and Fig. 2. Additionally, filtration efficiency is ensured by the mechanical force applied by the reverse screw in module 7, forming a dynamic plug that pressures the biomass.

As the pretreatment deconstructs the biomass, the viscosity of the matter rises, affecting the dynamic plug stability and efficiency filtration. As filtration is less effective, moisture content in extrudate is higher and the mineral matter of filtrate is lower. Consequently, dry matter of both outputs were monitored, and results are shown in Table 4 and iso-response curves (Fig. 3). The R 2 of the equation expresses the statistical accuracy of the fit. A significant R 2 means that the predicted values are close to the measured values. Statistical analysis performed on the extrudate composition demonstrates that the fit was accurate for IF, whereas for cellulose and hemicelluloses, R 2 of the models implies that only tendencies can be deducted.

Concerning lignin content in the extrudate, R² was not significant to be analyzed by the equation. For the IF, cellulose, and hemicelluloses, NaOH had a more significant impact than temperature as reflected by the p-values of their coefficients. Iso response curves of IF (Fig. 2(a)) display a strong influence of both parameters to diminish its value. When low NaOH and temperature are applied, IF fraction is essential. To disrupt the IF at 4% soda loading, a temperature higher than 150 °C is necessary. Nevertheless, the increasing soda loading disrupts IF, and its effect is reinforced by temperature.

Iso-response curves of H content (Fig. 2(b)) reveals a similar trend than IF, high H fraction is obtained with mild conditions and an increase in temperature and NaOH/SCC triggered its reduction. With temperatures below 90 °C, soda effect is only observed from 6%. As the temperature reaches 90 °C, NaOH solubilization of hemicelluloses is visible. Variations of H content due to NaOH are less noticeable at temperatures of 100-125 °C, with values between 45-52%. Values below 45% are attained with a temperature exceeding 125 °C and NaOH/SCC beyond 6%. To conclude, the temperature must surpass 50 °C to attain hemicelluloses solubilization. In the same way as IF, effect by soda is enhanced by an increase in temperature. Hence, reduction of IF is attributed to the partial solubilization of the hemicelluloses.

Considering the mass balance of the process, hemicelluloses removal (HR) rate from initial SCC was also analyzed (Fig. 2(d)). The HR is linear to the soda loading and the temperature attained. In the studied conditions, to reach more than 20% HR, both effects must be combined. To cite an instance, solubilization of hemicelluloses at 8% needs 50 °C to reach 20% removal, however, at 4% it needs more than 150 °C. Temperature needs to be attained to observe the NaOH effect on hemicelluloses. Although solubilization has occurred, a significant part of this fraction remained at the end of the process, considering that a minimum of 60% of the initial H is conserved, offering the possibility of improving ethanol yield by saccharification and co-fermentation of xylan.

The heterogeneous and amorphous structure of hemicelluloses makes it more susceptible to degradation by pretreatments than cellulose [START_REF] Saha | Hemicellulose bioconversion[END_REF]. By way of illustration, xylan losses were observed with Alkaline Twin Screw Extrusion pretreatment (ATSE process) by [START_REF] Liu | Alkaline twin-screw extrusion pretreatment for fermentable sugar production[END_REF] [START_REF] Liu | Alkaline twin-screw extrusion pretreatment for fermentable sugar production[END_REF]. During this process, biomass goes through a twin-screw extruder with 4 zones of compression as the alkaline solution is added. As soda loading was increased from 0.06 to 0.1 g/g of oven dry corn stover, the xylan loss raised from 21.7-39.3 %. Moreover, the study of soda concentration and barrel temperature using the present alkali-extrusion process by [START_REF] Duque | Optimization of integrated alkaline-extrusion pretreatment of barley straw for sugar production by enzymatic hydrolysis[END_REF] *indicate that term has a significant effect at 95% confidence interval.

Table 2

Extrudate composition, mass balance indicators and filtrate mineral matter for all samples analyzed. dry barley straw), hemicellulosic fraction breaks down. Nevertheless, the impact of temperature in the conditions tested was not significant [START_REF] Duque | Optimization of integrated alkaline-extrusion pretreatment of barley straw for sugar production by enzymatic hydrolysis[END_REF]. Effect of the temperature on the hemicelluloses recovery after alkaline pretreatment of corn stover was reported by [START_REF] Chen | Understanding of Alkaline Pretreatment Parameters for Corn Stover Enzymatic Saccharification[END_REF]. The statistical analysis, presented for the recovery of hemicelluloses, indicated that temperature plays a key role in the solubilization of hemicelluloses. Pretreatments in literature are difficult to compare due to the variations on the substrate, extrusion profiles and process conditions. However, it is generally agreed that the hemicellulosic structure is disrupted by soda attack, while the temperature effect depends on the conditions applied. Variation of the studied parameters demonstrates that hemicellulosic disruption by NaOH needs temperature, whereas lignin removal was detected at all conditions. At the mildest conditions (Sample 10), lignin concentration was already diminished to 4.7% instead of 6.1% in raw SCC. Nevertheless, delignification behavior is not regulated by the studied parameters, as the fitting had a low R 2 , implying more parameters could be involved in lignin breakdown. Lignin restructuration after hydrothermal pretreatment on corn stem [START_REF] Donohoe | Visualizing lignin coalescence and migration through maize cell walls following thermochemical pretreatment[END_REF] and steam explosion on poplar holocellulose [START_REF] Sannigrahi | Pseudo-lignin and pretreatment chemistry[END_REF] were studied. The results exhibit that lignin fraction can undergo physico-chemical modifications such as re-condensation and formation of pseudo-lignin with carbohydrates occasioned by the temperature and pressure employed. In extrusion pretreatment, the effect of reverse screw elements on corncobs produced lignin re-distribution at the fiber surface, thus, blocking pores [START_REF] Zheng | The effects of screw elements on enzymatic digestibility of corncobs after pretreatment in a twin-screw extruder[END_REF]. In the three studies mentioned, droplets of lignin on the biomass surface were observed by SEM. Although, in the present work lignin restructuration was not analyzed, these examples illustrate its sensibility for pretreatment conditions.

In regard to the cellulose contained in IF (Fig. 2(c)), the statistical analysis shows that its increase is attributed to the combined effect of soda and temperature, as expressed by the linear coefficients a 1 and a 2 of this response. Iso-response curves analysis indicate that NaOH loading lower than 6% at temperatures below 125 °C does not lead to modification of C content in the extrudate. C fraction increases as the NaOH loading is intensified and can attain more than 60% with temperatures higher than 130 °C and NaOH/SCC higher than 6%. The response surface of C content is inverted from the one of H content. The gain of the cellulose portion in the extrudate results from the extraction of hemicelluloses and lignin. Recovering a substrate with a higher cellulose portion than the untreated SCC proves to be an advantage for future saccharification step.

Additionally, cellulose yield (CY) was calculated from the mass balance of the process. A maximum yield of cellulose obtained for soda loadings between 5-7% NaOH/SCC and temperature of 75-125 °C, *indicates that the term has a significant effect at 95% confidence interval.

reached 90%. Globally, alkali extrusion preserves cellulose through the overall process and variations in different conditions were not significant. However, CY could be affected by the temperature, since at NaOH loading higher than 8% and temperature > 150 °C, the CY decreased. The results obtained imply the possibility of cellulose solubilization or charred as a result of the excessive severity of the pretreatment, as encountered by Silva et al., 2013 when using extrusion with ionic liquid at 180 °C [START_REF] Silva | Continuous pretreatment of sugarcane bagasse at high loading in an ionic liquid using a twin-screw extruder[END_REF].

As mentioned earlier, a statistical analysis of the extrudate and filtrate mineral was carried out to study the impact of parameters on filtration efficiency (Fig. 3). The fit is accurate for both responses (Table 4). Concerning extrudate dry matter, the coefficients exhibit that extrudate solid content is modified by the combination of alkali and temperature. It is negatively impacted by increasing reagent input and temperature augmentation (Fig. 3(a)). The liquid /solid ratio of base and acid were kept constant in the extrusion process. Therefore, variations of the solid content of the extrudate reflect the filtration performance. Study of IF and H content in extrudate demonstrates the partial lignocellulose solubilization takes place during pretreatment, thereby leading to variations of the biomass viscosity and the dynamic plug stability. Along with the produced IF solubilization, the dynamic plug is softened leading to a decrease in extrudate dry matter.

With respect to the filtrate mineral matter (Fig. 3(b)), the coefficients indicate that the soda loading mainly modifies it. As chemicals are added, the content of salts in the filtrate augment, however, as the temperature is increased at high chemical input, there is a decrease of the mineral matter, implying a lower filtration efficiency. Application of this process needs to take into consideration limitations due to the filtration step since higher IF breakdown leads to a weaker dynamic plug. Images of extrudates and filtrates (Figs. S 1 and S 2) depict the impact of the filtration efficiency in extrudate moisture content and solid content of the filtrate. It must be mentioned that filtration is critical in order to ensure the elimination of salts that inhibit fermentation step [START_REF] Casey | Effect of salts on the Co-fermentation of glucose and xylose by a genetically engineered strain of Saccharomyces cerevisiae[END_REF][START_REF] Deparis | Engineering tolerance to industrially relevant stress factors in yeast cell factories[END_REF], and to obtain high extrudate dry matter to easily manage L/S ratios in further phases.

Enzymatic hydrolysis

The enzyme efficiency to access polysaccharides for saccharification was mainly evaluated through glucose (GR) and xylose released (XR). These equations express the ability of enzymes to access and produce glucose and xylose from the extrudate. Simultaneously, the process was assessed by the analysis of the glucose and xylose yields. The results from enzymatic hydrolysis after 48 h were summarized in Table 5. Both sugar released and yields were studied statistically. Analysis of the XR was conducted by analysis of coefficients and the iso response curves (Table 6 and Fig. 4), however GR was studied through the values in Table 5 due to poor fitting with the model. General observation showed that the combination of twin-screw extrusion and thermo-chemical action upgraded the action of enzymes on SCC as GC and XC values of extrudates vary from 64% to 81% and 46%-81.7% which is an improvement of by 2-3 times compared to untreated biomass (Table 1). This is associated to the mechanical and chemical actions. The conveying and kneading screws have proven to increase porosity of biomass and reduce of particle size, promoting the enzyme action [START_REF] Zheng | The effects of screw elements on enzymatic digestibility of corncobs after pretreatment in a twin-screw extruder[END_REF]. Using alkaline extrusion, pretreated biomass exposes higher amount of non-crystalline cellulose at the surface, producing higher conversion of carbohydrates than sole alkali pretreated biomass [START_REF] Khatri | Determination of optimal biomass pretreatment strategies for biofuel production: investigation of relationships between surface-exposed polysaccharides and their enzymatic conversion using carbohydrate-binding modules[END_REF].

Increasing the NaOH load clearly leads to a modification of the structure of the lignocellulosic fraction in the extruded matter, primarily by the solvation of hemicelluloses and lignin. Since both components act as a physical barrier enveloping cellulose, its elimination from the substrate can lead to a higher production of glucose by EH. Mild pretreatment conditions lead to less destructuration of lignocellulose, thus lower GC values (Samples 1, 2, 10 and 11). Nonetheless, higher GR can be reached by increasing solubilization of hemicelluloses, either by NaOH increment or by temperature augmentation (Samples 3-5).

Concerning XR, according to the coefficients, the temperature and the NaOH/SCC impact its modification on the conditions studied (Table 6). The linear coefficients a 1 and a 2 are positive, meaning beneficial action of these parameters. However, quadratic coefficient for temperature is negative, putting into evidence a negative influence. As observed in the response surface, the best XR is reached with a combination of NaOH loading beyond 7% and temperature rise ˜125 °C.

Although, for all soda concentration, an increase in temperature has an adverse effect. These observations may be attributed to the lignin restructuration. As previously discussed, lignin is profoundly affected by temperature and can undergo structural modifications. In the previous study by [START_REF] Zheng | The effects of screw elements on enzymatic digestibility of corncobs after pretreatment in a twin-screw extruder[END_REF], the action of reverse screw induces high local temperatures that reached transition temperature of lignin. Although it was fixed at 75 °C, the outcome was lignin blocking pores and lower enzyme digestibility. Lignin content is not modified (around 4%) but could undergo structural changes due to temperature and pressure applied by the reverse screw placed after the filtration module (beginning of module 7), provoking the inactivation of enzymatic digestion.

Promoting enzyme efficiency to degrade carbohydrates through optimization of the pretreatment conditions is not the sole objective of this study. High production of fermentable sugar and the deduction of chemicals applied during the pretreatment is crucial to achieving a techno-economical viable process. Consequently, the evaluation of glucose (GY) and xylose yields (XY) is essential to determine the performance of the processes in the studied conditions. Fit for both responses did not reach significant values, in contrast, tendencies can be deduced from statistical analysis (Table 6). Iso-response curves of XY (Fig. 4(b)) reveal that for all alkali charges, temperatures between 90 and 125 °C allow the increment of XY beyond 60%. However, temperature and soda can also have adverse effects, for temperature above 150 °C and soda loading higher than 7%. Bringing into comparison the iso-response curves of XR and XY, it shows that XY is not correlated only with enzymatic accessibility. Although high XC is reached with augmented hemicelluloses solubilization, the XY is affected when the hemicelluloses yield decreased by the effect of soda and temperature.

Concerning GY (Fig. 4(c)), temperature boost GY between 6-7.5% NaOH/SCC, this effect is also limited for temperatures above 125 °C. This condition is consistent with high values of CY (˜90%) and GC (> 70%). These results illustrate the necessity to associate not only high enzymatic digestibility but significant cellulose conservation to achieve high GY. In addition, GY is negatively impacted by temperatures higher than 150 °C and NaOH/SCC ratios beyond 8%, likely due to cellulose deterioration as mentioned earlier. Furthermore, analysis of both iso-response curves allows to conclude that the conditions for optimal glucose and xylose yields are at 6-6.5% NaOH/dry SCC and internal temperature from 90 °C to 125 °C, where GY and XY could reach 25% each (500 g of sugars per 1 Kg of dry SCC).

Conclusions

In the present study, the effects of pretreatment alkali loading and temperature in the extruder were investigated. The results underline that several phenomena may be involved in biomass structural modification and the enzymatic hydrolysis. Extrudate insoluble fiber, cellulose and hemicelluloses are mainly affected by the employment of catalyst, and temperature enhances its effect. Under the test conditions, the higher NaOH charge and temperature resulted in the disintegration of hemicelluloses, producing a cellulose-rich substrate. Likewise, it improves the enzymes accessibility resulting from physical barrier elimination. Despite the positive effect of temperature and NaOH/SCC on enzymatic hydrolysis, there are limitations since an adverse effect is observed for cellulose conservation, xylose conversion, and sugar yields. Overall, sugar yield can be improved while decreasing the chemical content in the process through heat.

In terms of process, since there are not competitive uses for this material, the valorization of SCC into bioethanol is economically Table 6 Statistical analysis of the effect of studied parameters on enzymatic hydrolysis.

None.

interesting for the local market. Considering future steps (impregnation of enzymes, saccharification, and fermentation), the optimal conditions set by this work were at 6-7% NaOH/dry SCC with an internal temperature between 110-130 °C. Compared with conditions employed in earlier works, this study provides a decrease of 3% of sodium hydroxide and 2% of acid inputs. This contributes to the reduction of the pretreatment cost and its environmental impact. Advantageously, in these conditions, filtration is stable and solid content of extrudate is around 36%. This promotes working at high consistency for the future steps and the elimination of inhibitory molecules. The future work will focus on improved enzymatic impregnation for saccharification and co-fermentation at high consistency with additional feasibility studies at pilot scale.

Declarations of interest

Supplemental files

  the mass of carbohydrate in grams, the cellulose (C) and hemicelluloses (H) in extrudate correspond to each carbohydrate contained in 0.5 g of dry extrudate. The 1.1 value corresponds to the conversion factor of glucan into glucose (180/162) and hemicelluloses into xylose and arabinose (150/132). SCC hemicelluloses are composed of 85% of xylan and 15% arabinose (data not shown). Therefore, in this study, XR only accounts as an indicator of the xylose produced.

Fig. 1 .

 1 Fig. 1. Schematic diagram of the modular barrel of EV25 twin-screw extruders; N.C.: Non-controlled; TP: Temperature probe (to be viewed in color).

Fig. 2 .

 2 Fig. 2. Iso response curves indicating Insoluble Fiber (%DM) (a), hemicelluloses fraction (%IF) (b), cellulose fraction (%) (c) and hemicelluloses removal rate (%) (d) as function of the temperature (T) and NaOH/SCC ratio. The circles represent the experimental results. (To be viewed in color).

Fig. 3 .

 3 Fig. 3. Iso response curves indicating solid dry matter (DM) of extrudate (a) and mineral matter (MM) of the filtrate (b) as function of the temperature (T) and the NaOH/SCC ratio, the circles represent the experimental results. (To be viewed in color).

Fig. 4 .

 4 Fig. 4. Iso response curves indicating xylose released rate (%) (% of xylose released/g of hemicelluloses in extrudate)(a), xylose yield (XY)(% g of xylose released/g of dry SCC)(b) and glucose yield (GY)(% g of glucose released/g of dry SCC)(c). The circles represent the experimental results. (To be viewed in color).

Figure S 1 :

 1 Figure S 1: Photos of untreated biomass and extrudates obtained with part conditions evaluated in the present work

Table 1

 1 Raw sweet corn characteristics: solid content, composition and enzymatic hydrolysis results.

		Units	Values		S.D.
	Composition
	Dry matter	%	89.5	±	0.1
	Minerals	% DM	5.9	±	0.3
	Protein	% DM	8.5	±	0.22
	Insoluble fiber	%DM	75.2	±	0.5
	Cellulose	%IF (%DM)	40.8 (31.0)	±	0.8
	Hemicelluloses	%IF (%DM)	53.1 (39.0)	±
	Lignin	%IF (%DM)	6.1 (5.0)	±
	Enzymatic hydrolysis at 48 h
	GR	% (g/g cellulose in SCC)	28.7	±	1.7
	XR	% (g/g of hemicelluloses in SCC)	25.2	±	1.6
	GY	% (g/g of dry SCC)	9.8	±	0.6
	XY	% (g/g of dry SCC)	11.6	±	0.7

Note: DM: Dry Matter IF: Insoluble Fiber EH: Enzymatic Hydrolysis; GR: Glucose Released; XR: Xylose released; GY: Glucose Yield; XY: Xylose Yield.

  on barley straw, pinpointed that at the highest NaOH loading (i.e. 7.5% NaOH/g

	Sample	NaOH/SCC	Temperature	Extrudate composition				Mass balance indicators		Filtrate
				Solid content	IF	C	H	L	Q ext	CY	HR	MM
	Units	% (w/w)	°C	%DM	%DM	%IF (%DM)	%IF (%DM)	%IF (%DM)	Kg/h	%	%	%DM
	1	5.9	52	45.0	79.8	43.1 (34.4)	54.3 (40.5)	2.6 (4.9)	1.33	83.9	16.8	32.8
	2	5.9	92	40.0	79.0	43.4 (34.2)	53.9 (40.3)	2.8 (4.4)	1.34	83.3	16.8	33.0
	3	5.9	108	39.0	78.0	45.2 (35.3)	51.2 (37.9)	3.5 (4.8)	1.31	87.3	10.2	32.7
	4	5.9	124	38.0	76.0	46.4 (35.5)	50.0 (36.7)	3.5 (4.3)	1.27	87.4	11.2	33.0
	5	6.1	118	38.2	78.5	48.1 (37.7)	47.7 (37.0)	4.2 (3.7)	1.27	89.4	17.3	31.1
	6	6.0	129	38.0	76.5	49.7 (38.0)	45.3 (34.7)	4.9 (3.2)	1.27	87.5	22.8	31.2
	7	6.1	151	36.5	75.7	50.7 (35.3)	45.0 (34.3)	4.3 (3.6)	1.19	94.8	23.2	33.6
	8	6.1	168	36.3	73.3	48.2 (38.4)	47.3 (34.1)	4.4 (3.2)	1.22	95.4	28.9	33.0
	9	7.1	130	36.4	74.0	50.5 (37.4)	46.3 (34.3)	3.2 (2.4)	1.20	90.0	34.7	36.4
	10	3.9	95	45.1	81.6	43.1 (35.2)	52.2 (42.6)	4.7 (3.8)	1.21	85.2	31.6	31.5
	11	4.2	146	42.4	77.7	43.0 (33.5)	52.6 (40.9)	4.3 (3.4)	1.27	88.6	33.6	31.9
	12	8.0	92	38.3	72.8	49.4 (36.0)	45.7 (33.3)	4.9 (3.5)	1.22	86.4	34.7	39.2
	13	8.0	121	33.8	68.1	51.8 (35.2)	44.2 (30.1)	4.0 (2.7)	1.18	81.9	42.8	37.8
	14	8.0	139	32.9	67.8	50.7 (34.4)	44.9 (33.6)	4.3 (2.9)	1.14	77.0	44.2	37.0

Note: SCC: Sweet Corn Co-product; IF: Insoluble Fiber; DM: Dry Matter; Qext: Extrutate flow; CY: Cellulose Yield; HR: Hemicelluloses Removal; MM: Mineral Matter.

Table 3

 3 Statistical analysis of the effects of studied parameters on extrudate composition and mass balance indicators.

		Insoluble Fiber		Hemicelluloses		Cellulose		HR
	R 2	0.952		0.786		0.84		0.920
	Terms	Estimated	p value	Estimated	p value	Estimated	p value	Estimated	p value
	a 0	78.04	< 0.01 ***	49.50	< 0.01 ***	46.81	< 0.01 ***	21.175	< 0.01***
	a 1	-4.96	< 0.01 ***	-3.65	0.740 **	3.75	0.157 **	11.193	0.00477**
	a 11	-2.63	0.514 **	-0.51	70.7	-0.22	83.1	6.987	2.16*
	a 2	-3.29	0.0709 ***	-3.75	1.23 *	2.99	1.10 *	12.517	0.00541***
	a 22	-1.31	21.2	1.21	52.6	-1.02	49.2	-0.220	2.40*
	a 12	-0.89	52.8	-1.07	68.6	1.50	47.0	6.170	24.2

Table 4

 4 Statistical analysis of the effect of the normalized parameters (temperature and NaOH/SCC) on extrudate (a) and filtrate mineral matter (b).

		Extrudate Dry matter	Filtrate mineral matter
	R 2	0.981		0.900
	Terms	Estimated	p value	Estimated	p value
	a 0	38.76	< 0.01 ***	32.55	< 0.01***
	a 1	-4.06	< 0.01 ***	3.46	0.0128***
	a 11	1.01	2.83 *	2.31	0.686**
	a 2	-4.11	< 0.01 ***	-0.33	70.5
	a 22	1.94	0.646 **	0.55	56.0
	a 12	-2.07	2.36 *	-1.32	32

Table 5

 5 Enzymatic hydrolysis indicators. Results are based on glucose and xylose released after 48 h, and express the enzymatic performances to digest extrudate (GR and XR) and yield of the overall process (GY and XY).

	Sample NaOH/SCC Temperature Enzymatic hydrolysis
				GR	XR	GY	XY
	units	% (w/w)	°C	% g/g	% g/g of	% g/	% g/
				cellulose	hemicelluloses	g of	g of
				in	in extrudate	dry	dry
				extrudate		SCC	SCC
	1	5.9	52	67.7	50.9	23.6 19.4
	2	5.9	92	68.6	64.7	23.8 22.5
	3	5.9	108	70.3	73.0	24.7 26.7
	4	5.9	124	76.0	81.1	25.6 25.8
	5	6.1	118	71.3	82.3	25.0 27.6
	6	6.0	129	68.3	78.0	25.2 24.5
	7	6.1	151	72.4	73.0	23.8 21.0
	8	6.1	168	64.4	69.0	23.0 20.4
	9	7.1	130	81.0	79.1	27.4 22.8
	10	3.9	95	66.8	59.9	21.6 21.8
	11	4.2	146	67.5	66.6	23.0 21.9
	12	8.0	92	74.9	88.3	24.8 25.5
	13	8.0	121	70.9	89.3	22.0 22.4
	14	8.0	139	75.4	82.0	22.2 22.2
	Note: SCC: Sweet Corn Coproduct; GR: Glucose Released; XR: Xylose Released;
	GY: Glucose Yield; XY: Xylose Yield.

Table S 1: Normalized data using Z-score method

 S

			Experimental design
	Sample	NaOH/SCC	T	X1	X2
	units	% (w/w)	°C
	1	5.9%	52	-0.0500	-1.0000
	2	5.9%	92	-0.0500	-0.3103
	3	5.9%	108	-0.0500	-0.0345
	4	5.9%	124	-0.0500	0.2414
	5	6.1%	118	0.0500	0.1379
	6	6.0%	129	0.0000	0.3276
	7	6.1%	151	0.0500	0.7069
	8	6.1%	168	0.0500	1.0000
	9	7.1%	130	0.5500	0.3448
	10	3.9%	95	-1.0500	-0.2586
	11	4.2%	146	-0.9000	0.6207
	12	8.0%	92	1.0000	-0.3103
	13	8.0%	121	1.0000	0.1897
	14	8.0%	139	1.0000	0.5000

SCC: Sweet corn co-product; T: temperature probe 1

Acknowledgements

This work has been funded by the European Community (Horizon 2020) under Grant agreement reference n•654365 (BABET-REAL 5 PROJECT)). The authors would like to thank Dr. Franck Cerruti for the statistical analysis assistance.

Appendix A. Supplementary data

Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.indcrop.2019.03.044.


Aspect ratio effects in Rayleigh-Bénard convection of Herschel-Bulkley fluids

Keywords: Rayleigh-Bénard, Herschel-Bulkley, Aspect ratio, Finite element Research paper

Purpose

-The purpose of this paper is to analyze two dimensional steady state Rayleigh-Bénard convection within rectangular enclosures in different aspect ratios filled with yield stress fluids obeying the Herschel-Bulkley model. Design/methodology/approach -In this study, a numerical method based on the finite element has been developed for analyzing two dimensional natural convection of a Herschel-Bulkley fluid. The effects of Bingham number ݊ܤ and power law index ݊ on heat and momentum transport have been investigated for a nominal Rayleigh number range (5 * 10 ଷ < ܴܽ < 10 ହ ) , three different aspect ratios (ratio of enclosure length:height ܴܣ = 1, 2, 3 ) and a single representative value of nominal Prandtl number ݎܲ( = 10). Findings -Results show that the mean Nusselt number ݑܰ തതതത increases with increasing Rayleigh number due to strengthening of convective transport. However, with the same nominal value of ܴܽ the values of ݑܰ തതതത for shear thinning fluids ݊ < 1 are greater than shear thickening fluids ݊ > 1 . The values of ݑܰ തതതത decrease with Bingham number and for large values of ݊ܤ , ݑܰ തതതത rapidly approaches to unity which is indicated that heat transfer takes place principally by thermal conduction. The effects of aspect ratios have also been investigated and results show that ݑܰ തതതത increases with increasing ܴܣ due to stronger convection effects. Originality/value -Finite element analysis of viscoplastic fluids. Analysis of Rayleigh-Bernard flows involving Herschel-Bulkley fluids for a wide range of Rayleigh numbers, Bingham numbers and power law index. Study on the effects of aspect ratio on flow and heat transfer of Herschel-Bulkley fluids.

Introduction

Viscoplastic fluids are characterized by a yield stress, the existence of a "yield stress" is traditionally recognized to be responsible for the complicated transition between classical solid-like and liquid-like behavior. If the material is not sufficiently stressed, it will not deform and behaves like a solid, but once the yield stress is exceeded, it will deform and flow according to different constitutive relations.

Several systems and industrial processes are based on viscoplastic fluid behavior such as blood, drilling mud, mayonnaise, toothpaste, grease, and some lubricants. There are three ideal models based on viscoplastic model, the Bingham plastic, the Herschel-Bulkley model, and the Casson model [START_REF] Mitsoulis | Flows of viscoplastic materials: models and computations[END_REF]. Heat transfer analysis and especially convection had been one of the interesting parts of viscoplastic research works in last decades because of its applications, such as Nouar studies about combined forced and free convection heat transfer of a yield stress fluid [START_REF] Nouar | Thermal convection for a thermo-dependent yield stress fluid in an axisymmetric horizontal duct[END_REF] and thermal convection for an incompressible Herschel-Bulkley fluid along an annular duct [START_REF] Nouar | Numerical and experimental investigation of thermal convection for a thermodependent Herschel-Bulkley fluid in an annular duct with rotating inner cylinder[END_REF]. In this field instability of necking processes for an elasto-viscoplastic lithosphere has been studied by Regenauer-Lieb and his co-worker [START_REF] Regenauer-Lieb | Quasi-adiabatic instabilities associated with necking processes of an elasto-viscoplastic lithosphere[END_REF] and also, Soares et al. [START_REF] Soares | Heat transfer to viscoplastic materials flowing axially through concentric annuli[END_REF] have analyzed heat transfer in the entrance-region flow of Herschel-Bulkley fluids inside concentric annular spaces. Despite the number of researches carried out, information about natural convection process within a square enclosure containing viscoplastic material (especially for Herschel-Bulkley and Casson models) is still very scarce.

Rayleigh-Bénard convection as an important part of heat transfer has been investigated by some researchers. Actually Rayleigh-Bénard is a type of natural convection and occur in a fluid layer due to a horizontal temperature gradient, when the lower wall is heated and the upper wall cooled. First studies on this problem were motivated by the Bénard's experiments around 1900 [START_REF] Bénard | Les tourbillons cellulaires dans une nappe liquide[END_REF]) who considered the stability of a fluid layer heated from below. A linear stability analysis was proposed in first chapters of Chandrasekhar's book [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF] present the linear theory within the Boussinesq approximation. Non-linear approaches were reviewed in [START_REF] Busse | Non-linear properties of thermal convection[END_REF][START_REF] Newell | Order Parameter Equations for Patterns[END_REF]).

In the case of Newtonian and power law fluids, Aghighi et al. [START_REF] Aghighi | Non-incremental transient solution of the Rayleigh-Bénard convection model by using the PGD[END_REF]) solved transient and steady state Rayleigh-Bénard convection using the Proper Generalized Decomposition (PGD) method. They also analyzed the effects of nanoparticle's temperature and diameter on Rayleigh-Bénard convection of Nano-fluids [START_REF] Aghighi | Parametric solution of the Rayleigh-Benard convection model by using the PGD Application to nanofluids[END_REF]. Before that, the RB convection along the isochore of nitrogen has been numerically studied by Shen and his co-worker [START_REF] Shen | Rayleigh-Bénard convection in a supercritical fluid along its critical isochore in a shallow cavity[END_REF] using the SIMPLE algorithm. The other numerical study has been done by Kao and Yang [START_REF] Kao | Simulating oscillatory flows in Rayleigh-Bénard convection using the lattice Boltzmann method[END_REF]. They used the lattice Boltzmann method to simulate the oscillatory flows of the secondary instability in 2D Rayleigh-Bénard convection. The steady state condition of twodimensional Rayleigh-Bénard Convection was analyzed by [START_REF] Ouertatani | Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure[END_REF] using finite volume method. There are also some experimental works such as Maystrenko et al. [START_REF] Maystrenko | Structure of the thermal boundary layer for turbulent Rayleigh-Bénard convection of air in a long rectangular enclosure[END_REF]) research about measurements of the temperature distribution in the upper (cold) boundary layer of a rectangular Rayleigh-Bénard cell for air in large aspect ratios.

While for Newtonian and Power law fluids we can refer several research studies, in the case of viscoplastic material there are no many research works especially for Herschel-Bulkley and Casson models. In the case of Bingham fluids the effects of a fluid yield stress on the classical Rayleigh-Bénard instability were examined by Zhang and his co-workers [START_REF] Zhang | Yield stress effects on Rayleigh-Bénard convection[END_REF]) and Vikhansky [START_REF] Vikhansky | Thermal convection of a viscoplastic liquid with high Rayleigh and Bingham numbers[END_REF]) considered the effect of yield stress on the Rayleigh-Bénard convection of a viscoplastic material. A numerical analysis based on Fluent for Rayleigh-Bénard convection and Bingham model has been done by [START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF]. Recently, the effects of aspect ratio on steady state Rayleigh-Bénard convection of Bingham fluids have been analyzed by [START_REF] Yigit | Effects of aspect ratio on natural convection of Bingham fluids in rectangular enclosures with differentially heated horizontal walls heated from below[END_REF] using Fluent simulation. Finally in the field of experimental study, Darbouli and his co-workers [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF]) have investigated the influence of rheological and interfacial properties of yield stress fluids on the onset of the Rayleigh-Bénard convection. Besides developing a powerful numerical method based on the finite element which can be used for analyzing different types of viscoplastic fluids, this study also investigated the effects of aspect ratio ܴܣ( = 1, 2, 3) on steady-state natural convection of yield stress fluids obeying the Herschel-Bulkley model within rectangular enclosures over wide ranges of conditions as: Rayleigh number, 5 * 10 ଷ < ܴܽ < 10 ହ , Bingham number, 0 < ݊ܤ < ݊ܤ ௫ and power law index, 0.7 < ݊ < ݊ ௫ .

To the best of our knowledge, it has not been studied before and the results reported here are new.

Mathematical formulation

The case under consideration is the two dimensional steady natural convection of a Herschel-Bulkley fluid, in a rectangular cavity (Figure 1). No-slip conditions are considered at walls. The cavity is heated from below and cooled from above. It is assumed that the viscous dissipation terms are negligible and fluid has constant properties (except the density changes which produce buoyancy forces).

Assuming the Boussinesq approximation and using the characteristic scales ܪ for length, ݑ = )ܶ߂ܪߚ݃( ଵ ଶ ൗ for the velocity, ݐ = ܪ ݑ ൗ for the time, and  = ݑߩ ଶ for the pressure, the dimensionless continuity, momentum and energy equations for Rayleigh-Bénard convection could be written as:

ݑ߲ ݔ߲ + ݒ߲ ݕ߲ = 0 ݑ ݑ߲ ݔ߲ + ݒ ݑ߲ ݕ߲ = - ߲ ݔ߲ + ݎܲ ଶ ܴܽ ିଶ ଶ ቆ ߲߬ ௫௫ ݔ߲ + ߲߬ ௬௫ ݕ߲ ቇ ݑ ݒ߲ ݔ߲ + ݒ ݒ߲ ݕ߲ = - ߲ ݕ߲ + ݎܲ ଶ ܴܽ ିଶ ଶ ቆ ߲߬ ௬௬ ݕ߲ + ߲߬ ௫௬ ݔ߲ ቇ + ߠ ݑ ߲ߠ ݔ߲ + ݒ ߲ߠ ݕ߲ = (ܴܽ. )ݎܲ ି ଵ ଶ ൗ ቆ ߲ ଶ ߠ ݔ߲ ଶ + ߲ ଶ ߠ ݕ߲ ଶ ቇ (1)
The steady state boundary conditions for this problem would be:

ݑ = ݒ = 0 at all walls ߲ߠ ݔ߲ = 0 ݐܽ ݔ = 0 ܽ݊݀ ݔ = ܮ ܪ = 1, 2, 3 ߠ = 0.5 ݐܽ ݕ = 0 ߠ = -0.5 ݐܽ ݕ = 1 (2)
The Prandtl number, ,ݎܲ and Rayleigh number, ܴܽ, defined by:

ݎܲ = ܥߤ ݇ ܴܽ = ܪܶ∆ߚ݃ ଷ ߙߴ (3)
The various dimensional quantities above are defined as follows: ,ݑ ,ݒ ߠ and  are dimensionless horizontal velocity, vertical velocity, temperature and pressure respectively. ߤ is the dynamic viscosity, ܥ is specific heat capacity, ݇ is the thermal conductivity, ݃ is the acceleration due to gravity, ߚ is the coefficient of thermal expansion, α is the thermal diffusivity, ∆ܶ is the temperature difference between hot and cool walls, ߥ is the kinematic viscosity, ߬ is the viscous stress and ݊ is power law index.

The dimensionless temperature ߠ is defined by:

ߠ = ܶ -ܶ ܶ ு -ܶ (4)
Where ܶ and ܶ ு ( ܶ ு > ܶ ) are the temperatures enforced at the top and bottom cavity boundaries, respectively, and ܶ is a reference temperature:

ܶ = ( ܶ ு + ܶ ) 2 ⁄ .
The heat flux averaged over the hot wall is defined via the Nusselt number:

ݑܰ തതതത = -න ߲ߠ ݕ߲ ฬ ௬ୀ ݔ݀ (5)
ଵ

The stress-deformation behavior of viscoplastic materials based on Herschel-Bulkley model is given by:

߬ = ൬ߛሶ ିଵ + ݊ܤ ߛሶ ൰ ߛሶ ݎ݂ |߬| > ߬ ௬ ܽ݊݀ ߛ ሶ = 0 ݎ݂ |߬| < ߬ ௬ (6)
For ݅, ݆ = 1,2 with ( ݔ ଵ , ݔ ଶ ) = ( ,ݔ ݕ ), where ݊ is power law index, ݊ܤ Bingham number,

ߛሶ = ට ଵ ଶ ߛሶ ߛሶ and ߬ = ට ଵ ଶ ߬ ߬ .
The rate-of-strain tensor ߛሶ is defined by:

ߛሶ = ݑ߲ ݔ߲ + ݑ߲ ݔ߲ (7)
The dimensionless Bingham number ݊ܤ generally represents the ratio of yield stress to buoyancy stresses but based on equations ( 1) and ( 6) the value of ( .ݎܲ ܴܽ )

ି ଶ

ൗ should be added to it so:

݊ܤ = ( .ݎܲ ܴܽ ) ି ଶ ൗ ߬ ௬ ܪܶ∆݃ߚߩ (8)
Papanastasiou modifications method [START_REF] Papanastasiou | Flows of Materials with Yield[END_REF]) can be applied to Herschel-Bulkley model to avoid discontinuity between yielded and unyielded regions. He proposed an exponential regularization of stress equation. based on this model equation ( 6) would be rewritten as:

߬ = ൬ߛሶ ିଵ + 1(݊ܤ -exp (-݉ߛሶ ߛሶ ൰ γሶ (9)
where ݉ is regularization parameter and controls the exponential rise in the stress at low rates of strain.

While this model has been used by many researchers to analyze viscoplastic problems there are also other regularization techniques which circumvent the discontinuity inherent in these materials. The one which has gained wide acceptance in the literature is bi-viscous model. In this model, the unyielded material is assigned a very high value of viscosity. Comparison between the predictions of the yielded/unyielded regions using the two regularization techniques shows that the bi-viscous model is unable to identify the small unyielded regions as captured by the Papanastasiou regularization technique [START_REF] Sairamu | Natural convection from a circular cylinder in confined Bingham plastic fluids[END_REF]). However, these minor differences in the identification of yielded and unyielded regions by the two techniques have very little influence on the value of the Nusselt number [START_REF] Sairamu | Natural convection from a circular cylinder in confined Bingham plastic fluids[END_REF][START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF].

Numerical analysis

A numerical code based on the finite element method with quadrilateral, eight nodes elements was developed to solve the coupled conservation equations of mass, momentum and energy related to the two dimensional steady state Raleigh-Bénard convection within rectangular enclosures. When we proceed to the discretization of the weak form related to equation ( 1) some stability conditions must be ensured. One of them concerns the so-called LBB condition that restricts the free choice of pressure and velocity approximations. As in Aghighi et al. [START_REF] Aghighi | Non-incremental transient solution of the Rayleigh-Bénard convection model by using the PGD[END_REF] in what follows we consider a penalty formulation of the incompressibility constraint. Thus, the equation ( 1) can be rewritten as below:

ݑ߲ ݔ߲ + ݒ߲ ݕ߲ + 1 ߣ  = 0 ݑ ݑ߲ ݔ߲ + ݒ ݑ߲ ݕ߲ + ߲ ݔ߲ = ݎܲ ଶ ܴܽ ିଶ ଶ ቆ ߲߬ ௫௫ ݔ߲ + ߲߬ ௬௫ ݕ߲ ቇ ݑ ݒ߲ ݔ߲ + ݒ ݒ߲ ݕ߲ + ߲ ݕ߲ -ߠ = ݎܲ ଶ ܴܽ ିଶ ଶ ቆ ߲߬ ௬௬ ݕ߲ + ߲߬ ௫௬ ݔ߲ ቇ ݑ ߲ߠ ݔ߲ + ݒ ߲ߠ ݕ߲ -(ܴܽ. )ݎܲ ି ଵ ଶ ൗ ቆ ߲ ଶ ߠ ݔ߲ ଶ + ߲ ଶ ߠ ݕ߲ ଶ ቇ = 0 (10)
Where λ is a large enough constant. In this set of equations all variables are on the left side, but the stress terms on the right. The values of stress assumed to be known when the momentum equations analyzed. Initial values of velocity were considered to calculate shear rate and stress terms.

Then, using the obtained values of stress and solving equations ( 10) the new values of velocity were computed. These new values were used to recompute stress terms and the process repeated until the values of velocity and temperature converge. It should be indicated that the numerical analysis of equations ( 10) in this format is not stable so the Laplace values of velocity were added to both sides of momentum equations these values does not affect the results but make the equation stable. The CPU time shows that this convergence is very fast with an average of 400 second (Running on Mac 2.9

GHz Intel Core i5 and after about 1500 iteration). Using this method Rayleigh-Bénard convection has been analyzed successfully for viscoplastic fluids based on Herschel-Bulkley model in an enclosure with three different aspect ratios and results are presented below.

Results and discussion

Numerical method validation

Mesh convergence studies for each aspect ratio separately and based on the results, the meshes consist of 1825, 2133 and 2465 nodes have been used for AR=1, 2, 3 respectively.

The convergence of the solutions was checked by varying the penalty and regularization parameters.

Results show that the ݑܰ തതതത value converges within 0.1% by varying ߣ from 10 ଷ to 10 ସ and by varying ݉ from 10 ଶ to 10 ଷ . In the following, all results are obtained for ߣ = 10 ସ and ݉ = 10 ଷ .

To validate the numerical method Rayleigh-Bénard convection for Newtonian (݊ = 1, ݊ܤ = 0, ܴܣ = 1) and Bingham (݊ = 1, ܴܣ = 1) fluids was solved and the results compared with previous research works.

In the case of Newtonian fluids the results are compared with [START_REF] Ouertatani | Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure[END_REF]

Velocity and temperature

Due to the large number of results and to summarize, the results of velocity and temperature are only presented for mid-range Rayleigh number (ܴܽ = 5 * 10 ସ ).

One can see the effects of power law index in figure 3 where the variations of non-dimensional demonstrate the advantages of the present numerical method for solving the nonlinear complex models in critical conditions (critical power law index).

Comparing the results of figures 3-5 show that the magnitude of velocity increases with increasing ܴܣ due to stronger convection force there are also more maximum (and minimum) points of velocity because of more convection rolls in the enclosure where the numbers of convection rolls are equal to the value of aspect ratio. It should be noted that the vertical mid-plane of the enclosure passes through the center of convection roll when ܴܣ = 1 and 3, but from the boundary of them when ܴܣ = 2 so different distribution of the temperature can be seen.

One can see that the magnitude of velocity and non-linearity of temperature distribution decrease with increasing power law index (or Bingham number) indicating weakening of convective thermal transport due to the additional flow resistance and this trend strengthens with increasing power law index (Bingham number) cause flow tend toward the rest condition and temperature towards the linear variation with ܮ/ݔ which is indicated that the relative contribution of thermal conduction to the overall heat transfer increases with increasing power law index (Bingham number).

For high value of power law index (Bingham number), the viscous forces overcome the buoyancy force and as a result of this, the heat transfer is due to thermal conduction and there is no significant flow within the enclosure. One can find this result from figures 3-5 where the temperature distribution becomes linear and the magnitude of velocity converge to zero. For linear distribution of temperature, the mean Nusselt number approach to unity which means that the heat transfer is due to thermal conduction. The value of power law index or Bingham number at which ݑܰ തതതത approaches to ݑܰ തതതത = 1 is defined as critical power law index (݊ ௫ ) or critical Bingham number ( ݊ܤ ௫ ). Comparing the results from figures 3-5 show that the value of ݊ ௫ or ݊ܤ ௫ increases with increasing ܴܣ due to the stronger convection force which can overcome the flow resistance up to greater values of power law index (Bingham number). More details about heat transfer and mean Nusselt number are presented in the next section.

Heat transfer (Nusselt number)

The results are presented in two parts. First, the effects of power-law index ݊ are shown for three and ݊ܤ = 0.25. The black lines in these figures show the plug regions (TUR and AUR) where ߬ < ߬ ௬ . The truly unyielded regions (TUR) move with a plug velocity profile (no deformation), and the apparently unyielded regions (AUR) are in the corners, where the velocities are very small, so the area behaves as stagnant, and no deformation occurs [START_REF] Mitsoulis | Flows of viscoplastic materials: models and computations[END_REF]. One can see that the isothermal lines become more curved and plug regions decrease with decreasing ݊ because of the low viscosity effects. Furthermore, the number of convection rolls within the enclosure change with ܴܣ and thus the qualitative nature of isothermal line and velocity vectors have also changed with .ܴܣ It can be observed that increase in ܴܣ generally decrease plug regions due to stronger convection effects.

Conclusions

In this study a finite element numerical code has been developed for analyzing two dimensional steady state natural convection of viscoplastic fluids obeying the Herschel-Bulkley model within rectangular enclosures heated from below and cooled from above at constant temperature.

The effects of Bingham number ݊ܤ , power law index ݊ and aspect ratio ܴܣ

  for several Rayleigh numbers and in all cases very good agreement can be determined. Results of Bingham fluids are compared with (Turan et al. 2012) where the Fluent base solution of Rayleigh-Bénard convection for Bingham fluids could be found. Results of non-dimensional velocity ݑ , non-dimensional temperature ߠ and mean Nusselt number ݑܰ തതതത for ܴܽ = 5 * 10 ସ and ݎ = 10 are shown in figure 2 where the line shows the results of the present study and points show the Turan results.

velocity

  ݒ and temperature ߠ along the horizontal and vertical mid-planes of the cavity are shown for different values of Power law index ݊ at ܴܣ = 1, 2, 3 and ݊ܤ = 0.5. On the other hand, the effects of Bingham number are presented in figures 4 and 5 for shear thinning (݊ = 0.85) and shear thickening (݊ = 1.05) viscoplastic fluids respectively. It should be noted that while the power law index greater than one may not present real fluid behavior but the mathematical results are still considered to

selected

  Figures 10 and 11 show contours of non-dimensional temperature ߠ and velocity vectors for three different values of aspect ratio ܴܣ( = 1,2,3) and power law index (݊ = 0.85,1,1.05) at ܴܽ = 5 * 10 ସ

  (length: height) on heat and momentum transport have been investigated in the nominal Rayleigh number range ( 5 * 10 ଷ < ܴܽ < 10 ହ ) for a single representative value of nominal Prandtl number ݎܲ( = 10). Results show that the mean Nusselt number ݑܰ തതതത increases with increasing Rayleigh number due to stronger convection effects and decreases with increasing power law index because of strengthening viscous forces. For a given set of values of the Rayleigh number ܴܽ and power law index ݊ the mean Nusselt number ݑܰ തതതത is found to decreases with increasing values of the Bingham number. Results show that for a large value of power law index (݊ ≥ ݊ ௫ ) or Bingham number ݊ܤ( ≥ ݊ܤ ௫ ) the value of mean Nusselt number rapidly approaches unity ݑܰ( തതതത = 1) which means that the dominant mode of heat transfer in the enclosure is thermal conduction. The values of ݊ ௫ and ݊ܤ ௫ depends on ܴܽ , and increase with increasing ܴܽ due to stronger convection effects. The effects of aspect ratio on momentum and thermal transport have also been investigated in detail. The number of convection rolls and mean Nusselt numbers are found to increase with increasing ܴܣ for a given set of values of ܴܽ and power law index ݊. Results also show that for shear thinning fluids (݊ < 1) the effects of Bingham number on heat transfer decrease with decreasing Rayleigh number.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 12345 Figure 1: Schematic diagram of the physical model and coordinate system.

Figure 8 :

 8 Figure 8: Evolution of the mean Nusselt number with the Bingham number for three different values of aspect ratio and power law index at Ra = 5 ⇤ 10 3 (left) and Ra = 10 4 (right).

Figure 9 :

 9 Figure 9: Evolution of the mean Nusselt number with the Bingham number for three different values of aspect ratio and power law index at Ra = 5 ⇤ 10 4(left) and Ra = 10 5 (right).

Figure 10 :Figure 11 :

 1011 Figure 10: Contours of non-dimensional temperature ✓ and velocity vectors with plug regions (shown with black lines) for AR = 1 (on the left) and AR = 2 (on the right) and three different values of power law index (n = 0.85, 1, 1.05) at Bn = 0.25 and Ra = 5 ⇤ 10 4 .


Damir Bilic
email: damir.bilic@mdh.se

Etienne Brosse
email: etienne.brosse@softeam.fr

Andrey Sadovykh
email: andrey.sadovykh@softeam.fr

Dragos Truscan
email: dragos.truscan@abo.fi

Hugo Bruneliere
email: hugo.bruneliere@imt-atlantique.fr

Uwe Ryssel
email: uwe.ryssel@pure-systems.com

An Integrated Model-based Tool Chain for Managing Variability in Complex System Design

Keywords: Line Engineering, Model-based Systems Engineering, Integrated Tool Chain

HAL is

I. INTRODUCTION

Product Line Engineering (PLE) [START_REF] Van Der Linden | Development and Evolution of Software Architectures for Product Families[END_REF] is commonly used in the automotive industry [START_REF] Wozniak | How automotive engineering is taking product line engineering to the extreme[END_REF]. This approach allows to combine the economic benefits of mass production with customization for specific purposes or market segments. PLE is a relatively modern approach for the design and production of softwareintensive systems sharing a common platform or code base but dedicated to different systems -variants. Product variants are obtained by configuring and enabling various common and variable features in the product line. The approach has the potential to bring up to ten-fold productivity gains as shown in an extended survey [START_REF] Berger | A survey of variability modeling in industrial practice[END_REF].

In literature, the majority of PLE studies consider softwareonly systems [START_REF] Barros-Justo | Trends in software reuse research: A tertiary study[END_REF]. Complex systems, however, might include several domains, e.g., hardware and software. Model-based systems engineering (MBSE), [START_REF] Ramos | Model-based systems engineering: An emerging approach for modern systems[END_REF], is often applied for the development of such systems, e.g. [START_REF] Andersson | Experience from introducing unified modeling language/systems modeling language at saab aerosystems[END_REF]. In addition to the implementation of software, MBSE includes activities such as requirement engineering, system architecture specification, verification and validation, etc. In the MBSE world, the model is a first-class citizen of the engineering process. Therefore, we want to make it possible to manage variability in different engineering process phases, e.g., requirement engineering, system design, verification and validation, in addition to the variability in implementation. As models, in an MBSE environment, usually prescribe and coordinate development across different engineering domains, it is important to manage variability in the complete MBSE process.

Systems Modeling Language (SysML) [START_REF] Holt | SysML for systems engineering[END_REF] is a widely used modeling language in industrial practice [START_REF] Wolny | Thirteen years of sysml: a systematic mapping study[END_REF]. In this paper, we propose an approach to deal with SysML models in a PLE-based approach. We make it possible to add variability information from variability management tools into SysML models and generate individual system variants based on the variability information. The information between the variability management tool and system modeling tool is exchanged through the Variability Exchange Language (VEL) [START_REF] Schulze | Variability exchange language-a generic exchange format for variability data[END_REF].

We instantiate the approach in practice using an integrated tool chain based on tools from the MegaM@Rt2 [START_REF]The MegaM@Rt2 ECSEL Project: MegaModelling at Runtime -Scalable Model-Based Framework for Continuous Development and Runtime Validation of Complex Systems[END_REF] and REVaMP2 [START_REF] Sadovykh | REVAMP: Challenges and innovation roadmap for variability management in round-trip engineering of software-intensive systems[END_REF] projects.

The rest of the paper is structured as follows. Section II introduces the context in which Volvo Construction Equipment (Volvo CE) develops highly variable diesel engine systems and the motivation for migrating towards an MBSE-PLE approach. Section III provides a general overview of our proposed approach. The tool chain and its application to the Volvo CE example are described in Section IV. Preliminary results of the initial application suggest that PLE concepts in an MBSE environment have the potential to improve efficiency in the development of complex and variable systems. After the discussion from Section V, we briefly dive into the related work in Section VI. Finally, Section VII summarizes the paper.

II. MOTIVATING EXAMPLE

Volvo CE develops and manufactures off-road construction machinery such as articulated haulers, excavators, wheel loaders, as shown in Figure 1. All construction machinery share a common engine platform. The common platform is outsourced to a technology provider within the Volvo Group. The Engine Department of Volvo CE then extends the platform for Volvo CE specific applications. The required engine behavior and performance are vastly different between applications. For example, in the articulated hauler, the engine is connected to a transmission, while in an excavator, the engine is connected to a hydraulic pump. In addition to the variable requirements from product platforms (e.g. excavators, haulers), the engine system must fulfill legal requirements for exhaust emissions and diagnostic services in different legislation regions (US, EU, China, and non-regulated markets).

The technology provider develops both the base engine hardware (a bare-bones engine system) and engine control software. To adapt the base engine for Volvo CE specific applications, auxiliary components, such as sensors, oil filters, turbochargers, and cooling fans, are developed at the Engine Department and added to the base engine. Consequently, the software configuration parameters are adjusted to enable or disable software features that control the installed auxiliary components. The software then needs to be calibrated to meet engine performance requirements (e.g. max engine power and torque response) at the Engine Controls Department at Volvo CE. In order to cope with the increasing complexity of both the base engine platform and the auxiliary components, the Engine Department of Volvo CE is in the process of adopting an MBSE approach based on SysML. The aim is to bring together the two different engineering domains: hardware (physical engine and auxiliary components) and software (engine controls software). However, one of the major factors that have to be considered during systems engineering at the Engine Controls department is the fact that engine systems must support a large number of variants. At the moment, SysML does not provide a means to document variability within models.

As the models at the Engine Controls department include requirements, system structure and behavior, the system modeling approach must support variability on all SysML diagram types. Moreover, the approach must allow one to trace a feature at different abstraction levels, e.g. from requirements to behavior specifications to individual software parameters.

III. SOLUTION

In order to address the needs of organizations that are developing variable products, such as the described Volvo CE example, we present an approach that combines PLE and model-based development in a common process. Figure 2 illustrates the proposed approach. The left side in the figure depicts the problem domain where the variability in the system, the relations and the constraints between variable features are Fig. 2. A generic approach for combining variability modeling with system modeling described. In our approach, we propose the use of feature models, which are both documented extensively in literature ( [START_REF] Lee | Concepts and guidelines of feature modeling for product line software engineering[END_REF], [START_REF] Sinnema | Classifying variability modeling techniques[END_REF]) and applied in industrial practice [START_REF] Berger | A survey of variability modeling in industrial practice[END_REF].

Once the system variability is described in the variability model, the variability must also be implemented in the solution domain. The middle part of Figure 2 depicts the solution domain. This is where the system models of both common and variable features are developed.

To document variability in models we use 150% models ( [14], [START_REF] Heidenreich | Featuremapper: mapping features to models[END_REF], [START_REF] Mann | Dealing with variability in architecture descriptions to support automotive product lines: Specification and analysis methods[END_REF]) i.e. all system variants are modeled within a single over-specified model. The 150% approach to product lines has the potential to increase product quality, reduce amount of artifacts that need to be maintained, reduce leadtime, etc. [START_REF] Beuche | When 150% is too much: supporting product centric viewpoints in an industrial product line[END_REF]. The 150% model contains all system variations within a single model. It is important to consistently exchange variability information between the problem and solution domain. In our proposed tool chain, we use the Variability Exchange Language to integrate the design and variability management tools (cf. Section IV).

Variability models serve as a starting point for customers during the definition of a product variant. Customers select variable product features based on their needs. Constraints and relations between features restrict how features can be combined. The variability information then needs to be exchanged with the system modeling tool that contains the 150% model. The system modeling tool should then be able to generate a 100% system model of a single system variant. As illustrated in the right part of Figure 2, the derived model should include only model artifacts from the solution domain that were selected during feature selection. For example, these can be requirement models, system architecture, verification procedures, etc.

Since different organizations may have different systems engineering processes and use models at different levels of the process, we strive that our approach is generic enough to be applicable to any phase of their systems engineering process by making it possible to model variability on any diagram type, e.g. requirement models, structure, behavior, etc.

IV. TOOL CHAIN

This section describes a tool chain that supports the approach discussed in the previous section.

In the tool chain, we use the pure::variants tool from pure-systems [START_REF] Beuche | Variants and variability management with pure:: variants[END_REF] to manage variability in the problem domain and the Modelio [START_REF] Bagnato | High level methodologies in embedded system design[END_REF] modeling tool to model the system in the solution domain. In order to integrate the two tools, the Variability Exchange Language was used. It allows us to exchange variability information between the problem and solution domain. pure::variants, as a variability management tool, already supported VEL. However, Modelio had to be extended with mechanisms to support modeling of variability in system models and to exchange information via VEL.

In Figure 3, we illustrate the proposed tool chain workflow in more detail in the context of a Volvo CE use case. It covers several phases of an engine development process starting with requirements and continuing with system design and verification, to some extent.

In Step 1 , customers (product platforms) define their requirements on the system, e.g. needed torque, power, engine start/stop function, etc. In addition, depending on the legislation region in which the customer plans to use the engine, the exhaust emission and diagnostic requirements are varying.

In Step 2 , based on their requirements, customers select system features from the feature model with the pure::variants tool, as illustrated in Figure 4.

pure::variants is a variability management tool supporting VEL to facilitate product line development. As depicted in Figure 4, a feature model is specified to document all variability of the product line. A set of system features and the possible variants of the feature are shown. The example shows the variability of the Urea Dosing System (UDS) which is a subsystem of the engine. A single UDS consists of a Diesel Exhaust Fluid (DEF) 1 storage tank, a pump that pressurizes the DEF and and injection nozzle through which 1 DEF is a mixture of urea and water that chemically reacts with exhaust gases to reduce harmful emissions.

the DEF is dispersed into the exhaust stream. There are also a DEF quality sensor and a valve that controls engine coolant circulation throughout the system in order to defrost it at low temperatures. An actuator control unit (ACU) is used to control the pump in the system.

The constraints between features are defined manually by system engineers during the system design. The variability management tool only defines and manages the system feature model without any other knowledge about how the features are implemented. The feature model and feature constraints restrict the variability space of the system, reducing the risk of creating a variant with incompatible features. Some users might have requirements that are outside the scope of the current product family, i.e. a feature not present in the current set of features. In such a situation there are two possibilities: (i) extend the product line to include the new feature or (ii) implement the new requirement only for the individual variant. Both possibilities will be discussed later in this section. Once Modelio is an extensible UML/BPMN/SysML modeling tool. A module has been developed to extend the tool and allow variability information annotation on top of any kind of supported modeling language (e.g. UML/SysML in our context). As shown in Figure 5, a SysML Block Definition Diagram (BDD) of the UDS is annotated with variability information. It is a desensitized diagram that includes two variants of the UDS. The blocks of the two variants are annotated as 'Bosch' and 'Albonair' in Figure 5.

In order to exchange information between Modelio and pure::variants, Modelio first generates a serialized XML file in the VEL, as depicted in Figure 6. The XML file contains information of the 150% model. For simplicity purposes, we show one variation point, named 'Physical Architecture', with two possible variants, named 'Albonair' and 'Bosch', together with information of model elements that belong to either feature. The VEL file is then read by pure::variants and mapped to the feature diagram shown in Figure 4. Mapping is done through a naming convention based on the prototype implementation of the variability modeling functionality in Modelio. The features that the user did not select in Step 2 are then removed from the VEL file in pure::variants and a new VEL file is then generated.

This new file contains information with the feature set of a single variant, i.e., a 100% file. Modelio uses the new VEL file as an input for a model transformation where elements not contained in the VEL file of a single variant are removed from the 150% system model. As an example, in 3 all product family requirements are contained. Modelio, performs a model transformation based on the VEL file exchange to generate a 100% requirement model that is stored separately into 6 .

The requirements define the performance and driveability goals of the base engine system and auxiliary components. For example, the duration of the engine warm-up function when operating in cold conditions. The base engine software features and interfaces are configured in 4 . As all engine variants are built on a single engine platform, it is necessary to enable or disable software functions and system interfaces based on the application. The configuration is performed by enabling or disabling software parameters. The derivation of software configuration parameters is again automated based on the same feature selection from 2 and the single variant model elements are again stored into the new model repository, 6 .

Finally, in 5 it is necessary to select appropriate engine auxiliary components (sensors, oil tanks, hoses, radiators, etc.) and configure software functions that control these. Again, based on the feature selection from 2 , the annotated 150% model (containing the all auxiliary component design, software parameters, and tests) is transformed to remove all features that were not selected and generate a 100% model containing only the selected elements.

The variations in requirements, base engine, and auxiliary components are resolved and an individual model is generated and stored in a new model repository in 6 . This final system model then drives the engine assembly and activities such as performance calibration.

Now looking back at the case when customer requirements can not be mapped to any feature in the product line. The first possibility is to extend the product line and implement the required feature. That is, the feature model in 2 must be extended to include the new feature. Then, the requirements in 3 and system design in 4 and 5 must be extended to implement the new feature. The process from selecting variants in the feature model to product derivation needs to be performed again to derive the new variant. The second possibility is to select the closest possible variant in 2 , then once a product based on the feature selection is derived in 6 , it is modified to implement the new requirements. The first approach develops the new feature to be reusable in new projects in the future, while the second approach develops the feature for only one variant.

V. DISCUSSION

We first discuss the approach based on the quality criteria for product lines as defined in [START_REF] Apel | Feature-oriented software product lines[END_REF].

Pre-planning effort. Implementing product lines always requires some pre-planning, at least to define the scope of the product line. System models are often created iteratively and as it is with annotative approaches in general ( [START_REF] Apel | Feature-oriented software product lines[END_REF]), our approach enables the product line to be annotated as system models are evolving with little pre-planning. However, it is important to consider that no or little pre-planning might require significant effort to implement new features in later stages, as it might require large parts of the product to change.

Feature traceability evaluates support for mapping of features between the feature model and implementation models. The use of VEL for exchange of variability information allows Fig. 6. Generated VEL us to trace the features from the feature model directly to system model elements. VEL also allows us to trace a feature to multiple SysML models and model elements. That is, a feature can be mapped to any number of model elements that implement the given feature.

Separation of concerns, i.e. capability of the tool chain to support a separation between variable aspects of system features. Traditionally, in model-based development, views and viewpoints are used to express different concerns of the system design [START_REF] Cicchetti | Chess: a model-driven engineering tool environment for aiding the development of complex industrial systems[END_REF]. A view is a representation of the system from a perspective (viewpoint) of a stakeholder of a related concern. With SysML, a view can be considered as a set of model diagrams that allow the users to visualize, analyze and reason about a specific system concern or feature [START_REF]IEEE Recommended Practice for Architectural Description for Software-Intensive Systems[END_REF]. Separation of concerns mainly depends on the implementation of the system, i.e. whether the system is organized into cohesive components. Cross-cutting concerns, for example, the diagnostic function of the engine, however, span across multiple system components. In such cases, annotating model elements with variability information explicitly exposes the cross-cutting features.

Information hiding, i.e. separation of a component into its internal logic and external interfaces is inherently supported in SysML. Since our tool chain supports annotation on any model element, encapsulated model elements of any feature can be annotated with variability, thus enabling information hiding with regards to features.

Granularity, i.e. the level of detail at which variability can be implemented in system models. Our proposed annotative approach provides fine-grained granularity. That is, it allows annotation on all model element types and properties of such elements, as well as annotation on the relationship links between elements.

Uniformity states that variability in all product line artifacts should be encoded consistently. Our tool chain provides good uniformity as all product artifacts are SysML model elements and are annotated in the same fashion.

Based on the initial application, it was noted that the feature model and the constraints between features improve the consistency of system variants as the constraints restrict the possibility to create an invalid variant. The requirement elicitation time from customers (product platforms in the Volvo CE case) can be reduced since customers do not need to interact with domain experts to define a system variant except in the situation when customer requirements can not be mapped to any feature. This situation was described in Section IV. The feature model also received positive feedback from managers and non-technical stakeholders as it allowed them to understand the implications of certain system variants without needing extensive technical knowledge of the system or the domain.

VI. RELATED WORK

MBSE with SysML has been widely used in the industry in various domains, for example, aerospace [START_REF] Andersson | Experience from introducing unified modeling language/systems modeling language at saab aerosystems[END_REF] and automotive [START_REF] Góngora | Towards an architectural design framework for automotive systems development[END_REF]. SysML can be extended with profiles and stereotypes to support variability [START_REF] Weilkiens | Variant modeling with SysML[END_REF]. This approach does not require a separate variability management tool. However, compared to our approach, it does not support the modeling of variability in SysML behavior diagrams, which are an important part of the system design. Another approach for modeling variability with SysML is presented in [START_REF] Meacham | Sysml based design for variability enabling the reusability of legacy systems towards the support of diverse standard compliant implementations or standard updates: the case of ieee-802.15. 6 standard for e-health applications[END_REF]. But again, it is not clear how variability is captured in behavior models. As stated in an literature survey, [START_REF] Wolny | Thirteen years of sysml: a systematic mapping study[END_REF], activity and state machine diagrams are widely used in systems engineering. Thus, it is important to be able to define variants in behavior models as well.

The pure::variants tool was applied in a number of studies and industrial examples such as: [START_REF] Beuche | When 150% is too much: supporting product centric viewpoints in an industrial product line[END_REF] [START_REF] Rubin | Managing cloned variants: a framework and experience[END_REF]. The later being interesting to us as it uses VEL to add variability information on software code. Combining our approaches would allow to have feature traceability through the complete development cycle of a system, from requirements to the system design and then the actual system implementation.

Another use of the same variability management tool and VEL is reported in the experience report in [START_REF] Pohl | Variant management solution for large scale software product lines[END_REF]. It describes a methodical and tooling approach with similar concepts to ours. However, their approach uses a custom-built tool to support variability in the problem domain throughout the complete application lifecycle.

Feature modeling, of course, is not the only approach and many of them are compared in [START_REF] Sinnema | Classifying variability modeling techniques[END_REF]. For example, model views can also be used to handle variability in a modeling context [START_REF] Bruneliere | A Feature-based Survey of Model View Approaches[END_REF]. Another study suggests Orthogonal Variability Modeling for variability management in SysML models [START_REF] Bilic | Model-based product line engineering in an industrial automotive context: An exploratory case study[END_REF]. Nevertheless, feature-modeling seems to be the most popular one based on a survey with industrial practitioners [START_REF] Berger | A survey of variability modeling in industrial practice[END_REF].

VII. CONCLUSION

MBSE is a common approach in practice for the development of complex, software-intensive systems. However, as presented in the Volvo CE example, organizations face challenges where many similar products must be developed and maintained for different customers. To deal with the highly variable systems, PLE is often applied. To combine the benefits of both worlds, PLE and MBSE, we have presented an approach that allows us to integrate variability information into design models of complex systems. The approach is supported by our tool chain for which we extend the Modelio system modeling tool to support variability modeling in SysML models. In addition, we use the VEL to exchange variability information with the pure::variants tool for variability modeling and management. The approach is demonstrated through an example from Volvo CE. The initial application suggests a potential improvement in efficiency of the engine development process.

Fig. 1 .

 1 Fig. 1. Volvo CE Machinery

Fig. 3 .

 3 Fig. 3. Product Line Engineering Process with the Model-based Tool Chain

Fig. 4 .

 4 Fig. 4. Feature Model example

Fig. 5 .

 5 Fig. 5. Annotated SysML model

ACKNOWLEDGMENT

This work has received funding from the Electronic Component Systems for European Leadership Joint Undertaking under grant agreement No. 737494. This Joint Undertaking receives support from the European Union's Horizon 2020 research and innovation program and from Sweden, France, Spain, Italy, Finland and Czech Republic. This work was also partially funded by the ITEA3 15010 REVaMP2 project: FUI the le-de-France region and BPI in France.