Modelling and Simulation in Engineering

Volume 2013 (2013), Article ID 716383, 14 pages

http://dx.doi.org/10.1155/2013/716383

Research Article

An Implicit Algorithm Based on Continuous Moving Least Square to Simulate Material Mixing in Friction Stir Welding Process

Abdelaziz Timesli,1,2 Bouazza Braikat,¹ Hassane Lahmam,¹ and Hamid Zahrouni²

¹Laboratoire d’Ingénierie et Matériaux LIMAT, Faculte des Sciences Ben M’Sik, Université Hassan II Mohammedia-Casablanca, Sidi Othman, BP 7955, Casablanca, Morocco

²Laboratoire d’Etude des Microstructures et de Mécanique des Matériaux LEM3, CNRS UMR 7239, Université de Lorraine, Ile du Saulcy, 57045 Metz Cedex 01, France

Received 25 June 2013; Revised 16 September 2013; Accepted 30 September 2013

Academic Editor: Jean-Michel Bergheau

Copyright © 2013 Abdelaziz Timesli et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

An implicit iterative algorithm, based on the continuous moving least square (CMLS), is developed to simulate material mixing in Friction Stir Welding (FSW) process. Strong formulation is chosen for the modeling of the mechanical problem in Lagrangian framework to avoid the drawback of numerical integration. This algorithm is well adapted to large deformations in the mixing zone in the neighborhood of the welding tool. We limit ourselves to bidimensional viscoplastic problem to show the performance of the proposed implicit algorithm. The results show that the proposed algorithm can be employed to simulate FSW.

1. Introduction

Friction stir welding (FSW) [1] is a solid state welding process used to join two pieces (plates, tubes, spheres, etc.) of aluminum alloys. Joining two workpieces by FSW consists in heat generation mainly due to the shoulder and material mixing thanks to the pin that provokes both extremely high plastic deformation and also a high heat generation.

Numerical simulation of friction stir welding process is important in industry to control the product quality and optimize the fabrication cost and avoid the traditional tests. It has been investigated by several authors considering thermal or thermomechanical framework [2–8]. The difficulty in modeling this welding type resides to large deformations generated by the material mixing. Different formulations have been proposed in Eulerian, Lagrangian, or Arbitrary Lagrangian Eulerian frameworks. Eulerian formulation is appropriate to describe material flow for large deformations thanks to fixed mesh grid but this formulation cannot be used for problems involving free surfaces. Lagrangian formulation is well adapted for history dependence in solids and for simulation of material flow with free surfaces but mesh distortions require special procedure as remeshing. The ALE formulation has been proposed to benefit from both advantages of Eulerian and Lagrangian formulations [3]. In this technique, one has to describe the motion of the mesh and material particles separately with respect to a reference domain. ALE has been applied successfully in many fields and particularly in metal forming processes. Despite the intensive development of these techniques, material mixing in FSW process remains very difficult to achieve numerically and especially when using finite element method [2, 3].

An alternative solution for the simulation of this process is the use of meshless methods. Indeed, these techniques are developed for several decades to solve partial differential equations. They aim to avoid difficulties related to the mesh discretization and they have been intensively developed since the 90s. To our knowledge, smoothed particle hydrodynamics (SPH) was the first meshless method presented in 1977 by Lucia and Gingold and Monaghan for astrophysical problems [9, 10]. Another class of meshless methods called diffuse element method was proposed by Nayroles et al. [11, 12]. This method is based on local and moving least square approximation and a Galerkin discretization scheme. Belytschko et al. developed in 1994 element-free Galerkin method (EFG) [13].

The SPH method based on the notion of particles is best suited for large deformation problems [14] and the MLS approach [15] provides a better consistency for SPH. The first reason is that MLS approach uses monomial basis functions to ensure that a local polynomial expansion best fits the surrounding data points. The second reason is that there is less oscillation of the interpolating polynomial between particles with the MLS which is the source of the tensile instability [16]. In addition, the boundary integrals were dropped in the derivation of the particle equations for SPH, so corrections are needed near boundaries [17].

The aim of the proposed work is to implement an implicit algorithm to simulate material mixing in friction stir welding. This algorithm is based on the continuous moving least square method (CMLS) [18] which combines the advantages of both approaches SPH and MLS. CMLS is formulated here in the Lagrangian framework, using the strong form of partial differential equations. This framework allows one to avoid the drawback related to numerical integration. Thermal and mechanical fields are directly attached to the points used for the domain discretization by the MLS shape functions [19]. In this work, the boundary conditions are imposed with the boundary point collocation method [20], where collocation points are coincident with nodes along the boundary where the solutions are imposed. This allows one to easily simulate large deformations and material mixing around the welding tool.

The layout of this paper is as follows. In Section 2, we present the strong form of the governing equations. In Section 3, we detail the main idea of the CMLS method and we present space and time discretization procedures. In Section 4, we present the solution strategy adopted in this work. In Section 5, numerical application is proposed and a comparison with the industrial code Fluent is given to validate our model. This code is based on an Eulerian formulation and a finite volume discretization. Finally, some conclusions are presented in Section 6.

2. Governing Equations

Friction stir welding is a complex thermomechanical process that requires accurate knowledge of the relations existing between the main parameters of the process such as plunge depth, travel speed, rotation speed, thermomechanical properties, and tool geometry. In many contributions, the mixing zone is considered as a high viscosity incompressible fluid and the flow is obtained using computational fluid dynamics. In some other cases, the deformation is modeled using solid mechanics and numerical methods to solve the resulting nonlinear problem to compute the different variables as velocity, temperature, and so forth. In the present work, we consider a mechanical problem with a constitutive law that depends only on the equivalent strain rate. The resulting problem is then described by the conservation laws including the following.

Mass conservation:

Motion equation:

where is the material density, is the velocity vector, is the stress tensor, and the gravity is neglected. The stress tensor , under this assumption, is given by

where is the hydrostatic pressure, is the deviatoric stress tensor, and is the unity tensor. To satisfy the condition of incompressibility numerically in (1), we introduce, in general, a pressure term penalized by a large factor noted . This means that the incompressibility condition (1) is replaced by a viscous law [21]:

where is the strain rate tensor defined as

FSW is a process that generates large strains in which elastic strains are negligible. The advantage of this hypothesis is that the material can then be modeled as viscoplastic material:

where is the material viscosity. Taking into account (3), (4), and (6), the stress tensor is given by

For the viscosity, we choose a power law given by [22]

where is the temperature, , , and are the material properties, and is the equivalent strain rate given by

In this work, we focus our efforts on the mixing aspect that is difficult to achieve by the conventional methods such as finite element one. So, we assume that the temperature is constant and then the viscosity depends only on the equivalent strain rate . Finally, the governing equations can be summarized as follows:

Problem (10) will be completed by boundary and initial conditions.

3. Numerical Model

3.1. Principal of CMLS Technique

Moving least square (MLS) technique belongs to the family of meshless methods. Since the pioneering work about diffuse elements presented by Nayroles et al. in 1992 [12], several contributions using this technique have been developed in many mathematical and mechanical fields [23, 24]. Continuous moving least square method (CMLS) is a modified version of MLS as shown in [18]. In this method, the approximation of the function is defined as a polynomial of order but with nonconstant coefficients. The local approximation at the coordinate is given by

where is the basis function of the spatial coordinates and is the number of the basis functions. The vector coefficients are obtained by performing a weighted least square fit for the local approximation. This yields the continuous quadratic form

where represents the support domain of point , is the nodal parameter of at , and is a weighting function with a compact support . The weight function plays an important role in the performance of the MLS approximation. The exponential function and spline functions are often used in practice; for example, we present a weighting function in the form of cubic spline [25]:

where . The vector coefficients are determined in each point by minimizing the functional . This minimization leads to the following set of linear equations for :

where is the vector that collects the nodal parameters of field function for all the points in the support domain. is called the weighted moment matrix given by

Solving from (14) and substituting it into (11), the CMLS approximants can be defined as

The approximation defined by (16) is difficult to use due to its integral form. An equivalent discrete form of this equation constructed by a particle approximation is proposed. In this method, the approximation of (16) at an interpolated point of coordinate , in discrete notation, leads to the following approximation:

where is the number of points in the support domain, is the nodal parameter of at and is the volume associated with point .

The number of points , chosen in the support domain, should be sufficient to ensure that the matrix in (15) is invertible, so as to provide the interpolation stability. The choice of depends on the nodal distribution and the number of basis functions. In order to ensure the existence of and a well-conditioned , we usually let [26, 27]. Unfortunately, there is no theoretical best value of , and it has to be determined by numerical experiments.

Recalling the form of the approximation defined in (11), so

Finally, the CMLS shape functions are defined as

The CMLS shape function will be continuous in the entire global domain, as long as the weight functions are chosen properly. Field points included in this support domain are used to perform the CMLS approximation for the unknown function at this point.

In the case where the interpolated point is a particle, the SPH method is used for the calculation of the particle volume [18, 28] as follows:

where is the reference particle.

3.2. Algorithm for the Nearest Neighbors Search

In meshless methods, search for neighboring points in the influence domain is the most expensive procedure in terms of time consumption. The search of neighboring points must be performed for each point in the studied domain and many times during the computation. The naive technique consists in performing a double loop over all collocation points leading to a very expensive algorithm. Several methods are proposed in the literature to improve the performance of these algorithms [17, 29]. In this work, we have implemented an algorithm well adapted to the studied problem. Indeed, we have proposed to construct a fixed grid containing the domain subdivided into subdomains in the radial direction and subdomains along the circumferential direction. The size of each subdomain is noted , where is the subdomain length in the radial direction and is the subdomain length in the circumferential direction (see Figure 1(a)). The choice of and must depend on mesh density of collocation points in the domain and it can be considered as a user parameter. Note that the grid is fixed and each subdomain knows its neighboring subdomains. The procedure defining the neighbors of a given subdomain is executed only one time in the algorithm. For each subdomain, we have then the number of the subdomain and the list of its neighbors. The particular cases of the subdomains located at the domain edge are considered with a specific treatment. Each collocation point is defined by a given number of the point, its coordinates , and the number of the subdomain to which the considered point belongs (see Figure 1(b)). By using a simple mathematical operation, we can identify the number of the subdomain to which the considered point belongs. This operation consists in dividing the current position of the point by the subdomain size in radial or circumferential direction. Retaining only the integer part of the division, we obtain the position in the radial and circumferential directions of the subdomain to which the point belongs. All the points of this subdomain and of the neighboring subdomains are neighbors the considered point. The proposed algorithm does not require large computation time. At each end step, coordinates of collocation points and the corresponding subdomain are updated.

(a)

(b)

Figure 1: (a) Subdomain and its neighbors. (b) Illustration of the neighboring points of the point .

3.3. Space Discretization

The use of meshless methods to approximate the solution of the strong form of the mechanical problem leads to the following advantages: (i) the discretized equations can be obtained directly from the strong forms of governing equations of the problem; (ii) no numerical integration is required; (iii) no mesh is used.

To apply this technique to the problem (10), the shape functions presented in Section 3.1 are used to approximate the velocity vector , in the bidimensional case, at any point using a set of points in a local support domain of the considered point:

where is the matrix of shape functions of point and is the velocity vector of point . Using (5) and (21), the strain rate vector can be obtained using the approximated velocity vector:

where is the matrix of derivatives of shape functions of point . Substituting (21) and (22) in (10) and assembling, we obtain nonlinear system defined by

where is the global mass matrix and is a global matrix dependent of which is the vector that collects the point velocities.

3.4. Time Integration

Using an Euler implicit scheme, the time discretization of the previous problem leads to the following nonlinear system in terms of the new unknown :

where represents the value of the unknown at time .

4. Solution Strategy

Within the so-called iterative method [30], generally, a space discretization of this problem is carried out and the obtained discretized equations are solved iteratively. The nodal unknown is denoted by . In the first iteration, one assumes in the term the value of to start with. Therefore, the right hand side becomes known. Next one iterates until the convergence criterion is satisfied (see Figure 2). This process is controlled by a tolerance . The iterative process will continue until

Figure 2: Global sketch of the proposed implicit algorithm.

The solution of the problem (24) is the limit of the sequence which is the solution of the following recurrent algebraic systems:

5. Numerical Analysis and Applications

In this study, we present a numerical application to show the effectiveness of the proposed algorithm. The mixing zone which is subjected to high strain and strain rate is located in the neighborhood of the welding tool. That is why only a small and circular region is considered to model the mixing process leading to a reasonable computation time. The proposed algorithm is the first step to validate the material mixing procedure in a bidimensional analysis. It can be extended to a thermomechanical analysis by introducing the energy equation.

In these applications, we consider a circular plate of radius made of aluminum alloy with the viscoplastic material proprieties: , , and (see (8)). We perform a time integration by the implicit scheme and we choose , , and . The welding tool is considered rigid with radius . In the proposed work, we consider only the dwell phase and the following boundary and initial conditions:

where is the rotation speed. Results of our approach are compared with those obtained using Fluent software. This code is based on an Eulerian formulation and a finite volume discretization. To validate the results of our algorithm, we choose two equivalent configurations between the two formulations update Lagrangian and Eulerian. The CMLS calculation is performed until a time . In the Eulerian formulation, unsteady calculation with a time is performed. The two calculations use the same constitutive law (8).

5.1. Influence of the Distribution of Points

In this section, to verify the robustness of the proposed numerical algorithm, we propose to study at first some cases of simulation to test the algorithm with different distributions of points in the study domain. These tests allow us to perform a study on the sensitivity to the density of points. We have represented, in Figure 3, the domain configuration for four different distributions of points, where and is the interpoint distance. The results shown in these numerical tests shall be registered at time .

(a)

(b)

(c)

(d)

Figure 3: Different distributions of points for different interpoints distances (a) , (b) , (c) , and (d) .

Figures 4, 5, and 6 show the distribution of the strain rate components , , and obtained with the CMLS algorithm at different interpoints distances , , , and respectively.

(a)

(b)

(c)

(d)

Figure 4: Distribution of component of the strain rate given by CMLS algorithm at different interpoints distances , , , and , respectively.

(a)

(b)

(c)

(d)

Figure 5: Distribution of component of the strain rate given by CMLS algorithm at different interpoints distances , , , and , respectively.

(a)

(b)

(c)

(d)

Figure 6: Distribution of component of the strain rate given by CMLS algorithm at different interpoints distances , , , and , respectively.

We can conclude that the distribution of points has an influence on the distribution of the components of the strain rate as shown by the results. Analysis of the results shows that the algorithm converges from the interpoint distance , which is equivalent to 1780 points. In the following studies, the equivalent configuration to the interpoint distance is used in different simulations because this is the optimum choice.

5.2. Influence of the Support Size

The support size must be large enough to cover several collocation points and to avoid numerical drawback and it must be sufficiently small to avoid obtaining very smooth solution. In this simulation, we choose a fixed parameter . In the general case, this parameter should vary during the simulation to be adapted to the change of the point distribution in the domain. In this study, the variation of the interpoints distances is very small justifying the use of constant.

In this section, the goal is to determine the optimal size of the support (influence domain). For this, we chose a range of values of this size as shown in Table 1. The results show the existence of an optimum equal to for which the maximum relative difference, in comparison with Fluent results, is less than . In the following studies, the optimum choice of support size is used in different simulations.

Table 1: Influence of the parameter on the results of the proposed algorithm.

Support size

Maximum relative difference

5.3. Influence of the Weighting Function

The CMLS method based on strong formulation forces us to treat higher order derivatives; therefore, it is necessary to use higher order weighting function to ensure the numerical computation of these derivatives. We propose here to study the effect of these weighting functions on the effectiveness of the proposed algorithm. Table 2 shows that the quality of the solution obtained by the algorithm is better in using a weighting function of an order three (cubic spline) or more wherein the relative difference is less than .

Table 2: Influence of the weighting function on the results obtained by the proposed algorithm.

Weighting function
Quadratic

spline
Cubic

spline
Quartic

spline

Maximum relative difference
10%
5%
5%

5.4. Validation of the Implicit Algorithm Using the Industrial Code Fluent

In the present section, we propose a comparative and qualitative study of our algorithm with the industrial code Fluent using the optimum choices.

Figure 7 shows the geometry and the space distribution of the material points. Two colors are used to differentiate the two plates to be welded.

Figure 7: Initial distribution of the discretization points.

Figure 8 presents different configurations of material mixing at different times. One can observe the high deformation level around the welding tool. This result, which is difficult to achieve using finite element method, shows the potential of meshless techniques for modeling FSW process.

(a)

(b)

(c)

(d)

Figure 8: Material mixing around the tool at (a) , (b) , (c) , and (d) .

The distributions of and components of the velocity vector are shown in Figures 9 and 10, respectively. Figures 12 and 13 show the evolution of and along the horizontal and vertical sections and the comparison between CMLS and Fluent results (see Figure 11). One can observe that the maximum relative difference is less than confirming the validity of the proposed algorithm.

(a)

(b)

Figure 9: Distribution of component of the velocity vector given by CMLS (a) and Fluent (b).

(a)

(b)

Figure 10: Distribution of component of the velocity vector given by CMLS (a) and Fluent (b).

Figure 11: Horizontal cut view and vertical cut view .

Figure 12: Evolution of component of the velocity vector along the horizontal section .

Figure 13: Evolution of component of the velocity vector along the horizontal section .

Figures 14, 15, and 16 show the distribution of the strain rate components (, , ) with the comparison between MLS and Fluent results. A maximum relative difference of is obtained.

(a)

(b)

Figure 14: Distribution of component of the strain rate given by CMLS (a) and Fluent (b).

(a)

(b)

Figure 15: Distribution of component of the strain rate given by CMLS (a) and Fluent (b).

(a)

(b)

Figure 16: Distribution of component of the strain rate given by CMLS (a) and Fluent (b).

6. Conclusion

In the present work, we have proposed an implicit iterative algorithm to simulate material mixing observed in friction stir welding which is very difficult to achieve using finite element method. It is based on the combination of continuous moving least square spatial discretization, an Euler implicit scheme time discretization, and an iterative method. The performance of this implicit algorithm has been tested on the viscoplastic problem. Results of the proposed algorithm are compared with those obtained by the industrial code Fluent. The test carried out on the considered example establishes that the proposed algorithm is robust and efficient. In this work, we considered sticking conditions between tool and plates to be welded. The contact conditions between shoulder and plates are the next challenges to solve before this algorithm can be applied to FSW. We limited ourselves to bidimensional geometry but the work is in progress to extend the proposed algorithm to three-dimensional one with a complex geometry of the welding tool.

References

W. M. Thomas, E. D. Nicholas, J. C. Needham, M. G. Church, P. Templesmith, and C. Dawes, “Friction stir butt welding,” International Patent Application no. PCT/GB92/02203 and GB Patent Application no. 9125978.8, 1991

S. Guerdoux and L. Fourment, “A 3D numerical simulation of different phases of friction stir welding,” Modelling and Simulation in Materials Science and Engineering, vol. 17, no. 7, Article ID 075001, 2009.

O. Lorrain, J. Serri, V. Favier, H. Zahrouni, and M. E. Hadrouz, “A contribution to a critical review of friction stir welding numerical simulation,” Journal of Mechanics of Materials and Structures, vol. 4, no. 2, pp. 351–369, 2009.

D. Kim, H. Badarinarayan, I. Ryu et al., “Numerical simulation of friction stir welding process,” International Journal of Material Forming, vol. 2, no. 1, pp. 383–386, 2009.

D. Jacquin, B. de Meester, A. Simar, D. Deloison, F. Montheillet, and C. Desrayaud, “A simple Eulerian thermomechanical modeling of friction stir welding,” Journal of Materials Processing Technology, vol. 211, no. 1, pp. 57–65, 2011.

A. Bastier, M. H. Maitournam, K. Dang Van, and F. Roger, “Steady state thermomechanical modelling of friction stir welding,” Science and Technology of Welding and Joining, vol. 11, no. 3, pp. 278–288, 2006.

H. Schmidt and J. Hattel, “A local model for the thermomechanical conditions in friction stir welding,” Modelling and Simulation in Materials Science and Engineering, vol. 13, no. 1, pp. 77–93, 2005.

E. Feulvarch, Y. Gooroochurn, F. Boitout, and J.-M. Bergheau, “3D modelling of thermofluid flow in friction stir welding,” in Proceedings of the 7th International Conference on Trends in Welding Research, pp. 261–266, May 2005.

R. A. Gingold and J. J. Monaghan, “Smoothed particle hydrodynamics—theory and application to non-spherical stars,” Monthly Notices of the Royal Astronomical Society, vol. 181, pp. 375–389, 1977.

L. B. Lucy, “A numerical approach to the testing of the fission hypothesis,” Astronomical Journal, vol. 82, pp. 1013–1024, 1977.

B. Nayroles, G. Touzot, and P. Villon, “The diffuse approximation,” Comptes Rendus de L'Académie des Sciences, vol. 313, pp. 293–296, 1991.

B. Nayroles, G. Touzot, and P. Villon, “Generalizing the finite element method: diffuse approximation and diffuse elements,” Computational Mechanics, vol. 10, no. 5, pp. 307–318, 1992.

T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free Galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994.

A. Timesli, H. Zahrouni, B. Braikat, A. Moufki, and H. Lahmam, “Numerical model based on meshless method to simulate FSW,” in Particle-Based Methods: Fundamentals and Applications, vol. 2 of Computational Methods in Applied Sciences, pp. 651–662, 2011.

A. Timesli, B. Braikat, H. Zahrouni, A. Moufki, and H. Lahmam, “Toward friction stir welding simulation using moving least square technique,” in Proceding of the 2nd International Conference on Friction Stir Welding and Processing (FSWP '12), pp. 119–121, 2012.

G. A. Dilts, “Moving-Least-Squares-particle hydrodynamics—I. Consistency and stability,” International Journal for Numerical Methods in Engineering, vol. 44, no. 8, pp. 1115–1155, 1999.

T. Belytschko, Y. Krongauz, M. Fleming, D. Organ, and W. K. S. Liu, “Smoothing and accelerated computations in the element free Galerkin method,” Journal of Computational and Applied Mathematics, vol. 74, no. 1-2, pp. 111–126, 1996.

G. Shobeyri and M. H. Afshar, “Simulating free surface problems using discrete least squares meshless method,” Computers and Fluids, vol. 39, no. 3, pp. 461–470, 2010.

V. P. Nguyen, T. Rabczuk, S. Bordas, and M. Duflot, “Meshless methods: a review and computer implementation aspects,” Mathematics and Computers in Simulation, vol. 79, no. 3, pp. 763–813, 2008.

G. R. Liu and Y. T. Gu, An Introduction to Meshfree Methods and Their Programming, Springer, New York, NY, USA, 2005.

Y. Tillier, “Identification par analyse inverse du comportement mécanique des polyméres solides: application aux sollicitations multiaxiales et rapides,” Ecole Nationale Supérieure des Mines de Paris, 1998

G. Buffa, J. Hua, R. Shivpuri, and L. Fratini, “A continuum based fem model for friction stir welding—model development,” Materials Science and Engineering A, vol. 419, no. 1-2, pp. 389–396, 2006.

D. Chamoret, P. Saillard, A. Rassineux, and J.-M. Bergheau, “New smoothing procedures in contact mechanics,” Journal of Computational and Applied Mathematics, vol. 168, no. 1-2, pp. 107–116, 2004.

M. Oudjene, L. Ben-Ayed, A. Delamézière, and J.-L. Batoz, “Shape optimization of clinching tools using the response surface methodology with Moving Least-Square approximation,” Journal of Materials Processing Technology, vol. 209, no. 1, pp. 289–296, 2009.

S. Shao and E. Y. M. Lo, “Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface,” Advances in Water Resources, vol. 26, no. 7, pp. 787–800, 2003.

Y. Krongauz and T. Belytschko, “Enforcement of essential boundary conditions in meshless approximations using finite elements,” Computer Methods in Applied Mechanics and Engineering, vol. 131, no. 1-2, pp. 133–145, 1996.

G. R. Liu, Mesh Free Methods: Moving Beyond the Finite Element Method, CRC press, Boca Raton, Fla, USA, 2002.

G. R. Liu, Mesh Free Methods: Moving Beyond Finite Element Methods, Chemical Rubber, Boca Raton, Fla, USA, 2003.

S. Arya, D. M. Mount, N. S. Netanyahu, R. Silverman, and A. Y. Wu, “An optimal algorithm for approximate nearest neighbor searching in fixed dimensions,” Journal of the Association for Computing Machinery, vol. 45, no. 6, pp. 891–923, 1998.

T. R. Taha and M. I. Ablowitz, “Analytical and numerical aspects of certain nonlinear evolution equations. II. Numerical, nonlinear Schrödinger equation,” Journal of Computational Physics, vol. 55, no. 2, pp. 203–230, 1984.