The modeling of interfacial motion is an important task in many problems in science and engineering. In general,there are two categories for moving interface problems. The first category is the Lagrangian tracking method,in which the interface is explicitly represented by the Lagrangian particles (markers),and its dynamics is tracked by the motion of these particles. The second category is the capturing method,in which the interface is implicitly embedded in a scalar field function defined on a fixed mesh,such as a Cartesian grid. For more details,we refer to Ref.˙^[1].

In this paper,based on the idea of the method of the fundamental solution,an asymptotic approximation method is developed to find the minimizers of the Dirichlet integral

Click Browse formula

among all admissible configurations u(X) with the prescribed boundary value u | ∂Ω = 0,and constrained to remain in Ω above a prescribed obstacle ψ(X) with ψ(X) | ∂Ω≤0.

As shown in Fig. 1,the problem (1) describes the equilibrium position of a membrane that is “attached” at level 0 along ∂Ω and restricted to remain above the obstacle ψ(X) when the external force f is imposed. Hence,the problem is called an “obstacle problem”. In the standard Hilbert space H₀¹(Ω),after defining the closed convex set

Click Browse formula

Fig. 1.Deformation of equilibrium[3]

Figure options

we can state the obstacle problem for the membrane in the form of

Click Browse formula

Moreover,by performing the elementary computation,one finds that,if u solves (3),it must also be a solution of

Click Browse formula

The problem (4) is called an elliptic or variational inequality of first kind.

As for the questions like the existence,the uniqueness,and the regularity of the solution of the problem (4),we recommend Ref. ^[2]. We limit our work here only to proposing an approximate numerical method for solving this problem.

Note that some free boundary problems in mathematical physics can be written in the form of the elliptic variational inequality,such as the elastoplastic torsion of a cylindrical bar and a cavitation problem in hydrodynamic lubrication. Here,we make a brief introduction of these problems,and the reader can find more details in Ref. ^[3]. 1.1 Elastoplastic torsion of cylindrical bar

As shown in Fig. 2(a),an isotropic and homogeneous elastic cylinder,which has a simply connected cross section Ω ⊂ R²,is applied with a torsion at the end and with the lateral boundary stress to be free. By means of the so-called semi-inverse method and some assumptions,the cross section Ω is divided into two regions,i.e.,the plastic zone and the elastic zone. Based on the principle of minimum complementary work,similar to the membrane problem with a special obstacle,one can also describe the mathematical formulation of this problem as a variational inequality.

Fig. 2.Elastoplastic cylindrical bar (a) and cavitation problem in hydrodynamic lubrication (b) [3]

Figure options

The hydrodynamic lubrication of the complete journal bearing in the presence of cavitation is one of the first free boundary problems studied in the literature (see Fig. 2(b)),which can be regarded directly as an obstacle type problem. The physical problem consists of the study of the behaviour of a lubricant within the narrow clearance between two eccentric cylindrical bodies in relative motion. In a special case,the cavitation problem for the cylindrical journal bearing can be regarded as a natural consequence of a condition of physical equilibrium. 2 Numerical methods for obstacle problems

The problem (1) is nonlinear and difficult for the computation of approximate solutions. So far,the problem is usually solved by projection methods^[4],active set strategies^[5],multigrid and multilevel methods^[6],level set methods^[7],piecewise linear system methods^[8],discontinuous Galerkin methods^[9],moving obstacle to approach the contact methods^{[10, 11]},etc.

In unilateral obstacle problems,the domain Ω,shown in Fig. 3(a),is divided by the free interface Γ into a coincidence set Ω₁ and a non-coincidence set Ω₂. Let Ω ⊂ R² be strictly convex with a smooth boundary ∂Ω,u(X) = 0 on ∂Ω,and the obstacle function ψ be a strictly concave analytic obstacle satisfying > 0 and ψ < 0 on ∂Ω. Then,the free boundary Γ is an analytic closed curve,and the domain Ω₁ is closed while Ω₂ is open. In Ref. ^[12],a Uzawa algorithm was proposed to solve this kind of free boundary problems after the combination of the method of the fundamental solution (MFS) and the numerical integration.

Fig. 3.Domain Ω divided by free interface Γ into coincidence set Ω1 and non-coincidence set Ω2 (a) and multiply connected domain Ω2 with two exterior pseudo-boundaries Γ1 and Γ2 and collocated source points on them (b)

Figure options

Moreover,the problem (4) can be formulated as a Cauchy problem. In fact,the following relations hold: u(X) = ψ(X) when X ∈ Ω₁,and u(X) > ψ(X) and −∆u(X) = 0 when X ∈ Ω₂. That is,u is the solution of the boundary value problem in Ω as follows:

Click Browse formula

Note that the Cauchy problem (5) is ill-posed even for smooth initial data on the obstacle function ψ. Here,the double boundary condition on ψ holds to compensate the fact that it is not known a priori,and it can be used to solve a few very particular cases^[2].

Based on the problem (5),an asymptotic approximation method is proposed. The trick is first to turn the energy E(u) to be defined on the free boundary Γ as E(u(Γ)) by Green’s formula,and then use the MFS to solve the reduced Laplace equation in the non-coincidence set Ω₂. We minimize the energy E(u(Γ)) by the conjugate gradient method. The free boundary Γ is progressively updated along the opposite direction of the gradient ∇E. In light of the idea of inverse problems,we choose the Neuman boundary condition as the stop criterion. As for the step size α,we adopt the well-known Armijo rule.

Before presenting the algorithm,we introduce the MFS briefly.

The MFS was first introduced in Ref. ^[13] and was continuously developed by numerous mathematicians and scientists in the past 40 years. Now,the MFS has attracted great attention in the scientific and engineering community. It is extensively applied to problems in engineering (see ^{[13, 14, 15]} and the references therein).

We consider the Dirichlet problem for the Laplace equation as follows:

Click Browse formula

We denote by Φ the fundamental solution of the operator −∆,

Click Browse formula

where X ∈ Ω ⊂ R², ∈ R² \ Ω,and denotes the distance between the points X and . The MFS is an algorithm that provides us with an approximate solution of (6). It is described as follows^[16]:

(i) Choose N points ,which are called source points or charge points in the exterior of ∂Ω.

(ii) The approximate solution is sought in the following form:

Click Browse formula

where {l_k} (k = 1,2,· · · ,N) are constants to be determined.

(iii) Choose N points {X₁,X₂,· · · ,X_N},called collocation points,from the boundary ∂Ω.

(iv) Determine {l_k} (k = 1,2,· · · ,N) by the following N equations:

Click Browse formula

From (8),we obtain the following system of linear equations for the unknown coefficients l_k (k = 1,2,· · · ,N):

Click Browse formula

where A is the matrix with the entry . L = [l₁,l₂,· · · ,l_N]^T and b = [ψ(X₁),ψ(X₂),· · · ,ψ(X_N)]^T are the vectors of unknown coefficients l_k (k = 1,2,· · · ,N) and boundary values at the collocation points,respectively.

Since the condition number of the matrix A increases dramatically as the number of collocation points increases,direct approaches to the linear algebraic equation (9) would produce a highly unstable solution^[16]. Fortunately,several regularization procedures were developed to solve such ill-conditioned systems^[17]. In this paper,we select the well-known Tikhonov regularization method to solve the ill-conditioned system. In detail,for the system (9),we define its approximate solution Lκ=argmin,where ||·||₂ stands for the Euclidean norm,argmin{f} means the minimizer of the function f,and κ > 0 is called the regularization parameter and decided by the so-called “L-curve” method. As for the detailed description and implementation of the Tikhonov regularization,we refer to Refs. ^{[18, 19]}. 3 MFS for solving minimization problem

In light of the idea of the MFS and Refs. ^{[12, 20, 21]},as Fig. 3(b) shows,we discretize the interface Γ by collocating at M points (k = 1,2,· · · ,M) and choose M points as the interior source points on the pseudo-boundary Γ₁,where i is the imaginary unit,r_k = λR_k,and 0 < λ < 1 is a parameter denoting the distance between Γ and Γ₁.

As for the curve Γ,the work ^[22] presented a general characterization of the coincidence set Ω₁ without any geometric assumptions about Ω,i.e.,if the obstacle function ψ(X) ∈ C^m,α(Ω), then Γ is a closed C^m−1,α Jordan curve,where m ≥ 2,and 0 < α < 1. In this sense,the above process is available. Similarly,we collocate at M equidistant points on ∂Ω. Then,we choose a contour Γ₂ geometrically similar to the boundary ∂Ω with an appropriate distance d as the exterior pseudo-boundary and collocate M equidistant points on it as the exterior source points.

For simplicity,we introduce some notations as follows:

Click Browse formula

With M collocation points Xj on Γ and M collocation points X_j^′ on ∂Ω,respectively,we get a linear system with respect to the N = 2M expansion coefficients l_k as follows:

Click Browse formula

where b = [ψ₁,· · · ,ψ_M,0,· · · ,0]^T∈ R^N,L = [l₁,l₂,· · · ,l_N],and

Click Browse formula

It is noted that 2M source points locate on Γ₁ when k = 1,2,· · · ,M and on Γ₂ when k = M + 1,M + 2,· · · ,N. Once the coefficients {l_k,k = 1,2,· · · ,N} are determined,we can express the approximate solution by the formula as (7). 4 Asymptotic approximation method and numerical examples

Based on Green’s formula and the Dirichlet condition u = 0 on ∂Ω,we can express the energy E(u) as an integral on the boundary Γ as follows:

Click Browse formula

It is noted that the directional derivative in the normal direction at the points equals (k = 1,2,· · · ,M). Then,by locating the points (k = 1,2,· · · ,M) on Γ to discrete it,we can obtain the discretized energy

Click Browse formula

where R = [R₁,· · · ,R_M]^T ∈ R^M,θ = 2π/M ,and θk = 2π(k−1)/M (k=1,2,· · · ,M).

Moreover,the kth component of the gradient ofcan be written as

Click Browse formula

The terms and are defined,respectively,as

Click Browse formula

Recalling (7),we can find

Click Browse formula

Now,we develop an approximate method for solving the obstacle problem (5).

Given an initial guess of the coincidence set Ω₁,which is denoted by Ω⁽⁰⁾ and satisfies Ω⁽⁰⁾ ⊆ Ω₁,we solve the following boundary value problem in light of the MFS to get the numerical solution u^∗ in the region Ω/Ω⁽⁰⁾:

Click Browse formula

The boundary condition in the Cauchy problem (5) is taken as the strategy to determine whether the iteration terminates,namely,if ≤τ,where ||·||₂ is the normal 2-norm,and τ is a given tolerance,we stop the iteration and take Ω⁽⁰⁾ as Ω₁. Otherwise, we extend ∂Ω⁽⁰⁾ along the opposite direction of the gradient ,i.e.,let ∂Ω⁽¹⁾ = ∂Ω⁽⁰⁾− α,and repeat the above process in the region Ω⁽¹⁾,where ∂Ωⁱ (i=0,1) also denote the unit outer normal vectors of the corresponding boundaries ∂Ωⁱ (i=0,1). We use the Armijo rule to determine the step size α at each iteration. 4.1 Asymptotic approximation method for solving unilateral obstacle problem

Step 1 Choose a small region Ω⁽⁰⁾ ⊆ Ω₁.

Step 2 Solve the system (5) in Ω \ Ω⁽⁰⁾ in light of the MFS to get the solution u^∗. if≤τ,then stop and take Ω⁽⁰⁾ as Ω₁. Otherwise,go to the Step 3.

Step 3 Extend the boundary ∂Ω⁽⁰⁾ along the opposite direction of the gradient,i.e.,let

Click Browse formula

for obtaining the region Ω⁽¹⁾,and turn to Step 2.

We choose Example 4.2 in Ref. ^[7] as the example to show the efficiency of the present method. Let Ω = (−2,2) × (−2,2) and the obstacle function

Click Browse formula

With the Dirichlet boundary condition u(x,y)|_∂Ω = 0,the problem (1) with the above obstacle function ψ has an analytical solution in the form of

Click Browse formula

where ,R = 2,and r^∗ = 0.697 965 148 2 · · · satisfying

Click Browse formula

Three different initial guesses Ω⁽⁰⁾,i.e.,circle,ellipse,and square shown in Figs. 4-6,are considered to investigate the efficiency of the proposed algorithm. All experiments are performed on a ThinkPad workstation (2.53 GHZ Inter(R) Core(TM) i5 CPU) with Matlab codes.

Fig. 4. Numerical result (black solid line) and true free boundary (red solid line) when initial guess of Γ is circle x2 + y2 = 0.25 (a) and pointwise errors between numerical solution and analytical one (b) (order of maximum error is 10−9)

Figure options

Fig. 5. Numerical result (black solid line) and true free boundary (red solid line) when initial guess of Γ is ellipse x2/0.452 + y2/0.352 = 1 (a) and pointwise errors between numerical solution and analytical one (b) (order of maximum error is 10−7)

Figure options

Fig. 6. Numerical result (black solid line) and true free boundary (red solid line) when initial guess of Γ is square with length of side a = 0.4,and its center is located at origin point (a) and pointwise errors between numerical solution and analytical one (b) (order of maximum error is 10−4)

Figure options

First,we fix three parameters λ = 0.8,d = 2,and τ = 10⁻³ for all experiments,respectively. 255 equal angle points are collocated on the free boundary ∂Ω^(k) at each iteration step,i.e,the coordinates of collocation points are (R^j cos θ^j,R^j sin θ^j),where θ^j = 2π/255j (j = 0,1,· · · ,254), and R_i are the polar radii of the corresponding points on the free boundary. We first consider a special case,i.e.,∂Ω⁽⁰⁾ is taken as a circle with its radius R equal to 0.5. The numerical free boundary is shown in Fig. 4(a),and the pointwise errors between the numerical solution and the analytical one are shown in Fig. 4(b). We find that the asymptotic algorithm can efficiently track the free boundary associated problem (1). In fact,we can obtain a circle with a radius R^∗ = 0.697 9,i.e.,the numerical result of Ω₁ is the circle ∂Ω₁^∗ = {(x,y) |x² + y² = 0.697 9²}. After applying the MFS in the region Ω/Ω₁^∗ and defining the error criterion ,where n denotes the number of the nodes for comparison in Ω, we find e(u^∗,u) = 5.36 × 10⁻⁹.

Figure 5 shows the numerical results of the free boundary and the pointwise error between the analytical solution and the numerical one when the initial guess ∂Ω⁽⁰⁾ is taken as an ellipse {(x,y) ∈ R²|x²/0.45² + y²/0.35² = 1}. In this case,we obtain e(u^∗,u) = 1.87 × 10^−7.

A very different geometry of the true boundary Γ,i.e.,∂Ω⁽⁰⁾ = {(x,y) ∈ R²||x|≤0.4,|y|≤ 0.4},is also considered. In Fig. 6,the numerical results of the free boundary and the pointwise errors are shown. We find that the algorithm is also efficient,and e(u^∗,u) = 7.42 × 10⁻³. We also investigate different values of the parameters λ and d on the algorithm and find that it does not have significant impact on the numerical result. 5 Conclusions

Based on the form of an equivalent Cauchy problem,we propose an asymptotic algorithm to capture the corresponding free boundary associated with the particular unilateral obstacle problem. In this particular case,since the displacement function in the noncoincidence region satisfies the Laplace equation with the Dirichlet boundary condition,we can apply the MFS and obtain a high efficient numerical solution. In light of the idea of inverse problems,we take the Neuman boundary condition on the free boundary in the Cauchy problem as the stop criterion during the iteration. Numerical examples testify the feasibility of the proposed method.

Acknowledgements The authors would like to thank Prof. Xiao-liang CHENG (Zhejiang University) for helpful and stimulating discussion.