A scalar prototype problem of deep abrasive drilling with an infinite free boundary: an asymptotic modeling study

A vanishing at infinity solution of the three-dimensional Laplace equation is sought in an a priori unknown domain with overdetermined boundary conditions, the right-hand sides of which depend on the unit normal to the free boundary. A perturbation analysis of the nonstandard free boundary problem that models deep abrasive drilling has been performed under the assumption that the free boundary is close to the surface of a given semi-infinite cylinder, the longitudinal position of which depends on the boundary data, and the leading-order asymptotic solution has been developed.


Introduction
Free and moving boundary problems [1,2], where the problem domain is a priori unknown are encountered in different areas of science and engineering [3], including laser drilling of metals [4,5], crack propagation and delamination in composite materials [6,7], in Hele-Shaw flow [8,9], and elastic and plastic torsion [10,11].
An example of stationary free boundary problem for a harmonic function u can be formulated as an overdetermined problem with two boundary conditions: on a free boundary Γ . A number of approaches for numerical solution of stationary free boundary value problems have been developed, in particular, using trial methods [12,13], shape optimization [14], cost functional minimization [15], Newton's method [16], and level set techniques [17].
We note that the so-called Serrin-type overdetermined boundary value problems for Poisson's equation with the boundary conditions (1) with constant right-hand sides have been studied for both interior [18,19] and exterior [20,21] domains. It has been established that if g = const and h = const, then Γ is a sphere and u is radially symmetric about the center of the sphere. The case of a Bernoulli problem with non-constant gradient boundary constraint, that is (instead of the second boundary condition (1)) |∇u| Γ = h, (2) was studied in [22] under the assumption that its right-hand side depends on the unit normal to the unknown boundary.
In the present paper, we consider one nonstandard free boundary problem, where the right-hand sides of Eqs. (1) depend both on the unit normal vector to the surface Γ and its position. The problem was motivated by recent mathematical modeling studies [23,24] of deep drilling mechanics. Since the deep drilling problem considered in [23] is strongly nonlinear, we formulate its scalar prototype to highlight the novelty of the free boundary formulation. The term "abrasive" refers to the Neumann boundary condition (1) 2 , whereas the direct analogy with the deep drilling problem [23,24], which has been formulated in the framework of linear elasticity, implies the use of the Bernoulli boundary condition (2). We note that the Neumann boundary condition (1) 2 with a constant right-hand side models a nearly constant contact pressure distribution achieved in stationary abrasive wear according to Archard's law (see, e.g., [25,26]). Note also that confusion should be avoided with the term "free boundary" that is usually used to refer to traction-free boundary in solid mechanics.
An important feature of the deep drilling problem is that the free boundary Γ is assumed to be infinite and, in a sense, close to the surface of a semi-infinite cylinder, representing the drilling bit. To the best of the authors' knowledge, no published work, outside the works of Mikhailov and Namestnikova [23,24,27], has attempted to study an overdetermined boundary value problem with the Dirichlet boundary data defined as the sign distance to a given surface.
Another novelty of the deep abrasive problem is that the Dirichlet boundary condition, which is imposed on a proper subset of Γ , additionally contains a scalar parameter, δ, value of which depends on the Neumann boundary data via the unknown function u.
The rest of the paper is organized as follows. In Sect. 2, we formulate the scalar problem of deep abrasive drilling. Its perturbation analysis has been given in Sect. 3 under the assumption that the free boundary is close to the bit surface in some displaced position, corresponding to the level of imposed external loading. In Sect. 4, we linearize the free boundary problem and present the leading order asymptotic solution. A simple example, where some analytical results can be derived based on the known solutions, is considered in Sect. 5. Finally, in Sect. 6, we outline the discussion of the presented analysis.

Free boundary problem formulation
We consider stationary abrasive drilling of a semi-infinite bore-hole R 3 \Ω, spreading to x 3 = +∞ in an infinite space R 3 (see Fig. 1). Let the x 3 -axis of the Cartesian coordinate system x = (x 1 , x 2 , x 3 ) coincide with the borehole axis. For simplicity, we may assume that the bit is axially symmetric and, therefore, the domain outside the bore-hole Ω is axially symmetric as well. Let B denote the domain occupied by the bit in the unloaded state. We assume that the bit's boundary ∂ B is semi-infinite and is composed of two parts: a vertical cylindrical surface ∂ B and a curvilinear bottom surface ∂ B ∪ .
In the loaded state, the bit is moved into contact with the bore-hole surface. Let δ > 0 be the bit's vertical displacement in the opposite direction of the x 3 -axis, so that B δ and ∂ B δ = ∂ B δ ∪ ∂ B δ ∪ will denote the domain occupied by the displaced bit and its boundary (see Fig. 2). Further, in the loaded state, the bit's bottom surface ∂ B δ ∪ represents the rupture front ∂Ω δ ∪ , which is characterized by the boundary condition Here, n(x) is a unit outward (i.e., directed inward the bore-hole) normal vector, μ and p are positive constants. Note that at the center of coordinates (denoted as point O), we have n(O) = 1. For the sake of simplicity, we assume that in the loaded state, the vertical bore-hole surface ∂Ω = ∂Ω\∂Ω ∪ comes into full contact with the bit's vertical surface ∂ B δ , i.e., ∂Ω δ = ∂ B δ , and the contact is stress-and damagefree, so that In the loaded state, the bore-hole surface ∂Ω δ coincides with the bit surface ∂ B δ , whereas in the unloaded state, there is a gap between the surfaces ∂Ω and ∂ B. Hence, let g(x) denote the gap measured along the normal to the boundary ∂Ω in the inward direction with respect to the bore-hole. It is clear that the gap between the lateral surfaces ∂Ω and ∂ B is invariant with regard to any downward vertical translation of the cylindrical bit surface ∂ B , while the gap between the bottom surfaces ∂Ω ∪ and ∂ B δ ∪ will be, respectively, affected. So, let the gap between the surfaces ∂Ω ∪ and ∂ B δ ∪ be denoted by g δ (x). The unknown function u(x) is required to satisfy the Laplace equation inside the domain Ω, i.e., and to vanish at infinity, however, the problem domain itself is supposed to be unknown.
To complete the problem statement, we impose the second boundary condition on each of the boundary parts ∂Ω ∪ and ∂Ω as follows: Observe that by analogy with contact problems, one may additionally impose the so-called equilibrium condition with a constant P that, in view of Eq. (3), should be related to p. Indeed, the substitution of (3) into Eq. (9) leads to the relation where |ω ∪ | is the area of the projection ω ∪ of the surface ∂Ω ∪ on the plane x 3 = 0. Further, we note that there is some ambiguity in the definition of the bit surface ∂ B in the unloaded state. Here we will follow the convention usually used in contact problems that the surfaces ∂ B ∪ and ∂Ω ∪ coincide at a single point O lying on the vertical x 3 -axis. So, using the arbitrariness in the vertical position of the center of coordinates, we will assume that as the point O has been chosen as the coordinate center. Equation (11) immediately implies that where δ is the bit displacement. To simplify the calculations of the displaced gap function g δ under the assumption of relatively small displacement δ, the following approximation can be utilized: Finally, the partition of the bore-hole boundary ∂Ω into two parts ∂Ω ∪ and ∂Ω will be determined by the conditions so that the contribution to the total contact force P [see Eq. (9)] from the surface loading p, as it is specified by Eq. (3), will be positive on the entire rupture front ∂Ω ∪ . Equations (3)-(14) constitute a free boundary problem, where positions of the surfaces ∂Ω ∪ and ∂Ω must be determined in the process of solution.

Remark 1 Let H(x) be the Heaviside function, that is
Then, in view of (13) and (14), the boundary conditions (3), (4) and (7), (8) can be represented in the form (1) as follows: Here, Γ = ∂Ω is the entire boundary of the hole.

Perturbation analysis of the free boundary problem
We assume that the unknown bore-hole surface ∂Ω is close to the bit surface ∂ B ∪ , so that the gap g(x) is relatively small compared to the radius of the cylindrical surface ∂ B . Now, by applying a perturbation technique [28], we will move the boundary conditions from the surface ∂Ω to the surface ∂ B. With this aim, we introduce a local orthogonal curvilinear coordinate system on ∂ B, denoted as α = (α 1 , α 2 ) with the Lamé coefficients A 1 (α) and A 2 (α), and a unit outward (with respect to the bit domain B) vector ν(α) (it is assumed that ν = e 1 × e 2 ). Let also γ (α) denote the gap between the surfaces ∂ B and ∂Ω measured from the surface ∂ B to the surface ∂Ω in the direction of the vector ν(α). Then, in the local coordinate system (α 1 , α 2 , ν), the unknown surface ∂Ω can be parameterized as follows: By accounting only for the first-order perturbation effects, we can replace the boundary condition (8) with the following: Thus, by approximating the left-hand side of Eq. (18) as we linearize Eq. (18) as follows: Now, we can parameterize the surface ∂Ω by using the radius-vector r(α) of the surface ∂ B as follows: Differentiating the radius-vector R(α) with respect to the curvilinear coordinates and making use of the Rodrigues theorem, we obtain where R 1 and R 2 are the two principal curvature radii, e 1 and e 2 denote the unit coordinate vectors of the coordinate system (α 1 , α 2 ). Then, the gradient of a scalar field φ in the orthogonal coordinate system (α 1 , α 2 , ν) is defined by the formula Correspondingly, the Laplacian of a scalar field φ is defined as follows: Here, h 1 , h 2 , and h 3 = 1 are the Lamé coefficients of the coordinate system (α 1 , α 2 , ν), i.e., Observe also that formula (21) in the vicinity of the base surface ∂ B can be approximated as where ∇ α is the gradient operator in the curvilinear coordinate system on ∂ B.
Since the two vectors (20) (i = 1, 2) are tangent vectors to the surface ∂Ω , we can define the limit normal vector to the surface ∂Ω as follows: Recall that the normal vector n has been defined previously to be directed inward the bore-hole, i.e., in the opposite direction with respect to the normal vector ν. Substituting (20) (i = 1, 2) into the second equation (24) and making use of the assumptions |γ | min{R 1 , R 2 } and |∂γ /∂α i | 1 (i = 1, 2), we derive the second-order approximation Now, we are in a position to move the boundary condition (4) from the surface ∂Ω to the surface ∂ B . First, by definition, we have ∂u ∂n ∂Ω = n · ∇ x u r(α) + νν(α) where the gradient operator ∇ x is defined by (21). Second, by utilized formulas (23), we readily obtain and the subsequent application of the Maclaurin approximation yields Third, substituting the first-order approximations (25) and (27) into Eq. (26), we find Thus, the boundary condition (4) can be linearized as follows: In the same way, we replace the boundary condition (3), which is imposed on the surface ∂Ω ∪ , with the following one on the bit bottom surface: It is to emphasize here that p and μ are constants. Finally, let γ δ (α) denote the gap between the surfaces ∂ B δ ∪ and ∂Ω ∪ measured from the surface ∂ B δ ∪ in the direction of the vector ν(α). Then, in the framework of the first-order approximation, we replace the boundary condition (7) on the free boundary ∂Ω ∪ with the following, which is imposed on the bit bottom surface in the unloaded state: Now, the application of the perturbation method to Eq. (31) yields the linearized boundary condition Thus, on each part of the known boundary ∂ B, we again have two boundary conditions, one of which should be used for determination of the gap function.

Solution of the linearized free boundary problem
First of all, we emphasize that Eqs. (19), (29), (30), and (32) have been derived by neglecting second-order terms. Therefore, by making use of Eqs. (29) and (30), and staying within the limits of their accuracy, we can simplify Eqs. (19) and (32), respectively, as follows: The second leading-order solution of the linearized free boundary problem, which is composed by the Laplace equation and the boundary conditions (29), (30), (33), and (34), can be constructed by satisfying the limit forms of Eqs. (29) and (30), that are The boundary value problem (35), (36) has a unique solution satisfying the following asymptotic condition: At the same time, the leading-order approximations for the gap functions on the surfaces ∂ B and ∂ B ∪ , in light of Eqs. (33) and (34) are given by where the right-hand sides are known from the solution u 0 (x) of the Neumann problem (35)- (37). Some comment is necessary on the use of Eq. (39) instead of Eq. (34). Observe that Eqs. (35) and (36) are linear, and, therefore, the function u 0 (x) will be proportional to the dimensionless ratio p/μ. On the other hand, Eq. (38) has been derived under the assumption that γ 0 is relatively small (which turns out to be of the same order as u 0 ). Therefore, when considering the leading-order approximation, the second term on the left-hand side of Eq. (34) can be neglected.
As can be seen from Eq. (39), there arises a difficulty in determining the gap function γ 0 on the bit bottom surface ∂ B ∪ . This problem can be solved as follows.
First, by recollecting Eqs. (11) and (12), we can write Therefore, Eqs. (39) and (40) yield where δ 0 is the leading-order approximation for the bit displacement. Let r(α 1 , α 2 ) be the radius-vector of the surface ∂ B ∪ . Then, the displaced surface ∂ B δ ∪ can be defined by the radius-vector r(α 1 , α 2 ) − δi 3 , where i 3 is the x 3 -axis unit vector. It is to emphasize that, because ∂ B δ ∪ is obtained from ∂ B ∪ by translation along the vertical, the normal vector ν(α 1 , α 2 ) of the surface ∂ B ∪ will serve as the normal vector for the surface ∂ B δ ∪ as well. Therefore, the unknown surface ∂Ω ∪ can be approximated by the radius vector where δ 0 and γ δ 0 (α) are given by Eqs. (41) and (39), respectively. Thus, since the parametrization (42) is known, the gap γ δ 0 (α) between the surface ∂ B ∪ and the surface defined by Eq. (42), which is measured from the surface ∂ B ∪ in the direction of the vector ν(α) can be determined by direct calculations.
Finally, the first-order approximation can be obtained by solving the boundary value problem and using the equations where the zero index refers to the leading-order approximations.

Example of cylindrical bit
Let B be a semi-infinite circular cylinder of unit radius, so that the surface ∂ B ∪ coincide with the circular part 0 ≤ r ≤ 1 of the plane z = 0. For this bit configuration, the leading-order approximation problem (35)- (37) represents a special case of the external Neumann problem for a semi-infinite cylinder studied in [29], i.e., where f 0 is a constant, The function U 0 (r, z) is represented as follows [29]: Here, J 0 (x) and H (1) 0 (x) are Bessel functions of the first and third kind. In turn, the functions A 0 (λ) and B 0 (λ) depend on the solution ω 0 (x) of the integral equation via the function as follows: Here, I 1 (x) and K 1 (x) are modified Bessel functions of the first and second kind, Re and Im denote the real part and the imaginary part. Let us rewrite Eq. (52) in the form where we have introduced the notation Then, the use of Taylor and asymptotic expansions of the modified Bessel functions yields Observe that the solvability of singular integral equations (55) of type (56) was studied in detail in [30,31]. To solve Eq. (55), we use the Bubnov-Galerkin method and look for the solution in the form of the expansion where and L n (x) are the Laguerre polynomials defined as The substitution of (57) into Eq. (55) results in the infinite system of linear algebraic equations x + y dy.
In view of (56), all integrals (59) are convergent. The infinite system (58) is then solved by truncation. Observe that by means of formula (3.353.5) from [32] the iterated improper integral (59) 2 can be evaluated as follows: Here, Ei(x) is the exponential integral. We note also that the substitution of (57) into (53) leads to the integrals ∞ 0 x m e − px (λ 2 + x 2 ) −1 dx, which can be evaluated in terms of the sine and cosine integrals Si(x) and Ci(x). In the problem under consideration, the following quantities are of interest: Here, according to Eqs. (53) and (54), we have Note that in the derivation of formula (63), we made use of the known relations H (2) The results for the leading-order approximations for the relative displaced gap functionγ δ = (μ/ p)γ δ and the relative gap functionγ = (μ/ p)γ under the bit's base are shown in Fig. 3. It is interesting to observe that, since n = e 3 and ν = −e 3 , the bottom surface ∂Ω ∪ of the bore-hole is concave (i.e., curving inward).

Discussion
Let us comment on the derivation of the linearized free boundary problem (35)- (37), and in particular to answer the question of why we have chosen the boundary conditions (36). Consider, first, the lateral boundary ∂ B , where we have got two equations (29) and (33). By excluding the unknown function γ , we arrive at the equation Now, by accounting for Eq. (5) and the formula where H = R −1 1 + R −1 2 /2 is the mean curvature, and we can rewrite Eq. (64) as follows: Thus, it is readily seen that the second equation (36) has been obtained as a linearization of Eq. (65). Further, the derivation of the nonlinear boundary condition from Eqs. (30) and (32) has been complicated by the presence of γ δ along with γ . For shallow surfaces, one can make use of the approximation γ δ γ + δν · i 3 and the equations γ δ = u and γ δ (O) = −δ to obtain γ u + u(O)ν 3 . The latter equation allows us to rewrite Eq. (30) in the form Again, the linearization of Eq. (66) leads to the corresponding boundary condition (see Eq. (36) 1 ) for the leadingorder approximation.
By analogy with the contact problem of drilling, which was formulated in [27], we can consider the following rupture front equation instead of Eq. (3): Here, ∇ s u is the gradient of the function u in the local orthogonal curvilinear coordinate system s = (s 1 , s 2 ) on the surface ∂Ω ∪ , σ c is a positive constant. The boundary condition (67) should be supplemented by the inequality which ensures the load transfer from the bit to the bore-hole surface. Hence, in light of (68), Eq. (67) can be rewritten as follows: The boundary condition (69) can be moved from the surface ∂Ω ∪ to the surface ∂ B ∪ , and in the light of leadingorder approximation one can arrive at the following result: Observe that the nonlinearity of Eq. (70) is due to the intrinsic nonlinearity of the boundary condition (67). Note that the solvability of the so-called exterior gravitational problem has been studied in [33,34].
Another modification of the considered problem formulation that can be introduced by analogy with the deep drilling problem [27] is the assumption that only the bore-hole bottom surface ∂Ω ∪ is a free boundary, whereas the lateral surface ∂Ω is assumed to be cylindrical with the radius determined by the homogeneous Neumann boundary condition (4) and the continuity with the adjacent part of the surface ∂Ω ∪ . In this case, it can be shown that the linearized problem equations (35)-(38), (40), and (41) still hold. Note also that, by employing the terminology used in [35], such free boundary problem can be termed as partially overdetermined.
Finally, in the same way as it was done in Sect. 3, it can be shown that the boundary conditions from the free boundary Γ can be moved onto any nearby surface Γ j . This simple idea, together with a approximate rule for iteration of the boundary, forms the basis of numerical methods for free boundary problems. Assuming that u j (x) solves the Neumann problem with the boundary condition −μ ∂u j ∂n (x) = pH n j 3 (x) , x ∈ Γ j , the new surface Γ j+1 can be found (see, e.g., [36]) by moving Γ j in its normal direction n j so that u j x j+1 ≈ u j x j + ∂u j ∂n x j d j x j , where x j+1 = x j + n j d j (x j ), and the distance d j (x j ) is determined by the requirement that x j+1 equals the