© T. Sekine, Published by EDP Sciences, 2021

1 Introduction

Surface topography is one of the vital factors in product features. The manufacturing technologies are being continuously required to create a value-added surface in many industries such as automobile, aerospace, and electronic device. Against the background of industrial demands, a lot of contributions have been made from various perspectives [1–4]. Among them, a machined surface in milling is attracting persistent attention from industrial society and researchers [5–11]. The major factors affecting a machined surface has been gradually revealed and considered mostly as cutting parameters, thermal parameters, dynamic parameters, machine tool parameters, tool properties, and workpiece properties [12]. As the technical findings, controlling surface topography in milling enables us to add the surface to a functional property [13,14]. It is also well known that optimizing a path interval in tool path generation can improve not only surface topography but also a balance between cost and product quality in milling [15,16]. The methodologies to determine a single optimum point have been mostly developed based on an algebraic equation derived from geometrical analysis of milling process.

There exist two path intervals on a surface machined by milling. The one is a path interval along the feed direction of a tool, and the other is a path interval along the cross-feed direction of a tool. Path interval determination have mostly focused on the latter [17,18]. Scallop height acts as a dominant factor in the prediction of path interval. Although there are various kinds of tools for milling in these days, typical tool tip geometries have been studied so far in path interval determination with various machining situations [19–23]. Among them, multi-axis flat and filleted end milling requires the elaborate treatments according to a machining situation. There are four general classifications of the situations in multi-axis filleted end milling [24]. The derivation of geometrical relationships tends to be unachievable with the two-dimensional expression for the cutting edge geometry, whereas that of three-dimensional (3D) expression entails mathematical complexities according to intersection problems of 3D geometries [25]. A novel procedure to overcome the complexities was reported with respect to an accurate prediction of path interval in filleted end milling with a tool inclination [26]. The results represented that introducing a reference point was important to estimate a path interval correctly in the 3D consideration.

It is crucial to optimize a milled surface topography in high level through selecting suitable machining parameters in milling. There are a variety of methodologies aimed at a point-based optimization of machining parameters in milling [27–31]. Cutting force is one of the high-interest topics to enhance finish surface quality and to avoid tool failure. Lazoglu et al. proposed a feedrate scheduling technique based on their force model [32]. Habibi et al. also reported a computationally-inexpensive approach minimizing flute engagement to adjust tool orientation for optimizing surface errors in five axis ball end milling [33]. Moreover, cutter/workpiece engagement is an important factor to figure out suitable processing parameters. Zhang et al. investigated an optimization of tool orientation in 5-axis ball-end milling through a model of cutter/workpiece engagement [34]. Chip thickness is theoretically modeled to estimate the other machining factors. Lotfi et al. gave a model of instantaneous undeformed chip thickness [35]. Residual stress in machined surface is also a measure to optimize processing parameters; in addition, Masmiati and Sarhan revealed an effect of cutting parameters on residual stress in machined surfaces [36]. Although there are dominant processing parameters' optimizations from various perspectives, a few comprehensive strategies have been proposed to determine machining conditions in filleted end milling. The one based on path interval determination has been scarcely reported so far in practice.

The purpose of this study is to demonstrate theoretical approaches useful for practical determination of machining conditions affecting machined surface topography in filleted end milling. Tool orientation is intensively investigated in particular. The remainder of this paper is organized as follows. Section 2 describes geometrical description of multi-axis filleted end milling. Then, multi-layer approach will be explained to determine a suitable path interval in Section 3. The other theoretical approaches for several dominant processing parameters will also be given in Section 4. The demonstration and discussion are made in Section 5. Finally, conclusion and future work are presented in Section 6.

2 Geometrical description of filleted end milling

Several coordinate systems is firstly provided to express a machining state of filleted end milling. Henceforth, an axis of a coordinate system is invariably denoted by a normalized vector. This study introduces three coordinate systems as shown in Figure 1. These right-handed coordinate systems are labelled as G, M, and T. G coordinate system comprises X, Y, Z axes, whose components are (1, 0, 0), (0, 1, 0), and (0, 0, 1), respectively. M coordinate system is also defined based on a surface to be machined. MZ is a unit surface normal on a workpiece. MY is orthogonal to MX, and it is oriented to the scanning direction of a tool; moreover, MY can be obtained from the cross product of MX and MZ.

Two sections are instantaneously set to consider the other coordinate system and inclination angles in filleted end milling. The one is a section based on the scanning direction of a tool (i.e., MZ–MX plane), while the other is a section based on the cross-feed direction perpendicular to the scanning direction (i.e., MZ–MY plane). Let these planes be scanning section and cross-feed section, respectively. T coordinate system can be expressed using the two sections. The position of TZ corresponds to that of the rotational axis of a tool, and the direction is oriented to the shank side of tool. As illustrated in Figure 2, tool inclination angle θ denotes the angle between MZ and TZ. TX is geometrically calculated as a vector perpendicular to TZ. The direction depends on the directional relationship between TZ and MZ. Specifically, the angle between TX and MX–MY plane is inevitably set to be θ. Then, TY can be obtained from the cross product of TX and TZ.

The angles θ_s and θ_c are additionally introduced to derive θ. As shown in Figure 3, θ_s is the inclination angle on scanning section, designating the angle between MZ and the orthogonal projection of TZ onto scanning section. Likewise, θ_c is the inclination angle on cross-feed section, indicating the angle between MZ and the orthogonal projection of TZ onto cross-feed section. Counter clockwise rotation is the positive rotational direction of θ_s and θ_c.

Given that M coordinate system coincides completely with G coordinate system for simple modeling of filleted end milling, θ_s and θ_c can be calculated as follows:(1) (2)where TZ_x, TZ_y, and TZ_z are three components of TZ in each axis of G coordinate system. Since TZ is given as a normalized vector,(3)

Hence, equations (1) and (2) can lead to the following components of TZ:(4) (5)

Substituting equations (4) and (5) to equation (3) provides the following formula:(6)

Accordingly, the angle θ can be offered as an inverse trigonometric function with TZ_z:(7)

Figure 4 illustrates the positional relationship of two cutting edges in filleted end mill. Although there exist various kinds of cutting edge geometry, this study focuses on a filleted end mill with two flutes and without a helix angle of cutting edges. With reference to Figure 4, a cutting point on a cutting edge in T coordinate system TPe can be calculated as follows: when d_ae ≥ R − R_cr,(8)and when d_ae < R − R_cr,(9)

It should be careful that clockwise rotation is the positive rotational direction of ξ and ζ. Common end mills have a number of cutting edges, so that ξ depends on the number of cutting edges n_ce. In the case that the angular position of a cutting edge overlaps TX axis, ξ can be derived as a simple form:(10)where m is an assigned number of a cutting edge. From the above explanation, a cutting point on a cutting edge Pe can be calculated through the following formula:(11)where R_T is coordinate transformation matrix, providing a transformation from T to O coordinate system; moreover, Ptc is an instantaneous tool center shown in Figure 4a, depending on a tool path up until a present position. The geometrical description explained in this section assumes a filleted end mill without a helix angle of cutting edges, while it can effortlessly apply to the common, commercial ones through dividing a tool into small disk elements.

Fig. 1 Three coordinate systems in filleted end milling.

Fig. 2 Two instantaneous sections and tool inclination angle θ between MZ and TZ.

Fig. 3 Inclination angles θs and θc.

Fig. 4 Typical tool profile and cutting edges' position.

3 Multi-layer approach

This section describes a novel procedure for theoretical estimation of machined surface topography. Multi-layer approach is introduced as a careful, effectual selection procedure with path interval determination. There are two path intervals associated with a machined surface topography after milling. The one is a path interval in cross-feed direction, and the other is a path interval in feed direction. Let each path interval be path interval L_c and feed interval L_f, respectively. From these definitions, scallop heights can also be given as scallop height h_c and feed mark's height h_f, respectively.

Considering a variation of the tool orientation for determining a path interval L_c, a variety of the tool orientation derived from a tool geometry projected onto an instantaneous section can be classified into four cases according to θ_s and θ_c [24]. The possible situations in 3D geometry can be considered thoroughly based on a tool inclination angle θ. This section mainly deals with filleted end milling under the case that θ ≠ 0.

3.1 Multi-layer concept derived from path interval determination

Two path intervals commonly express a distance between adjacent tool centers in a tool trajectory. Although a path interval in cross-feed direction L_c can be generally calculated as a point-to-point distance, this study focuses mainly on L_c/2. This variable indicates a distance from a tool center to an expediential section located with a predetermined scallop height h_c. Figure 5 illustrates several radii, i.e. R, R_b, and R_cr, of filleted end mill with a tool inclination angle θ. A torus is used to express machining states in the cutting edge geometry of filleted end mill.

A machining situation of cutting edge geometry is given in Figure 6, and a torus is used as a model of cutting edge. The coordinate systems in Figure 6 are identical to the ones in Figure 1. The direction of YT axis is the same as that of Y axis, and tool feed direction coincides with the direction of X axis. There are several planes, i.e. a designed surface, workpiece's surface, and two pseudo planes regarding a scallop height h_c and a feed mark's height h_f. The distance of two planes is situationally changeable in accordance with a variation of designed surface topography. n_s is given as a surface normal in each surface and plane. The direction of n_s changes depending on the orientation of reference surface or plane.

Fig. 5 An assigned torus and several radii.

Fig. 6 An assigned torus and several radii.

3.2 Path interval determination in a cross-feed direction

Figure 7 shows the pseudo-code of path interval determination in filleted end milling [26]. The algorithm focuses on h_c pseudo plane as an exemplified explanation, whereas the computational process is similarly available for workpiece's surface. In the algorithm, a torus has inclination θ and contacts at a point on a designed surface. When a section as an exact circle is cut out from the torus, there exist three positional relationships between a torus section and h_c pseudo plane. The one is an intersection between a torus section and h_c pseudo plane, while there exists a positional relationship without an intersection. The other situation is a single contact point between a torus section and h_c pseudo plane.

The notation P, t, and u indicate a positional vector in 3D space, a tangent vector at each torus section's center, and a directional vector. Each subscript of these vectors is mainly associated with the positional relationships between a torus section and h_c pseudo plane. The subscript base is used for P_base and t_base which express vectors at torus section's center with a contact point between the torus and a designed surface. Moreover, the distance between P_base and a designed surface is completely equal to R_cr. The subscript ap is used to express an arbitrary position. γ is an angular parameter for determining an arbitrary position on a circle with R_b. Rearranging the following formula can provide the initial value of γ designated as γ_ap in Figure 7:(12)

A search range in iterative calculation is introduced along a circle with R_b. P_ap is an arbitrary position on the circle. P_start and P_end are temporarily provided as the starting and ending point of search range, respectively. P_start is a position vector when γ = 0, while P_end is a position vector when γ = 0.5π. A position vector of P_ap can be easily derived through rotating P_base around ZT axis. Moreover, a tangent vector t_ap can be also calculated using t_base in the same manner. Figure 8 illustrates some variables in a torus section at P_ap. P_hp indicates a positional vector on h_c pseudo plane. A direction cosine η between t_ap and h_c pseudo plane can be obtained by the following equation:(13)

A distance d_s between P_ap and P_hp can be given as follows:(14)

The subscript scp means a single contact point between a torus section and h pseudo plane. P_scp and t_scp are obtained through updating P_ap and t_ap in iterative calculation, which is the first do-while statement in Figure 7. ε is set as the convergence condition of iterative calculation. The second iterative calculation in Figure 7 can ascertain an intersection between a torus section and h_c pseudo plane. Through finding out the farthest intersection from the tool center point, the algorithm can provide a suitable path interval L_c/2.

Fig. 7 Computational algorithm for determining Lc/2 [26].

Fig. 8 Some variables associated with a torus section.

3.3 Path interval determination in a feed direction

Feed interval L_f and feed mark's height h_f have been scarcely studied so far. The interval is closely associated with feed per tooth f_t. A tool moves along with X axis, and the direction of TX–TZ plane includes a contact point between a torus and a designed surface. Here, let ω be the rotational angle between X axis and TX–TZ plane. The rotation is about Z axis. Then, the following formula can be made to calculate a feed interval:(15)where f_t is feed per tooth. Moreover, feed mark's height h_f can be mathematically expressed in the same manner to path interval determination in ball end milling.

4 Estimation of the other dominant processing parameters

A milled surface topography can be directly and indirectly affected by the other dominant processing parameters, so that the effectual selection procedures for determining these parameters are also indispensable to obtain a decent topography. With the aim of optimal parameters' selection, this section provides the some theoretical procedures in filleted end milling. Note that the following calculations are unalterably given for one cutting edge during one rotation of a tool.

4.1 Maximum frictional distance

Frictional distance is a critical factor to predict tool wear at the cutting edge. It also has an influence on a machined surface feature. The following formula can provide the maximum length:(16)

In the above formula, α_fd is obtained from an arccosine calculation for a scalar product of u_scp and −TZ. Moreover, a rotational angle around TZ axis β_fd is utilized as the one between TX–TZ plane and a center of torus section having a contact point with a workpiece's surface. Each angle can be extracted as partial results arising inevitably and functionally in the execution of computational algorithm.

4.2 Maximum contact arc length

Contact arc length is defined based on a cutting-related part in a cutting edge, so that it is also an important factor in considering tool surface damage and a machined surface feature. The maximum length can be expressed:(17)

where α_cal is easily identified from an angle calculation using R_cr and the depth of cut d_doc which means a distance between a designed surface and workpiece's surface principally. In contrast, α_hf can be computed from an angle calculation based on a feed mark's height h_f. The execution of computational algorithm including a partial result arising inevitably and functionally can also offer the angle.

4.3 Maximum cutting speed

It is well known that cutting speed has actual impact on both tool surface damage and cutting force. The following expression can be made to estimate the maximum cutting speed:(18)

where s is a spindle speed in cutting process, and α_cs is easily identified using the depth of cut d_doc. The execution of computational algorithm with partial functions can instantaneously compute the angle.

4.4 Average uncut chip thickness

Uncut chip thickness is one of the vital factors having great impact on cutting force. It is undoubted that uncut chip thickness affects tool surface damage in practice since cutting force directly depends on material removal rate in machining process. The following formula enables us to calculate the average value:(19)

where α_uct indicates an angular parameter for determining an arbitrary position on a cutting edge. In addition, β_uct designates an angle between TX axis and GY–GZ plane, and the rotation is about TZ axis. The above formula can estimate an uncut chip thickness at any position on a cutting edge. As a brief estimation, an average uncut chip thickness t_ave was computed in a cutting edge's rotational position where the instantaneous cutting load and area were the largest in both measures, and three angles of α_uct were considered for the calculation. The first angle was obtained based on the depth of cut d_doc. The second angle was set to 5 deg. a priori. Finally, their median was also set to the third angle, and the average calculation was conducted with t_uct obtained using these angles.

5 Demonstration and discussion

This section describes a demonstration of the estimation approaches proposed above. Especially, we investigated influence of tool orientation on each dominant processing parameter in filleted end milling. The discussion will be given after visualizing characteristics of these parameters. The machining conditions used in the demonstration are shown in Table 1. The tool diameter of end mill was set on the basis of a commonly-used size, and the tool had two straight cutting edges. The unit of angle was set as deg. to aid an intuitive understanding.

The influence of tool orientation on θ is shown in Figure 9. In the figure, each curve for respective θ_s in graph legends represented increasing tendency with increasing θ_c in any case. The increasing tendency gradually diminished with increasing θ_s. Moreover, the differences between each curve became smaller with increasing θ_c. The results indicated that the variation of θ depended largely on the larger one of two angles θ_s and θ_c. As a numerical example, when θ_s = 20 deg. and θ_c = 5 deg., the value of θ is 20.52 deg. It was obvious that the effect of θ_c on the value of θ was extremely small.

The influence of tool orientation on L_c is given in Figure 10. In the figure, each curve for respective θ_s in graph legends showed decreasing tendency with increasing θ_c in any case. The decreasing tendency gradually reduced with increasing θ_s, and an asymptotical behaviour could be observed in any curve. Moreover, the differences between each curve became drastically smaller with increasing θ_c. The results clearly expressed that the value of L_c was independent of tool orientation in the case of θ having an angle more than 10 deg. To take a numerical example, when θ_c = 20 deg., the values of L_c were 0.22 mm for θ_s with 5 deg., 0.24 mm for θ_s with 10 deg., and 0.26 mm for θ_s with 20 deg. From the numerical values, the differences in each condition were small. What should be careful here is that these difference would vary according to a pre-determined condition of h_c.

The influence of tool orientation on L_f is displayed in Figure 11. In the figure, each curve for respective θ_s in graph legends represented decreasing tendency with increasing θ_c in any case. Unlike the decreasing tendency of L_c, that of L_f mildly weakened with increasing θ_s. The differences between each curve became larger with increasing θ_c, whereas it can be presumed that there existed the limit of L_f in the case of θ having an angle more than 20 deg. The results revealed that the changing tendency of L_f was completely different from that of L_c. In contrast, it was distinctive that L_f became completely equal in any case without θ_c. In other words, the fact means that these intersection condition between a torus section and h_f pseudo plane coincided perfectly despite the difference in θ_s.

The influence of tool orientation on sl_fd is shown in Figure 12. In the figure, each curve for respective θ_s in graph legends represented decreasing tendency with increasing θ_c in any case. They can be seemingly observed as straight decline. Moreover, the differences between each curve became gradually smaller with increasing θ_c. The results denoted that the values of sl_fd appeared to converge towards a certain value with increasing θ_c. As an example of numerical results, when θ_c = 20 deg., the values of sl_fd were 4.82 mm for θ_s with 5 deg., 4.71 mm for θ_s with 10 deg., and 4.32 mm for θ_s with 20 deg. Numerically, the values of sl_fd are likely to become convergent in the case of θ_c with more than 20 deg.

The influence of tool orientation on sl_cal is given in Figure 13. In the figure, each curve for respective θ_s in graph legends represented decreasing tendency with increasing θ_c in any case. The differences between each curve became larger with increasing θ_c, whereas it could be presumed that there existed the limit of sl_cal in the case of θ_c having an angle more than 20 deg. The results indicated that contact arc length was largely unaltered in any case. They also implied that the intersection conditions between a torus section and workpiece's surface were almost identical despite the difference in θ. In the pre-determined condition of d_doc, the values of sl_cal were 0.95 mm in any case when θ_c = 0 deg. In addition, when θ_c = 20 deg., they were 0.91 mm for θ_s with 5 deg., 0.92 mm for θ_s with 10 deg., and 0.93 mm for θ_s with 20 deg. From the numerical values, there was little difference in each tool orientation.

The influence of tool orientation on v_cs is represented in Figure 14. In the figure, each curve for respective θ_s in graph legends showed increasing tendency with increasing θ_c in any case. The increasing tendency gradually diminished with increasing θ_s. Moreover, the differences between each curve became smaller with increasing θ_c. As a numerical example, when θ_s = 20 deg. and θ_c = 5 deg., the value of v_cs was 47.92 m/min. Likewise, when θ_s = 5 deg. and θ_c = 20 deg., the value of v_cs was 47.92 m/min. In the two machining states, θ was 20.52 deg. It was obvious from the results that the value of v_cs was completely identical under the same θ.

The influence of tool orientation on t_ave is displayed in Figure 15. In the figure, each curve for respective θ_s in graph legends indicated decreasing tendency with increasing θ_c in any case. The differences between each curve became larger with increasing θ_c, whereas it can be presumed that there existed the limit of t_ave in the case of θ_c having an angle more than 20 deg. The results implied that an instantaneous cutting force acting on a cutting edge decreased with increasing θ. With reference to the variation of t_ave, the decreasing rate would be especially prominent in the small value of θ_s.

Wojciechowski et al. reported several relationships between three average forces in cutting process and average uncut chip thickness under some inclination angles in ball end milling [37]. From their results, tangential and radial average forces increased with increasing average uncut chip thickness, whereas axial one decreased in the same condition. In addition, the changing tendencies of three average forces became moderate with increasing an inclination angle. They also provided that there were precipitous variations of average forces within an inclination angle less than 15 deg. In contrast, the variations had a little change when an inclination angle was more than 15 deg. It was also denoted that a width of flank wear depended on an amount of average forces.

Budak and Ozlu investigated some relationships between cutting forces and feed rates in machining process [38]. The results showed that a cutting force increased with increasing a feed rate. In addition, an amount of cutting force depended on a depth of cut.

Bouzakis et al. studied some relationships between tool orientation and surface roughness in ball end milling [39]. Their results revealed that oblique plunge up and down milling were extremely desirable tool orientations in terms of surface roughness and cutting force. Especially, oblique plunge up milling was recommended from the experimental results obtained using some materials. Moreover, surface roughness around 10 deg. could stay as the smallest value.

As one possible conclusion from the results and findings described above, optimal angle of θ is in a range of 15–20 deg. in filleted end milling without θ_c. This condition can achieve both high production efficiency and decent surface feature since two path intervals L_c and L_f are theoretically large. Oblique plunge up milling is also recommended highly under the condition of a tool inclination along a tool feed direction.

Fig. 9 Influence of θs and θc on θ.

Fig. 10 Influence of θs and θc on Lc.

Fig. 11 Influence of θs and θc on Lf.

Fig. 12 Influence of θs and θc on slfd.

Fig. 13 Influence of θs and θc on Lcal.

Fig. 14 Influence of θs and θc on vcs.

Fig. 15 Influence of θs and θc on tave.

6 Conclusions

In this study, theoretical approaches were demonstrated to determine machining conditions affecting machined surface topography in filleted end milling. After geometrical description was explained to model multi-axis filleted end milling, multi-layer approach and the other theoretical approaches were proposed to obtain decent surface topography generated in filleted end milling. The analytical example focusing on tool orientation was given with discussion. As a result, some characteristics of theoretical approaches were revealed with visual evidences. The findings led to one possible conclusion that optimal angle of θ was in a range of 15–20 deg. without θ_c. Oblique plunge up milling is also recommended highly under the condition of a tool inclination along a tool feed direction.

As a future work, the further detailed analysis will be conducted in a wide variety of conditions, and the experimental verification will be made to evaluate validity and applicability of theoretical approaches proposed in this study.

Nomenclature

X, Y, Z : global, stationary coordinate system (O coordinate system)

MX, MY, MZ : M coordinate system

TX, TY, TZ : T coordinate system

θ: tool inclination angle [rad]

θ_s: tool inclination angle in feed direction [rad]

θ_c: tool inclination angle in cross-feed direction

R: tool radius

R: corner radius of cutting edge

R_b: major radius of torus

P_tc : tool center

P_e : a cutting point on a cutting edge

TP_e : a cutting point on a cutting edge in T coordinate system

R_T : transformation matrix

ψ: positional angle in a filleted part on a cutting edge [rad]

ξ: initial angle between TY and each cutting edge [rad]

ζ: rotational angle of a cutting edge [rad]

n_ce: the number of cutting edges

L_c: path interval in cross-feed direction [mm]

L_f: path interval in feed direction [mm]

h_c: scallop height [mm]

h_f: feed mark's height [mm]

P_base, P_start, P_end : a positional vector

P_ap, P_scp, P_hp : in each applicable part

t_base, t_ap, t_scp : a tangent vector in each applicable part

n_s: surface normal

m_scp, u_scp: a directional vector in each applicable part

γ: an angular parameter for determining an arbitrary position on a circle with R_b [rad]

η: a direction cosine between t_ap and an applicable surface or plane

d_s: a distance between P_ap and P_hp [mm]

f_t: feed per tooth [mm/tooth]

ω: the rotational angle between X axis and TX − TZ plane [rad]

Sl_fd: maximum frictional distance [mm]

Sl_cal: maximum contact arc length [mm]

v_cs: maximum cutting speed [m/min]

t_uct: uncut chip thickness [mm]

t_ave: average uncut chip thickness [mm]

α_fd, β_fd, α_cal, α_hf: an angular parameter

α_cs, α_uct, β_uct: at each applicable part [rad]

Acknowledgments

The authors would like to thank the financial support provided by OSG Fund, Shotoku Science Foundation, and the research grant from Faculty of Science and Technology, Seikei University.

References

All Tables

All Figures

	Fig. 1 Three coordinate systems in filleted end milling.
In the text

	Fig. 2 Two instantaneous sections and tool inclination angle θ between MZ and TZ.
In the text

	Fig. 3 Inclination angles θs and θc.
In the text

	Fig. 4 Typical tool profile and cutting edges' position.
In the text

	Fig. 5 An assigned torus and several radii.
In the text

	Fig. 6 An assigned torus and several radii.
In the text

	Fig. 7 Computational algorithm for determining Lc/2 [26].
In the text

	Fig. 8 Some variables associated with a torus section.
In the text

	Fig. 9 Influence of θs and θc on θ.
In the text

	Fig. 10 Influence of θs and θc on Lc.
In the text

	Fig. 11 Influence of θs and θc on Lf.
In the text

	Fig. 12 Influence of θs and θc on slfd.
In the text

	Fig. 13 Influence of θs and θc on Lcal.
In the text

	Fig. 14 Influence of θs and θc on vcs.
In the text

	Fig. 15 Influence of θs and θc on tave.
In the text