A Simple FEM for Solving Two-Dimensional Diffusion Equation with Nonlinear Interface Jump Conditions

: In this paper, we propose a numerical method for solving parabolic interface problems with nonhomogeneous ﬂux jump condition and nonlinear jump condition. The main idea is to use traditional ﬁnite element method on semi-Cartesian mesh coupled with Newton’s method to handle nonlinearity. It is easy to implement even though variable coefﬁcients are used in the jump condition instead of constant in previous work for elliptic interface problem. Numerical experiments show that our method is about second order accurate in the L ∞ norm.


Introduction
Parabolic problems with nonhomogeneous flux jump condition and nonlinear jump condition are often encountered in scientific computing and engineering. For example, in Radu et al. [Radu, Meir and Bakker (2004); Hetzer and Meir (2007); Bakker and Meir (2003)], a numerical model is introduced to describe the concentration u of an ion I in aqueous solution Ω aq and adjoining polymeric membrane Ω org . Since the diffusion coefficients of the aqueous solution and of the membrane are different, this an interface problem with interface Γ at the membrane. The equations for the problem are as follows: where Ω = Ω aq Ω org . The diffusion coefficient is given by The ion flux continuity across the interface gives a flux jump condition k aq u aq,x − k org u org,x = 0.
The analyte ion I exchanges with ion J of the same valency at the interface, therefore u R,org u aq − K pot I,J u J,aq u org = u aq u org , which can be rewritten as the nonlinear jump condition σ aq u aq − σ org u org = σu aq u org . ( Given proper boundary condition and initial condition, the above mentioned problem is a one-dimensional nonlinear interface problem with homogeneous flux jump condition. The numerical solution of interface problems has attracted much attention in recent years. The numerical method using non-body-fitted grids for elliptic interface problem was first proposed by Peskin [Peskin (1977)]. After that, various numerical methods are developed, including the immersed interface method (IIM) [LeVeque and Li (1994); Li and Ito (2006); Li (1998)], the immersed finite element method(IFEM) [Li, Lin and Wu (2003); He, Lin and Lin (2011); Gong, Li and Li (2008)] and its partially penalized version [Lin, Lin and Zhang (2015)], the immersed finite volume element method(IFVE) [Ewing, Li, Lin et al. (1999)], the matched interface and boundary method (MIB) [Xia, Zhan and Wei (2011);Yu, Zhou and Wei (2007); Zhou, Zhao, Feig et al. (2006)], the non-traditional finite element method(NTFEM) [Hou, Wang and Wang (2010); Hou, Song, Hou, Li, Wang et al. (2012)], the weak Galerkin method [Mu, Wang, Wei et al. (2013)], the kernelfree boundary integral (KFBI) method [Ying and Henriquez (2007)], and so on. Gradient recovery technique Yang (2018, 2017)] is applied to improve the accuracy of the gradient.
In this paper, we propose a numerical method for a much extended scope of problems, which are two-dimensional parabolic interface problems with nonhomogeneous flux jump condition and nonlinear jump condition. Motivated by the numerical model described by Eqs.
(1)-(3), and inspired by the traditional finite element method(TFEM) for solving elliptic interface problem we proposed in Wang et al. [Wang and Shi (2014)], we propose a new numerical method for solving the two-dimensional parabolic interface problem with nonhomogeneous flux jump condition and nonlinear jump condition. We use traditional finite element method coupled with Newton's method to deal with these two special kinds of jump conditions and use Crank-Nicolson scheme to deal with the time derivative. The reason we did not use the NTFEM is because the nonlinear term of this problem is special. It comes from the jump condition. Normally it is effective to use Newton's method to linearize the nonlinear system. However, when the nonlinearity comes from the jump condition along the interface, we need to solve the nonlinear system at the local cell for a time step if we use the NTFEM, and the result of this system will be a linear combination of the unknown solution located at the vertices of the element cell. Newton's method need to have an initial entry at every iteration step. However, the NTFEM can only provide an unknown linear combination, which will cause a serious problem at the linearization part. That is why we choose to use the TFEM. Compared with our recent work on elliptic interface problem with nonlinear jump condition in [Wang, Hou and Shi (2017)], the new contributions include: first, for the parabolic interface problem we need to use Crank-Nicolson's method to handle the time derivative term. This needs to be done carefully with presence of interface. Second, in our previous work, the coefficients σ, σ + and σ − are all constants for simplicity. In this paper, all these are assumed to be variables depending on both time and space, which makes the nonlinear scheme more difficult to derive. Numerical examples show that our method is very effective and stable, it is efficient for matrix coefficient and sharp-edged interface, and can achieve about second order accuracy in the L ∞ norm.

Mathematical formulation
In this paper, we consider a parabolic interface problem defined on an open bounded domain Ω ⊂ R 2 cross with a time interval [0, T ]. Let Γ be an interface that divides Ω into disjoint open subdomains, Ω − and Ω + , hence Ω = Ω − Ω + Γ. Assume that the boundary ∂Ω and the boundary of each subdomain ∂Ω ± are Lipschitz continuous as submanifolds. Since ∂Ω ± are Lipschitz continuous, so is Γ. Therefore a unit normal vector of Γ can be defined a.e. on Γ, see Fig. 1. We seek solutions of the variable coefficient parabolic problem away from the interface Γ given by in which x = (x, y) denotes the spatial variables and is the gradient operator. The coefficient β( x, t) is assumed to be a 2 × 2 matrix that is uniformly elliptic on each disjoint subdomain, Ω − and Ω + , and its components are continuously differentiable on each disjoint subdomain, but they may be discontinuous across the interface Γ. The righthand side f ( x, t) is assumed to lie in L 2 (Ω). Note that f is piecewise defined on two sides of the interface. The flux jump condition along the interface Γ is given by and the nonlinear jump condition along the interface Γ is prescribed by The " ± " superscripts refer to limits taken from within the subdomains Ω ± . Finally, we prescribe initial and boundary conditions for a given function g on the boundary ∂Ω and u 0 on the domain Ω.
In order to derive the weak formulations of the mathematical model defined by Eqs. (4)- (7), we need to define some notations as in Chen et al. [Chen and Zou (1998); Huang and Zou (2002)]. For a Banach space B, define where u H m (0,T ;B) is the norm of H m (0, T ; B) and is given by Now we are ready to derive the weak formulation of Eq. (4). Using the traditional finite element method, with the standard procedure of multiplying by a test function ψ and integrating by parts, we deduce the weak solution: where ψ is in H 1 0 (Ω). Below the weak formulation is employed to solve the parabolic interface problem. First divide the time interval [0, T ] into n t equally spaced subintervals [t n−1 , t n ], n = 1, 2, · · · , n t , with t n = n∆t, ∆t = T /n t . According to the Crank-Nicolson method in Wang et al. [Wang and Shi (2015)], the following semidiscrete problem of Eq. (4) is derived: Further, the weak formulation of the semidiscrete problem reads 3 Numerical method In this section, a numerical method is implemented to solve the nonlinear system described in Eqs. (4)- (7). The grid used in this paper is the semi-Cartesian grid introduced in Wang et al. [Wang and Shi (2014); Wang, Hou and Shi (2017)  For any cell T o,p,q , o, p, q ∈ {1, 2, · · · , N }, where N is the number of grid points, we define two sets bellow: Then the body-fitting triangulation T h of the domain Ω is obtained as and Γ h is used to define the interface segment on the interface cell T . The set T R of regular cells is defined by Note that for body-fitted grid the interface does not cut through the triangle, it only touches one vertex or one side of the triangle. Therefore there are only two kinds of interface cells, see Fig. 3. It is this desirable property that ensures our method straightforward to implement.
(a) (b) Figure 3: Interface cells Now we introduce the canonical piecewise linear finite element space X h on the triangulation T h : where P 1 (τ ) denotes the space of linear functions on τ . We shall construct an approximate solution to the interface problem by using the space X h , but taking into account the jump conditions. First note that the flux jump condition [(β∇u) · n] = b along the interface Γ is already absorbed into the weak formulation. Hence, the most important and challenging part is the nonlinear jump condition prescribed in Eq. (6). The approximate solution is piecewise linear on both subdomains Ω − and Ω + , but discontinuous along the interface Γ.
For grid points that are involved in two different solutions, like point o in Fig. 3(a) or points o, q in Fig. 3(b), we have the following definitions: Next we have notations defined as below [He, Lin and Lin (2011); Wang, Hou and Shi (2017)]: With these notations, a globally piecewise linear approximation u h can be defined as Before deriving the discrete formulation for the weak Eq. (10), for simplicity consideration, let Then we have the following method: Although we can get this discrete Eq. (12), it is a nonlinear equation. In this paper, we use the Newton's method to linearize and solve the problem. Unlike other nonlinear problems having a nonlinear term with specific definition in the governing equation defined on the whole domain, the nonlinear part of this problem comes from a local place: it located at the interface. Therefore we need to separate the interface points from the interior points. Then the nonlinear term needs to be coupled in the governing equation such that we can use Newton's method to linearize the system. Since σ is nonzero in this paper, the solution jump is nonhomogeneous, that means there will be two solutions u ± and two governing equations defined on the point located at the interface. When we are calculating the Jacobian matrix of the nonlinear system, we need to figure out it is with respect to which equations and which unknowns. This is a tricky issue for the linearization procedure. Based on the above discussion, and take the nonlinear jump condition into consideration, define Hence the iteration method can be written as where l refers to the iteration step of the Newton's method, u l is the unknown vector, such that u l+1 = u l + δu l , and ∂R(u l ) ∂u l is the Jacobian matrix. Let with the notations defined above, and from Eqs. (13) and (14), we have the following method: , vol.119, no.1, pp.73-90, 2019 Method 3.2. Find a solution u h such that for all ψ j ∈ X h , we have

CMES
∇u,bo +I n u,in\oi + I n u,oi + I n u,bo − ∆t 2 I n ∇u,in\oi + I n ∇u,oi + I n ∇u,bo −∆tI b + ∆tI f .
Remark 3.1. In our implementation, the integrals are computed with Gaussian quadrature rules. For each triangular cell, the midpoint of each edge is denoted by p ij . In numerical computation, the average of three f (p ij ) in each cell is utilized.

Numerical experiments
In this section, we use three specific examples to testify the efficiency of our method (Method 3.2) for solving the two-dimensional diffusion equations with nonlinear interface jump conditions. In all numerical experiments below, the level-set function φ, the coefficients β ± , σ ± , σ, and the solution u + in Ω + are given. Hence Since the solutions are given, g and u 0 are obtained as a proper Dirichlet boundary condition and an initial condition. The right hand side functions of the PDEs are evaluated using the left hand side and the true solution. The termination condition for the Newton iteration method in this paper is min(∆x,∆y) 10 9 . All errors in solutions are measured in the L ∞ norm in the whole domain Ω.

Conclusion
We developed a simple, efficient numerical method for matrix-valued coefficient parabolic equations coupled with nonhomogeneous flux jump condition and nonlinear jump condition. We used TFEM coupled with Newton's method to deal with these two special and complicated jump conditions and used Crank-Nicolson's scheme for time derivative. The coefficients in the jump conditions are variable instead of constant in our previous work. We found that our previous work on NTFEM is not suitable treating such problem  due to the linearization part, and this is not a trivial extension of our previous work on TFEM for elliptic interface problem due to the time derivative and the coefficients not being constant in the jump conditions. This is not a trivial extension of the previous work on the one dimensional nonlinear interface problem either, because the geometry in two dimensions is much more complicated than one dimension. As far as we know, a solver for this problem has not been reported in the literature. Extensive numerical experiments indicate that the method is about second-order accurate in the L ∞ -norm, and it is stable even for problems with very complicated interfaces.