Numerical study via total Lagrangian smoothed particle hydrodynamics on chip formation in micro cutting

Numerical simulation is an effective approach in studying cutting mechanism. The widely used methods for cutting simulation include finite element analysis and molecular dynamics. However, there exist some intrinsic shortcomings when using a mesh-based formulation, and the capable scale of molecular dynamics is extremely small. In contrast, smoothed particle hydrodynamics (SPH) is a candidate to combine the advantages of them. It is a particle method which is suitable for simulating the large deformation process, and is formulated based on continuum mechanics so that large scale problems can be handled in principle. As a result, SPH has also become a main way for the cutting simulation. Since some issues arise while using conventional SPH to handle solid materials, the total Lagrangian smoothed particle hydrodynamics (TLSPH) is developed. But instabilities would still occur during the cutting, which is a critical issue to resolve. This paper studies the effects of TLSPH settings and cutting model parameters on the numerical instability, as well as the chip formation process. Plastic deformation, stress field and cutting forces are analyzed as well. It shows that the hourglass coefficient, critical pairwise deformation and time step are three important parameters to control the stability of the simulation, and a strategy on how to adjust them is provided.


Introduction
Cutting, as a fundamental mechanical process, has been studied for over one hundred years. With the advancement of machine and cutting tools, the controllability of material removal is being enhanced step by step, leading to higher form accuracy and surface quality of components. From conventional to precision machining, the thickness of the undeformed chip decreases from macro to micro/nanoscale, and the material mechanism becomes more complex because of the size effect [1].
To optimize process parameters and improve production efficiency, comprehensive knowledge about cutting mechanism, which is always obtainable through numerical and experimental studies, is essential. The merits of numerical simulation are that the cutting conditions can be easily improved over a wide range with a high repeatability, and that the deformation/stress field in the material can be easily obtained with a high resolution.
In general, methods for cutting simulation can be categorized into mesh-based and particle-based. One example of a mesh-based approach is the finite element method (FEM). As reported in literature, however, there are some issues in cutting simulations performed using FEM. Firstly, the grids near the tool tip undergo a large distortion, which reduces the accuracy of the FEM algorithm. As a result, mesh reconstruction is usually required to modify the topology of the nodes, which costs a high amount of time [2]. Secondly, there are difficulties in handling the chip separation. For example, the geometrical or material criterion involves the deletion of node connections or mass elements as the strain or stress reaches critical values, which can induce arbitrariness and error. There is also no uniform standard to rely on [2,3]. This issue can be avoided using the Eulerian approach, but the chip shape has to be pre-defined, and investigation on residual stress is limited [4]. The arbitrary Lagrangian-Eulerian (ALE) formulation, which combines the merits of the two FEM methods, has always been employed [5], but this method still suffers from the large deformation [6,7].
By contrast, these disadvantages can be easily avoided using a particle-based method. For example, molecular dynamics (MD) has been widely used in the 1990s to study nanometric cutting science [8]. In MD simulation, each particle represents an atom in the material. With the use of a potential function, forces on the particles are calculated, and the system evolves following Newton's law. However, because the model is built at the atomic level, large amounts of computation resources are required, and hence the simulation scale in which this method is feasible is extremely limited. Although some efforts have been made to improve the use of computational resources [9,10], it is still difficult for a common workstation to handle the issues at the micrometer level, like in a precision cutting process.
Smoothed particle hydrodynamics (SPH) is a numerical method that, like MD, it is particle-based, and that, like FEM, it is capable of simulating the process from micro to macroscale. SPH was firstly developed for astrophysical simulation [11] and is now widely used in engineering, especially computational fluid dynamics and fluid-structure interactions [7]. The key concept of SPH is to represent a function and its derivative using two steps: kernel approximation and particle approximation [12]. Basic equations of continuum mechanics can then be converted into ''particle-form''. In SPH, a particle represents a material element with a finite volume, which would deform and move under a stress field. However, conventional SPH exhibits some shortcomings, such as tensile instability, lack of consistency, zero-energy mode, and difficulty in enforcing essential boundary conditions [6,13]. To solve these problems, various kinds of methods to improve the accuracy and robustness of the algorithm have been developed. These methods include kernel and integral correction [14], mass density and kernel gradient correction [15], and adaptive density and variable smooth length [16]. Artificial stress is also introduced to dampen the excessive motion of particles.
However, although conventional SPH is a Lagrangian method, the kernel function is Eulerian, which causes the instability. Thus, the total Lagrangian SPH (TLSPH) is being developed [17], and ''hourglass control'' is employed to eliminate zero-energy modes [18]. In cutting simulations performed using conventional SPH, sources of instability, such as non-uniform particle distribution and low cutting speed, have been reported [19]. Influences of both the material model and the SPH parameters on simulation result have also been discussed [20,21]. Thus far, there have been many SPH studies on macro/micro metal cutting, grinding of brittle material [22], and laser-or vibration-assisted machining [23,24]. On the other hand, TLSPH, which is more suitable for solid material simulations, such as for fracture and machining, has been attracting more research interest, but only in recent years [25,26]. For example, a series of studies on the scratching process by a spherical indenter has been conducted [27][28][29][30][31]. Surface topology and scratching forces under various load conditions have been compared. Effects of process temperature, material model, scratching velocity, and ploughing phenomenon have also been investigated. However, compared to studies using SPH, fewer studies on the use of TLSPH to simulate cutting have been reported.
Unlike the MD method, where a model is constructed exactly, atom by atom, and only atomic interaction is considered, numerical instability can also arise in a cutting simulation performed using TLSPH, because of many factors. These factors can be classified into two groups: TLSPH algorithm setting and cutting model parameter. The effects of these two groups are not yet well understood. Therefore, a systematic numerical study is necessary to improve the stability and reliability of the simulation, which is the objective of this paper.
In the next section, the principles of conventional and total Lagrangian SPH are briefly reviewed. Details of the cutting model and simulation setup are then presented. Numerical instability, chip formation process, stress field, and cutting force under different simulation conditions are comprehensively discussed in Sect. 3. The result shows that the instability can be eliminated by properly choosing the hourglass coefficient, and the criterion for updating the reference configuration and time step. A strategy for improving these parameters is also provided. 2 Methods and model 2.1 Principles of conventional and total Lagrangian SPH SPH can be thought as a reformulation of basic equations in continuum mechanics, such as conservation of mass, momentum, and energy, which are a group of differential equations. Therefore, the topic of this subsection is to represent a function and its derivative in a special form, shown as follows To obtain the function value at a special position x i , values at n neighbor points f(x j ) weighted by a kernel function W are used. In addition, a mass m and a density q are attached to each of the points, which constitute the material meaning of SPH particles. The derivative of function f is then represented as a derivative of the kernel W. From Eqs. (1) and (2), continuum mechanics equations can be rewritten to perform a numerical simulation. SPH particles move under the forces on them, and as a result, their neighbor lists would be updated according to the new particle distribution, as in an MD simulation. This is the Eulerian feature of the kernel function. The corresponding representations in TLSPH are shown as follows The most important difference between conventional and total Lagrangian SPH is that, in TLSPH, the material coordinate X is used instead of the spatial coordinate x, which is just the Lagrangian description in continuum mechanics. The material coordinate X is used to distinguish each particle in the initial state of the system, i.e., reference configuration, and would not change with subsequent material deformations in a current configuration. Because the kernel W is now a function defined in reference configuration, the particle's neighbors would not change with system evolution. The superscript ''0'' indicates that the values and derivative are taken with respect to the material coordinate.
As reported in Ref. [32], internal force, which is decided by the constitutive model, is required to solve the equations of motion, wherein differences occur for solids and fluids, resulting in the two formulations above. For an elastic solid, the internal force is a function of strain energy, so the deformation gradient relating the initial and current positions of a particle is necessary, and thus the material coordinate has to be selected (TLSPH). Meanwhile, for a fluid, the internal force is a function of density, so it is more convenient to use the spatial coordinate (conventional or updated SPH). These, in turn, lead to different forms of physical objects. For example, the first Piola-Kirchhoff stress and Cauchy stress are used in TLSPH and conventional SPH, respectively. It is also noted that the physical laws driving the system evolution are the same, irrespective of which description is used.
In this work, TLSPH simulation is conducted based on a customized package SMD in LAMMPS [33]. A two-dimensional spiky kernel function, which reduces the clustering effect caused by a spline kernel [34], is employed where h is the smoothing length and R is the distance between two points. The hourglass control algorithm has been embedded into the code to improve the simulation stability and convergence through the addition of a correction term to the nodal force. As shown in Ref. [18], the main part of the nodal force on particle i can be expressed as Eq. (6), where j denotes the neighbors in the domain S of kernel function W, V * 0 the particle volume in the reference configuration, P * the first Piola-Kirchhoff stress, and L * -1 the reverse of a correction matrix (see Eq. (7)) to fulfill first-order completeness. The subscript of the kernel function indicates which particle the derivative is with respect to.
To obtain the force f i , material behavior law, i.e., the constitutive model, which will be discussed in the next section, is required. The hourglass term of the nodal force is where c is the hourglass control coefficient and E * is Young's modulus. d i ij is a fix term for particle separation to minimize the error that arises from the mapping between the reference and deformed configurations where F i is the deformation gradient of particle i.
Numerical study via total Lagrangian…

Cutting model and simulation setup
A two-dimensional cutting model is used as shown in Fig. 1. The workpiece and tool are constituted by SPH particles located on the simple cubic lattice. The lattice constant is a, which is set to either 25 nm, 50 nm or 100 nm for different particle resolutions. The left and bottom sides of the workpiece are fixed during the simulation, and the tool moves along the -x direction with a speed of v, which is either 10 m/s, 20 m/s or 50 m/s here. Absolutely rigid and elasto-plastic tools are used, with rake and clearance angles of 0°and 12°, respectively. Geometrically, the tool edge radius is zero. In SPH, however, there is a volume attached onto each particle, so the edge radius can be thought of as the distance between SPH particles (25 nm or 50 nm, in most cases). The rigid tool, whose particles do not interact with each other and do not change relative positions, is used in most of the simulation tests. While the stress distribution in the tool is being studied, the elasto-plastic body is used. Finally, the cutting depth d is set to either 1 lm or 5 lm. The effects of five TLSPH algorithm parameters are considered in this study. The first is the hourglass control coefficient (c), which ranges from 1 to 20. The second is the criterion for updating the reference configuration. As mentioned above, TLSPH represents a function and its derivative in the reference configuration. For a process with localized strain and high strain rates, such as cutting, particles would fly away, and some non-physical force pulses would occur (see Fig. 2) if the reference configuration is unchanged during the simulation. Therefore, the reference configuration should be periodically substituted by the current state. For this to be done, a pairwise deformation is defined If g exceeds the critical pairwise deformation (g c ), the reference configuration will be updated. It is natural that the smaller the g c , the more frequent the updates. The value of g c is set in the range of 0.1-5. The third parameter is the smoothing length (h), which decides the range of neighbor summation in Eqs. (3) and (4), just like the cut-off radius in MD. In this study, h starts from a and increases to 4a, where a is the initial SPH distance.
The fourth parameter is the scale factor s of time step Dt, which is determined by the CFL criterion where l is the characteristic distance between particles, and c 0 is the speed of information propagation and is decided by the bulk modulus K, shear modulus G, and mass densityq for a solid The time step used in this study is 0.5-3 ps. The last parameter is the artificial viscosity (j), which is widely used in conventional SPH simulations on fluid and solid to dampen excessive particle vibration. Stress can be decomposed into isotropic and deviatoric parts and is used to calculate the nodal force, where p is the pressure, I the unit matrix, and r d the deviatoric tensor. The pressure is decided by the equation of state, where K is bulk modulus, and q and q 0 are mass densities in the deformed and reference configurations, respectively. The deviatoric stress of workpiece is described by the Johnson-Cook model, which determines the effective yield stress, as follows where A is the initial yield stress, B the hardening constant, e p the effective plastic strain, n the hardening exponent, C the strain rate constant, _ e p and _ e p;0 the strain rate and its reference value, respectively, m the thermal softening exponent, and T, T 0 and T m current temperature, room temperature, and melting point, respectively. The workpiece material is a kind of stainless steel (Stavax ESR) widely used in injection molding. The chemical Fig. 1 Two-dimensional TLSPH cutting model components are the following: C (0.38%), Si (0.9%), Mn (0.5%), Cr (13.6%), V (0.3%), and Fe (majority). The heat treatment includes preheating at 600-850°C and quenching at 1 020-1 030°C. The material parameters are obtained from Ref. [35] and the supplier. For the non-rigid tool, the linear elastic/ideal plastic model is used where G is the shear modulus, and e d and _ e d are the deviatoric parts of the strain and strain rate tensors, respectively. If the second invariant of r d trial is larger than the yield stress r yield , plastic deformation would take place. In this study, the material parameters for the non-rigid tool are those of a diamond. It should be mentioned that in practice, a diamond tool cannot be used in the machining of ferrous metals without extra manipulation such as ultrasonic vibration, because chemical tool wear will otherwise occur. Meanwhile, the objective here is just to make a tool with a finite strength and investigate its response during cutting. Interaction between the tool and the workpiece is described by the Hertzian contact model where E contact is the contact stiffness, which is equal to the elastic modulus of the tool; R i and R j are the contact radii of two particles, each of which is half of the initial distance between SPH particles, i.e., 0.5a; and r ij is their mutual distance. All of the model information and simulation settings are summarized in Table 1.

Results and discussion
In this section, the effects of various TLSPH parameters are firstly discussed. Cutting under different process conditions is then simulated. Based on the results, guidance on how to reduce the instability in the cutting simulation by TLSPH, as well as the recommended parameters, are obtained. Figure 3a shows the deformation and von Mises stress distributions under different hourglass intensities. Without hourglass control and artificial stress (c = 0, j = 0), strong numerical instability occurs. The chip profile is not well-defined, and the machined surface is fractured, which should not happen in the micro cutting of ductile metal. The low plastic strain in the chip is also a non-physical phenomenon, because the material experiences severe shear deformation during the chip formation. As c increases, the boundaries of both the deformation and the stress field in the primary shear zone become clearer, and a well-defined curled chip can be observed. In addition, the arrangement of SPH particles in the chip become more ordered. As mentioned in Ref. [18], this improvement in the particles topology is a result of hourglass control and benefits having a stable solution.

Effects of TLSPH settings
It should be emphasized that no information about material lattice deformation (e.g., dislocation and phase Critical pairwise deformation for updating the reference configuration (g c ) 0.1-5 Smoothing length of the kernel function (h) a-4a Initial distance between SPH particles or lattice constant (a)/nm 25-100 Artificial viscosity coefficient (j) 0.5-5 Time step factor (s) 0.01-0.5 transition) can be obtained from the particle arrangement, because SPH particles represent continuum material elements rather than atoms. However, a too large c also triggers instability originating from the tool tip and the machined surface. A potential interpretation may be that excessive hourglass control restrains the particle disorder near the tool edge. In other words, an improper c interferes with the large plastic deformation during cutting. The instability on the machined surface is interesting. The pit evolves from the flat surface when the tool has already passed by, so it is not caused by the tool-workpiece interaction. The reason is thought to come from the nature of the TLSPH method. As the neighbor list of an SPH particle is built in the reference configuration, where the particle's coordinates would not change during deformation, ghost forces may arise between particles separated by long distances in the current configuration but are near to each other in the reference configuration. To solve this problem, updating the reference configuration is required, which would be discussed later. As shown in Fig. 3b, the principal cutting force exhibits an asymptotic convergence with the increase in c. For the largest c, which is 20, the frequency domain analysis of cutting forces reveals some high-frequency peaks (63 MHz and 112 MHz for the principal cutting force, 50 MHz for the thrust force), which implies an excessive increase in material stiffness.
Further tests are conducted using different criteria for updating the reference configuration. As shown in Fig. 4a, reducing the g c , which means to update the reference configuration if a smaller change in the pairwise distance occurs, indeed improves the numerical stability. By contrast, the simulations are all interrupted with an increase in g c , and the stress field profile changes are accompanied by an intensification of the elastic wave. The chip morphology also has a strong dependence on the g c and c. In general, larger g c and c depress the formation of localized shear band as well as the chip curling. The chip formation even becomes difficult when g c is 5 and c is larger than 10. On the other hand, the localized shear band also undergoes blurring as c decreases to 1. It should be carefully noted that the serrated chip (g c = 0.1 and c = 1) is just a result of numerical instability. Because chip morphology is related Numerical study via total Lagrangian… to the cutting mechanism, analysis on the deformation and stress fields is necessary to avoid being misled by the results. The system state will be kept after the configuration is updated, so there is no obvious influence on the cutting forces as long as the simulation is stable (see Fig. 4b). Finally, computing times under different g c are compared in Fig. 4c. Theoretically, a smaller g c would result in a higher updating frequency and a longer time, but this effect is slight in the current test.
The smoothing length of the kernel function decides the size of the neighbor list for an SPH particle. It also influences the smoothness of the numerical solution. As shown in Fig. 5a, no material removal and deformation take place when this length is equal to the initial particle distance a. With increasing h, more neighbors are accounted while the function value is calculated according to Eqs. (3) and (4). As a result, plastic deformation and stress fields become smoother, and some details, such as localized shear bands, are swept out. The cutting force converges, and the vibration is reduced (see Fig. 5b). For example, the thrust force component at 0.  MHz is reduced by half as the h increases from 3a to 4a. In practice, the smoothing length should be properly set to balance the convergence and the spatial resolution. Computing time should also be considered, because its increase with the neighbor size is nonlinear, as shown in Fig. 5c. In this study, it is verified that an h of 3a can meet the requirement above.
The initial lattice constant a determines the material element volume represented by one SPH particle. The smaller the a, the higher the model resolution, and the g c should be reduced to avoid instability. As shown in Fig. 6a, fine structures occur in the chip, and the residual stress layer is clear with an increase in the particle resolution. Because more shear bands are formed, the cutting energy introduced by the tool is further released by plastic deformation, which results in lower cutting forces (see Fig. 6b). The force vibration is also obviously reduced when a is 25 nm, so the residual stress layer is thinner than that in a coarse particle model (a = 100 nm).
However, a new type of numerical instability happens at the end of the simulation with the smallest a, which is 25 nm, as the chip collapses, as shown in Fig. 7a. Based on the discussion above, a large hourglass control or pairwise deformation criterion may cause the instability, but reducing c and g c does not mitigate the problem. Although a smaller c eliminates the collapse, deformation in the chip is unreasonably low, which has also been mentioned. This issue is finally solved by reducing the time step factor s in Eq. (11). In principle, any numerical method discretizes the time domain in the simulation of a dynamic process. If the time step is too large, the discretization error would induce instability. As shown in Fig. 7b, if s increases to 0.5, many SPH particles fly away just as the tool contacts the workpiece, while if s decreases to 0.05, the simulation becomes stable. Furthermore, reducing the time step has no influence on the stress field and cutting forces, which is to be expected (see Fig. 7c).

Numerical study via total Lagrangian…
Finally, the influence of artificial viscosity is studied (see Fig. 8). Similar to the hourglass control, increasing j to some extent improves the numerical stability. It smoothens the residual stress field and reduces the highfrequency vibration of the cutting forces. However, artificial viscosity may excessively impose fluidal behavior on the workpiece when j is larger than 2, and thus it is not used in the following tests.

Effects of cutting parameters
Cutting-speed-induced instability is firstly studied. As shown in Fig. 9a, the lower the speed, the stronger the instability. Furthermore, the plastic deformation in the chip is abnormally low when v is 10 m/s. Once again, the instability is eliminated only after the time step is reduced from 0.1 to 0.05. However, the low deformation state still exists, so the hourglass control coefficient should be increased. This, in turn, triggers instability (c = 30), and the time step has to be further reduced (s = 0.01). In general, a larger c and a smaller s are required for a slower cutting speed, but this, in turn, reduces the simulation efficiency drastically. The mechanism of the low-speedinduced instability will be further investigated in the following work. Figures 9b and c reveal that as v increases, the chip becomes less curled, and the principal cutting force increases, which is a result of strain rate hardening.
To study the influence of cutting depth, the ratio of the spatial resolution of SPH settings to the model size is kept. In other words, in the simulation where d is 5 lm, the SPH particle distance, kernel, and contact radii, as well as the  Numerical study via total Lagrangian… length and height of the workpiece, are all five times those where d is 1 lm. The cutting speed, on the other hand, is unchanged. This test could, therefore, give empirical knowledge about the dependence of numerical stability on the scaling of the simulation. As shown in Fig. 10a, instability occurs in the scale-up model. There is neither low plastic deformation in the chip, nor collapse at the tool edge or the machined surface, so reducing the time step is the only way to resolve the instability. Chip morphology and stress distribution are similar for two cutting depths (see Fig. 10b). The principal cutting force increases with the cutting depth under the same scale factor of five. Violent vibration of the thrust force occurs where d is 5 lm, which implies a more severe flank wear (see Fig. 10c). Figure 11 shows the results of simulations with rigid and non-rigid tools. When the material strength of the tool is considered, a high stress concentration occurs at the tool edge, and the chip morphology is not influenced (see Fig. 11a). Cutting forces encounter the most change. As the tool just touches the material, violent force vibration takes place, which results in large instantaneous force values (see Fig. 11b). This result corresponds to the occurrence of an intensive elastic wave in the tool (see Fig. 11c). The wave originates from the tool edge, then propagates and reflects through the tool body. In the following 1.5 lm distance, the cutting evolves to a steady state, where the forces converge to the values in the simulation involving a rigid tool. The stress field in the tool becomes almost time-invariant and is concentrated only at the edge, chip-tool contact face, and flank face. Although there is no plastic deformation of the tool, owing to its high stiffness and yield point, the high-frequency force peak is a potential risk that may fracture the tool edge, especially in the cutting of hard or brittle materials, so tool-workpiece impact should be avoided in practice.

Strategy for parameter adjustment
As discussed in Sects. 3.1 and 3.2, many items, including the TLSPH settings and process parameters, influence the simulation results. In this section, these results are summarized and further analyzed, in order to obtain a strategy about how to conduct a stable cutting simulation using TLSPH.
Firstly, the hourglass coefficient, reference configuration updates, and time step are three major parameters that are used to control the stability. Their effects are different, and if the instability is caused by one improper parameter, it is difficult to remedy by adjusting the others. As shown in Fig. 12, plastic deformation in the chip can be modified   Secondly, instabilities caused by cutting parameters can also be classified into three types, i.e., low plastic deformation, micro hole formation at the tool edge and machined surface, and chip collapse. As a result, it may be necessary to adjust the three TLSPH parameters to stabilize the simulation. In detail, a low cutting speed results in both a low plastic deformation and the chip collapse, so c should be increased, and the s should be reduced. The chip collapse also occurs in a scale-up model (for a simulation using larger cutting depth) and when the SPH particles resolution is improved. Therefore, it is necessary to reduce  Table 2. In addition, the parameter settings in Table 3 are recommended as an initial choice when the cutting process is simulated using TLSPH.

Conclusions
In this paper, micro cutting process is studied using the TLSPH method. The influences of various parameters, including the algorithm settings and the cutting model, on the numerical instability and the chip formation are comprehensively discussed. A methodology for improving the simulation stability and a group of recommended settings are provided. The results give further knowledge about the behavior of TLSPH when used to study micro cutting, and the main conclusions can be summarized as follows.
(i) Hourglass control is an efficient approach to modify the plastic deformation state in the chip. The larger this coefficient, the stiffer the material. It should not be too large to avoid instability and difficulty in the chip formation. In most of the cases in this study, a value of 5 or 10 is used. (ii) Updating the reference configuration is also necessary in the simulation of large-strain or strain-rate process, such as cutting. If the instability comes from the cutting edge or the machined surface, the update frequency should be increased. (iii) The chip collapse is difficult to remedy by adjusting the hourglass control and reference configuration update, and thus reducing the time step is required. (iv) Influences of various parameters of cutting model are also investigated. The issue of low-cuttingspeed-induced instability is preliminarily resolved. Simulation using a non-rigid tool reveals the violent cutting force vibration due to the toolworkpiece collision.
Numerical study on micro cutting using TLSPH is a research angle that is currently just emerging, and some issues still need to be further addressed. For example, in practice, the cutting speed is usually less than 10 m/s. To perform an experimental validation, an efficient approach without severe increase in simulation time is required to deal with the low-speed-induced instability. In addition, investigation on the effect of tool edge radius is also important in micro machining, and tests on various materials could lead to greater understanding about the universality of TLSPH as a tool for cutting simulation. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.