Hydration and Durability of Concrete Containing Supplementary Cementitious Materials
View this Special IssueResearch Article  Open Access
Wei Wang, Xiaozu Su, "Algorithm of DRM with Kinetic Damping for Finite Element Static Solution of StrainSoftening Structures", Advances in Materials Science and Engineering, vol. 2017, Article ID 8060987, 7 pages, 2017. https://doi.org/10.1155/2017/8060987
Algorithm of DRM with Kinetic Damping for Finite Element Static Solution of StrainSoftening Structures
Abstract
In order to deal with the divergence and instability due to the illposedness of the nonlinear finite element (FE) model of strainsoftening structure in implicit static analysis, the dynamic relaxation method (DRM) was used with kinetic damping to solve the static increments in the incremental solution procedure so that the problem becomes wellposed. Moreover, in DRM there is no need to assemble and inverse the stiffness matrix as in implicit static analysis such that the associated computational cost is avoided. The ascending branch of static equilibrium path was solved by load increments, while the peak point and the descending branch were solved by displacement increments. Two numerical examples illustrated the effectiveness of such application of DRM in the FE analysis of static equilibrium path of strainsoftening structures.
1. Introduction
The static equilibrium path analysis is important in structural analysis. The finite element (FE) method is an important means for structural analysis. In nonlinear FE analysis of structures, the incremental method is the key [1]. The static equilibrium path analysis of a strainsoftening structure is a tough problem due to the existence of descending branch after peak load point. Concrete is a wellknown strainsoftening material widely used in civil engineering.
Implicit static algorithm (ISA) [2], explicit quasistatic algorithm (EQSA) [2], and dynamic relaxation method (DRM) [3] are the three main methods used for structural static analysis.
According to the increment control type, ISA can be divided into loadcontrol, displacement control, and arclength control. Firstly, for the loadcontrol method, it was proven mathematically that ISA such as Newton method, modified Newton method, and quasiNewton method cannot snap the load peak, search the descending branch of static equilibrium path, and handle local instability because of the nonpositive stiffness matrix or the instability of node displacement [4]. Secondly, comparing to the loadcontrol method, usually a better convergence rate as well as crossing the load peak point can be obtained by adopting the displacecontrol method [5]. However, not only do dissymmetry and disbandment of raise the compute cost [6], but also it is difficult to select the unstable displacement components as displacement control ones [7]. Otherwise, it was reported that displacement control method cannot trace the snapback equilibrium path [8] and cannot be fitting for analyzing strainsoftening structure [5]. Thirdly, the arclength control method [9], which is usually used in analysis of structure static problem, can snap load peak point and track descending branch caused by geometrical nonlinearity [8] but cannot deal with local instability problem caused by strainsoftening or local buckling and cannot snap the bifurcation point [2, 10, 11]. For these reasons, Duan proposed local arclength method which decomposes a structure into damaged zone, failure zone, and dysfunction zone and when doing structure analysis by the local arclength method, the constraint equation only contains the displacement increment of nodes in failure zone and dysfunction zone [11]. From the numerical case studies that have been reported, local arclength method can conduct analysis of strainsoftening structure, but these studies mainly deal with strainsoftening under tension [11, 12]; moreover, the constraint equation needs to be updated during incrementaliterative procedure according to the changing of damaged zone and failure zone.
The literature [2] suggests employing explicit dynamic algorithm (EDA) to do static analysis of structures characterized with local instability; and, in this way, static problem is treated as quasistatic problem in which the inertia effect can be ignored. And we define the EDA as EQSA when it is used to solve quasistatic problem. Some applications of EQSA can be seen in literature [13, 14]. Being different from implicit algorithm, EQSA uses explicit time integral formula to update displacement vector, which can avoid the construction and inversion of [5]. However, there are difficulties in using EQSA. First of all, EQSA is a conditionally stable algorithm; it is that the stable solution of static problems depends on the minimum size of finite elements under the hypothesis that every element has the same physical properties [2]. This condition increases the computational cost of the discrete model with smalldimension elements. Secondly, it is required to reduce the loading rate if a reasonable pseudo static solution is prospected, but it will increase the number of time steps (computational cost). Increasing the loading rate will introduce noise in the loading response (nonnegligible inertia effect) [2] and will even make the computation result completely deviate from the test result [15]. In order to reduce the computational cost, the node mass scaling may also introduce noise into the load response [2]. The determination of viscous damping coefficient has a strong empirical dependence [2].
The explicit time integration algorithm is used to update the vector, and the static solution is obtained by the kinetic energy dissipation of virtual dynamic process which includes viscous damping type and kinetic damping type [16]. DRM can use the virtual mass and kinetic damping, so as to avoid the situation that computational cost of EQSA depends on the small size of the model, and the viscous damping coefficient depends on the analysis experience. The application of DRM in the static analysis of structures began in the 50th and 60th of the last century [17–19], and its application scope gradually expanded in recent years [20–22], and recently it was often used in the analysis of tension structures [3, 23, 24]. However, it is less common in the application and analysis of strainsoftening structures.
In the paper, kinetic damping DRM is used to calculate the static equilibrium path of strainsoftening structure. Firstly, the paper gives the problem description and then describes the process of structural static problem based on DRM in two levels. And two examples and the corresponding results are described. Finally, the discussion and conclusions are given.
2. Description of the Problem
The problem to be solved is a boundary problem in static solid analysis. The condition for material softening is given by [25, 26]where is the stress rate vector and is the strain rate vector. The material strainsoftening makes the equation of the boundary problem no longer elliptic [27], and accordingly the solution lost wellposedness [28]. That introduced time into the problem can make the equation become wellposed [29], but the original boundary problem changes to initial boundary value problem. In order to solve the original boundary problem, we can firstly divide the whole problem solving process into several incremental loads or displacement steps; secondly we attach masses to the nodes to form the diagonal mass matrix and then use DRM to solve the incremental steps. The essence of DRM is that the original solid problem is turned into the one of wellposed damped vibration problem: the dissipation of mechanical energy caused by damping will make the quasistatic displacement solution of the vibration problem approaches the solution of the original problem as . And at any time point , if the displacement vector satisfies the convergence criterion, then it is also the solution of the original static problem.
For the FE discrete model of the solid, all the degrees of freedom (DOFs) can be divided into those with loading and those without loading. Here, we define the DOFs whose nodes are either with load increment or with displace increment as loading DOFs and the remaining DOFs as nonloading DOFs. In the nonlinear FE static analysis of the structure, it is assumed that the boundary of the solid is sufficiently restrained such that it becomes geometrically stable system with no rigid body displacement. The set of nonlinear equilibrium equations to be solved corresponding to all DOFs iswhere is the residual nodal load vector, is the external nodal load vector, is the nodal displacement vector, is internal nodal load vector, and is the total number of DOFs in the system.
Let the index set , the index set , and the index set , where denotes the total DOFs, denote the loading DOFs, denotes the nonloading DOFs, and . For an arbitrary vector with elements, define and . Also define the constitutive status of Gaussian points of elements as the set in which is strain vector, is stress vector, and is status vector of stressstrain relation. The norm of the vector is defined as .
During loading increment, is the unknown and is the known; in this case, the elements in are all not zero, while . During displacement increment, is the known and and are the unknown, while . The vector is computed according to [1], where is geometry matrix, is stress vector, is the total number of elements in , and is the space occupied by the solid.
3. Solution of the Problem
In this paper, the incremental method is used to solve the static problem (2), and DRM is used to solve the increment. In the solution process of DRM, there are several processes in which the kinetic damping is applied. And the method for applying kinetic damping in this paper is to set the nodal velocity vector corresponding to maximum kinetic energy of the system to zero.
3.1. Incremental Solution Method and Solution Process for the Static Problem
Either load increment or displacement increment is used in the global solution of static problem (2) using incremental method. When load increment is used, the load increment ; when displacement increment is used, the displacement increment , where is a reasonably assumed quantity such as the displacement mode of a structure observed during a strong earthquake.
The solution process for increment step is as follows.
At the beginning of the increment : the known quantities are , , and , where is the nodal incremental load vector, is the nodal incremental displacement vector, and the subscripts denote the sequence number of increment steps; here no superscripts are used because the quantities are the converged ones; now either loading increment or displacement increment is used depending on the computation results of increment step (refer to Section 3.3); if , use load increment, , , and , where , , and represent the status of the discrete system before increment Step .
When load increment is used, the known quantities are , , and , and the unknown are and . When displacement increment is used, the known quantities are , , and , and the unknown are , , and .
The unknown given above are all solved by using DRM, so we get , , and .
If the convergence criterion is satisfied, go to the solution process for increment step .
3.2. The Increment Step Solution Method and Solution Process for DRM
Before solving the increment step by using DRM, the time domain is discretized with fixed interval to get the points ih along the time axis. During the solution, the discrete time domain is sequentially divided into series in which is th motion segment with application of kinetic damping. In detail, includes an undamped motion process and a final time point at which the kinetic damping is applied. The beginning point of is the ending point h of , and the ending point of is the beginning point of . Through the series , the dynamic responses , , (the superscripts represent the time point sequence number) approach the static solution , , .
The control process based on dynamic response for solving is as follows:
Determining the status quantities at the beginning time point of h of .
If , then ; the status at this time point is as follows. Firstly, the nodal velocity vector ; secondly, the vector if load increment is used, or if displacement increment is used, in which and ; thirdly, the status vector ; fourthly, the vector is derived from according to (7); fifthly, the vector if load increment is used, or and if displacement increment is used.
If , then the status at this time point is as follows. The vectors , , and are known; secondly, the vector is derived from according to (7); thirdly, the vector if load increment is used, or , if displacement increment is used.
Using the status quantities (, , , , and ) at the beginning time point of as the initial status, the status quantities at can be obtained by solving the undamped motion equationin which is the diagonal mass matrix, and the method for determining and h is shown in Section 3.3. The status quantities , , , at ih time point are solved by using explicit time integration algorithm.
We compute and and the control based on .
The central difference method formulas , based on velocity and acceleration are used to construct explicit time integration algorithm to solve (3). The formulas for are [30]where .
We compute from , , .
If the system kinetic energies , at time points , , respectively, obtained by difference calculation from , , and that had been calculated satisfy the criterion then it is taken that the maximum kinetic energy of the system happens at time point , and, thereby, the kinetic damping is applied, and the nodal velocity vector at time point is simultaneously set to zero (), at the same time, the process ends, which gives the status quantities at time point : , , . After that, the process begins.
If satisfies the divergence criterion when loading increment is used, where is a limiting value, then it is taken that is near or surpasses the peak point of static equilibrium path, and the computation for increment step is restarted and the displacement increment is used.
We compute and and the control based on and .
(221) If load increment is used, ; if displacement is used, and ; the latter guarantees that (see (4)).
(222) The vectorin which the subscripts for increment step sequential number and the superscripts for time point sequential number are both omitted [5]. In (7), e is the element number; is the total number of elements; , , , and are the internal nodal load, geometry matrix, stress vector, and occupied space of element , respectively; , , and are , , and matrix in which is the DOFs number of elements ; is the transformation matrix for element with the size where is the total DOFs of the system.
If and satisfy convergence criterionwhere is the limiting value, then the dynamic responses , , and at the time point are taken as the solution of increment step ; that is, and increment step is started. If displacement increment is used in increment step and if at convergence, then the current point is at the hardening segment of the static equilibrium path, and therefore load increment will be used in increment step .
If the computed and do not satisfy the convergence criterion, then , and repeat steps “ and ”.
3.3. Parameter Setting
The explicit time integration algorithm is conditionally stable. For discrete linear system, the condition for stability is , where is the highest circular eigenfrequency of the system. The computation of is costly [10]. It is pointed out that [23] if then the algorithm is stable, where is nodal mass and is nodal linear displacement stiffness. The nonlinear discrete system can be approximated as linear discrete system within one incremental step. The DRM with kinetic damping solves for the static solution of increment step through virtual vibration process. Before starting increment step , the value of can be chosen so that the nodal mass can be adjusted for all nodes to satisfy the stability condition. We takewhere is a coefficient greater than 1 and is the maximum linear stiffness among the four displacement directions at each node in the directions of the orthogonal coordinate system in the 2dimensional case. This stiffness is estimated by applying infinitesimal displacement node by node to calculate the corresponding nodal load. The term is lower bound value for the nodal mass, whose function is to avoid the formation of local negative stiffness associated with negative mass. Thereby the global diagonal mass matrix can be constructed.
4. Numerical Case Studies
4.1. Model I: Plane Truss
Model I is a plane truss constructed by 3 tension bars (Figure 1). The section area of each bar is 1 mm^{2}. Each tension bar is divided into 4 line elements, and linear polynomial is used as the shape function. By calculating the axial stress of line element, can be computed precisely by using direct forcebalance method. And the axial stress of the element which is used to calculate can be written as , where is the undeformed axial length of the element, and is the deformed axial length.
The artificial onedimensional stressstrain constitutive relation used in Model I is so constructed that it is characteristic of strainsoftening and history dependence. And it is defined as (Figure 1)where and are maximum tension strain and maximum tension stress during the strain history of material, respectively.
In Figure 1 the loaded nodes and loaded degrees of freedom are also shown. Some of the analysis parameters are valued as s, , , and kg, and the values of other analysis parameters are shown in Table 1.

If there is an element whose in each of the three bars, then the computation is finished.
4.2. Model II: A Concrete Column Subjected to Axial Compression
Model II is a concrete column subjected to axial compression (Figure 2), whose thickness is 150 mm. The meshing is illustrated in Figure 2. The meshingstyle in direction is denoted as IIA, IIB, and IIC whose are , , and , respectively, where is the number of elements and is the length of each element in direction. The element of Model II is rectangular type with 4 nodes, and the quadratic polynomial is chosen as the shape function; can be computed approximately by using 4 Gauss integration points.
A scalar damage constitutive model which was put forward by Mazars and PijaudierCabot [31, 32] is used in Model II. In the constitutive model, the material is supposed to behave elastically and to remain isotropic. The constitutive equation can be written asin which is tress tensor, is damage scalar, is initial stiffness tensor, is strain tensor, and the symbol “” indicates the tensorial product contracted on two indices. can be constructed by elasticity modulus and Poisson ratio . The damage scalar can be calculated by the following:in which is the maximum damage value during the strain history of material and the value of the expression can be regarded as realtime damage value. In (14), the weights and are functions of the strain state; and are defined as uniaxial tension and compression damage value [31]. The expressions of and arewhere is the initial damage threshold, , , , and are characteristics of the material, and is called the notion of equivalent strain, whose definition and expression were given in [31]. In Model II, the material parameters are as follows: N/mm^{2}, ; , , = 11,500, , and = 2,000. And by following these material parameters, we got the uniaxial compression strength N/mm^{2} and corresponding strain . The illustration of the uniaxial compressive stressstrain relation is shown in Figure 2.
In Figure 2 the loaded nodes and loaded degrees of freedom are also shown. Some of the analysis parameters are chosen as s, , , and kg, and the values of other analysis parameters are shown in Table 1.
When the expression is satisfied, then the computation is finished.
5. Results and Discussion
For Model I, the loaddeflection curve of node 1 in direction is given in Figure 3, which shows not only the hardeningstages, softeningstages, and load peaks, but also the loadbreaks caused by softening of elements.
The loaddisplacement curves of node 3 in direction of Models IIA, IIB, and IIC are presented in Figure 4. These curves show the hardeningstage and softeningstage of Model II under the axial compression. It is shown in Figure 4 that the different meshes not only present the same load peak value but also give the same equilibrium path before this load peak. The displacement value and load value at peak point of the equilibrium path are close to their estimated values mm, and kN with the assumption that and . After the load peak, different static equilibrium paths are obtained, which is in the accordance with the nonuniqueness of static solution for strainsoftening solid [1].
In Model I, an artificial, onedimensional strainsoftening constitutive model is used such that the static equilibrium path for a simple plane truss is solved. In spite of the simple structure, the solution process incorporates the concept associated with the structure of the universal first law of thermodynamics that is basis of the kinetic damping DRM for the solution of the increment steps. That is, it reflects the law of dissipation of kinetic energy and conversion among the external potential energy, internal potential energy, and kinetic energy. Therefore, the application of the method is in essence not limited by the constitutive relation of material and the type of elements. Thereby, further application of the method can be expected in which constitutive model can be calibrated against experimental data so that it is suitable for use in static analysis of engineering structure.
In Model II, the Mazars elastic damage constitutive model with one scalar damage variable is used. Although the resulting stressstrain relation is simple and smooth enough, it is relatively poor in reflecting the constitutive relation under multiaxial state of stress [33]. In future studies, the elastic and/or plastic damage constitutive model with 2 scalar damage variables should be used in the static analysis of different concrete members in order to truly reflect the tensioncompression difference of concrete material under complicated condition [34].
6. Conclusions
The process for the solution of static problems based on nonlinear finite element DRM method is clarified with three depths of global, increment step, and the motion process with application of kinetic damping. The method of DRM makes the original illposed static problem become wellposed and solvable by introducing the time domain for the virtual dynamic process such that the static solution is obtained by way of dissipation of mechanical energy.
The global analysis including the ascending branch, the peak point, and the descending branch of the equilibrium path of the strainsoftening structure is realized by application of kinetic damping to the motion process, by introducing the rule for conversion between load increment and displacement increment and by introducing the divergence criterion.
Two instances of structures with strainsoftening and the associated illposedness, one model with geometrical nonlinearity and the other model with mesh size sensitivity, are successfully solved using DRM.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was financially supported by the National Natural Science Foundation of China (51178328).
References
 O. C. Zienkiewicz, R. L. Taylor, and D. D. Fox, The finite element method for solid and structural mechanics, Elsevier/Butterworth Heinemann, Amsterdam, Seventh edition, 2014. View at: Publisher Site  MathSciNet
 ABAQUS 6.14 analysis users guide (volume II): Damping in dynamic analysis, Dassault systemes simulia corp,.
 J. Garcia, G. RIO, and JM. Cadou, Numerical study of dynamic relaxation methods: Contribution to the modeling of inflatable lifejackets, LAP LAMBERT Academic Publishing Gmbh Co, Saarbrücken, Germany, 2012.
 H. Kardestuncer and DH. Norrie, Finite element handbook, Science Press, Beijing, China, 1996.
 R. D. Borst, M. A. Crisfield, and J. J. C. Remmers, Nonlinear Finite Element Analysis of Solids and Structures, JOHN WILEY SONS, Chester, UK, 2012.
 J. J. Jiang, X. Z. Lu, and L. P. Ye, Finite Element Analysis of Concrete Structures, Tsinghua University Press, Beijing, China, 2005.
 Z. He and J. P. Ou, Nolinear Finite Element Analysis of Reinforced Concrete Structures, Harbin Institute of Technology Press,Harbin, China, 2007.
 M. A. Crisfield, NonLinear Finite Element Analysis of Solids and Structures, Vol 1: Essentials, JOHN WILEY SONS, Chichester, UK, 1991.
 G. A. Wempner, “Discrete approximations related to nonlinear theories of solids,” International Journal of Solids and Structures, vol. 7, no. 11, pp. 1581–1599, 1971. View at: Publisher Site  Google Scholar
 Abaqus 6.14 theory guide, Dassault systemes simulia corp,.
 Y. L. Duan, “Loacal arclength methoda solution procedure for nonlinear finite element method,” ActaMechanicaSinica, vol. 29, no. 1, pp. 116–122, 1997. View at: Google Scholar
 Z. Yang and J. Chen, “Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arclength methods,” International Journal of Solids and Structures, vol. 41, no. 34, pp. 801–826, 2004. View at: Publisher Site  Google Scholar
 S. B. Cai, P. S. Shen, and Z. A. Hou, “Pseudotransient Finite Element Method for Geometric and Material Nonlinear Static Analysis of RC Shear Wall,” Journal Guizhou University(Natural Science, vol. no.3, pp. 144–152, 1997. View at: Google Scholar
 L. Z. Liu, M. Y. Zhang, and X. H. Liu, “Finite Element Simulation of Rolling Process by Static Implicit Method and Dynamic Explicit Method,” Journal of Plasticity Engineering, vol. 8, no. 4, pp. 81–83, 2001. View at: Google Scholar
 P. Natário, N. Silvestre, and D. Camotim, “Web crippling failure using quasistatic FE models,” ThinWalled Structures, vol. 84, pp. 34–49, 2014. View at: Publisher Site  Google Scholar
 L. A. Zhang, T. X. Yu, and R. Wang, “maDR Method and its applications , Computational structural mechanics and applications,” maDR Method and its applications , Computational structural mechanics and applications, vol. 7, no. 3, pp. 25–33, 1990. View at: Google Scholar
 N. A. Newmark, “A Method of Computation for Structural Dynamics,” Journal of the Engineering Mechanics Division, vol. 85, no. 3, pp. 67–94, 1959. View at: Google Scholar
 J. R. H. Otter, “Computations for prestressed concrete reactor pressure vessels using dynamic relaxation,” Nuclear Engineering and Design, vol. 1, no. 1, pp. 61–75, 1965. View at: Publisher Site  Google Scholar
 J. R. Otter, A. C. Cassell, and R. E. Hobbs, “Dynamic relaxation,” Proceedings of the Institution of Civil Engineers, vol. 35, no. 4, pp. 633–656, 1966. View at: Publisher Site  Google Scholar
 J. Harding, “The elastoplastic analysis of imperfect cylinders,” Proceedings of the Institution of Civil Engineers, vol. 65, no. 4, pp. 875–892, 1978. View at: Publisher Site  Google Scholar
 S. Webb and P. Dowling, “Largedeflexion elastoplastic behaviour of discretely stiffened plates,” Proceedings of the Institution of Civil Engineers, vol. 69, no. 2, pp. 375–401, 1980. View at: Publisher Site  Google Scholar
 G. Ramesh and C. S. Krishnamoorthy, “Geometrically nonlinear analysis of plates and shallow shells by dynamic relaxation,” Computer Methods in Applied Mechanics and Engineering, vol. 123, no. 14, pp. 15–32, 1995. View at: Publisher Site  Google Scholar
 S.E. Han and K.S. Lee, “A study of the stabilizing process of unstable structures by dynamic relaxation method,” Computers and Structures, vol. 81, no. 17, pp. 1677–1688, 2003. View at: Publisher Site  Google Scholar
 J. Rodriguez, G. Rio, J. M. Cadou, and J. Troufflard, “Numerical study of dynamic relaxation with kinetic damping applied to inflatable fabric structures with extensions for 3D solid element and nonlinear behavior,” ThinWalled Structures, vol. 49, no. 11, pp. 1468–1474, 2011. View at: Publisher Site  Google Scholar
 R. Hill, “A general theory of uniqueness and stability in elasticplastic solids,” Journal of the Mechanics and Physics of Solids, vol. 6, no. 3, pp. 236–249, 1958. View at: Google Scholar
 D. C. Drucker, “A definition of stable inelastic material,” vol. 26, pp. 101–106, 1959. View at: Google Scholar  MathSciNet
 V. I. Nechaev, “Ivan Matveevich Vinogradov,” Ministerstvo Obrazovaniya Rossi\u\i sko\u\i \ Federatsii. Matematika v Shkole, no. 3, pp. 7980, 1983. View at: Google Scholar  MathSciNet
 R. de Borst, L. J. Sluys, H.B. Muhlhaus, and J. Pamin, “Fundamental issues in finite element analyses of localization of deformation,” Engineering Computations, vol. 10, no. 2, pp. 99–121, 1993. View at: Publisher Site  Google Scholar
 Vinogradov, I. Matveevich, S. I. Adian et al., Математическая энциклопедия, vol. 3, Science Press, Beijing, China, 1997.
 O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The finite element method: its basis and fundamentals, Elsevier/Butterworth Heinemann, Amsterdam, Seventh edition, 2013. View at: Publisher Site  MathSciNet
 J. Mazars, “A description of micro and macroscale damage of concrete structures,” Engineering Fracture Mechanics, vol. 25, no. 56, pp. 729–737, 1986. View at: Publisher Site  Google Scholar
 J. Mazars and G. PijaudierCabot, “Continuum damage theory: application to concrete,” Journal of Engineering Mechanics, vol. 115, no. 2, pp. 345–365, 1989. View at: Publisher Site  Google Scholar
 J. Li and X. D. Ren, “A review on the constitutive model for static and dynamic damage of concrete,” Advances in mechanics, vol. 40, no. 3, pp. 284–297, 2010. View at: Google Scholar
 C. Comi and U. Perego, “Fracture energy based bidissipative damage model for concrete,” International Journal of Solids and Structures, vol. 38, no. 3637, pp. 6427–6454, 2001. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2017 Wei Wang and Xiaozu Su. 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.