Application of two-grid interpolation to enhance average gradient method for solving partial differential equations

Previously presented method of calculating local average gradients for solving partial differential equations (PDEs) is enhanced by using interpolating gridpoints and triangular grids. The interpolating mesh provides finer computational grid, which is then used for solving the PDE. The combined use of the finer interpolating grid together with the original sparser grid is a two-grid method. By comparing the previous application of rectilinear grid for diffusion from initial point concentration to the new triangular two grid method, it was found that the application of triangular two-grid method improves stability of the solution and it provides more rapid convergence to the correct analytical solution.


Introduction
The material microstructure is in many case decisively responsible for the macroscale properties of the material.For this reason, it is very important to obtain the capability of modelling the microstructure evolution in detail.Although the mean field models, such as the ones used in [1,2,3,4], are useful in rapid modelling and optimization of the relative quantities related to the microstructure, they lack the capability of including local effects, which can play important role in the microstructure evolution.For example in [2], it was recently observed that in order to understand and to model stabilization of austenite regions in steel processing, it is necessary to consider local carbon enrichment and the effect of local mechanical strains.Such effects can be considered by using full field models that explicitly simulate the microstructure evolution.
In the physical science based numerical full field modelling of mesoscale phenomena, it is often necessary to solve partial differential equations describing the physical phenomena such as phase transformations [5,6,7], deformation [8,7], diffusion [5,7], fluid flow [9,10], recrystallization [11,12] etc.During processing of materials, deformation of the material is often necessary in order to obtain desired shape.In addition, the deformation can have several beneficial effects to the material properties by affecting the internal microstructure and its development.Also, microscopic phenomena, such as phase transformations, can introduce deformations to the material microstructure [13,8,7] due to the strains caused by the transformations.For these reasons, it is desirable to apply numerical solvers that can handle both regular and deformed numerical grids so that the physical phenomena can be simulated during and after the material deformation.The perhaphs simplest and easiest method to apply for solving PDEs is the finite difference method.Unfortunately, it is difficult to handle complex deformation in the finite difference methods, since the numerical grid in the standard formulation cannot change arbitrarily.To solve the problem in arbitrary geometries, several sophisticated numerical methods for solving partial differential equations exist: the finite element, finite volume, spectral methods, etc.Also, it is worth mentioning in this connection, that I recently [14] used another kind of method for calculating the movement of iso-contours of the field, which can be used for obtaining a solution to a PDE.However, the sophisticated numerical methods usually are also more complex, which requires more effort and time for their implementation.For this reason, it is worthwhile to obtain a method, that is as simple to implement as the standard finite difference method, and which is capable of performing the solution procedure in deformed computational grids, where the gridpoints follow the deformation.These goals were achieved and presented in the previous study [15], where rectilinear grid was used for obtaining the solution by approximating the average local gradients with a plane equation.Unfortunately, it was found that numerical instabilities were introduced to the solution near maximas and minimas of the field.While this could be corrected by using interpolation points near the maximas and minimas, the correction procedure introduced unnecessary complication, which turned out to be difficult to implement in a stable way in the case of more complex differentiation.Particularly, difficulties were encountered in simulating the plastic deformation of a two-phase region where the elastic constants of the phases were different.For these reasons, in the current study, a two-grid interpolation of triangular grid was tested.The triangular grid was chosen, since it allows for a unique linear interpolation between all neighbouring gridpoints in the mesh.It was found that the presented method avoided the numerical instabilities shown in the previous study, and it provided a better and more rapid convergence to an analytical solution for a diffusion from a initial concentration point than the previous study with single quadrilateral numerical grid.

Theory
The fundamental idea of approximating the average local gradient with a plane was presented in the previous article [15].In the current study the method is enhanced by using triangular grid and interpolating points that essentially make the method more stable.
The plane equation is fitted to the neighbouring gridpoints shown in Fig. 1 (ordered grid) and Fig. 2 (disordered grid).The local average gradient is then calculated as the weighted average of the gradients of the neigbouring regions R i , where the areas A i of the regions are used as the weights.
Previously this procedure was applied using a rectilinear numerical grid.In the current study, triangular grid is applied, which makes it possible to construct a finer interpolating grid by line interpolation between each neighbouring gridpoint, as shown in Fig. 3.The triangular grid enables the unique linear interpolation point between each pair of neighbouring gridpoints.Once the finer numerical grid is obtained, it is used for solving the PDE for the next time step.After the solution is obtained, it contains also the values in the original sparse grid.
For completeness, the theory is presented from the first principles, although the original idea of approximating the gradient with a plane equation is the same as in the previous study [15].
As described in the previous study [15], the plane Eq. ( 1) is used to approximate the local average gradient, where u(x, y) is the function value at point (x, y), (x 0 , y 0 ) is the local origin of the plane, where the function value u 0 = u(x 0 , y 0 ).

IC-MSQUARE
The linear coefficients then approximate the average local partial derivatives in each region R i (see Figs. 1 and 2), as described by Eq. (3).
The local average gradient at (x 0 , y 0 ) is evaluated as a weighted average of the gradients of the regions R i , where the area of the region is used as the weight, as described by Eq. ( 4).
The areas A i can be obtained using the cross product of the of the vectors originating from p 0 = (x 0 , y 0 ) and directed to the neighbouring gridpoints p i = (x i , y i ) and p i+1 = (x i+1 , y i+1 ), as described by Eq. (5).
Once the first order partial derivatives of the functions are calculated, the second order partial derivatives can be obtained by differentiating the functions ∂ x u and ∂ y u, as described by Eq. (6).
The numerical solution was tested in a similar way as in the previous study [15].For a test case the diffusion equation ( 7) was solved.Since there exist an analytical solution for a diffusion from an initial point concentration [16], Eq. ( 8), the numerical solution was compared to it.The values M = 0.1 and D = 1 were used.
Similar tests for ordered and disordered grids were made as in the previous study to see the enhancement the triangular two-grid method provided.In the tests, a regular initial grid shown in Fig. 4 and a deformed grid shown in Fig. 5 were used.The number of initial gridpoints was 42 in the x-and y-directions in both cases.

Results and Discussion
It was found that the current approach provided the correct solution to the diffusion from a point source.A surface plot of a two dimensional diffusion simulation is shown in Fig. 6.Especially for sparse grids it was evident that the current triangular two-grid method provided better solution than the previous use of single rectilinear grid.Also, the method provided accurate solution, even if the initial field included a sudden jump at a local maxima.This finding is depicted in Fig. 7, where a) shows the initial field and b) shows the calculated field at time t = 25 compared to the analytical solution provided by Eq. ( 8).The solution is plotted along the line which is oriented along the x direction and passes the origin To confirm that the current method provides correct solution also for deformed grids, the similar checks were made as in the previous study.It was found that the solution was capable of providing correct solution also for the deformed grids.The results for defomed grids are shown in Figs. 8 and 9.

Conclusions and outlook
A triangular two grid method was implemented.The method provides unique interpolation points between two neighbouring initial gridpoints.The solution procedure is performed using the finer grid and then projected back to the original grid.It was found that the tested method provided better stability than the previous use of the single rectilinear grid.The method can be used effectively in regular and undeformed grids.

Figure 1 :
Figure 1: The points of a local regular triangular grid surrounding a gridpoint located in the middle.

Figure 2 :Figure 3 :
Figure 2: The points of a local triangular deformed grid surrounding a gridpoint located in the middle. a

Figure 4 :Figure 5 :
Figure 4: Regular triangular grid which was used in the diffusion calculation.

Figure 7 :
Figure 7: a) Initial condition of the field, line plot along x direction passing through origin.b) Field at time instant t = 25 of the diffusion calculation, line plot along x direction.Regular grid.D = 1, M = 0.1.