**3. Content of Uzbekistan cotton germplasm collection**

#### **3.1. Main collections**

germplasm resources, which are being preserved and maintained for cotton genetics and breeding research that target current and future needs of the cotton improvement for different

Aspects of Uzbekistan cotton germplasm resources including the content, distribution, descriptions, characterizations, utilization, genetic and molecular diversity, maintenance, and ongoing and prospective research previously have been highlighted in several documents [2; 3; 7; 8; 9; 10]. In this chapter, we provide a detailed inventory of the Uzbekistan cotton germplasm collection, review previous reports and add updated information including the

The past century of cotton production in Uzbekistan has developed well-established cotton research programs and distinguished cotton scientists who initiated and devoted themselves to collecting important materials for cotton research. As highlighted in previous reports [3; 8] the cotton germplasm collection initiative was began by Drs. N. I. Vavilov and F. M. Mauer in 1930 in the former Soviet Union. Subsequently, Uzbekistan cotton germplasm founder and leader Dr. A. Abdullaev and his group expanded this initiative and established a collection of Uzbekistan germplasm materials through the (1) coordination of scientific efforts of continu‐ ous selection of cultivated cotton varieties, (2) continuation of collecting and preserving wild cotton species and landraces from centres of origin during many scientific expeditions, and (3)

According to Abdullaev et al. [8], several expeditions to Central Asia, Afghanistan, China, India, Turkey, Iran, Korea and Japan to obtain germplasm were made during 1920-1930, Dr. N. I. Vavilov, P. M. Jukovsky and Dr. F. M. Mauer, and in 1950s Dr. D. V. Ter-Avanesyan. In later periods from 1974 to 1998, Dr. A. Abdullaev visited to Mexico, Peru, China, India and Sri-Lanka, Australia and Pakistan and obtained germplasm. These historic scientific expedi‐ tions enriched Uzbekistan collection with Old World (Afro-Asian and Indian) diploid cottons (*G. herbaceum* and *G. arboreum*), and a number of wild, exotic and cultivated tetraploid and

The Uzbekistan collection has been periodically enriched as a result of germplasm exchanges with collections worldwide. In the most recent exchanges within the framework of USDA-Uzbekistan Cooperation programs, approximately ~ 1000 *G. hirsutum* exotic and varietal accessions were exchanged with the US cotton germplasm collection [11; 12]. Annually, 100-120 cotton accessions from the collection of Uzbek Research Institute of Cotton Breeding and Seed Production (UzSRICBSP) are exchanged with world centres. During the period of 2001-2003, the Institute received 990 accessions from and sent 260 cotton accessions to the US

soil-climatic regions of Uzbekistan [1; 2; 3; 8].

290 World Cotton Germplasm Resources

germplasm exchange worldwide.

diploid cottons around the world [8].

cotton germplasm collection [7].

development and characterization of novel germplasm resources.

**2. History and development of Uzbekistan cotton collection**

The main cotton germplasm collections are being historically preserved at the research centres and institutions of the Academy of Sciences of Uzbekistan (ASUz), Ministry of Agriculture and Water Resources of Uzbekistan (MAWR), and the biology department of the National University (NU) of Uzbekistan. Table 1 summarizes and highlights the general content and description of cotton germplasm resources of these main collections. These collections maintain cultivars, wild and primitive, predomesticated landraces, hybrids breeding and genetic stocks, cytogenetic and mutant lines of widely cultivated allotetra‐ ploids(*G. hirsutum* and *G. barbadense*) representing the primary gene pool, and two cultivated Asian diploids (*G. herbaceum* and *G. arboreum*) representing the secondary gene pool. Although some redundancy of accessions maintained by collections could be possible and is a subject for future comparative inventory work, each collection has its own specifici‐ ties and has been formed according to distinctive research efforts conducted by the scientists of these institutions for past decades.

collection, also referred to as "wild collection", was formed as a result of scientific expeditions to the centres of origin of *Gossypium* species and a 50 years of research efforts by IG&PEB scientists and research programs. IG&PEB cotton collection maintains more than 40 wild A-to G and K-genome *Gossypium* species [2; 3; 8] and targets the study of the taxonomy, phylogeny, evolution, hybridization compatibility, and introgression of wild germplasm for breeding purposes, all of which make the collection unique compared to others. The detailed descrip‐ tion, content, geographic coverage and history of *G. hirsutum* and *G. barbadense* germplasm development in Uzbekistan were discussed by Abdullaev et al., 2013, where authors report the representation of 4 continents, ~33 geographic regions and ~103 countries by the IG&PEB collection [8].There are a large number di-and tri-genomic hybrids and their diploid and allotetraploid progenies recovered from tedious sexual crosses within and/or between wild, primitive and domesticated genotypes within different gene pools [13; 14; Dr. Rizaeva, personal communication]. Some of examples such as tri-genomic hybrids were highlighted in

Cotton Germplasm Collection of Uzbekistan http://dx.doi.org/10.5772/58590 293

MAWR has two distinctive cotton collections: one is preserved at the UzSRICBSP, anoth‐ er is in the Uzbek Research Institute of Plant Industry (UzRIPI), which was a Central Asian branch of All Union scientific-research institute after N.I. Vavilov. The UzSRICBSP collection preserves more than 12,000 cotton accessions from primary and secondary gene pools and refereed to as "breeding" germplasm resources that resulted from continuous breeding and selection efforts of the institute's scientists as well as cotton germplasm exchange efforts [16]. Geographically, this collection represents 107 countries of origin for cotton acces‐ sions [7]. The uniqueness of UzSRICBSP collection is its wide representation and cover‐ age of cultivar germplasm developed and collected over the past century from worldwide breeding efforts. This collection also maintains synthetic tetraploid, pentaploid and octoploid hybrids [7; 16]. There are small differences in germplasm accession numbers reported here (Table 1) and by Ibragimov et al. [7]. However, our inventory is based on the latest information obtained from this collection (Dr. H. Saydaliev, a germplasm curator

The UzRIPI cotton collection has contents similar to those of the IG&PEB collection and preserves a total of over 6,000 accessions (Table 1) from primary and secondary gene pools as well as accessions of other gene pools of wild species. Among all the collections, UzRIPI cotton collection is the richest resource for primitive, and pre-domesticated landrace stocks for all *Gossypium* gene pools in the country. However, there is a need to conduct comparative reinventory between UzRIPI and IG&PEB collection to identify the distinctive versus redundant germplasm accessions. Because UZRIPI was the Central Asian branch of All Union scientificresearch institute after N.I. Vavilov some level of redundancies to Russian VIR collection is

The NU collection is tasked with maintaining of a total about 1200 germplasm resources that include 771 genetic stocks and recombinant inbred lines, formed during study of key cotton traits and mutations [3; 17; 18]. Additionally, the NU collection has a unique set of over 400 cytogenetics stocks primarily derived from radio-mutagenesis of a single genotype of *G.*

previous reports [14; 15].

of the UzSRICBSP, personal communication).

expected that requires future study.

**Table 1.** Summary of the content of Uzbekistan cotton germplasm collections

Cotton germplasm collection of the Institute of Genetics and Plant Experimental Biology (IG&PEB) of the ASUz, founded and headed by academician and Prof. Abdumavlyan Abdullaev, preserves a total of 7,500 cotton accessions. The collection specifically aims to gather, maintain, and study wild cotton species, primitive, pre-domesticated landraces and domesticated genotypes from entire *Gossypium* genus. The IG&PEB cotton germplasm collection, also referred to as "wild collection", was formed as a result of scientific expeditions to the centres of origin of *Gossypium* species and a 50 years of research efforts by IG&PEB scientists and research programs. IG&PEB cotton collection maintains more than 40 wild A-to G and K-genome *Gossypium* species [2; 3; 8] and targets the study of the taxonomy, phylogeny, evolution, hybridization compatibility, and introgression of wild germplasm for breeding purposes, all of which make the collection unique compared to others. The detailed descrip‐ tion, content, geographic coverage and history of *G. hirsutum* and *G. barbadense* germplasm development in Uzbekistan were discussed by Abdullaev et al., 2013, where authors report the representation of 4 continents, ~33 geographic regions and ~103 countries by the IG&PEB collection [8].There are a large number di-and tri-genomic hybrids and their diploid and allotetraploid progenies recovered from tedious sexual crosses within and/or between wild, primitive and domesticated genotypes within different gene pools [13; 14; Dr. Rizaeva, personal communication]. Some of examples such as tri-genomic hybrids were highlighted in previous reports [14; 15].

 The collection specifically aims to gather, maintain, and study wild cotton species, primitive, pre-domesticated landraces and domesticated genotypes from entire *Gossypium* genus. The IG&PEB cotton germplasm

**Table 1.** Summary of the content of Uzbekistan cotton germplasm collections

771 - - - - 771

986 286 - - - 1272

20 - - - - 20

1852 1852

24 24

**Other species**

**Total**

0

MAWR has two distinctive cotton collections: one is preserved at the UzSRICBSP, anoth‐ er is in the Uzbek Research Institute of Plant Industry (UzRIPI), which was a Central Asian branch of All Union scientific-research institute after N.I. Vavilov. The UzSRICBSP collection preserves more than 12,000 cotton accessions from primary and secondary gene pools and refereed to as "breeding" germplasm resources that resulted from continuous breeding and selection efforts of the institute's scientists as well as cotton germplasm exchange efforts [16]. Geographically, this collection represents 107 countries of origin for cotton acces‐ sions [7]. The uniqueness of UzSRICBSP collection is its wide representation and cover‐ age of cultivar germplasm developed and collected over the past century from worldwide breeding efforts. This collection also maintains synthetic tetraploid, pentaploid and octoploid hybrids [7; 16]. There are small differences in germplasm accession numbers reported here (Table 1) and by Ibragimov et al. [7]. However, our inventory is based on the latest information obtained from this collection (Dr. H. Saydaliev, a germplasm curator of the UzSRICBSP, personal communication).

The UzRIPI cotton collection has contents similar to those of the IG&PEB collection and preserves a total of over 6,000 accessions (Table 1) from primary and secondary gene pools as well as accessions of other gene pools of wild species. Among all the collections, UzRIPI cotton collection is the richest resource for primitive, and pre-domesticated landrace stocks for all *Gossypium* gene pools in the country. However, there is a need to conduct comparative reinventory between UzRIPI and IG&PEB collection to identify the distinctive versus redundant germplasm accessions. Because UZRIPI was the Central Asian branch of All Union scientificresearch institute after N.I. Vavilov some level of redundancies to Russian VIR collection is expected that requires future study.

The NU collection is tasked with maintaining of a total about 1200 germplasm resources that include 771 genetic stocks and recombinant inbred lines, formed during study of key cotton traits and mutations [3; 17; 18]. Additionally, the NU collection has a unique set of over 400 cytogenetics stocks primarily derived from radio-mutagenesis of a single genotype of *G.* *hirsutum* line L-458 [19; 20]. Readers can find detailed description of Uzbekistan's cytogenetic cotton collection in this book.

**4. Storage, maintenance and funding**

al. [3] and Abdullaev et al. [8].

traits [3; 8].

regenerated lines.

The above-mentioned main collections and novel germplasm resources are stored and maintained in each institution and managed separately by its scientists. The IG&PEB, UzS‐ RICBSP, and NU collections are stored under room temperature conditions (20-25°C) and there is no facility available for cold storage of germplasm accessions [3; 8]. In contrast, UzRIPI [33] and CGB collections are stored in short term (under+4°C) cold room facilities that were established as a result of government and international funding (in the case of UzRIPI, [33]. No long term cold storage (-20 or-80 °C) facilities, requiring attention and investment, are available for any of the germplasm collections in Uzbekistan, as highlighted by Campbell et

Cotton Germplasm Collection of Uzbekistan http://dx.doi.org/10.5772/58590 295

Germplasm accession seeds are preserved according to commonly practised procedures used over the decades of germplasm maintenance efforts in each collection. For instance, according to previous reports [8] germplasm seeds are ginned and put into paper bags with a label of catalogue number, accession name, year of collection and origin. Paper bags also contain "*the weight of seeds (50 g or 100 g individual or total pick respectively*") and "*bags are stored in special metal boxes (30 x 11 cm) and boxes are placed in wooden-cases"* [8]. Other collections follow similar storage procedure with some modifications in types of storage boxes and variations in labelling of bags. Since 2003, after the reconstruction of the building for germplasm resources at UZRIPI,

Consequently, germplasm resources, in particular those without short term cold room facility, are scheduled for seed renewal every 8-10 years under forced self-pollination in the open field conditions [3]. Each organization has its own, but very similar protocols, schemes of planting, growing and evaluating germplasm (see [8] for detailed protocol for IG&PEB cotton collec‐ tion), government research grants and field extension stations with up to 8-10 staff working on germplasm maintenance. During a seed renewal year, accessions of re-grown germplasm are phenotypically evaluated for major agronomic and morpho-biological and fibre quality

Germplasm evaluation records from each collection are maintained as a hard copy catalogue book that contains all descriptions about accessions (origin, year, collector, collected sites, seed renewal, etc), and data from the past 50 years of evaluations [e.g. 34] In UzRIPI, there is a "Unified Council for Mutual Economic Assistance (COMECON) list of descriptors for the species *Gossypium L* [35]. IG&PEB has started using modified "Cotton descriptor" of International Plant Genetic Resource Institution [see below for detail description; 3; 8]. Moreover, data records on germplasm accession characteristics and description are being currently entered into electronic Microsoft database formats (personal communications with germplasm heads of all collections). All cotton collections also have their own green‐ house facilities to vegetatively maintain [7; 8] wild and primitive accessions, unique multigenomic hybrids, mutants, and cytogenetics stocks as well as transformed and somatically

cotton germplasm seed have been stored in plastic containers [33].

#### **3.2. Novel resources**

Efforts focused on genetic mapping of important traits, application of marker-assisted breeding as an aid for contemporary breeding, and the development of cotton tissue culturebased transgenomics programs and their integration into conventional cotton improve‐ ment efforts have resulted in the creation and collection of novel germplasm resources in Uzbekistan. These novel germplasm resources were developed in the past decade by scientists of Centre of Genomics and Bioinformatics (CGB), ASUz, MAWR, and "UzcottonIn‐ dustry association" within the framework of international collaborations and government funding [1; 21; 22]. The CGB collection with over 5,000 germplasm resources (Table 1) includes (i) panels of association mapping and nested-association mapping populations [22; 23; 24], (ii) hybrids and recombinant inbred lines (F1-4) derived from the combination of sexual top crosses between 9 commercialized Uzbek cotton cultivars and 16 different chromosome substituted lines (CS-B) [25; 26; 27; 28; 29], (iii) germplasm resources devel‐ oped through marker-assisted selection (MAS) programs that bear novel quantitative trait loci (QTL) loci mobilized from unique donors to the genetic background of commercial Uzbek cotton cultivars via molecular markers, and (iv) tissue culture-derived, genetically modified (GM) germplasm and their hybrids to local cultivars that bear RNA interference (RNAi), synthetic hairpin oligonucleotides, anti-sense, or overexpression genetic con‐ structs for *de novo* characterized genes and sequence signatures in the CGB laboratories [1; 22; 30; 31; 32].

It is noteworthy to mention that association mapping individuals (Table 1) were selected from the IG&PEB collection and re-grown at the Mexico Winter Nursery of USDA-ARS by Drs. Russel Kohel and John Yu, (cotton germplasm Unit of USDA-ARS at College Sta‐ tion,Texas) for phenotypic evaluations and seed increase. Increased seeds grown at the Mexican environment kindly were sent back to Uzbekistan by Dr. Richard Percy (USDA cotton germplasm curator) and currently backed-up at the CGB and IG&PEB collections. Additionally, chromosome substituted (CS-B) lines were received within the frame of USDA-Uzbekistan Cooperation Programs, kindly provided by Dr. David Stelly (Texas A&M University), Dr. Sukumar Saha and Dr. Johnie Jenkins, USDA-ARS, Starkville, Mississip‐ pi, and now are preserved in both CGB and IGPEB collections. CGB scientists in collabora‐ tion with USDA partners are developing CS-B specific chromosome substituted recombinant inbred lines (CSRILs) in the background of important Uzbekistan cultivars. Further, development of cotton tissue culture and trangenomics efforts [21], and the mobilization of useful genetic constructs into commercialized cultivars has created novel germplasm resources, useful for cotton improvement and helpful to address many problems associat‐ ed with improving and boosting yield and quality [1; 22; 32].


2.5.2. Theoretical Analysis

In order to confirm the validity of the numerically predicted values of all performance characteristics of the thermoelectric module, the same performance characteristics were evaluated with a theoretical approach. The voltage at maximum power indicates optimum voltage and is calculated from Equation (6) [1,31,32]. Here, α is the Seebeck coefficient, which is measured as the amount of thermoelectric voltage induced in any material due to temperature difference across that material. Each thermoelectric material has its own unique value of Seebeck coefficient. The Seebeck coefficient is also recognized as thermopower, thermoelectric power, and thermoelectric sensitivity:

$$V\_{\text{theoretical},opt} = \alpha \left( T\_H - T\_{\text{C}} \right) \tag{6}$$

Current depends on voltage, external load resistance, and internal resistance—Equation (7) [10]. *R* = *RL* is the condition to obtain optimal current wherein power is maximum [4,31]:

$$I\_{\text{theoretical},opt} = \frac{V\_{\text{theoretical,opt}}}{R + R\_L} = \frac{\alpha \left(T\_H - T\_{\odot}\right)}{2R} \tag{7}$$

By manipulating the equations of *I* = α (*TH*−*TC*) *R*+*RL* and *P* = *I* 2 (*R* + *RL*), the power equation was obtained. To obtain maximum power, constraint of ∂*P* ∂*RL* = 0 was applied on the power equation, which gave the maximum power condition of *R* = *RL*. When *R* = *RL* is used in the power equation, Equation (8) is derived, which calculates the theoretical maximum power of the thermoelectric module [4,33]:

$$P\_{\text{theoretical,max}} = \frac{\alpha^2}{4R}(T\_H - T\_C)^2\tag{8}$$

The theoretical maximum efficiency was calculated using Equation (9) under the same condition of *R* = *RL* in the theoretical efficiency equation of η*theoretical* = ∆*T*·*RL TH*(*R*+*RL*)− ∆*T*·*R* 2 + (*R*+*RL* ) 2 *ZR* [4].

$$
\eta\_{\text{theoretical},\text{max}} = \Delta T \frac{\sqrt{1 + Z\overline{T}} - 1}{T\_H \sqrt{1 + Z\overline{T}} + T\_{\text{C}}} \tag{9}
$$

Resistance for square prism and cylindrical legs was calculated with Equation (10) [1].

$$R = \frac{(\rho\_p + \rho\_n)L}{A} \tag{10}$$

Resistance of the trapezoidal legs cannot be calculated using Equation (10) because of the non-uniform cross-sectional area, which varies along the height of the thermoelectric legs. In the trapezoidal legs, an unevenly distributed area along the centerline of leg in height was converted to an evenly distributed area along the same centerline of the leg such that volume of the leg was same in both cases. Using the uniform area approach, the resistance for the trapezoidal legs was calculated with Equation (11) [11].

$$R = \frac{\sigma\_p + \sigma\_n}{2\sigma\_p \sigma\_n \frac{A\_0}{L} \left(\frac{R\_A - 1}{R\_A + 1}\right)} \text{ } \ln(R\_A) \tag{11}$$

where *RA* is the ratio of the top side area of the trapezoidal leg (*AT*) to the bottom side area of the trapezoidal leg (*AB*) from Equation (12). *A*0 is the equivalent uniform cross-section area of the trapezoidal leg:

$$R\_A = \frac{A\_T}{A\_B} \tag{12}$$

In order to evaluate the maximum power and maximum efficiency of the thermoelectric module with the two-stage arrangement, it is necessary to derive the intermediate temperature between the two stages. This intermediate temperature can be obtained using voltage and temperature difference correlation equation of *V* = α (*TH* − *Ti*) or α (*Ti* − *TC*). In the present study, the voltages applied to both the stages were different; hence, the equations were modified to *VTop* = α (*TH* − *Ti*1) and *VBottom* = α (*Ti*2 − *TC*). *Ti*1 is the cold side temperature for the first stage, whereas *Ti*2 is the hot side temperature of the second stage. The difference of *Ti*1 − *Ti*2 shows the temperature loss due to the thermal resistance of the intermediate plate. The voltage boundary conditions applied to the high potential sides of both the stages were inserted in the corresponding voltage equations, obtaining *Ti*1 and *Ti*2 temperatures. These temperatures are used in Equations (8) and (9) in order to calculate the maximum power and maximum efficiency of the two-stage arrangement of the thermoelectric module. The segmented thermoelectric module consists of a combination of two materials in a single thermoelectric leg. Hence, the concept of equivalent Seebeck coefficient, equivalent resistivity, equivalent thermal conductivity, and equivalent Figure of Merit was used. Using this equivalent parameter concept, maximum power and maximum efficiency of the segmented thermoelectric module were calculated theoretically from the following equations [22,34]:

$$P\_{\text{theoretical},\max} = \frac{\alpha\_{\text{equivalent}}^2}{R\_{\text{equivalent}}} (T\_H - T\_\text{C})^2 \tag{13}$$

$$
\alpha\_{equivalent} = \frac{\alpha\_2 k\_1 + \alpha\_1 k\_2}{k\_1 + k\_2} \tag{14}
$$

$$\eta\_{\text{theoretical},\text{max}} = \Delta T \frac{\sqrt{1 + Z\_{\text{equivalent}}\overline{T}} - 1}{T\_H \sqrt{1 + Z\_{\text{equivalent}}\overline{T}} + T\_\mathbb{C}} \tag{15}$$

$$Z\_{equivalent} = \frac{\alpha\_{equivalent}^2}{(\rho\_1 + \rho\_2) \* \left(\frac{k\_1 k\_2}{k\_1 + k\_2}\right)}\tag{16}$$

Thermal stress induced in the thermoelectric modules due to various temperature conditions, as well as due to mismatch in thermal properties of different materials, were calculated in the form of von Mises stress. This is the equivalent of the second deviatory stresses and computed using the following equation [1,3]:

$$\sigma\_{\upsilon} = \sqrt{\frac{\left(\sigma\_{11} - \sigma\_{22}\right)^2 + \left(\sigma\_{22} - \sigma\_{33}\right)^2 + \left(\sigma\_{33} - \sigma\_{11}\right)^2 + 6\left(\sigma\_{12}^2 + \sigma\_{23}^2 + \sigma\_{31}^2\right)}{2}} \tag{17}$$

Normal stress in the longitudinal direction can be calculated using the following equation. In this equation, normal strain ε*n* = 0 and curvature K = 0 for the present study [3] are:

$$
\sigma\_{\rm xx} = \frac{E}{1 - v^2} \left( \varepsilon\_{\rm II} + Z\_{\rm c} K \right) - \frac{E \alpha\_{\rm T}}{1 - v} (T\_H - T\_{\rm C}) = -\frac{E \alpha\_{\rm T}}{1 - v} (T\_H - T\_{\rm C}) \tag{18}
$$

Thermal stress in the direction of temperature gradient can be calculated using the following equations, which are based on temperature difference and thermal properties.

$$
\sigma\_T = \alpha\_T \mathbb{E}(T\_H - T\_\mathbb{C}) \tag{19}
$$

$$
\sigma\_T = \alpha\_T E \Delta T \tag{20}
$$

Thermal strain is the ratio of thermal stress to the Young modulus, which is represented by Equation (21).

$$
\varepsilon = \frac{\sigma\_T}{E} \tag{21}
$$


**5. Biomarkers in Onco-Immunotherapy**

As we have already described, immunotherapies are one of the most promising approaches to treat cancer patients, but despite the demonstrated success in a variety of malignant tumours, the responses only occur in a minority of patients. Furthermore, these treatments involve inflammatory toxicity and a high cost. Therefore, determining which patients would derive clinical benefit from immunotherapy is an important goal. This requires the identification and validation of prognostic biomarkers. The integration of multiple tumour and immune response parameters, such as protein expression, genomics and transcriptomics, may be necessary to accurately predict clinical benefit.


*2.1. International Tourism and Income*

Previous research suggested that income, prices, and the currency rate have a considerable effect on tourist inflow, but there is significant variation in the degree of influence between variables. Many factors can explain this. According to Peng et al. (2015), there is a large spread in the income elasticity depending on the destination, country of origin, season, and type of holiday. Income in the origin country seems to be a dominant variable in explaining the level of international tourism. Sanchez-Rivero and Pulido-Fernández (2020) suggested that the average income elasticity of visitors crossing international borders is substantially higher than 1.0 (their estimate was just over 2.0). If this is the case, an increase of one percent will cause the demand for tourist travel to rise far beyond one percent. Economic theory refers to this as a luxury good. There may be wide gaps for the same destination depending on the country of origin, as well as large differences within the same country of origin depending on the destination. Due to a lower income level, income elasticity tends to be higher for visitors from countries with low GDP per capita compared to countries with high GDP per capita. There is limited holiday time, and the choice of destination can be sensitive to changes in some important variables. Crouch (1996) reported an income elasticity of 1.5 for international tourism. Recent research reports suggest that income elasticity is between 1.0 and 2.0, but with significant variations Kumar and Kumar (2020); Ongan et al. (2017); Sanchez-Rivero and Pulido-Fernández (2020); however, Jensen (1998) pointed out that the income elasticity of foreign visitors is considerably higher than that of domestic tourist visitors.


# **Women in Artificial intelligence (AI)**

Aida Valls and Karina Gibert Printed Edition of the Special Issue Published in *Applied Sciences*

Edited by

www.mdpi.com/journal/applsci

## **Women in Artificial Intelligence (AI)**


**1. Introduction**

As a reaction wheel in spacecraft, a ball bearing wheel, a magnetic bearing wheel, and an active magnetic bearing wheel have mainly been used. First, a ball bearing wheel uses a ball bearing to maintain the separation between the bearing races.  Ball bearings tend to have a lower load capacity than other kinds of rolling element bearings mainly due to the small contact area between balls and races. Also, ball bearings should be lubricated periodically with a lubricant such as oil and grease for ball bearings to operate properly [1]. Next, a magnetic bearing is used in a magnetic bearing wheel. A magnetic bearing supports a load by the magnetic levitation principle. In magnetic bearing wheels, permanent magnets are used to carry a wheel, a control system is used to hold a wheel stable, and power is used when a levitated wheel deviates from its target position. A magnetic bearing wheel also requires a back‐up bearing in case of control system or power failure and

© 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

during initial start‐up conditions. A magnetic bearing has two kinds of instabilities. One is that attractive magnets provide an unstable static force that decreases at distant distances and increases at close distances. The other is that magnetism gives rise to oscillations that may cause loss of suspension if driving forces are present [1]. Finally, in an active magnetic bearing wheel, a rotating shaft is levitated by the principle of electromagnetic suspension. A wheel is supported in an active magnetic bearing wheel without physical contact. The contactless operation of the active magnetic bearing wheels eliminates the need of lubrication of the bear‐ ing components, which allows them to operate cleanly. Moreover, it can accommodate irregu‐ larities in the mass distribution automatically, which allows it to spin around its center of mass with very low vibration, and can suppress the nutation and precession of the rotor effectively. The components of an active magnetic bearing wheel are an active magnetic bearing, a wheel, a control system, an electromagnet assembly, power amplifiers, and gap sensors. This bias current is mediated by a control system that offsets the bias current by equal but opposite perturbations of current as the rotor deviates by a small amount from its center position [1].

The active magnetic bearing wheel exhibits very lower vibration than ball bearing wheels and magnetic bearing wheels. Thus, it is a desirable reaction wheel for the spacecraft attitude control since vibration is the critical factor for the high precision spacecraft attitude control. The active magnetic bearing is the very important component among components of active magnetic bearing wheels. Due to this importance, various kinds of active magnetic bearings have been developed and their control methods have been studied (e.g., [2–10]).

In this chapter, the author presents an electromagnetic levitation system for active mag‐ netic bearing wheels. This system consists of a rotor, a shaft, a cover, and a base. Also, this system does not include a mechanism for spinning the rotor around its rotating axis. The author derives a meaningful electromagnetic force by using the singular value decom‐ position [11]. The proportional‐integral‐derivative (PID) controller is used to control the position of the rotor and regulate the two gimbal angles of the rotor. The author gives the numerical simulation and experimental results on the control of the electromagnetic levita‐ tion system.

## **2. Electromagnetic levitation system**

In **Figure 1**, the schematic of the electromagnetic levitation system developed in this chapter is given. The cover protects the rotor, and the base supports the rotor, shaft, and cover. This system can levitate the rotor up to 0.8 mm from the ground in the *z*‐axis, rotate the rotor, and gimbal the rotor within a small angle of ±0.2°.

The dynamic equations of motion of the electromagnetic levitation system are given as follows:

$$m\mathcal{Z} = -F\_z \tag{1}$$

$$I\Phi \quad = \quad T\_x \tag{2}$$

$$I\ddot{\Theta} = -T\_y \tag{3}$$

In Eqs. (1)–(3), *φ* and *θ* are the gimbal angles of the rotor in the *x*‐ and *y*‐axes, respectively, *z* is the displacement of the rotor in the *z*‐axis, *F*z is the control force in the *z*‐axis, *T*x and *T*y are the control torques applied to the rotor in the *x*‐ and *y*‐axes, respectively, *m* = 0.72 kg is the mass of the rotor, and *I* = 877.367 × 10‐6 kg m2 is the inertia of the rotor for the *x*‐ and *y*‐axes.

Let us consider the four pairs of electromagnets shown in **Figure 1**. Then, the control inputs *F*z in Eq. (1), *T*x in Eq. (2), and *T*y in Eq. (3) can be represented as follows:

$$\sum\_{t=1}^{4} F\_{st} = \quad \quad F\_z \quad + \text{ mg} \tag{4}$$

$$D\_{\epsilon} \begin{pmatrix} F\_{\epsilon 2} \ - \ \ F\_{\epsilon 4} \end{pmatrix} \quad = \quad T\_x \tag{5}$$

$$D\_{\epsilon} \begin{pmatrix} F\_{\epsilon 3} \ - \ F\_{\epsilon 1} \end{pmatrix} = \begin{array}{c} T\_y \end{array} \tag{6}$$

In Eqs. (4)–(6), *F*ei is the electromagnetic force generated by the *i*th pair of electromagnets, and *g* = 9.8 m/s2 is the acceleration of gravity. Then, Eqs. (4)–(6) can be written as follows:

$$A \, F\_e = \begin{array}{c} \boldsymbol{u} \ \text{ } + \ \boldsymbol{u}\_g \end{array} \tag{7}$$

where

during initial start‐up conditions. A magnetic bearing has two kinds of instabilities. One is that attractive magnets provide an unstable static force that decreases at distant distances and increases at close distances. The other is that magnetism gives rise to oscillations that may cause loss of suspension if driving forces are present [1]. Finally, in an active magnetic bearing wheel, a rotating shaft is levitated by the principle of electromagnetic suspension. A wheel is supported in an active magnetic bearing wheel without physical contact. The contactless operation of the active magnetic bearing wheels eliminates the need of lubrication of the bear‐ ing components, which allows them to operate cleanly. Moreover, it can accommodate irregu‐ larities in the mass distribution automatically, which allows it to spin around its center of mass with very low vibration, and can suppress the nutation and precession of the rotor effectively. The components of an active magnetic bearing wheel are an active magnetic bearing, a wheel, a control system, an electromagnet assembly, power amplifiers, and gap sensors. This bias current is mediated by a control system that offsets the bias current by equal but opposite perturbations of current as the rotor deviates by a small amount from its center position [1].

The active magnetic bearing wheel exhibits very lower vibration than ball bearing wheels and magnetic bearing wheels. Thus, it is a desirable reaction wheel for the spacecraft attitude control since vibration is the critical factor for the high precision spacecraft attitude control. The active magnetic bearing is the very important component among components of active magnetic bearing wheels. Due to this importance, various kinds of active magnetic bearings

In this chapter, the author presents an electromagnetic levitation system for active mag‐ netic bearing wheels. This system consists of a rotor, a shaft, a cover, and a base. Also, this system does not include a mechanism for spinning the rotor around its rotating axis. The author derives a meaningful electromagnetic force by using the singular value decom‐ position [11]. The proportional‐integral‐derivative (PID) controller is used to control the position of the rotor and regulate the two gimbal angles of the rotor. The author gives the numerical simulation and experimental results on the control of the electromagnetic levita‐

In **Figure 1**, the schematic of the electromagnetic levitation system developed in this chapter is given. The cover protects the rotor, and the base supports the rotor, shaft, and cover. This system can levitate the rotor up to 0.8 mm from the ground in the *z*‐axis, rotate the rotor, and

The dynamic equations of motion of the electromagnetic levitation system are given as follows:

 *mz*¨ = *Fz* (1) *I*ɸ¨ = *Tx* (2) *I* θ¨ = *Ty* (3)

have been developed and their control methods have been studied (e.g., [2–10]).

tion system.

28 Bearing Technology

**2. Electromagnetic levitation system**

gimbal the rotor within a small angle of ±0.2°.

$$A \triangleq \begin{bmatrix} 1 & 1 & 1 & 1 \\ 0 & D\_\epsilon & 0 & -D\_\epsilon \\ -D\_\epsilon & 0 & D\_\epsilon & 0 \end{bmatrix} \tag{8}$$

*F*e ≜ [ *F*e1 *F*e2 *F*e3 *F*e4] T , *u* ≜ [*F*z  *T*x  *T*y ] T , and *u*g ≜ [*mg* 0 0]T . After designing the control inputs *F*z , *T*x , and *T*y , the four electromagnetic forces *F*ei, *i* = 1,…,4 have to be determined by Eq. (7). Among solutions for Eq. (7), the minimal norm solution is derived by using the singular value decomposition [11]. Let the singular value decomposition of the matrix *A* ∈ *R*3×4 in Eq. (8) be *U*Σ*V*T and define

$$A^\* \triangleq \,^\mathsf{V}\Sigma^\* \mathcal{U}^\mathsf{T} \tag{9}$$

**Figure 1.** Schematic representation of the electromagnetic levitation system.

where *U*∈ *R*3×3 and *V* ∈ *R*4×4 are orthogonal matrices, Σ ∈ *R*3×4, and *A*+ ∈ *R*4×3 and Σ+ ∈ *R*4×3 denote the pseudoinverse matrices of the matrices A and Σ, respectively. With some calculations, we obtain the following by the singular value decomposition of the matrix *A* in Eq. (8)

$$V = \begin{bmatrix} \frac{1}{2} & 0 & -\sqrt{2} \ -1 \\ 1 & \sqrt{2} & 0 & 1 \\ 1 & 0 & \sqrt{2} & -1 \\ 1 & -\sqrt{2} & 0 & 1 \end{bmatrix} \tag{10}$$

$$
\Sigma^{+} = \frac{1}{2} \begin{bmatrix} 1 & 0 & 0 \\ 0 \frac{\sqrt{2}}{D\_{\cdot}} & 0 \\ 0 & \frac{\sqrt{2}}{D\_{\cdot}} \\ 0 & 0 & 0 \end{bmatrix} \triangleq \begin{bmatrix} \Sigma\_{1}^{-1} \\ 0\_{1 \times 3} \end{bmatrix} \tag{11}
$$

$$\begin{array}{rcl} \mathcal{U} & = & \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{array} \tag{12}$$

In Eq. (11), 01×3 implies the 1 × 3 zero matrix. Then, with *A*+ of Eq. (9), we obtain

$$F\_{\varepsilon} = A^{\star} \begin{pmatrix} u + u\_{\varepsilon} \\ \end{pmatrix} = \begin{bmatrix} \frac{1}{4} \begin{pmatrix} F\_{z} + mg \end{pmatrix} - \frac{1}{2} \frac{T\_{y}}{D\_{\varepsilon}}\\\frac{1}{4} \begin{pmatrix} F\_{z} + mg \end{pmatrix} + \frac{1}{2} \frac{T\_{z}}{D\_{\varepsilon}}\\\frac{1}{4} \begin{pmatrix} F\_{z} + mg \end{pmatrix} + \frac{1}{2} \frac{T\_{y}}{D\_{\varepsilon}}\\\frac{1}{4} \begin{pmatrix} F\_{z} + mg \end{pmatrix} + \frac{1}{2} \frac{T\_{z}}{D\_{\varepsilon}}\\\frac{1}{4} \begin{pmatrix} F\_{z} + mg \end{pmatrix} - \frac{1}{2} \frac{T\_{z}}{D\_{\varepsilon}} \end{pmatrix} \tag{13}$$

The following condition holds for any other solution *F* ^ *e* to Eq. (7) [12]

$$\|F\_{\epsilon}\|\_{2} \ll \|\hat{F}\_{\epsilon}\|\_{2} \tag{14}$$

where '‖ ⋅ ‖ 2 ' denotes the Euclidean norm of a vector (i.e., for a vector *x* ∈ *R*n, ‖*x*‖2 ≜ √ \_\_\_\_\_ ∑ *i*=1 *n xi* 2 ).

Since the two gimbal angles, φ and θ, are very small, we can approximate sin(*φ*) ≅ *φ* and sin(*θ*) ≅ *θ* by the small‐angle approximation.  (15) to (18), the system state *z*, *φ,* and *θ* can be calculated as follows:

$$\begin{array}{rcl} z & = & \frac{1}{4} \Big( \sum\_{i=1}^{4} L\_{g^i} \, - \quad \sum\_{i=1}^{4} l\_{g^i} \Big) \end{array} \tag{19}$$

$$\phi\_{\pm} = \frac{1}{4D\_{\pm}} \left[ \left( l\_{\pm^3} + l\_{\pm^4} - l\_{\pm^1} - l\_{\pm^2} \right) - \left( L\_{\pm^3} + L\_{\pm^4} - L\_{\pm^1} - L\_{\pm^2} \right) \right] \tag{20}$$

$$\Theta = \frac{1}{4D\_{\rm g}} \left[ \begin{pmatrix} l\_{g1} + \ \ l\_{g4} \ \ - \ \ l\_{g2} \ - \ \ l\_{g3} \end{pmatrix} - \begin{pmatrix} L\_{g1} + \ \ L\_{g4} \ - \ \ L\_{g2} \ - \ \ l\_{g3} \end{pmatrix} \right] \tag{21}$$

Similarly, the displacement from the bottom surface of the *i*th pair of electromagnets to the top surface of the rotor in the *z*‐axis can be calculated as follows:

$$L\_{\iota1} = \begin{array}{c} L\_{\iota1} \ + \end{array} \begin{array}{c} \begin{array}{c} \end{array} L\_{\iota1} \ = \begin{array}{c} L\_{\iota1} \ - \end{array} \begin{array}{c} z \ + \end{array} D\_{\iota} \sin(\theta) \ \not\equiv \begin{array}{c} L\_{\iota1} \ - \end{array} z \ + \end{array} D\_{\iota} \theta \tag{22}$$

$$L\_{\alpha} = \begin{array}{c} L\_{\alpha} \ + \end{array} \begin{array}{c} \delta L\_{\alpha} \ \end{array} = \begin{array}{c} L\_{\alpha} \ - \end{array} - z \ - \begin{array}{c} D\_{\epsilon} \sin \{\phi\} \end{array} \not\equiv \begin{array}{c} L\_{\alpha} \ - \end{array} - z \ - \begin{array}{c} D\_{\epsilon} \phi \end{array} \tag{23}$$

$$L\_{\odot} = \ \ L\_{\odot} + \delta L\_{\odot} = \ \ L\_{\odot} - z - \ \ D\_{\epsilon} \sin(\theta) \ \ \text{"} \ \ L\_{\odot} - z - \ \ D\_{\epsilon} \theta \tag{24}$$

$$L\_{o4} = \begin{array}{c} L\_{o4} \ + \ \delta L\_{o4} \end{array} = \begin{array}{c} L\_{o4} \ - \ \end{array} + \begin{array}{c} D\_{\epsilon} \sin \{\phi\} \ \nightharpoonup \ L\_{o4} \ - \ \end{array} + \begin{array}{c} D\_{\epsilon} \phi \end{array} \tag{25}$$

By the Maxwell's equation [13], the following equation is obtained for the control currents supplied to the coils of the four pairs of electromagnets

$$\dot{\mathbf{q}}\_{i} = \frac{2}{\pi} \frac{l\_{\dot{n}}}{\sqrt{\mu\_{\dot{0}} \mathbf{G}'}} \sqrt{\frac{F\_{\dot{n}}}{\mu\_{\dot{0}} \mathbf{G}'}}, \ \dot{\mathbf{n}} = \mathbf{1}, \ldots, \mathbf{4} \tag{26}$$

where *μ*0 = 4π × 10‐7 N/A2 is the permeability constant of free space, *n* = 240 is the number of coil turn, *i* i is the control current of the *i*th pair of electromagnets, and *G* = 50.265 × 10‐6 m2 is the cross‐sectional area of a pair of electromagnets. The author limits each control current in Eq. (26) by 1 A.

#### **3. Numerical simulation and experimental results**

where *U*∈ *R*3×3 and *V* ∈ *R*4×4 are orthogonal matrices, Σ ∈ *R*3×4, and *A*+ ∈ *R*4×3 and Σ+ ∈ *R*4×3 denote the pseudoinverse matrices of the matrices A and Σ, respectively. With some calculations, we

> 0 √ \_ 2 0 − √ \_ 2

> > 0 0 √ \_ \_ 2 *De* 0

⎤

⎥ ⎦

1 0 0 0 1 0

⎡

\_1 4

> 4

\_1 4

\_1 4

⎢

⎣

‖2 < ‖ *F* ^ *e*

Since the two gimbal angles, φ and θ, are very small, we can approximate sin(*φ*) ≅ *φ* and sin(*θ*) ≅ *θ* by the small‐angle approximation. Thus, the displacement from the bottom surface of the *i*th gap

' denotes the Euclidean norm of a vector (i.e., for a vector *x* ∈ *R*n, ‖*x*‖2 ≜ √

*g*1 = *Lg*1 + δ *Lg*1 = *Lg*1 − *z* + *Dg* sin(*θ*) − *Dg* sin(*φ*) ≅  *Lg*1 − *z* + *Dg θ* − *Dg φ* (15)

*g*2 = *Lg*2 + δ *Lg*2 = *Lg*2 − *z* − *Dg* sin(*θ*) − *Dg* sin(*φ*) ≅  *Lg*2 − *z* − *Dg θ* − *Dg φ* (16)

*g*3 = *Lg*3 + δ *Lg*3 = *Lg*3 − *z* − *Dg* sin(*θ*) + *Dg* sin(*φ*) ≅  *Lg*2 − *z* − *Dg θ* + *Dg φ* (17)

*g*4 = *Lg*4 + δ *Lg*4 = *Lg*4 − *z* + *Dg* sin(*θ*) + *Dg* sin(*φ*) ≅  *Lg*4 − *z* + *Dg θ* + *Dg φ* (18)

Then, from Eqs. (15) to (18), the system state *z*, *φ,* and *θ* can be calculated as follows:

*g*1 − *l*

≜ [ Σ1 −1

(*Fz* + *mg*) − \_1

(*Fz* + *mg*) + \_1

(*Fz* + *mg*) + \_1

(*Fz* + *mg*) − \_1

^

\_1

− √ \_ 2 0 √ \_ 2 0

⎤

⎥ ⎦ (10)

(13)

\_\_\_\_\_ ∑ *i*=1 *n xi* 2 ).

0 1×3 ] (11)

0 0 1 ] (12)

of Eq. (9), we obtain

⎤

⎥

⎦

 ‖2 (14)

2 *T* \_*y De*

2 *T*\_*x De*

2 *T* \_*y De*

2 *T*\_*x De*

*e* to Eq. (7) [12]

*gi*) (19)

*g*2) − (*Lg*3 + *Lg*4 − *Lg*1 − *Lg*2)] (20)

obtain the following by the singular value decomposition of the matrix *A* in Eq. (8)

⎡

⎢ ⎣

2

⎡

0 √ \_ \_ 2 *De* 0 0

[

⎢ ⎣

2

(*u* + *ug*) =

sensor to the top surface of the rotor in the *z*‐axis can be calculated as follows:

*g*4 − *l*

*V* = \_\_1

Σ+ = \_\_1

*U* =

*Fe* = *A*+

where '‖ ⋅ ‖

30 Bearing Technology

*l*

*l*

*l*

*l*

2

*z* = \_\_1

*φ* = \_\_\_1

In Eq. (11), 01×3 implies the 1 × 3 zero matrix. Then, with *A*+

The following condition holds for any other solution *F*

‖*Fe*

4 (∑ *i*=1 4

> 4 *Dg* [(*l g*3 + *l*

 *Lgi* − ∑ *i*=1 4 *l* The author gives the numerical simulation and experimental results on the control of the electromagnetic levitation system in this section.

The author uses the following discretized PID controller to control the displacement of the rotor in Eq. (19) and the two gimbal angles of the rotor in Eqs. (20) and (21).

$$\begin{array}{rcl} \stackrel{\circ}{\text{rotor in Eq. (19) and the two gimbal angles of the rotor in Eqs. (20) and (21)}.}\\ \iota(t) &=& K\_p \underline{x}\_i(t) + K\_i T\_s \sum\_{i=0}^t \underline{x}\_i(i) + K\_d \frac{\left(\underline{x}\_i(t) - \underline{x}\_i(t-1)\right)}{T\_s} \end{array} \tag{27}$$

where *T*s is the sampling time and given by *T*s = 1 ms, *x*e ≜ *x*t − *x* is the error between the sys‐ tem state *x* ≜ [*z φ θ*]T and the target system state *x*t ≜ [*z*t  *φ*t  *θ*t ]T and *K*p ≜ diag[*K*pz, *K*p*φ*, *K*pθ], *K*i ≜ diag[*K*iz, *K*i*φ*, *K*iθ], and *K*d ≜ diag[*K*dz, *K*d*φ*, *K*dθ] denote the 3 × 3 diagonal positive definite matrices.

The target position and target gimbal angles of the rotor are set to be *z*t = 0.3 mm and *φ*t = 0° and *θ*t = 0°, respectively. The author initially decides the feedback gains of the PID controller in Eq. (27) that can achieve the control objective by adopting the well‐known Ziegler‐Nichols method [14] and then finely tunes the feedback gains of the PID controller in Eq. (27) by an experiment. Also, the antiwindup compensator is used to make the overshot as small as pos‐ sible. According to the Ziegler‐Nichols method [14], first, *K*i = *K*d = diag[0, 0, 0] are set, and the proportional gain *K*p is then increased until the system just oscillates. The proportional gain is then multiplied by 0.6, and the integral and derivative gains are calculated as *K*p = 0.6 *K*m, *K*i = *K*p(*ω*m/π), and *K*d = *K*p(0.25 π/*ω*m) where *K*m ≜ diag[*K*mz, *K*m*φ*, *K*mθ] denotes the 3 × 3 diagonal positive definite matrix with the gain elements at which the proportional system oscillates, and *ω*m ≜ diag[*ω*mz, *ω*m*φ*, *ω*mθ] denotes the 3 × 3 diagonal positive definite matrix with the oscil‐ lation frequency elements. As a result, the feedback gains of the PID controller in Eq. (27) are chosen as follows: *K*pz = 1000, *K*p*φ* = *K*pθ = 20, *K*iz = 5000, *K*i*φ* = *K*iθ = 8.33, *K*dz = 100, and *K*d*φ* = *K*dθ = 0.013. Consequently, we see that the dominant feedback gains in this PID controller are *K*p and *K*i, and thus, one can obtain *K*mz ≅ 1666.67, *K*m*φ* = *K*mθ ≅ 33.33, *ω*m ≅ 15.71 rad/s and *ω*m*φ* = *ω*mθ ≅ 1.31 rad/s.

The control flow diagram of the system is shown in **Figure 2**. After we measure the displace‐ ments from the bottom surfaces of the four gap sensors to the top surface of the rotor, we cal‐ culate the displacement of the rotor and the two gimbal angles of the rotor by Eqs. (19)–(21), respectively. The control input *u* ≜ [*F*z  *T*x  *T*y ] T is made by the PID controller in Eq. (27). Then, we calculate the four electromagnetic forces and the displacements from the bottom surfaces of the four pairs of electromagnets to the top surface of the rotor by Eq. (13) and the Eq. (22)–(25), respectively. After we calculate the control currents by Eq. (26), they pass through the current limiters and are supplied to the coils of the four pairs of electromagnets by the power electronics.

**Figure 2.** Control flow diagram of the electromagnetic levitation system.

In the numerical simulation, it will be demonstrated that the electromagnetic force *F*e of Eq. (13) satisfies the condition of Eq. (14) with respect to the following electromagnetic force *F* ^ *e* , which is another solution to Eq. (7)

$$\begin{aligned} \hat{F}\_{\epsilon} \text{ which is another solution to Eq. (7)}\\ \hat{F}\_{\epsilon} &= \begin{bmatrix} mg\\ \frac{1}{2} F\_z + \frac{1}{2} \frac{\left(T\_z - T\_y\right)}{D\_\epsilon} - \frac{1}{2} mg\\ \frac{T\_y}{D\_\epsilon} + mg\\ \frac{1}{2} F\_z - \frac{1}{2} \frac{\left(T\_z + T\_y\right)}{D\_\epsilon} - \frac{1}{2} mg \end{bmatrix} \triangleq \begin{bmatrix} \hat{F}\_{\epsilon 1} & \hat{F}\_{\epsilon 2} & \hat{F}\_{\epsilon 3} \end{bmatrix}^{\mathrm{T}} \tag{28}$$

With the initial system state given by *z*ini = 0 mm, *z*˙*ini* = 0 mm/s, the values are *φ*ini = *θ*ini = 0° and *φ*˙ *ini* = *θ*˙ *ini* = 0°/s. In **Figures 3**–**5**, the numerical simulation results on the control of the electromagnetic levitation system using *F*e of Eq. (13) and *F* ^ *e* of Eq. (28) are shown. As shown in **Figure 3**, the rotor reaches the target position *z*t = 0.3 mm, and the control forces using *F*e of Eq. (13) and *F* ^ *e* of Eq. (28) show the same behaviors, respectively. In **Figure 4**, we see that each control current using *F*e of Eq. (13) reaches the same value, the control currents of *i* 1 and *i* 3 using *F* ^ *e*1 and *F* ^ *e*3 of Eq. (28), respectively, reach the same value, and the control currents of *i* 2 and *i* 4 using *F* ^ *e*2 and *F* ^ *e*4 of Eq. (28), respectively, reach the same value. In **Figure 5**, we see that the Euclidean norm of *F*e of Eq. (13) is about 3.162 times smaller than that of *F* ^ *e* of Eq. (28) in the steady‐state region. Therefore, the numerical simulation results shown in **Figure 5** illustrate that *F*e of Eq. (13) satisfies the condition of Eq. (14) with respect to *F* ^ *e* of Eq. (28).

and *ω*m ≜ diag[*ω*mz, *ω*m*φ*, *ω*mθ] denotes the 3 × 3 diagonal positive definite matrix with the oscil‐ lation frequency elements. As a result, the feedback gains of the PID controller in Eq. (27) are chosen as follows: *K*pz = 1000, *K*p*φ* = *K*pθ = 20, *K*iz = 5000, *K*i*φ* = *K*iθ = 8.33, *K*dz = 100, and *K*d*φ* = *K*dθ = 0.013. Consequently, we see that the dominant feedback gains in this PID controller are *K*p

The control flow diagram of the system is shown in **Figure 2**. After we measure the displace‐ ments from the bottom surfaces of the four gap sensors to the top surface of the rotor, we cal‐ culate the displacement of the rotor and the two gimbal angles of the rotor by Eqs. (19)–(21),

we calculate the four electromagnetic forces and the displacements from the bottom surfaces of the four pairs of electromagnets to the top surface of the rotor by Eq. (13) and the Eq. (22)–(25), respectively. After we calculate the control currents by Eq. (26), they pass through the current limiters and are supplied to the coils of the four pairs of electromagnets by the power electronics.

In the numerical simulation, it will be demonstrated that the electromagnetic force *F*e

 – \_1 2 *mg*

 – \_1 2 *mg*

+*mg*

*mg*

(*Tx* − *T* \_*y*)

*De*

(*Tx* + *T* \_*y*)

*De*

\_1

*T* \_*y De*

Eq. (13) satisfies the condition of Eq. (14) with respect to the following electromagnetic force

⎤

⎥

≜ [*F* ^ *e*1 *F* ^ *e*2 *F* ^ *e*3 *F* ^ *e*4] T

⎦

 *T*x  *T*y ]

and thus, one can obtain *K*mz ≅ 1666.67, *K*m*φ* = *K*mθ ≅ 33.33, *ω*m ≅ 15.71 rad/s and *ω*m*φ* =

T is made by the PID controller in Eq. (27). Then,

of

(28)

and *K*i,

32 Bearing Technology

*F* ^ *e*

*F*

*ω*mθ ≅ 1.31 rad/s.

respectively. The control input *u* ≜ [*F*z

, which is another solution to Eq. (7)

⎡

**Figure 2.** Control flow diagram of the electromagnetic levitation system.

⎢

2 *Fz* + \_1 2

⎣

\_1 2 *Fz* – \_1 2

^ *e* = The experimental results on the control of the system using *F*ei, *i* = 1,…,4 of Eq. (13) are shown in **Figures 6**–**9**. As shown in **Figures 6** and **7**, the PID controller successfully levitates the rotor at the target position *z*t = 0.3 mm with well regulating the two gimbal angles. The trajectories of control currents are shown in **Figure 8**. And the Euclidean norm of *F*e of Eq. (13) is shown in **Figure 9**. As shown in **Figures 4**, **5**, **8**, and **9**, the trajectories of control currents and Euclidean norm of *F*e of Eq. (13) obtained by the experiment move around the values obtained by the numerical simulation.

**Figure 3.** Time histories of the position of the rotor and the control force, which are obtained by the numerical simulation.

**Figure 4.** Time histories of the control currents, which are obtained by the numerical simulation.

**Figure 5.** Time histories of the Euclidean norm of the electromagnetic force, which are obtained by the numerical simulation.

**Figure 6.** Time histories of the system state of the electromagnetic levitation system using Fe of Eq. (13), which are obtained by the experiment.

**Figure 7.** Time histories of the control inputs of the electromagnetic levitation system using Fe of Eq. (13), which are obtained by the experiment.

**Figure 5.** Time histories of the Euclidean norm of the electromagnetic force, which are obtained by the numerical

**Figure 4.** Time histories of the control currents, which are obtained by the numerical simulation.

simulation.

34 Bearing Technology

**Figure 8.** Time histories of the control currents of the electromagnetic levitation system using Fe of Eq. (13), which are obtained by the experiment.

**Figure 9.** Time histories of the Euclidean norm of Fe of Eq. (13), which are obtained by the experiment.

In order to simulate an external disturbance, a human hand presses down hard on the rotor to the ground at about 20 s, and it is removed from the rotor momentarily. In **Figure 6**, we see that, after we remove the external disturbance from the rotor, the system state becomes a steady state within 2 s. Also, in **Figure 8**, we see that, as the external disturbance applied to the rotor increases, each control current increases to resist the external disturbance. It should be remarked that the operating parameters like an applied load and a speed of the rotor may influence on the design of electromagnetic levitation system because these parameters make an impact on the system dynamics.


**2. Finite element simulation on homogeneous SP model**

### **2.1. Homogeneous SP model**

### *2.1.1. Introduction of SOLID164 element*

The explicit dynamic analysis of ANSYS/LS-DYNA program provides a rich element library, including 3D bar element (LINK160), 3D beam element (BEAM161), thin shell element (SHELL163), solid element (SOLID164), spring damping element (COMBI165), mass element (MASS166), cable element (Link167), and ten-node tetrahedral element (Tet-Solid168). Each explicit dynamic element has a corresponded linear displacement function, and the default is set to a single point of integration (one of the reduced integrals). It has been proved that the explicit dynamic element with linear displacement function and the element integration can be used for solving the nonlinear problems effectively, such as the cases of large deformation and material's failure [34]. The element utilized in simulation of SP is the SOLID164 element, which is a 3D explicit solid element consisting of eight nodes. Each node in SOLID164 has the degrees of freedom in *x*, *y*, *z* direction of translation, velocity, and acceleration. The geometry, node locations, and the coordinate system for SOLID164 element are shown in **Figure 1** [35]. This element is only used in the dynamic explicit analysis, which supports all licensed nonlinear characteristics. By default, SOLID164 is the single point integral with the viscous hourglass control to speed up the calculation.

points. So, using the simplified integration of elements can save the data storage capacity and reduce the number of calculations greatly, for instance, the single point integral in SOLID164 element, that is, an element has only one integration point, which is in the center of element. The utilization of single points can save a lot of calculation time, but may lead to the hourglass phenomenon. Mesh deformation with an hourglass effect is called the hourglass phenomenon, and in this situation, the typical feature is that the nonrigid element grid is distorted irregularly, resulting in a mathematically stable state but a physically impossible state. The presence of hourglass phenomenon will distort the results of the solution and even the solution cannot be carried out. Therefore, when using the simplified integration of elements, the hourglass phenomenon should be controlled. In ANSYS/LS-DYNA finite element analysis process, if the result of hourglass energy calculated by the model is less than 10% of the total

Finite Element Dynamic Analysis on Residual Stress Distribution of Titanium Alloy and Titanium…

The choice of material model is not only related to the success of simulation, but also directly related to the rationality and reliability of calculated results. The process of SP causes a high plastic deformation in the surface layer, so Cowper-Symbols in ANSYS is adopted to achieve 3D finite element dynamic analysis [34]. Cowper-Symbols model is a piecewise linear plasticity model, and the yield stress can be obtained via the strain rate, which is shown in Eq. (1).

> ) = *σy*(*εeff p* )[1 + ( *ε*̇ *eff p* \_\_\_ *C* ) \_\_1 *p*

original yield stress while the strain rate has not been considered. In this model during simulation, the input data include totally: the density of material *ρ*, the elastic modulus of material *E*, Poisson's ratio *ν*, the yield stress, the tangent modulus, the strain rate parameters *P* and *C*, and the true stress–strain curve. If the load curve is used, the yield stress and tangent modulus

Before establishing SP model for (TiB+TiC)/Ti-6Al-4V, the 3D model for homogeneous matrix Ti-6Al-4V should be established firstly. For reducing the number of elements and the calculation time, the symmetry of SP sample and shot balls should be considered, and a 1/2 model can be set up. The 3D homogeneous model established in this work is shown in **Figure 2** [36], including the top four-layer of shot balls, and the bottom of peened target. The dimension of

of the intensive impacting in the near surface layer, a mesh refinement is adopted and each mesh depth is 0.02 mm. SOLID164 dynamic analysis element is chosen for meshing element,

, in which *R* is the average radius of shot balls. Because

*ε*̇ represents the effective strain rate, *P* and *C* are the parameters for strain rate, and *σy*(*εeff*

] (1)

http://dx.doi.org/10.5772/intechopen.73120

27

*p* ) is the

energy, this simulation result and the model can be identified as credible.

*p* , *ε*̇ *eff p*

are ignored. If *P* and *C* are set to 0, the strain rate effect is omitted.

*2.1.2. Piecewise linear plasticity model*

*σy*(*εeff*

*2.1.3. Establishment of homogeneous SP model*

peened target is 12*R* × 6*R* × 2.1 mm3

and total mesh number is 120,000.

In the process of explicit dynamic analysis, the processing of element integral is the most time-consumption, and the processing time by CPU is proportional to the number of integral

**Figure 1.** Geometry properties of SOLID164 [35].

points. So, using the simplified integration of elements can save the data storage capacity and reduce the number of calculations greatly, for instance, the single point integral in SOLID164 element, that is, an element has only one integration point, which is in the center of element. The utilization of single points can save a lot of calculation time, but may lead to the hourglass phenomenon. Mesh deformation with an hourglass effect is called the hourglass phenomenon, and in this situation, the typical feature is that the nonrigid element grid is distorted irregularly, resulting in a mathematically stable state but a physically impossible state. The presence of hourglass phenomenon will distort the results of the solution and even the solution cannot be carried out. Therefore, when using the simplified integration of elements, the hourglass phenomenon should be controlled. In ANSYS/LS-DYNA finite element analysis process, if the result of hourglass energy calculated by the model is less than 10% of the total energy, this simulation result and the model can be identified as credible.

### *2.1.2. Piecewise linear plasticity model*

**2. Finite element simulation on homogeneous SP model**

26 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

The explicit dynamic analysis of ANSYS/LS-DYNA program provides a rich element library, including 3D bar element (LINK160), 3D beam element (BEAM161), thin shell element (SHELL163), solid element (SOLID164), spring damping element (COMBI165), mass element (MASS166), cable element (Link167), and ten-node tetrahedral element (Tet-Solid168). Each explicit dynamic element has a corresponded linear displacement function, and the default is set to a single point of integration (one of the reduced integrals). It has been proved that the explicit dynamic element with linear displacement function and the element integration can be used for solving the nonlinear problems effectively, such as the cases of large deformation and material's failure [34]. The element utilized in simulation of SP is the SOLID164 element, which is a 3D explicit solid element consisting of eight nodes. Each node in SOLID164 has the degrees of freedom in *x*, *y*, *z* direction of translation, velocity, and acceleration. The geometry, node locations, and the coordinate system for SOLID164 element are shown in **Figure 1** [35]. This element is only used in the dynamic explicit analysis, which supports all licensed nonlinear characteristics. By default, SOLID164 is the single point integral with the viscous hourglass control to speed up the

In the process of explicit dynamic analysis, the processing of element integral is the most time-consumption, and the processing time by CPU is proportional to the number of integral

**2.1. Homogeneous SP model**

calculation.

*2.1.1. Introduction of SOLID164 element*

**Figure 1.** Geometry properties of SOLID164 [35].

The choice of material model is not only related to the success of simulation, but also directly related to the rationality and reliability of calculated results. The process of SP causes a high plastic deformation in the surface layer, so Cowper-Symbols in ANSYS is adopted to achieve 3D finite element dynamic analysis [34]. Cowper-Symbols model is a piecewise linear plasticity model, and the yield stress can be obtained via the strain rate, which is shown in Eq. (1).

$$
\sigma\_y \left( \varepsilon\_{gt'}^p, \dot{\varepsilon}\_{gt}^p \right) = \sigma\_y \{ \varepsilon\_{gt}^p \} \left[ 1 + \left( \frac{\dot{\varepsilon}\_{gt}^p}{C} \right)^{\frac{1}{p}} \right] \tag{1}
$$

*ε*̇ represents the effective strain rate, *P* and *C* are the parameters for strain rate, and *σy*(*εeff p* ) is the original yield stress while the strain rate has not been considered. In this model during simulation, the input data include totally: the density of material *ρ*, the elastic modulus of material *E*, Poisson's ratio *ν*, the yield stress, the tangent modulus, the strain rate parameters *P* and *C*, and the true stress–strain curve. If the load curve is used, the yield stress and tangent modulus are ignored. If *P* and *C* are set to 0, the strain rate effect is omitted.

### *2.1.3. Establishment of homogeneous SP model*

Before establishing SP model for (TiB+TiC)/Ti-6Al-4V, the 3D model for homogeneous matrix Ti-6Al-4V should be established firstly. For reducing the number of elements and the calculation time, the symmetry of SP sample and shot balls should be considered, and a 1/2 model can be set up. The 3D homogeneous model established in this work is shown in **Figure 2** [36], including the top four-layer of shot balls, and the bottom of peened target. The dimension of peened target is 12*R* × 6*R* × 2.1 mm3 , in which *R* is the average radius of shot balls. Because of the intensive impacting in the near surface layer, a mesh refinement is adopted and each mesh depth is 0.02 mm. SOLID164 dynamic analysis element is chosen for meshing element, and total mesh number is 120,000.

*σt* = *σε*(1 + *εε*) (2)

*εt* = ln(1 + *εε*) (3)

[37]. The true stress–strain curve is used to provide the corresponding deformation parameters for the shot peened material in the finite element model. Because of the existence of reinforcements in composite and the influences of which on the matrix's residual stress distribution, the strain–stress curves with a single strain rate (10−3 s−1) are utilized during simulation. The dynamic stress–strain curves become flatter with increasing strain rate, which reduces the material strain hardening [38]. Thus, considering the variation of strain rates, based on the references

The true stress–strain curves of Ti-6Al-4V and (TiB+TiC)/Ti-6Al-4V can be obtained from the tensile tests. The dimension of tensile specimens with a gauge section of 4 × 1.8 × 18 mm is shown in **Figure 4(a)**, and the axes direction is parallel to the hot-forging direction. The tensile tests are performed using a Zwick T1-FR020TN test machine in air, and the initial strain rates is 10−3 s−1. The true stress–strain curves are shown in **Figure 4(b)**, which are implemented in the simulation process as the elastic–plastic deformation curves. In addition, the typical

In order to show the deformation of material obviously after SP, the coverage rate of 300% is chosen and the shot velocity is 100 m/s. The simulation results of *x*-direction principal stress (σxx) are shown in **Figure 5**. From the **Figure 5(a)**, it can be found that after SP, there is a compressive stress field with a certain depth in the plastic deformation area below the shot balls, and the

mechanical parameters of shot balls, target, and reinforcements are shown in **Table 1**.

is the true strain, and *εε*

Finite Element Dynamic Analysis on Residual Stress Distribution of Titanium Alloy and Titanium…

is the engineering strain

29

http://dx.doi.org/10.5772/intechopen.73120

is the engineering stress, *εt*

[39, 40] and Eq. (1), the parameters *C* and *P* are set as *C* = 1300 and *P* = 5.

*σt*

is the true stress, *σε*

**2.2. Simulated results and discussion**

*2.2.1. Residual stress distribution on surface*

**Figure 3.** Coverage rate of SP with different impacts [36].

**Figure 2.** 3D SP model for homogeneous materials [36].

For improving the computational efficiency, the mesh size gradually increases at the model boundary and the lower half. In order to avoid the influence of reflected stress wave within the target on the distribution of residual stress during the process of impacting, nonreflective boundary conditions are implemented on the bottom and flank of the target. Because the XOZ plane is symmetrical, the symmetrical boundary conditions are applied. Because a small size model (12*R* × 6*R* × 2.1 mm3 ) may lead to the nonreal oscillation while simulating the process of impacting, in order to eliminate the effect of nonreal oscillation, the alpha damping constraint is applied to the model in the dynamic analysis. The alpha damping is a damping coefficient proportional to the mass and is very effective for the low frequency oscillations.

In 3D SP model, the top four-layers are shot balls made by case steel, and the hardness and strength of matrix is smaller than that of the case steel balls. The deformability of shot balls is very weak, which hardly affects the results of residual stress distribution.  During the explicit dynamics analysis, the degrees of freedom of all nodes in the rigid body are coupled to the mass center of rigid body, thus greatly reducing the computational time of the explicit analysis. Moreover, the corresponding mechanical parameters are given to the rigid body center to describe the dynamic characteristics. In addition, the coverage rate is very important to the distribution of residual stresses after SP. In order to keep the simulation accuracy, the coverage rate is defined and the schematic is shown in **Figure 3** [36].  On a single impact, the coverage rate is about 25%, and after four impacts, it approaches 100%. If coverage rate increases to 200%, eight impacts should be carried out in turn.

Because the piecewise linear model is adopted during simulation, the true stress–strain curve is needed. Usually, the stress–strain curve obtained by experiment is the engineering stress– strain curve, which is needed to be transformed into the true stress–strain curve, and the formula of transformation is shown in Eqs. (2) and (3).

Finite Element Dynamic Analysis on Residual Stress Distribution of Titanium Alloy and Titanium… http://dx.doi.org/10.5772/intechopen.73120 29

$$
\sigma\_i = \sigma\_i (1 + \varepsilon\_i) \tag{2}
$$

$$
\varepsilon\_t = \ln(1 + \varepsilon\_t) \tag{3}
$$

*σt* is the true stress, *σε* is the engineering stress, *εt* is the true strain, and *εε* is the engineering strain [37]. The true stress–strain curve is used to provide the corresponding deformation parameters for the shot peened material in the finite element model. Because of the existence of reinforcements in composite and the influences of which on the matrix's residual stress distribution, the strain–stress curves with a single strain rate (10−3 s−1) are utilized during simulation. The dynamic stress–strain curves become flatter with increasing strain rate, which reduces the material strain hardening [38]. Thus, considering the variation of strain rates, based on the references [39, 40] and Eq. (1), the parameters *C* and *P* are set as *C* = 1300 and *P* = 5.

The true stress–strain curves of Ti-6Al-4V and (TiB+TiC)/Ti-6Al-4V can be obtained from the tensile tests. The dimension of tensile specimens with a gauge section of 4 × 1.8 × 18 mm is shown in **Figure 4(a)**, and the axes direction is parallel to the hot-forging direction. The tensile tests are performed using a Zwick T1-FR020TN test machine in air, and the initial strain rates is 10−3 s−1. The true stress–strain curves are shown in **Figure 4(b)**, which are implemented in the simulation process as the elastic–plastic deformation curves. In addition, the typical mechanical parameters of shot balls, target, and reinforcements are shown in **Table 1**.

### **2.2. Simulated results and discussion**

For improving the computational efficiency, the mesh size gradually increases at the model boundary and the lower half. In order to avoid the influence of reflected stress wave within the target on the distribution of residual stress during the process of impacting, nonreflective boundary conditions are implemented on the bottom and flank of the target. Because the XOZ plane is symmetrical, the symmetrical boundary conditions are applied. Because a small size

impacting, in order to eliminate the effect of nonreal oscillation, the alpha damping constraint is applied to the model in the dynamic analysis. The alpha damping is a damping coefficient

In 3D SP model, the top four-layers are shot balls made by case steel, and the hardness and strength of matrix is smaller than that of the case steel balls. The deformability of shot balls is very weak, which hardly affects the results of residual stress distribution.  During the explicit dynamics analysis, the degrees of freedom of all nodes in the rigid body are coupled to the mass center of rigid body, thus greatly reducing the computational time of the explicit analysis. Moreover, the corresponding mechanical parameters are given to the rigid body center to describe the dynamic characteristics. In addition, the coverage rate is very important to the distribution of residual stresses after SP. In order to keep the simulation accuracy, the coverage rate is defined and the schematic is shown in **Figure 3** [36].  On a single impact, the coverage rate is about 25%, and after four impacts, it approaches 100%. If coverage

Because the piecewise linear model is adopted during simulation, the true stress–strain curve is needed. Usually, the stress–strain curve obtained by experiment is the engineering stress– strain curve, which is needed to be transformed into the true stress–strain curve, and the

proportional to the mass and is very effective for the low frequency oscillations.

rate increases to 200%, eight impacts should be carried out in turn.

formula of transformation is shown in Eqs. (2) and (3).

) may lead to the nonreal oscillation while simulating the process of

model (12*R* × 6*R* × 2.1 mm3

**Figure 2.** 3D SP model for homogeneous materials [36].

28 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

### *2.2.1. Residual stress distribution on surface*

In order to show the deformation of material obviously after SP, the coverage rate of 300% is chosen and the shot velocity is 100 m/s. The simulation results of *x*-direction principal stress (σxx) are shown in **Figure 5**. From the **Figure 5(a)**, it can be found that after SP, there is a compressive stress field with a certain depth in the plastic deformation area below the shot balls, and the

**Figure 3.** Coverage rate of SP with different impacts [36].

**Figure 4.** (a) Specimen size for room temperature tests (unit: mm); (b) stress-strain curves of the matrix and composite (strain rate: 10−3 s−1).

on the simulation results. So, the tensile stress field of surface boundary has no effect on the residual stress distribution in the plastic deformation zone. Moreover, in the actual SP experiment, the whole surface of material is subjected to SP process, and the uniform CRS can be

Finite Element Dynamic Analysis on Residual Stress Distribution of Titanium Alloy and Titanium…

http://dx.doi.org/10.5772/intechopen.73120

31

The effect of SP on the residual stress distribution is mainly described by the residual stress variation with the increase of depth and the four characteristic parameters (the surface residual stress, the max residual stress, the depth of residual stress layer, and the depth of max residual stress). Usually, in experiment, after SP, the surface layer of samples are subjected to electrochemical etching and removed layer by layer (each layer is about 15–20 μm), and then the stress value of each layer can be measured by XRD method, and the curves of residual stress variation along the depth are obtained. In order to avoid measurement error, the irradia-

simulation, for obtaining the distribution of residual stress along the depth, the stress value of all nodes at a certain depth along Z direction are selected and averaged, which can represent the average stress in a certain depth. Using this method to average the stress values of all nodes, the curve of residual stress distribution with the increase of depth can be obtained.

After SP, the surface coverage rate refers to the ratio of the area occupied by the shot craters to the area of surface required for SP. During the process of SP, the coverage rate is usually required for reaching or exceeding 100%. Moreover, in the experiment, the coverage rate more than 100% can be expressed as a multiple of the time required for a full coverage rate of 100%, for example, the coverage rate of 200% means that the SP time is two times of full coverage.

, so actually, the measured stress by

. In this numerical

obtained. The surface tensile stress on the boundary disappears at the same time.

**Figure 5.** Simulation results of residual stress field of SP (σxx), (a) 3D result; (b) 2D result of surface.

experiment is the average stress under the statistical of surface area in 1 mm2

tion area of X-ray on the surface is generally about 1 mm2

*2.2.2. Influence of coverage rate on residual stress distribution*

**Table 1.** Typical mechanical parameters of shot balls and target.

distribution of CRS in the plastic deformation zone is uniform, which is about −900 to −1000 MPa. However, in the 2D result of surface in **Figure 5(b)**, the distribution of residual stress is not uniform, and the max CRS is about −1453 MPa, which results from the stress concentration on the impact point of shot ball. In addition, there is tensile residual stress in some area between two craters. As the impacting time of shot ball on material's surface is different, the size of crater on the surface is not consistent, and the area of craters formed firstly is significantly reduced by the extrusion of follow-up craters. In these extruded area, the tensile residual stress is formed, but after repeated SP, the distribution of residual stress is gradually tended to uniform.

Due to the repeated impact and extrusion of shot balls, the plastic deformation appears and many craters are formed on surface. The formation of craters results in the increase of surface area and the appearance of tensile residual stress, which is also the reason that the max CRS does not appear on surface. On surface, the residual stress includes both the compressive stress formed by the plastic deformation and the tensile stress formed by the increase of surface area. With the increase of depth, the tensile residual stress decreases gradually but the CRS is gradually increased, then the max CRS is reached at a certain depth. At last, due to the reduction of plastic deformation, the CRS is gradually reduced to the level of stress before SP. In **Figure 5(b)**, it can be found that there is tensile stress on the surface boundary. The formation of such tensile stress field is mainly due to the limitation of finite element model size (only 12*R* × 6*R*), and the SP area only 8*R* × 4*R* in order to avoid the effect of model boundary Finite Element Dynamic Analysis on Residual Stress Distribution of Titanium Alloy and Titanium… http://dx.doi.org/10.5772/intechopen.73120 31

**Figure 5.** Simulation results of residual stress field of SP (σxx), (a) 3D result; (b) 2D result of surface.

on the simulation results. So, the tensile stress field of surface boundary has no effect on the residual stress distribution in the plastic deformation zone. Moreover, in the actual SP experiment, the whole surface of material is subjected to SP process, and the uniform CRS can be obtained. The surface tensile stress on the boundary disappears at the same time.

The effect of SP on the residual stress distribution is mainly described by the residual stress variation with the increase of depth and the four characteristic parameters (the surface residual stress, the max residual stress, the depth of residual stress layer, and the depth of max residual stress). Usually, in experiment, after SP, the surface layer of samples are subjected to electrochemical etching and removed layer by layer (each layer is about 15–20 μm), and then the stress value of each layer can be measured by XRD method, and the curves of residual stress variation along the depth are obtained. In order to avoid measurement error, the irradiation area of X-ray on the surface is generally about 1 mm2 , so actually, the measured stress by experiment is the average stress under the statistical of surface area in 1 mm2 . In this numerical simulation, for obtaining the distribution of residual stress along the depth, the stress value of all nodes at a certain depth along Z direction are selected and averaged, which can represent the average stress in a certain depth. Using this method to average the stress values of all nodes, the curve of residual stress distribution with the increase of depth can be obtained.

### *2.2.2. Influence of coverage rate on residual stress distribution*

distribution of CRS in the plastic deformation zone is uniform, which is about −900 to −1000 MPa. However, in the 2D result of surface in **Figure 5(b)**, the distribution of residual stress is not uniform, and the max CRS is about −1453 MPa, which results from the stress concentration on the impact point of shot ball. In addition, there is tensile residual stress in some area between two craters. As the impacting time of shot ball on material's surface is different, the size of crater on the surface is not consistent, and the area of craters formed firstly is significantly reduced by the extrusion of follow-up craters. In these extruded area, the tensile residual stress is formed, but

**Figure 4.** (a) Specimen size for room temperature tests (unit: mm); (b) stress-strain curves of the matrix and composite

**Shot balls Ti-6Al-4V TiB TiC**

Due to the repeated impact and extrusion of shot balls, the plastic deformation appears and many craters are formed on surface. The formation of craters results in the increase of surface area and the appearance of tensile residual stress, which is also the reason that the max CRS does not appear on surface. On surface, the residual stress includes both the compressive stress formed by the plastic deformation and the tensile stress formed by the increase of surface area. With the increase of depth, the tensile residual stress decreases gradually but the CRS is gradually increased, then the max CRS is reached at a certain depth. At last, due to the reduction of plastic deformation, the CRS is gradually reduced to the level of stress before SP. In **Figure 5(b)**, it can be found that there is tensile stress on the surface boundary. The formation of such tensile stress field is mainly due to the limitation of finite element model size (only 12*R* × 6*R*), and the SP area only 8*R* × 4*R* in order to avoid the effect of model boundary

after repeated SP, the distribution of residual stress is gradually tended to uniform.

Density ρ (Kg·m−3) 7800 4500 4520 4930 Poisson's ratio ν 0.30 0.34 0.14–0.15 0.18–0.2 Young's modulus *E* (GPa) 210 112 371–485 400–460 Constraint of the part Rigid Elastic–plastic Elastic–plastic Elastic–plastic

**Table 1.**

After SP, the surface coverage rate refers to the ratio of the area occupied by the shot craters to the area of surface required for SP. During the process of SP, the coverage rate is usually required for reaching or exceeding 100%. Moreover, in the experiment, the coverage rate more than 100% can be expressed as a multiple of the time required for a full coverage rate of 100%, for example, the coverage rate of 200% means that the SP time is two times of full coverage. In order to simplify the model, the collision between two shot balls are not considered, and SP position of each ball can be precisely controlled. Therefore, the different coverage rate can be simulated by the multi-layer shot balls and the coverage increases linearity with the number of layers. In this work, the characteristics of 100% coverage rate (4-layer shot balls), 200% coverage rate (8-layer shot balls), and 300% coverage rate (12-layer shot balls) are simulated. The shot velocity are set as 50 and 100 m/s based on experiment (mentioned in Section 2.2.4), and the average radius of shot balls are chosen as 0.15, 0.3, and 0.6 mm. The influence of coverage rate, radius of shot balls (*r*), and shot velocity (*v*) on residual stress distribution are investigated and discussed.

**Figure 6** shows the distribution of residual stress along depth with different coverage rates. While the radius and speed of shot balls are constant, the similar variation trend of residual stress can be obtained. With the increase of coverage rate, the CRS and the depth of stress layer are improved. In addition, the depth of max CRS decreases with the increase of coverage rate, which is more obvious in **Figure 6(e)** and **(f)**. Comparing the results under coverage rate of 200 and 300%, it can be found that the increment of surface residual stresses are not obvious with increasing the coverage rate, since the surface of almost all area is covered by craters and the stress field reaches saturation. Though there is a little bit increment of the max CRS while coverage rate increases from 200 to 300%, there is not obvious comparing the coverage variation from 100 to 200%.

The simulated pictures about residual stress distribution under different coverage rates are shown in **Figure 7**, while *r* = 0.3 mm and *v* = 100 m/s. With increasing coverage, the number of craters on surface increases, and the surface residual stresses becomes more uniform. As well, the uniformity of residual stress distribution on subsurface is also improved while viewed from the cross section, and slight tensile stresses in the deep surface decrease gradually and disappear with the increase of coverage rate.

### *2.2.3. Influence of shot balls' radius on residual stress distribution*

According to the formula of *m* = \_\_4 3 *R*3 , the mass of shot ball is not only linearly related to the density, but also is proportional to *R*3 . In most of the industrial production, the criterion for choosing the size of shot ball is based on the surface roughness after SP. If it requires high quality of surface (small roughness), the small radius of shot ball is chosen. If there is no high quality requirement of surface, the big radius of shot ball is determined, because the cost of shot balls with small radius are much higher than that of big shot balls. Besides the influence of shot ball's radius on surface roughness, the shot ball also influences the residual stress distribution. So the influence of three kinds of different radius (*r* = 0.15, 0.3, and 0.6 mm) on the residual stress distribution of Ti-6Al-4V is simulated by the model, while the shot velocity is 100 m/s. The results of the residual stress versus depth are shown in **Figure 8**. From the simulation results, it can be found that the surface CRS and max CRS are higher while using small shot balls, but the layer depth of CRS is smaller and the CRS decreases rapidly with increasing depth. When using the big shot balls, both the surface and max CRS are smaller, however, the layer depth of CRS is deeper and CRS decreases slowly with increasing depth.

Some work have shown that when the coverage rate is 100%, the depth of CRS layer has the following relationship with the crater diameter and the shot ball's diameter.

*Z*0 = *a* × *D* − *c* × *R* (4)

**Figure 6.** Residual stress distribution in depth with different coverage rates, (a) *r* = 0.15 mm, *v* = 50 m/s; (b) *r* = 0.15 mm, *v* = 100 m/s; (c) *r* = 0.3 mm, *v* = 50 m/s; (d) *r* = 0.3 mm, *v* = 100 m/s; (e) *r* = 0.6 mm, *v* = 50 m/s; (f) *r* = 0.6 mm, *v* = 100 m/s.

Finite Element Dynamic Analysis on Residual Stress Distribution of Titanium Alloy and Titanium…

http://dx.doi.org/10.5772/intechopen.73120

33

constant coefficients. In general, the value of *a* is between 1 and 1.5, and the value of *c* is in

is the depth of CRS layer, R is the shot ball's diameter, and *a, c* are

*D* is the crater diameter, Z0

Finite Element Dynamic Analysis on Residual Stress Distribution of Titanium Alloy and Titanium… http://dx.doi.org/10.5772/intechopen.73120 33

In order to simplify the model, the collision between two shot balls are not considered, and SP position of each ball can be precisely controlled. Therefore, the different coverage rate can be simulated by the multi-layer shot balls and the coverage increases linearity with the number of layers. In this work, the characteristics of 100% coverage rate (4-layer shot balls), 200% coverage rate (8-layer shot balls), and 300% coverage rate (12-layer shot balls) are simulated. The shot velocity are set as 50 and 100 m/s based on experiment (mentioned in Section 2.2.4), and the average radius of shot balls are chosen as 0.15, 0.3, and 0.6 mm. The influence of coverage rate, radius of shot balls (*r*), and shot velocity (*v*) on residual stress distribution are

**Figure 6** shows the distribution of residual stress along depth with different coverage rates. While the radius and speed of shot balls are constant, the similar variation trend of residual stress can be obtained. With the increase of coverage rate, the CRS and the depth of stress layer are improved. In addition, the depth of max CRS decreases with the increase of coverage rate, which is more obvious in **Figure 6(e)** and **(f)**. Comparing the results under coverage rate of 200 and 300%, it can be found that the increment of surface residual stresses are not obvious with increasing the coverage rate, since the surface of almost all area is covered by craters and the stress field reaches saturation. Though there is a little bit increment of the max CRS while coverage rate increases from 200 to 300%, there is not obvious comparing the coverage variation from 100 to 200%.

The simulated pictures about residual stress distribution under different coverage rates are shown in **Figure 7**, while *r* = 0.3 mm and *v* = 100 m/s. With increasing coverage, the number of craters on surface increases, and the surface residual stresses becomes more uniform. As well, the uniformity of residual stress distribution on subsurface is also improved while viewed from the cross section, and slight tensile stresses in the deep surface decrease gradually and

choosing the size of shot ball is based on the surface roughness after SP. If it requires high quality of surface (small roughness), the small radius of shot ball is chosen. If there is no high quality requirement of surface, the big radius of shot ball is determined, because the cost of shot balls with small radius are much higher than that of big shot balls. Besides the influence of shot ball's radius on surface roughness, the shot ball also influences the residual stress distribution. So the influence of three kinds of different radius (*r* = 0.15, 0.3, and 0.6 mm) on the residual stress distribution of Ti-6Al-4V is simulated by the model, while the shot velocity is 100 m/s. The results of the residual stress versus depth are shown in **Figure 8**. From the simulation results, it can be found that the surface CRS and max CRS are higher while using small shot balls, but the layer depth of CRS is smaller and the CRS decreases rapidly with increasing depth. When using the big shot balls, both the surface and max CRS are smaller, however, the

Some work have shown that when the coverage rate is 100%, the depth of CRS layer has the

, the mass of shot ball is not only linearly related to the

. In most of the industrial production, the criterion for

investigated and discussed.

disappear with the increase of coverage rate.

According to the formula of *m* = \_\_4

density, but also is proportional to *R*3

*2.2.3. Influence of shot balls' radius on residual stress distribution*

32 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

3 *R*3

layer depth of CRS is deeper and CRS decreases slowly with increasing depth.

following relationship with the crater diameter and the shot ball's diameter.

**Figure 6.** Residual stress distribution in depth with different coverage rates, (a) *r* = 0.15 mm, *v* = 50 m/s; (b) *r* = 0.15 mm, *v* = 100 m/s; (c) *r* = 0.3 mm, *v* = 50 m/s; (d) *r* = 0.3 mm, *v* = 100 m/s; (e) *r* = 0.6 mm, *v* = 50 m/s; (f) *r* = 0.6 mm, *v* = 100 m/s.

$$Z\_{\mathbf{o}} = \mathbf{a} \times \mathbf{D} - \mathbf{c} \times \mathbf{R} \tag{4}$$

*D* is the crater diameter, Z0 is the depth of CRS layer, R is the shot ball's diameter, and *a, c* are constant coefficients. In general, the value of *a* is between 1 and 1.5, and the value of *c* is in

the range from 0 to 0.1. In order to obtain the values of *a* and *c*, the SP experiment should be carried out. Firstly, the same material as the workpiece is chosen and then the surface of this material is subjected to SP treatment. Secondly, the diameter of crater and the depth of CRS are measured. And at last, the values of *a* and *c* are calculated by linear fitting. Based on the values of *a* and *c*, the depth of CRS layer can be estimated. At the same shot velocity, the diameter of crater increases with increasing shot diameter, and thus the depth of CRS layer increases.

Finite Element Dynamic Analysis on Residual Stress Distribution of Titanium Alloy and Titanium…

SP is the process of consuming shot balls' kinetic energy and transfer the kinetic energy to the deformation energy of target material. So, after SP, the elastic and plastic deformations are introduced in the surface layer of target material. The shot balls' mass and velocity directly affect the value of kinetic energy. When the material of shot balls is same, the kinetic energy increases with the improvement of shot velocity. In SP experiment by using an air blast machine, the shot velocity can be varied and obtained by adjusting the air pressure. During the flight of shot balls, the velocity will be decreased because of the collision between them and the effect of air resistance, and the attenuation is related to the distance between the nozzle and the material. The smaller diameter of shot ball, the velocity attenuation is more obvi-

shot distances are shown in **Figure 9** [41]. From this figure, it can be found that when the shot distance is less than 2 m, the attenuation rates are proportional to the distance increment.

The shot velocity is also affected by the shot angle in addition to the attenuation with the distance. When the shot balls impact on the surface of workpiece at a certain angle, the velocity can be decomposed into two directions. One is perpendicular to the surface (normal velocity) and the other is parallel to the surface (tangential velocity). The former velocity contributes to the plastic deformation of the surface layer, but the latter only promotes the friction effect. Based on the above analysis, in the experiment of this work, the distance between nozzle and samples is 100 mm (0.1 m),

) with different diameters at different

http://dx.doi.org/10.5772/intechopen.73120

35

*2.2.4. Influence of shot velocity on residual stress distribution*

ous.  Relationship between decrement of shot velocity and shooting distance [41].

**Figure 7.** Residual stress distribution after SP with different coverage rates: (a) 100%; (b) 200%; (c) 300%. [36].

**Figure 8.** Influence of shot radius on depth distribution of residual stress (*v* = 100 m/s and coverage = 100%).

the range from 0 to 0.1. In order to obtain the values of *a* and *c*, the SP experiment should be carried out. Firstly, the same material as the workpiece is chosen and then the surface of this material is subjected to SP treatment. Secondly, the diameter of crater and the depth of CRS are measured. And at last, the values of *a* and *c* are calculated by linear fitting. Based on the values of *a* and *c*, the depth of CRS layer can be estimated. At the same shot velocity, the diameter of crater increases with increasing shot diameter, and thus the depth of CRS layer increases.