**4. The Re.N.Di.S. platform for planning**

Flood damage quantification and assessment is not only a (optional) requirement for the implementation of EFD, but it has recently become a requirement for the evaluation of the projects of countermeasures and for the eligibility of public financing.

Actually, with the recent 2015 decree of the President of the Ministry Council (2015DPCM in the following), a new discipline to evaluate the priority of public financing for flood protection projects came into force in Italy. To set priorities is fundamental when financial resources are scarcer than necessities. Actually, at the moment, the total needs for protection (by considering all the kinds of natural hazards, which are floods, landslides, flash floods, debris flows) amount to about 30G€, by far higher than the resources available at the moment (**Table 2**).

The model for prioritization proposed by the 2015DPCM is based on a score approach.

According to the 2015DPCM, flood defense projects have to be collected into the so-called Re.N.Di.S. procedure, which is a national repository of projects for soil defense. For a given project, either the assessments by public administration, or the level reached by design, or the effectiveness of the designed countermeasures are scored. The Re.N.Di.S. procedure allows to associate a total score to each proposed project, included in the Re.N.Di.S. procedure. The criteria are listed in **Table 3**.

For the sake of simplicity, it can be said that, for a given project of flood protection measures included in the Re.N.Di.S procedure, the priority to financing is based on its total score: the higher is its score, the higher is its ranking and therefore the probability for public financing. As shown in **Table 3**, flood damage evaluation score is relevant with respect to other criteria: if the damage evaluation is given, there is a 10-point score.

Without a uniform methodology for damage assessment and quantification, the damage assessment and the cost/benefit analysis (upon which also the prioritization of the public administration is based) is meaningless, either for the aims of EFD or for the Re.N.Di.S. pro-

**Table 3.** List of the evaluation criteria of the projects, with the respective scores (table taken from the 2015DPCM).

**Criteria contained in the DPCM Description of the criteria Min-max score**

considered

Design level (low/medium/high) Italian laws consider three levels of design. The

Directly endangered people Higher score for a higher number of protected population

Frequency Higher score if the project allows the protection

*reduction*

*involved*

environment

Compensation and mitigation measures Presence of compensation and mitigation

The higher is the total score, the higher is the ranking of the project for public financing.

administration and expressed in three levels. The cost/effectiveness rate should be

Index of Proportional Risk (IRP) Flood-Risk Assessment Model and Comparison to Collected Data

higher the level, the higher the score

Higher score for projects that allow the

*The criteria allow a higher score for projects which quantify the effectiveness in terms of damage*

*The criteria allow a higher score for projects which allow the reduction of the total number of people*

measures, which can mitigate the effects on

complementary to other projects

protection of goods at risk

to more frequent floods

0–20

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

23

3.3–10

0–10

0–60

0–30

4.2–30

0–30

0–5

*0—no quantification 10—damage quantification*

Priority by public administration The priority level is expressed by public

Completion Higher score for projects that are

Therefore, a uniform methodology for risk evaluation is required. It should be homogeneous at the national scale and based on databases available and preferably free to public and

In order to overcome the constraints of the present methodologies for risk mapping, Regione Piemonte public administration in collaboration with the Politecnico of Turin and the University developed a methodology for risk assessment and quantification, which is based

cedure. Heterogeneity would undermine the usability of the analysis.

on the quantification of the *Index of Proportional Risk* (IRP).

updated.

**5. The IRP model**

Goods at risk (properties, communication lines, etc.)

*damage reduction*

involved people

*Effectiveness of the project in terms of*

Effectiveness of the project in terms of the reduction of the total number of

**Table 2.** Total amount of proposed interventions for risk reduction and financing necessities (taken from [12]).

**Table 3.** List of the evaluation criteria of the projects, with the respective scores (table taken from the 2015DPCM).

Without a uniform methodology for damage assessment and quantification, the damage assessment and the cost/benefit analysis (upon which also the prioritization of the public administration is based) is meaningless, either for the aims of EFD or for the Re.N.Di.S. procedure. Heterogeneity would undermine the usability of the analysis.

Therefore, a uniform methodology for risk evaluation is required. It should be homogeneous at the national scale and based on databases available and preferably free to public and updated.


**Advances in Transboundary Aquifer Assessment**

Editors

**Sharon B. Megdal Anne-Marie Matherne**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Sharon B. Megdal Water Resources Research Center University of Arizona Tucson United States

Anne-Marie Matherne U.S. Geological Survey New Mexico Water Science Center Albuquerque United States

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Water* (ISSN 2073-4441) (available at: www.mdpi.com/journal/water/special\_issues/transboundary\_ aquifer).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-8451-5 (Hbk) ISBN 978-3-0365-8450-8 (PDF)**

Cover image courtesy of Andrew Robertson

© 2023 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.


*2.3. PSO Optimization*

The PSO algorithm will make the population evolve more intelligently after each iteration and can accumulate search knowledge, which is called an evolutionary algorithm [29–31]. The PSO algorithm does not use the survival of the fittest but uses a mechanism in which each individual in the population competes with the others to generate the global optimal solution. It generates the optimal solution through information sharing and a mechanism of cooperation between the individuals in the population [32,33].

Suppose the PSO consists of multiple particles. In the *D*-dimensional search space, the particle swarm contains *n* particles. The position of the *n*th particle in the *D*-dimensional space is defined as *xi*:

$$\mathbf{x}\_{i} = (\mathbf{x}\_{i1}, \mathbf{x}\_{i2}, \cdots, \mathbf{x}\_{iD}), \quad i = 1, 2, \cdots, n. \tag{5}$$

Suppose the current velocity of particle *xi* and its individual optimal historical position are *vi* and *pi*, respectively, as follows:

$$\begin{array}{l}v\_{i} = (v\_{i1}, v\_{i2}, \dots, v\_{iD})\\p\_{i} = (p\_{i1}, p\_{i2}, \dots, p\_{iD}).\end{array} \tag{6}$$

Then, for the entire particle swarm, the global optimal position is *Pg*:

$$p\_{\mathbb{X}} = (p\_{\mathbb{X}^1}, p\_{\mathbb{X}^2}, \dots, p\_{\mathbb{X}^D}).\tag{7}$$

At the *t*th moment, the velocity update formula of the *d*th dimension of particle *xi* is:

$$
\omega\_{\rm id}(t+1) = \omega v\_{\rm id}(t) + c\_1 r\_1 (p\_{\rm id}(t) - \mathbf{x}\_{\rm id}(t)) + c\_2 r\_2 (p\_{\mathcal{K}}(t) - \mathbf{x}\_{\rm id}(t)),
\tag{8}
$$

where *ω* is the inertia weight and is in the range [0, 1]. *c*1 and *c*2 are the acceleration coefficients. *r*1 and *r*2 are random coefficients, both of which are in [0, 1], which determine the motion of semirandom particles affected by the single and global optimal solutions.

The particle velocity update process has three main parts: the current initial velocity, the self-motion trajectory and self-trajectory correction. The influence of the current speed on the particle update speed can be adjusted by the inertia weight. The influence of the particle's own trajectory on the particle update speed can be adjusted by the acceleration coefficients and the random coefficients. When the trajectory is inaccurate, it needs to be corrected with the help of global optimization.

For particle *xi*, we can update the position *xid* of the *d*th dimension according to the velocity:

$$
\dot{x}\_{id}(t+1) = \dot{x}\_{id}(t) + \upsilon\_{id}(t+1). \tag{9}
$$


*Edited by Faieza Abdul Aziz*

This book attempts to bring together selected recent advances, tools, application and new ideas in manufacturing systems. Manufacturing system comprise of equipment, products, people, information, control and support functions for the competitive development to satisfy market needs. It provides a comprehensive collection of papers on the latest fundamental and applied industrial research. The book will be of great interest to those involved in manufacturing engineering, systems and management and those involved in manufacturing research.

ISBN 978-953-51-0530-5

Manufacturing System

Photo by efks / iStock


**Meet the editor**

Dr Soha Maad is the founder and manager of IGROW EEIG (Integrated Technologies and Services for Sustainable Growth European Economic Interest Grouping) at INVENT Dublin City University, Ireland. Dr Maad obtained her PhD from the University of Warwick, UK. She has 12 years of research experience in the application of Information and Communication Technologies,

including Grid and High Performance Computing to various domains. Dr Maad is an Associate Honorary Member of HRB Centre for Primary Care Research, Royal College of Surgeons, Ireland. She worked as a research fellow at Trinity College Dublin and University College Cork, Ireland; INRIA, France, and was an ERCIM fellow at Fraunhofer Institute of Media Communication, Germany.

Contents

**Preface IX**

**Section 1 Advances in Grid Computing - Workflow 1**

Takashige Hoshiai and Yoshiyori Urano

**Resources Management 47**

**Intensive Tasks on Grids 49**

Imran Ahmad and Shikharesh Majumdar

**Discovery in Grid Computing 71**  Leyli Mohammad Khanli, Saeed Kargar

Chapter 5 **Task Scheduling in Grid Environment Using**

Chapter 6 **Efficient Parallel Application Execution**

**Simulated Annealing and Genetic Algorithm 89**  Wael Abdulal, Ahmad Jabas, S. Ramachandram

Alfredo Goldman, Fabio Costa and Raphael Camargo

**Section 3 Advances in Grid Computing - Parallel Execution 111**

**on Opportunistic Desktop Grids 113**  Francisco Silva, Fabio Kon, Daniel Batista,

**Section 2 Advances in Grid Computing -**

Chapter 3 **Resource Management for Data**

Chapter 4 **A New Approach to Resource**

and Ali Kazemi Niari

and Omar Al Jadaan

Chapter 1 **w-TG: A Combined Algorithm to Optimize the Runtime**

Chapter 2 **On the Effect of Applying the Task Clustering for Identical Processor Utilization to Heterogeneous Systems 29**  Hidehiro Kanemitsu, Gilhyon Lee, Hidenori Nakazato,

**of the Grid-Based Workflow Within an SLA Context 3**  Dang Minh Quan, Joern Altmann and Laurence T. Yang


**4. Discussion**

The antibacterial agents used in the poultry for treatment or prophylaxis are implicated for the development of bacterial resistance [29]. Treatment in large groups of chickens is often done by oral administration [30]. Some studies reported that antibacterial agents or their metabolites are excreted in manure and residue can therefore be found in the environment [16,31,32].


**3. Results**

In this section, an example of a biplex network is analyzed in detail in order to highlight the possibilities offered by disposing of a centrality measure that establishes a classification of the nodes in order of importance.

For this purpose, let us consider a graph of 20 nodes, where each node represents a physical person; specifically, a player from a football team. Around this group the example is developed.

With these 20 nodes, let us proceed to construct a network with two layers that relate the nodes in a different way, taking two datasets (one for each of the layers). The calculation of the APABI centrality on this biplex network will allow us to determine the importance of each member of the team and obtain a classification of them in order of importance.

First, a layer *l*1 is constructed with the 20 nodes where the relationships between the members of the team are analyzed from the point of view of social or virtual relationships. Thus, an undirected graph is constructed with an adjacency matrix *A*1 in which two nodes are joined by an edge if they are related or linked through a social network. The graph of social relations between the nodes is shown in Figure 2. The data that are considered in this layer associated with each node are related to the number of messages that each person receives from their teammates in a period of time.

**Figure 2.** Graphs of the first layer (**left**) and the second one (**right**).

Secondly, a second layer *l*2 is constituted from the same 20 nodes of layer *l*1, but analyzing in *l*2 how the players relate to each other within the game, that is, if they associate with each other in the game. Thus, two players who combine with each other or pass the ball with some assiduity during a match are connected by an edge. From these links, it is possible to build a new adjacency matrix associated with this layer, which we denote by *A*2. In most team sports with a ball, each player occupies a specific position in the field of play, covering a certain area. Players (nodes) that occupy closer positions associate or relate more easily with each other than with those who are further away. For example, a defender is associated more with a midfielder than with a forward. The graph of game relations between the nodes is shown in Figure 2.

Both the links and the data associated to each of the 20 nodes of the graph are summarized in Table 1. There, the second column specifies the social links of every node while in the four column the game links between them are detailed. So, for example, the table shows that node 1 (player labeled as 1) has social interactions with the nodes {2, 5, 7, 9, 16, 17, 19, 20} while presents strong interactions within the game with the nodes {2, 4, 5, 6, 9, 12, 13, 14, 18, 19}.

**Table 1.** Data associated to the biplex network constructed from the team.

Datasets about the number of messages received during one day through social networks and the number of games that have played along the season are detailed in columns three and five, respectively. So, node 1 has received 15 messages in a day and has played 33 games in the season.

In this example, one of the advantages of working with biplex networks becomes manifest, such as the possibility of studying different relationships between the same set of nodes, analyzing the correlation between them. This example shows the advantage offered by adding what we can call a data layer in each of the multilayers of the network. We can assign the data that we consider appropriate to the specific relationships that we are representing by means of the corresponding graph. Thus, as can be seen in this example, in the layer where the social relations between the nodes are considered, we introduce the data corresponding to the number of received messages. However, in the second layer where game relations are represented, data are completely different, since now the number of games played are considered. Thus, each layer allows us to introduce one or more data sets related to the relationships of the nodes. This leads us to affirm that the inclusion of data in each of the layers enriches the nature of the problems that can be analyzed.

The objective in this example is to determine the most important or influential players within the team. For this task, two different aspects may be evaluated; on the one hand, the importance of the nodes from the point of view of the social relations that are established between them through messages, social networks, or any other virtual means. The nodes that have an intense social activity in the group create a very important union within the group, being very influential for other nodes. On the other hand, the importance of the nodes from the point of view of the game may be also evaluated, that is, which players are more important in the game, for their participation or quality. In other words, it is decisive to look for the leaders of the group, analyzing their importance from the social and technical point of view.

In order to determine the importance of the nodes of the biplex network of this example, the APABI centrality described in Section 2.3 has been calculated. Algorithm 2 has been executed using the information shown in Table 1. The numerical results shown in Table 2 are graphically displayed in Figure 3. This calculation gives us the importance of the nodes relating both layers. Figure 4 shows the final graph considering the information of two layers and the final result of the APABI centrality for each node, representing the size of each node according to its importance.

**Table 2.** Adapted PageBreak algorithm (APA) centrality for layers *l*1 and *l*2 and APA biplex (APABI) centrality for biplex network.

**Figure 3.** Centralities shown in Table 2.

**Figure 4.** Adapted PageRank algorithm biplex (APABI) centrality for the example object of this study.

The APABI centrality shows that the nodes that can be classified as the most important, the leaders within the group, are nodes 20 and 18, respectively. Note that node 20, which is the most important, is not the node that receives the most messages from its colleagues, being the node 18 is the one that receives the most messages from his teammates, but is second in the ranking.

Finally, it should be noted that, as discussed in Section 2.4, the use of the power method to calculate the stationary vector of the Markov chain that forms the stochastic matrix *MBI* provides the numerical stability needed in the implemented algorithm.

We have performed several tests with randomly generated adjacency matrices of different sizes in a range from 10 to 10,000 and we have obtained stable results. Matrices of dimension exceeding 105 cannot be stored in the central memory of most computers, except for sparse matrices. Consequently, the only matrix-arithmetic operation that is easily performed is a matrix-vector product. This makes possible to use this centrality algorithm for large matrices. In the scope of our research with urban network matrices, the sizes of the case studies are relatively large, with 2000–5000 nodes, approximately.


**1. Introduction**

Myocarditis is a clinical syndrome characterized by inflammation of myocardium. It can be produced by a myriad of etiologies including infectious, autoimmune, myocardial toxins, hypersensitivity reactions and physical agents. Virtually any infectious agent can produce myocardial inflammation and injury. Human myocarditis is most frequently caused by viral infection. Ongoing viral infection, myocardial destruction, and adverse remodeling can lead to persistent ventricular dysfunction and dilated cardiomyopathy. The modern molecular techniques have facilitated new insights into inflammatory autoimmune processes that affect the myocardium and ultimately result in acute or chronic dilated cardiomyopathy.

The clinical manifestations are highly variable, ranging from asymptomatic electrocardiographic or echocardiographic abnormalities to acute myocardial infarction-like syndrome, overt congestive heart failure, fulminant condition with new atrial or ventricular arrhythmias or profound cardiogenic shock and death. Myocarditis is occasionally the unrecognized culprit in cases of sudden cardiac death. Autopsy series have reported much higher rates of myocarditis than expected with overt clinical manifestation from different etiological agents. The prospective postmortem data have implicated myocarditis in sudden cardiac death of young adults at rates of 8.6 percent to 12 percent (Doolan et al., 2004; Fabre & Sheppard, 2005). Furthermore, it has been identified as a cause of dilated cardiomyopathy in 9 percent of cases in a large prospective series (Felker et al., 1999).

The clinical history in patients presented with myocarditis remains essential to envelop a wide variety of etiologies in the clinical scenarios, many of which are infectious (Brodison & Swann, 1998). In the past 10 years, however, viruses, including adenovirus, parvovirus B19, hepatitis C, and herpes virus 6, have emerged as significant pathogens (Mahrholdt et al., 2006). The geographical distribution can be of relevance for some forms of myocarditis. In selected countries, Chagas disease, Lyme myocarditis, acute rheumatic fever and disorders associated with advanced human immune deficiency virus infection are significant causes. Other important infrequent clinicopathologic variants in the etiological spectrum are systemic disorders like giant cell myocarditis, cardiac sarcoidosis and eosinophilic myocarditis. Additionally, drugs, vaccinations, toxins, physical agents like radiation, heat stroke and hypothermia can be the key point for some rare clinical diagnoses. The physical examination in patients with myocarditis might be normal, but more severe cases frequently evident for significant physical findings. Although histological findings remains the gold standard for establishing the diagnosis of myocarditis, low risk patients are often given a

Clinical Presentation 5

transplant free 11 years following initial biopsy, compared with only 45 percent in those

Classically, patients with acute myocarditis presents with a less distinct onset of illness with nonspecific symptoms related to the heart. Viral prodrome occurs between 20 and 80 percent of the cases, can be readily missed by the patient, and thus cannot be relied upon for diagnosis. They present with an established ventricular dysfunction and may respond to immunosuppressive therapy or their condition may progress to dilated cardiomyopathy. In a series of 245 patients with clinically suspected myocarditis, the most common symptoms include fatigue (82 percent); dyspnea on exertion (81 percent); arrhythmias (55 percent, both supraventricular and ventricular); palpitations (49 percent); and chest pain at rest (26 percent), (Kuhl et al., 2005). The presentation can mimic acute coronary syndromes in view of troponin release, ST segment elevation on electrocardiogram, and segmental wall motion

This group of patients with chronic active myocarditis represents the majority of older adult with myocarditis. They are also presents with a less distinct onset of illness and often insidious with symptoms compatible with moderate ventricular dysfunction such as fatigue and dyspnea. Affected patients may initially respond to immunosuppressive therapy but often have clinical and histologic relapses and develop ventricular dysfunction associated with chronic inflammatory changes, and mild to moderate fibrosis on histological study

This group of patients with chronic persistent myocarditis, whom also presents with a less distinct onset of illness, is characterized by a persistent histological infiltrate, often with foci of myocyte necrosis but without ventricular dysfunction, despite other cardiovascular

The previously depicted four forms are still used to describe the clinical presentation and progression of myocarditis, particularly in the absence of ongoing histological evaluation. These categories may also provide some prognostic information and may suggest which

A new diagnostic criteria derived from limited data was proposed in 2009. The Lake Louise Consensus Criteria utilizes the cardiac magnetic resonance imaging (CMR) for the diagnosis of myocarditis (Friedrich et al., 2009). The CMR enhances the ability to detect myocardial inflammation through noninvasive means, as well as to improve diagnostic accuracy. In these criteria, four major domains are considered when making the diagnosis including, clinical presentation compatible with myocarditis, evidence of new or recent onset myocardial damage, increased T2 signal or delayed enhancement on CMR (compatible with myocardial edema and inflammation), and endomyocardial biopsy evidence of myocardial inflammation. The use of CMR appears suitable to identify patients with significant ongoing inflammation, which may be especially important for patients with recurrent or persisting symptoms and in patients with new onset heart failure. The awareness came out that the

with more classic forms of acute myocarditis.

**2.2 Acute myocarditis**

abnormalities on echocardiogram.

**2.3 Chronic active myocarditis**

**2.4 Chronic persistent myocarditis**

symptoms such as chest pain or palpitation.

patients can or cannot benefit from immunosuppressive therapy.

including giant cells.

presumptive diagnosis if imaging studies and a compatible clinical scenario suggest new onset cardiomyopathy. The clinical manifestations of suspected patients with myocarditis will be discussed in details.


*4.7. Cloning of MdSAPs and qRT-PCR Analysis*

We extracted total RNA from previously frozen apple tissues according to the CTAB method and from *Arabidopsis* leaves, using Trizol reagent (Thermo Fisher Scientific-CN; Shanghai, China; https://www.thermofisher.com/cn/zh/home.html) [53]. Two micrograms of total RNA were collected for synthesizing first-strand cDNA. For cloning *MdSAP*, complete open reading frames were obtained via RT-PCR from fully expanded leaves of "Golden Delicious" apple, using specific primers listed in Table 5. The 50 - and 30 -untranslated regions (UTRs) were obtained with a Rapid Amplification for cDNA Ends kit (TaKaRa, Dalian, China). For the qRT-PCR assays, reverse transcription was performed with 1 µg of total RNA from each sample, followed by PCR-amplification of 1 µL of the product. We conducted the qRT-PCR assays in 20-µL reaction mixtures that contained 10 µL of SYBR® Premix Ex Taq™ (TaKaRa; Beijing, China; http://www.takarabiomed.com.cn), and used an iQ5 instrument (Bio-Rad, Hercules, CA, USA) as described before [34]. Thermal cycling included an initial 3 min at 95 ◦C; then 40 cycles of 10 s at 95 ◦C, 30 s at 58 ◦C, and 15 s at 72 ◦C; followed by 3 min at 72 ◦C and then 81 cycles of 7 s each, increasing by an increment of 0.5 ◦C from 55 ◦C to 95 ◦C. Three biological replicates were tested in each assay, and ∆Ct values were calculated by using *MdMDH* as our endogenous control [54]. Relative quantification was calculated according to the 2 −∆∆*C*t method [55], and dissociation curve analysis was performed for determining the specificity of the amplifications.


*2.5. Machining*

Finally, the gear was machined. After the roughing and finishing operations, the resultant gear is shown in Figure 4. Machined strategies were R-I (cavity mill) for teeth 1–25 and R-II (variable contour) after R-I and for teeth 1–25. For the finishing strategies: F-I for teeth 1–12, and F-II for teeth 13–25.

**Figure 4.** Roughing strategies result (**a**); Finishing strategies result (**b**).

#### **3. Predictive Model of Topography on Gear Flank**

The developed model estimates the gear teeth surface topography depending on machining parameters such as tool inclination and orientation, cutting geometric parameters and tool feed and speed values. The model was tested for two different finishing operations; a 5 continuous axes machining operation and a 3 + 2 axes machining operation in order to determine the influence of the machining number of axes on the surface finish.

The model follows the following steps:

1. **Interdental gap points and trajectories representation (Figure 5).** In this first step, for each interpolation point the model obtains from the machining program (CL data): the tool tip point position (*xj*, *yj*, *zj*) and the tool axis orientation defined by a direction vector (*uj*, *vj*, *wj*) in the workpiece reference system *XYZ*.

Milling trajectories can be obtained by means of the coordinates of successive tool tip points. Figure 5a shows the tool tip point positions obtained from the machining program of an interdental gap. In Figure 5b, the milling trajectories are represented. In the next steps, milling trajectories are evaluated after the elimination of initial and final noncutting movements, as shown in Figure 6a.

**Figure 5.** Interdental gap points (**a**) and trajectories (**b**) representation.

2. **Simulation area definition (Figure 6).** Secondly, as can be seen in Figure 6a, five simulation areas are defined (in black) in the feed direction. Each of these simulation areas is defined by 6 points in the feed direction and about 30 points in the direction perpendicular to the feed direction. The dimension of the simulation area depends on the programmer's criteria. The reference trajectory for the study is defined in magenta, corresponding to the one located in the middle of the simulation area. As an example, Figure 6b shows the milling trajectories followed by the tool tip in one of the simulation areas considered in Figure 6a. Black points in each milling trajectory represent the interpolation points given by the machining program for the selected simulation area. Next, the surface topography generated in each simulation area is predicted. In order to achieve this, a local reference system *OWXWYWZW* is defined for each simulation area (Step 4). First, the positions of the tool center point (*C*) in the selected simulation area are deduced in Step 3.

**Figure 6.** Simulation area (**a**) and reference trajectory points (**b**) for the study definition.

3. **Tool center position determination (Figure 7).** In order to obtain the coordinates of the tool center point (*C*) (Figure 7f), tool tip point coordinates (*xj*, *yj*, *zj*) and tool axis direction vector (*uj*, *vj*, *wj*) given in the machining program and the ball end mill, radius *R* needs to be taken into account. First, the line passing through the point (*xj*, *yj*, *zj*) and parallel to the vector (*uj*, *vj*, *wj*) is considered by means of the following equation:

$$\frac{\mathbf{x} - \mathbf{x}\_j}{u\_j} = \frac{y - y\_j}{v\_j} = \frac{z - z\_j}{w\_j} \tag{1}$$

Next, taking into account that the tool center point (*C*) is located at a distance equal to the tool radius *R* from the tool tip point of coordinates (*xj*, *yj*, *zj*), the coordinates (*x*, *y*, *z*) of the tool center point must fulfill this equation:

$$\sqrt{\left(\mathbf{x} - \mathbf{x}\_{j}\right)^{2} + \left(y - y\_{j}\right)^{2} + \left(z - z\_{j}\right)^{2}} = R$$

Therefore, in order to obtain the coordinates (*x*, *y*, *z*) of the tool center point when the tool tip point is located at a point of coordinates (*xj*, *yj*, *zj*) and tool axis vector is (*uj*, *vj*, *wj*), Equations (1) and (2) must be solved for *x*, *y* and *z*.

The red points shown in Figure 7a represent the positions of the tool center point (*C*) for the simulation area shown in Figure 6b. The surface generated by the positions of the tool center point approximates to a surface parallel to the tool tip point positions.

**Figure 7.** Tool center points (red) in relation to tool tip points (black) and to tool axis orientation and local reference system *OWXWYWZW* attached to workpiece (gear tooth). (**a**) Tool center points position; (**b**) M point definition; (**c**) Reference trajectory; (**d**) Tool tip and center point positions; (**e**) Reference for the prediction of the topography generated in the surface of gear teeth; (**f**) Ball end mill geometry.

4. **Definition of workpiece local coordinate system.** In order to predict the surface topography generated in each simulation area, a local coordinate system *OWXWYWZW* attached to the workpiece (gear tooth surface) is defined. The definition of this local coordinate system for each simulation area is based on the tool tip point positions and the tool axis orientations given in the machining program for the reference milling trajectory shown in magenta in Figures 6 and 7. The origin *OW* and the axis *ZW* are the elements of this system that are obtained first. In order to achieve this, a set of auxiliary elements (a point *M* and two unit vectors **u**ˆ and **q**ˆ) is considered. Firstly, two interpolation points, named *A* and *B,* which are located in the center of the reference trajectory (Figure 7b) are selected. The middle point *D* of the linear interpolation between points *A* and *B* is considered. The coordinates (*xD*, *yD*, *zD*) of point *D* are calculated as a function of the coordinates (*xA*, *yA*, *zA*) and (*xB*, *yB*, *zB*) of points *A* and *B*:

$$x\_D = \frac{x\_A + x\_B}{2}; y\_D = \frac{y\_A + y\_B}{2}; z\_D = \frac{z\_A + z\_B}{2} \tag{3}$$

Taking into account the tool axis orientation when the tool tip goes from point *A* to point *B*, the position of the tool center point when the tool tip point is located at point *D* can be calculated through Equations (1) and (2). Therefore, this position defined in this paper as point *M* (Figure 7b), with coordinates (*xM*, *yM*, *zM*), can be obtained from the resolution of the following equations:

$$\frac{\mathbf{x}\_M - \mathbf{x}\_D}{u\_{AB}} = \frac{y\_M - y\_D}{v\_{AB}} = \frac{z\_M - z\_D}{w\_{AB}} \tag{4}$$

$$\sqrt{\left(\mathbf{x}\_M - \mathbf{x}\_D\right)^2 + \left(y\_M - y\_D\right)^2 + \left(z\_M - z\_D\right)^2} = R \tag{5}$$

where (*uAB*, *vAB*, *wAB*) is the direction vector of tool axis, given by the machining program, when the tool tip goes from point *A* to point *B*. The point *M* is calculated for each simulation area and employed for the definition of the local system attached to each simulation area. In this paper, it is assumed that axis *ZW* of the local system passes through this point *M*.

Next, the vector defining the direction of axis *ZW* is calculated by considering two unit vectors. A unit vector **u**ˆ = (*ux*, *uy*, *uz*) parallel to the tool linear motion direction between points *A* and *B* is defined as:

$$\hat{\mathbf{u}} = \frac{\mathbf{p}}{|\mathbf{p}|} \text{ where } \mathbf{p} = (\mathbf{x}\_B - \mathbf{x}\_{A\prime}, y\_B - y\_{A\prime}, z\_B - z\_A) \tag{6}$$

A second unit vector **q**ˆ = (*qx*, *qy*, *qz*) that takes into account the direction at point *M* of the middle points of the milling trajectories selected in the simulation area is also defined (Figure 7c). Once the unit vectors **u**ˆ and **q**ˆ are obtained, a unit vector **w**ˆ = (*wx*, *wy*, *wz*) perpendicular to **u**ˆ and **q**ˆ is defined as follows:

$$
\hat{\mathfrak{u}} := \hat{\mathfrak{u}} \times \hat{\mathfrak{q}}\tag{7}
$$

This vector **w**ˆ is assumed to coincide with the direction of axis *ZW* (Figure 7d). In addition, in order to define the position of the origin *OW*, it is assumed that this point is located at a distance *R* from the point *M*. Therefore, taking into account that point *OW* is also located on a line passing through point *M* and parallel to vector **w**ˆ , the coordinates (*xOw*, *yOw*, *zOw*) of point *OW* are obtained by solving these equations:

$$\frac{x\_{\text{Ow}} - x\_M}{w\_x} = \frac{y\_{\text{Ow}} - y\_M}{w\_y} = \frac{z\_{\text{Ow}} - z\_M}{w\_z} \tag{8}$$

$$\sqrt{(\mathbf{x}\_{\text{Ow}} - \mathbf{x}\_{M})^2 + (y\_{\text{Ow}} - y\_{M})^2 + (z\_{\text{Ow}} - z\_{M})^2} = \mathbb{R} \tag{9}$$

Once the origin *OW* and the axis *ZW* have been obtained, axes *XW* and *YW* are defined. The axis *XW* is assumed to have the same direction as the unit vector **u**ˆ (Figure 7d). Therefore, the axis *XW* coincides with the feed direction of the tool between points *A* and *B*. Finally, the axis *YW* is perpendicular to axes *XW* and *ZW*. The direction of axis *YW* is defined by a unit vector ˆ**v** = (*vx*, *vy*, *vz*) calculated as the cross product between vectors **w**ˆ and **u**ˆ:

$$
\hat{\mathbf{v}} = \hat{\mathbf{w}} \times \hat{\mathbf{u}}\tag{10}
$$

Figure 7e shows an equivalent representation of Figure 7d. In Figure 7e, the local system *OWXWYWZW* is taken as a reference for the prediction of the topography generated in the surface of gear teeth. The surface topography will be simulated in a rectangular area defined along axes *XW* and *YW*, as shown in Figure 7e. In order to model the surface roughness, the equations expressing the trajectories of tool cutting edges in this local system are deduced in Step 6 as a function of the ball end mill geometry (Step 5), the tool axis orientation and the milling trajectories followed by the tool in the simulation area. Firstly, the coordinates of points defining the milling trajectories and the direction vector of tool axis orientations expressed in the workpiece system *XYZ* (given by the machining program and obtained in Step 1) must be transformed into the local system *OWXWYWZW*.

In order to express the coordinates of a point in the local system *OWXWYWZW* from its coordinates in the workpiece system *XYZ*, a homogeneous transformation matrix **T** is defined as a function of the unit vectors **u**ˆ, **v**ˆ and **w**ˆ and the coordinates of the origin *OW* in the system *XYZ*. The coordinates (*xj*, *yj*, *zj*) of the tool tip point positions given in the machining program can be expressed in the local system *OWXWYWZW*, (*xj W*, *yj W*, *zj W*), as:

$$
\begin{bmatrix} x\_j \\ y\_j \\ z\_j \\ 1 \end{bmatrix} = \begin{bmatrix} u\_x & v\_x & w\_x & x\_{Ow} \\ u\_y & v\_y & w\_y & y\_{Ow} \\ u\_z & v\_z & w\_z & z\_{Ow} \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x\_j^W \\ y\_j^W \\ z\_j^W \\ 1 \end{bmatrix} = \mathbf{T} \begin{bmatrix} x\_j^W \\ y\_j^W \\ z\_j^W \\ 1 \end{bmatrix} \rightarrow \begin{bmatrix} x\_j^W \\ y\_j^W \\ z\_j^W \\ 1 \end{bmatrix} = \mathbf{T}^{-1} \begin{bmatrix} x\_j \\ y\_j \\ z\_j \\ 1 \end{bmatrix} \tag{11}
$$

Similarly, a direction vector (*uj*, *vj*, *wj*) given in the machining program can be expressed in the local system *OWXWYWZW*, (*uj W*, *vj W*, *wj W*) by means of a matrix **T**1 depending on the unit vectors **u**ˆ, **v**ˆ and **w**ˆ

$$
\begin{bmatrix} u\_j \\ v\_j \\ w\_j \\ 1 \end{bmatrix} = \begin{bmatrix} u\_x & v\_x & w\_x & 0 \\ u\_y & v\_y & w\_y & 0 \\ u\_z & v\_z & w\_z & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} u\_j^W \\ v\_j^W \\ w\_j^W \\ 1 \end{bmatrix} = \mathbf{T}\_1 \begin{bmatrix} u\_j^W \\ v\_j^W \\ w\_j^W \\ 1 \end{bmatrix} \rightarrow \begin{bmatrix} u\_j^W \\ v\_j^W \\ w\_j^W \\ 1 \end{bmatrix} = \mathbf{T}\_1^{-1} \begin{bmatrix} u\_j \\ v\_j \\ w\_j \\ 1 \end{bmatrix} \tag{12}
$$

Once the tool tip point positions and the tool axis orientations are expressed in the local system *OWXWYWZW*, the equations for the cutting edge trajectory are deduced. In order to achieve this, the geometry of the tool cutting edges is first modeled.

	- The reference system origin *OT* is located on the tool tip being coincident with the tool axis
	- The axis *ZT* corresponds to the tool axis
	- The axis *XT* is radial and tangent to edge 1 projection in the plane containing point *OT* and perpendicular to *ZT* axis
	- The axis *YT* is perpendicular to axes *XT* and *YT* forming a right-handed system.

The position angle *ϕk* of cutting edge *k* with respect to axis *XT* can be expressed as:

$$
\varphi\_k = \frac{2\pi}{N\_t}(k-1) \tag{13}
$$

a) b)

**Figure 8.** Ball end mill geometry (**a**) and *k* edge geometry definition (**b**).

The position of a cutting edge point *P*(*i*, *k*), located at a height *zi* on the edge *k*, can be written in the tool reference system *OTXTYTZT* as a function of (a) the radius *Ri* and the position angle *βi* of cutting edge point at height *zi* and (b) the position angle *ϕk* of edge *k* as follows:

$$\mathbf{x}\_{P(i,k)}^T = \boldsymbol{R}\_i \cdot \cos(\boldsymbol{\beta}\_i + \boldsymbol{\varphi}\_k) \tag{14}$$

$$\boldsymbol{y}\_{P(i,k)}^T = \boldsymbol{R}\_i \cdot \sin(\boldsymbol{\beta}\_i + \boldsymbol{\varphi}\_k) \tag{15}$$

$$z\_{P(i,k)}^T = z\_i \tag{16}$$

where, from ball end mill geometry, the radius *Ri* and the position angle *βi* are:

$$R\_i = \sqrt{2Rz\_i - \left(z\_i\right)^2} \tag{17}$$

$$\beta\_i = z\_i \cdot \tan(i\_0) / R \tag{18}$$

Next, the equations of trajectories followed by cutting edge points are deduced.

6. **Tool axis orientation and points trajectories determination (Figures 9 and 10).** In this step, the trajectory followed by any cutting edge point in a five-axis milling operation is expressed as a function of cutting parameters, tool axis orientation and milling trajectories defined in the machining program. The cutting parameters are the feed value (F) in mm/min and the spindle speed (S) in rpm. Therefore, the tool feed in mm per revolution can be calculated as:

$$f\left(\frac{\text{mm}}{\text{rev}}\right) = \frac{F(\text{mm}/\text{min})}{S(\text{rpm})}\tag{19}$$

In order to define the cutting trajectories, as an example, the linear motion shown in Figure 7e, when the tool tip point goes from a point *O* to a point *I* in a milling trajectory of the simulation area, is considered. The procedure presented below is carried out in every section of the milling trajectories located inside the simulation area.

From the machining program, the coordinates (*xO*, *yO*, *zO*) and (*xI*, *yI*, *zI*) of points *O* and *I* and the direction vector (*uOI*, *vOI*, *wOI*) of tool axis orientation during this linear motion can be known. By means of the matrices given by Equations (11) and (12), the coordinates (*xOW*, *yOW*, *zOW*) and (*xI W*, *yI W*, *zI W*) of these points and the direction vector (*uOIW*, *vOIW*, *wOIW*) can be expressed in the local system *OWXWYWZW* attached to each simulation area.

Taking into account the feed direction of the ball end mill along the linear motion between points *O* and *I*, feed values *fx*, *fy* and *fz* along axes *XW*, *YW* and *ZW* can be defined. In order to achieve this, a unit vector **r**ˆ = (*rx*, *ry*, *rz*) that considers the tool feed direction between points *O* and *I*, is defined as a function of their coordinates:

$$\hat{\mathbf{r}} = \frac{\mathbf{q}}{|\mathbf{q}|} \text{ where } \mathbf{q} = \begin{pmatrix} x\_I^W - x\_O^W \ \ y\_I^W - y\_O^W \ \ z\_I^W - z\_O^W \end{pmatrix} \tag{20}$$

The components *fx*, *fy* and *fz* of the tool feed can be expressed as:

$$f\_{\mathbf{x}} = f r\_{\mathbf{x}} \tag{21}$$

$$f\_y = f r\_y \tag{22}$$

$$f\_z = f r\_z \tag{23}$$

By means of feed values *fx*, *fy* and *fz*, the position of the tool tip point along the linear displacement between points *O* and *I* is expressed as a function of the tool rotation angle *α* in radians. In Figure 9, it is assumed that the tool tip point is located at a point *O*1 whose coordinates are:

$$\mathbf{x}\_{O1}^{W} = \mathbf{x}\_{O}^{W} + f\_{\mathbf{x}} \frac{\mathbf{a}}{2\pi} \tag{24}$$

$$y\_{O1}^{W} = \left. y\_O^{W} + f\_{\text{y}} \frac{\alpha}{2\pi} \right. \tag{25}$$

$$z\_{O1}^{W} = \left. z\_{O}^{W} + f\_{z} \frac{a}{2\pi} \right. \tag{26}$$

where the tool rotation angle *α* goes from 0 to *αI-αO*, being *αI* and *αO*, the tool rotation angles when the tool tip point is located at points *I* and *O*, respectively.

**Figure 9.** 5-axis milling scheme with tool tilt angle *β* and lead angle *γ*.

In order to simplify the deduction of cutting edge trajectories, the tool axis orientation is defined in this step by means of two angles, which are represented in Figure 9: a tilt angle *β* and a lead angle *γ*.

These angles *β* and *γ* can be written as function of the tool axis direction vector (*uOIW*, *vOIW*, *wOIW*):

$$\beta = \arccos(w\_{OI}^W) \tag{27}$$

$$\gamma = a \tan(-\frac{v\_{OL}^W}{u\_{OL}^W}) \tag{28}$$

In order to obtain the equations expressing the trajectory of cutting edge points in the local system *OWXWYWZW* from their coordinates in the tool system *OTXTYTZT*, three auxiliary systems *O*1*X*1*Y*1*Z*1, *O*2*X*2*Y*2*Z*2 and *O*3*X*3*Y*3*Z*3, shown in Figure 10, are defined as a function of:

• The current position of the tool tip point given by the coordinates of point *O*1 (Figure 9). A translation of system *O*1*X*1*Y*1*Z*1 with respect to system *OWXWYWZW* (Figure 10a) is considered. The system *O*1*X*1*Y*1*Z*1 is shifted in *XW*, *YW* and *ZW* by distances *xO*1 *W*, *yO*1 *W* and *zO*1 *W* respectively. The homogeneous transformation matrix **T***W*1(*fx*, *fy*, *fz*, *α*) of this translation in *XW*, *YW* and *ZW* is:

$$\mathbf{T}\_{W1}(f\_{x'}, f\_{y'}, f\_{z'}, \mathbf{a}) = \begin{bmatrix} 1 & 0 & 0 & \mathbf{x}\_{O1}^W \\ 0 & 1 & 0 & y\_{O1}^W \\ 0 & 0 & 1 & \mathbf{z}\_{O1}^W \\ 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & \mathbf{x}\_O^W + f\_x \mathbf{a} / (2\pi) \\ 0 & 1 & 0 & y\_O^W + f\_y \mathbf{a} / (2\pi) \\ 0 & 0 & 1 & \mathbf{z}\_O^W + f\_z \mathbf{a} / (2\pi) \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{29}$$

• The lead angle *γ* (Figure 10b). The system *O*2*X*2*Y*2*Z*2 is rotated by an angle *γ* about axis *Z*1 in clockwise direction. The homogeneous transformation matrix **T**12(*γ*) of this rotation about *Z*1-axis is:

$$\mathbf{T}\_{12}(\gamma) = \begin{bmatrix} \cos(\gamma) & \sin(\gamma) & 0 & 0 \\ -\sin(\gamma) & \cos(\gamma) & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{30}$$

• The tilt angle *β* (Figure 10c). The system *O*3*X*3*Y*3*Z*3 is rotated by an angle *β* about axis *Y*2 in a clockwise direction. The homogeneous transformation matrix **T**23(*β*) of this rotation about *Y*2-axis is:

$$\mathbf{T}\_{23}(\beta) = \begin{bmatrix} \cos(\beta) & 0 & \sin(\beta) & 0 \\ 0 & 1 & 0 & 0 \\ -\sin(\beta) & 0 & \cos(\beta) & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{31}$$

• The rotation angle *α* (Figure 10d). The tool system *OTXTYTZT* is rotated by an angle *α* about axis *Z*3 in clockwise direction. The homogeneous transformation matrix **T**3*T*(*α*) of this rotation about *Z*3-axis is:

$$\mathbf{T}\_{\text{3T}}(\boldsymbol{a}) = \begin{bmatrix} \cos(\boldsymbol{a} + \boldsymbol{a}\_{O}) & \sin(\boldsymbol{a} + \boldsymbol{a}\_{O}) & 0 & 0 \\ -\sin(\boldsymbol{a} + \boldsymbol{a}\_{O}) & \cos(\boldsymbol{a} + \boldsymbol{a}\_{O}) & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{32}$$

where *αO* is the tool rotation angle when the tool tip point is located at point *O*.

Finally, the trajectory of any cutting edge point *P(i*, *k)* located at a height *zi* on the edge *k*, as shown in Figure 8, is expressed in the system *OWXWYWZW* of the gear tooth (Figure 9) as a function of the previous transformation matrices **T***W*1(*fx*, *fy*, *fz*, *α*), **T**12(*γ*), **T**23(*β*) and **T**3*T*(*α*) and the coordinates of the cutting edge point *P(i*, *k)* in the tool system *OTXTYTZT* (Equations (14)–(16) in Step 5):

$$
\begin{bmatrix} \mathbf{x}\_{P(i,k)}^W\\\mathbf{y}\_{P(i,k)}^W\\\mathbf{z}\_{P(i,k)}^W\\\mathbf{1} \end{bmatrix} = \mathbf{T}\_{W1} \begin{bmatrix} f\_{x'} \ f\_{y'} \ f\_{z'} \ \mathbf{a} \end{bmatrix} \cdot \mathbf{T}\_{12}(\gamma) \cdot \mathbf{T}\_{23}(\boldsymbol{\beta}) \cdot \mathbf{T}\_{37}(\mathbf{a}) \cdot \begin{bmatrix} \mathbf{x}\_{P(i,k)}^T\\\mathbf{y}\_{P(i,k)}^T\\\mathbf{z}\_{P(i,k)}^T\\\mathbf{1} \end{bmatrix} \tag{33}
$$

Once the cutting edge trajectories are obtained, surface topographies generated in each simulation area are obtained through the procedure described in Step 7.

**Figure 10.** Definition of auxiliary systems *O*1*X*1*Y*1*Z*1 (**a**,**b**; *O*2*X*2*Y*2*Z*2 (**c**) and *O*3*X*3*Y*3*Z*3 (**d**).

7. **Surface topographies determination (Figure 11).** The surface topography in the gear teeth is obtained from successive point positions of cutting edge points located on the *Nt* edges with the tool rotating movement. The surface topography is simulated in a rectangular area defined along axes *XW* and *YW*, as shown in Figure 7e. This area is divided into a discrete number of planes perpendicular to axes *XW* and *YW*. For each plane, taking into account the trajectories of cutting edge points given by Equation (33), the model predicts the area swept by the tool cutting edges during the tool rotation and feed motion. The profile generated at each plane is obtained from the lowest positions of marks left by cutting edge points. By considering the profiles generated in those planes, the 3D surface topography in the rectangular area can be predicted. This allows the profiles generated along axes *XW* and *YW* to be analyzed. As an example, Figure 11 shows the surface topography predicted for the simulation area selected in Figure 6b. In Figure 11a, a 3D representation of the predicted surface topography is shown. When the 2D profiles generated at two planes (*YW* = 0 and *XW* = 0) are considered (Figure 11b,c), it can be observed that profiles are not only composed of roughness marks left by tool cutting edges. In Figure 11b, the straight sections defining the tool trajectory (red line) can be observed. In Figure 11c, the form of the gear tooth surface can also be observed. In order to analyze the surface roughness, the effect of tool trajectory in each milling pass is removed for each profile predicted along *YW*-axis. As a consequence, roughness profiles (in black) are obtained. In Figure 11c, the form of the gear surface is also removed from the predicted profile and the roughness profile in black is obtained. In Figure 11d, the predicted surface roughness without the influence of the milling trajectory is shown. It can be observed that the step over between milling passes has a significant influence on the topography and the roughness peak-to-valley values. However, for this case, the effect of tool feed on roughness is less.

**Figure 11.** 3D simulated surface topography (**a**); roughness profiles (**b**,**c**) and surface roughness (**d**): ball end mill with 3 mm diameter and *Nt* = 3 cutting edges and feed value *f* = 0.18 mm/rev.

### **4. Model Validation Results**

In order to validate the surface roughness predictive model explained in the previous section, measured roughness after gear machining, and model predicted roughness, were compared.

#### *4.1. Gear Roughness Measurement*

For measuring gear surface roughness (*Ra* and *Rz*) (Table 4), confocal tridimensional Leica® DMC 3D and contact profilometer Taylor Hobson® Form Taylorsurf were used for the 3D and 2D roughness measurement, respectively.

**Table 4.** Roughness measured results for gear teeth after finishing strategies.

**Table 4.** *Cont.*

As it is mentioned in Section 3, the developed model predicts surface roughness in 5 different gear tooth zones. Nevertheless, surface roughness is only measured in two zones, always avoiding tooth edges (zones 1 and 5) and coinciding with the intermediate tooth zone. This fact must be taken into consideration when comparing predicted and measured roughness values. Therefore, for model validation, measured roughness values were compared to predicted values from zones 2–4.

#### *4.2. Gear Roughness Prediction*

The developed roughness model analyses 5 gear tooth zones and predicts surface roughness and generates surface topography and roughness profiles in each of these 5 zones. Thanks to this, different roughness profiles are predicted along tool feed direction *XW*. They are simulated for different tool edges number (*Nt*) and feed values (*f*) as can be seen in Figures 10 and 11.

In Table 5, roughness values for each gear tooth flank depending on machining type, cutting patterns and programmed scallop height are shown. These values also depend on tooth flank (concave or convex), and on the analyzed zone (1–5). In Table 5, the first of the rows corresponds to the concave flanks and the second to convex ones.

In this case, taking into account that the selected gear design method corresponds to the Gleason method, the roughness values (*Ra* and *Rz*) are slightly different in the 5 analyzed zones. This is a consequence of the Gleason method, which generates a gear with a variable gearing height. For this reason, in the areas closest to the outer diameter, the obtained roughness values are higher.

The shape of the roughness profiles is a consequence of the straight sections between two points of the machining program (interpolation points) that define tool trajectory. The distance between the different interpolation points within the same path is approximately 0.56 mm.

In the left-side figures (Figures 12 and 13), the straight sections defined by machining G1 (linear movement) are represented by a red line. The blue color represents the roughness profile. As can be seen, the roughness profile depends on the tool edge number (*Nt*) and programmed tool feed value (*f*). In the right-side figures (Figures 12 and 13), the form associated to workpiece component is eliminated, and for each roughness profile, both the arithmetic mean roughness parameter (*Ra*) and the average roughness parameter (*Rz*) are obtained.

**Figure 12.** Predicted roughness profiles for different feed values and for the same tool.

**Table 5.** Roughness predicted results for gear teeth.

**Figure 13.** Predicted roughness profiles for tools with different edge numbers and working with different feed values.

#### *4.3. Gear Roughness Predictive Model Validation*

Figure 14 shows a comparison, for the same gear tooth, between the profile generated after gear measurement (after the machining process) and the profile generated by the developed roughness predictive model. The variations shown in some roughness profiles peaks correspond to the trajectory when interpolating between the different interpolation points that make up the same machining path.

**Figure 14.** Roughness comparison of predictive model vs. measured roughness.

#### **5. Discussion**

The obtained roughness results, which determine gear quality, agree with cited works related to gear manufacturing in CNC machines, and, indicate improved surface quality and process versatility. Moreover, more detailed results are added in this work for optimal machining strategies and predicting roughness.

There are several aspects to mention in relation to machining strategies, cutting parameters and roughness results. On one hand, for3+2 axis machining, for both zig and zig-zag strategies, the obtained *Ra* and *Rz* values are slightly different for the zig strategy. The tool works down milling and machining time is increased due to non-cutting movements. Machining time is almost doubled. For 5 continuous axis machining, obtained results show differences that need to be taken into account. For zig cutting patterns, obtained roughness values for *Ra* and *Rz* are considerably lower than the ones obtained for zig-zag cutting patterns, about 50% lower. If both, 3 + 2 axis and 5 continuous axis machining are compared, better surface roughness values were obtained for 5 continuous axis machining with zig cutting pattern. On the contrary, the worst surface roughness results are obtained for 5 continuous axis machining with zig-zag cutting pattern. In relation to programmed stepovers, when programmed stepovers are reduced and the tool works up milling, the tool tends to go back to already machined path increasing roughness values. This is called 'rail effect' and is a consequence of tool deflection, which could explain why surface roughness values are higher for those machining trajectories with zig-zag cutting patterns.

The surface roughness prediction model was also used to analyze 5 axes and 3 + 2 axes machining strategies behavior. In the3+2 axes machining case, the roughness values obtained from concave and convex gear flanks were similar. Instead, in the 5 continuous axis machining case, obtained surface roughness values were higher in convex gear flanks. Thanks to the developed roughness model, this problem can be corrected by changing the tool axis attack angle in the machining surface in order to improve surface roughness. Indirectly, the model also detects the so-called 'rail effect' previously mentioned, that is, the obtained roughness values almost doubled. The presented model also determined the influence of cutting parameters such as tool feed values (mm/rev) and tool edge number. Thus, it can be concluded, that for a higher tool edge number (maintaining tool feed values), the obtained roughness results were significantly lower. Clearly there is a direct link between programmed feed values and obtained roughness results. The results showed that a small increase in the feed value has a significant effect on the final surface quality obtained. After roughness values analysis in the tool feed direction and comparing these with those generated in the direction perpendicular to it, it is clear that the roughness due to the application of the different stepovers, for the same feed value programmed, is always higher. This makes it a more restrictive parameter, which is why it is measured in each of the different machined flanks.

#### **6. Conclusions**

In recent years, technological advances have enabled gear manufacturing in general purpose machines. This is a feasible process that is less limited to gear size and geometry than traditional technologies. The machining of gears in multitasking machines is presented here as a real application for this type of technology due to its flexibility, size and the variety of geometries that can be machined in this type of machine. It can be concluded that the presented roughness predictive model fulfills its function, making it possible to predict and control the cutting strategies and parameters to obtain the required surface finish. The use of standard tools for the machining of gears requires less time for tool supply, provides a greater variety range and competitive time and costs. In addition, it is possible to exchange tools between different machines. Specific gear manufacturing machines require 'blank' material for gear manufacturing. On the contrary, general purpose machines with multitasking technology enable the entire gear machining in a single machine and take a single reference in the workpiece. This ensures tight execution times and the required quality consecution. Therefore, the presented study has demonstrated that although spiral bevel gears geometry is a complex geometry, it is possible to machine it with 3 + 2 machining axes kinematics. The rotary axis can be positioned and fixed from the beginning for a more robust process. The obtained surface roughness values are acceptable in this case. Moreover, 3 + 2 machining is cheaper because less interpolated axes are required. Also, 3 + 2 axes programming is not as challenging as 5 axes programming. This work validates the developed surface topography model for a ball mill. There is concordance between experimentally obtained values and theoretical values obtained by the model. Moreover, the obtained values for different programmed stepovers maintain the same tendency. Topographical simulation becomes an essential tool after the programming of each of the finishing strategies because it optimizes machining results without resorting to the trial and error method. Thus, costs and time is reduced, which is desirable under current market conditions.

**Author Contributions:** Á.Á. and A.C. designed and performed the experiments. Additionally, M.A. developed roughness model. A.C. wrote the paper. H.G. analyzed roughness measurements. Finally, L.N.L.d.L. contributed to the resources (machine, tools, material, etc.) and supervised all of the work carried out in this research.

**Funding:** This research received no external funding.

**Acknowledgments:** Thank you to the Department of Education, and to the Universities and Research program of the Basque Government for their financial support, by means of the ZABALDUZ program. We also thank the UFI in Mechanical Engineering department of the UPV/EHU for its support of this project.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; and in the decision to publish the results.