Abstract
The contact is a typical non-linear problem that exists in various projects. For traditional three-node triangular mesh and four-node quadrilateral mesh, the accuracy and convergence of the calculation results are affected by the quality of the mesh. The test space and trial space in the virtual element method (VEM) do not need to be accurately calculated, avoiding mesh dependence. In this paper, the formulation of linear elasticity and the formulation of the frictionless node-to-segment (NTS) contact model via VEM are shown. There are four numerical simulations. The sensitivity of the virtual element method to mesh distortion is studied in the first numerical simulation. The exactness and convergence of the algorithm are investigated by the second numerical example. The second numerical example simultaneously explores the penalty factor’s effect on the results. The third example investigated the impact of mesh shape and number of Voronoi mesh elements on the results by comparing normal contact stresses. The fourth numerical example studies the application of the method to engineering. Those numerical examples show that the virtual element method is insensitive to mesh distortion and could solve the joint contact in engineering.
1 Introduction
The contact problem is a typical nonlinear problem, widespread in actual engineering such as geotechnical engineering, building structure, water project, and machinery engineering. In the past 2 decades, with the development of electronic computers and the rise and development of various numerical methods, there are powerful means to handle the contact problems such as finite element method (FEM), numerical manifold method (NMM) (Yang et al., 2019; Zheng et al., 2019; Yang et al., 2020a; Yang et al., 2020b; Yang et al., 2021a; Yang et al., 2021b), boundary element method (BEM).
For those numerical methods, FEM has been the greatest broadly utilized. Hughes (Hughes et al., 1976) and Francavilla (Padmanabhan and Laursen, 2001) are considered pioneers who solve the contact problem by using FEM. In order to improve the classical contact discretization, there are some methods that have been studied, which are based on constraint element node enforcement. (Hautefeuille et al., 2012; Khoei et al., 2006; Wriggers et al., 2001). More recently, some researchers (Wriggers et al., 2001; da Veiga et al., 2014; Sheng and Yuan, 2012; Krstulovic-Opara et al., 2002; Liu et al., 2007; Flemisch et al., 2005) prefer the mortar method for the discretization of contact constraints. The stable interpolation condition for contact constraints is provided by those methods because of the weak form based on the mortar method. What is known to us is that the gradient of the FEM with a standard degree of freedom is not continuous on internal element edges. When FEM is employed, it has been proved that results are highly sensitive to mesh quality (Lee and Bathe, 1993; Liu et al., 2007; Yang et al., 2014). For contact problems, the results are influenced by the quality of the meshes where they are in the possible contact area. When the contact boundary between two contact bodies is irregular, the calculation results using FEM are greatly affected.
Because of the insufficient FEM, the VEM is proposed by Brezzi (Beirão Da Veiga et al., 2013). Since the birth of the VEM, it has been applied in many aspects by many scholars. The two-dimensional Poisson problem that is discretized by polygonal discretization is solved by VEM (Sutton, 2017). In the VEM framework, the maximum entropy basis function is employed to settle the Poisson problem and linear elastic problem by Ortiz-Bernardin (Ortiz-Bernardin et al., 2017). Sun and Lin (Sun et al., 2020) studied the stability of stony soil slope under excavation using VEM. However, there are few papers that solve the contact problem by VEM.
The contact constraint is commonly handled by the Lagrange multiplier method (LMM) (Béchet et al., 2010; Hautefeuille et al., 2012), the penalty method (PM) (Liu and Borja, 2008; Liu and Borja, 2010a; Liu and Borja, 2010b), and the augmented Lagrange method (ALM). The PM can convert the non-linear contact problem into material nonlinearity. The advantage of the PM is that the global system is not extended when the contact conditions are introduced. The disadvantage of PM is that the contact constraints can be satisfied approximately. The contact constraint can be satisfied accurately. However, the global system needs to be extended by an additional variable for the Lagrange multiplier. To evade the disadvantages of PM and LMM, the ALM was proposed. However, there are sub-iterations in each calculation step in the ALM, which is its primary deficiency. This paper aims to find the solution to the contact problem in engineering, and the result of PM can fully meet the needs of engineering. In summary, the contact constraint is handled by PM in this paper.
The rest of this article is composed of the following parts. The application of the lowest order VEM for the linear elasticity problem is presented in Section 2 and Section 3 shows the NTS contact model. In Section 4, the numerical examples are given for performances of the VEM in different situations, including the response to mesh quality, the performance of normal contact pressure in Hertzian contact, the influence of varying mesh shapes for the normal contact pressure in the horizontal contact interface and the application of algorithms in engineering. The discussion and concluding remarks are presented in Section 5 and Section 6, respectively.
2 Linear Elasticity
2.1 The Model
The elastic body is composed of an open domain . The Dirichlet boundary and Neumann boundary are represented by , . The displacement field of the elastic body can be represented by . The Dirichlet boundary conditions are . The linear elastic boundary-value problem can be expressed that discovering satisfies the following conditions (Ortiz-Bernardin et al., 2017):Where represents the stress tensor, is the body force, is the unit normal of boundary, the is the external traction. The equivalent Calerkin variational formulation can be expressed that finding satisfies the equationWhere and represent the displacement test and trial space.
2.2 Discrete Bilinear Form
The domain is divided into non-overlapping polygonal elements that make up the area .
The virtual function space is defined to beWhere is the local space on the element sustains some properties (Sutton, 2017).
The basis function of space can be expressed as .
For any .
Therefore
The basis function of polynomial space can be expressed as (Nguyen-Thanh et al., 2018). is a subspace of .
We can define a projection from the virtual function space to the polynomial space . The defined projection needs to satisfy the following equation:
The defined vectors form isWhere is the constant term projection operator (Sutton, 2017).
2.2 Element Stiffness Matrix
Due to , The equation can be obtainedWhere represents the expression of the matrix under the basis function .
Additionally, on account of , The equation can be obtained
From the Eq. 2.9, the equation is obtained
The equations of Eq. 2.11 and Eq. 2.10 are brought into Eq. 2.2, and the following equation is obtained
The matrixes of and are computed as follows
Eq. 2.10 can be written as a matrix expression as
The element stiffness matrix is (Beirão da Veiga et al., 2014)
The meaning of the Eq. 2.15 can refer to the articles (Chen, 2015; Ortiz-Bernardin et al., 2017).
3 The Contact Problem for VEM
3.1 Description of Different Contact Condition
Figure 1 presents a two-dimensional frictionless contact model. For domain and , the possible contact boundaries are and , g represents the contact gap. let and are the contact force vectors in the contact region.
FIGURE 1
Contact problem description.
Consider and are the infinitesimal region where is coming into contact. The virtual work done by the contact traction is
In the contact area, every point should satisfy the equilibrium equation
Thus, we could consider the integral in Eq. 3.1 along the contact line or .
3.2 NST Contact Model
In this paper, the NST contact model is employed. A mapping is defined in Figure 2.
FIGURE 2
NTS contact model.
In the Figure 2, the tangent vector and normal vector of the master element establish a local coordinate system . The can be written in matrix form as
The isWhere stands for the length of node to and is coordinates of node .
The projection node coordinates of node is on the surface composed of nodes and is
The relative displacement from node to the corresponding segment could be calculated as
The indicates the normal gap, and it can be computed by the following equation
Bring the Eq. 3.5 into Eq. 3.3, The contact integral can be calculated
The is replaced by and is replaced by . The contact force can be split into
In the frictionless contact problem, the tangential part disappears . Thus, the integral of Eq. 3.7 can be written as
3.3 Penalty Method
When the PM is employed for the frictionless contact problem, contact traction can be written as
Thus, the Eq. 3.9 can be written as
Using numerical methods to catch the ball nonlinear problems, like contact problems, solved by iterative methods. The method requires the derivatives of the weak form . This method is called linearization.
Differentiating Eq. 4–13 with respect to time , we getWithout considering the rotation of the contact body, the contact matrix is
4 Numerical Examples
There are three numerical examples to be carried out in this section. In example 1, the cantilever beam with free end applied to bending moment can be employed to study the VEM and FEM mesh distortion tolerance, respectively. The accuracy and convergence of the algorithm are explained by Hertz contact in example 2. The influence of the mesh shape and amount of element is presented in example 2. The application of algorithmic reengineering is shown in example 3. The results of the FEM are obtained by Abaqus.
4.1 Example 1: Cantilever Beam With Free End Imposed to Bending Moment
In this part, the geometry of the cantilever beam with free end imposed to bending moment is presented in Figure 3A. The geometric parameters (Yang et al., 2014), which are , are applied to computation. In the Figure 3B, two plane strain quadrilateral elements are used to discretize the model. The distortion of the element is expressed by the distortion parameter (Zhang and Rajendran, 2008; Remacle et al., 2012; Stavroulakis, 2013), which is always considered in articles that explore the response of calculation method to mesh distortion.
FIGURE 3
The cantilever beam with imposed to bending moment. (A) Discretion for cantilever beam. (B) Geometric parameters of cantilever beam (Yang et al., 2014).
In this case, the value which is calculated point is compared with the exact value which is to illustrate the sensitivity of different numerical methods to the quality of the mesh. In Figure 4, the following conclusions can be obtained:
FIGURE 4
Curve of the ratio of calculated value to real value with twist parameter.
For the twist factor is zero, the performance of FEM is better than VEM. With distortion parameters increasing, the accuracy of the FEM decreases faster than the VEM. The result of FEM is equal to VEM when the distortion parameter is near 1. When the distortion parameter exceeds 1, the accuracy of VEM is superior to that of FEM.
4.2 Example 2: Hertzian Problem
The second simulation is to testify the convergence rate and exactness of the algorithm. The reason for choosing the Hertzian contact as the second numerical example is an analytical solution to the Hertzian contact. The accuracy of the method is illustrated by comparing the computational and analytical solutions.
The model geometry is shown in Figure 5. The disc of the radius is loaded by a pressure . The and is selected as material parameters. Rectangular at the bottom of the model is taken as . The quadrilateral meshes are used. Because of the model’s symmetry, half of the models are selected for research.
FIGURE 5
Hertz contact model.
It is known that the penalty factor has an impact on the result. When the number of mesh is 600, the number of iterations and the maximum contact force under different penalty functions are listed in the following Table1. From Table 1, the penalty factor is finally selected as . The model is discretized by 1,155 quadrilateral elements, as exhibited in Figure 6A. The contour of normal pressure is exhibited in Figure 6B.
TABLE 1
| Penalty iterations | |||
|---|---|---|---|
| 1 | −227.11 | −227.10 | −227.11 |
| 2 | −98.20 | −98.01 | −102.62 |
| 3 | −84.71 | −84.57 | −87.65 |
| 4 | −80.59 | −80.59 | −82.29 |
| 5 | −79.86 |
The convergence process of penalty function for Hertzian contact.
FIGURE 6
The finite mesh and contour of normal stress is the normal stress in y-direction.
In Figure 7, the analytical solution, numerical solutions obtained by the VEM and the FEM for contact force are shown, respectively. Some conclusions can be drawn:
FIGURE 7
The normal contact stress distribution along the contact zone.
In the contact region where is less than 1.4, the difference between the normal stress obtained by the VEM and the analytical solution is smaller than that obtained by the FEM. The maximum stresses of analytical solution, virtual element method, and FEM are 83.60, 82.21, and 80.57. When is more significant than 1.4, The curve of the VEM coincides with the FEM.
In practical problems, we are more concerned about the maximum normal contact traction. The contact stress with 389, 600, and 1,155 meshes are shown in
Figure 8. There are some conclusions reached from
Figure 9.
1) When the same mesh discretizes the structure, the maximum stress from the VEM is closer to the analytical solution.
2) The convergence rate of the VEM is higher than FEM for the number of the element from 389 to 600.
3) As the number of mesh elements increased from 600 to 1,155, the convergence rate of the VEM was the same as the FEM.
FIGURE 8
Maximum normal contact stress under different numbers of elements.
FIGURE 9
Model geometry with straight lines in the contact interface.
4.3 Example 3: A Horizontal Interface Under Uniform Compression
The third numerical model has been researched by Hirmand (Hirmand et al., 2015). This example is to compare the influence of different mesh shapes and the number of elements for contact stress. The geometry is shown in Figure 9. A displacement loads the upper rectangle in the y-direction . The displacements of the bottom of the model are both fixed. The Young’s modulus is . The Poisson’s ratio is .
The advantage of VEM is to it calculates arbitrary mesh shapes. The Voronoi mesh (Talischi et al., 2012) is used to discrete the model. The normal stress along the contact interface of the different shape mesh with Hirmand is shown in Figure 10A. From Figure 10A, it can be concluded that the maximum normal contact stress is basically the same, with slight differences at both ends. Therefore, it can be noticed that the normal contact stress is slightly affected by the mesh shape. In this example, the influence of the different amounts of elements for the normal contact stress is studied. The 50, 100, 200, 300, and 500 Voronoi elements are employed to discrete the model. In Figure 10B, mesh 1, mesh 2, mesh 3, mesh 4, and mesh 5 correspond to 50, 100, 200, 300, and 500 Voronoi elements. Figure 10B presents the normal contact stress for different number elements. The following conclusions are obtained from Figure 10B: When the number of elements is 200, 300, and 500, the normal contact stress curves remain coincident.
FIGURE 10
The distribution of normal contact traction in horizontal contact interface. (A) The distribution of normal contact traction for different mesh shapes. (B) The distribution of normal contact traction for different number element.
Figure 11A shows the contour of vertical displacement obtained by VEM under the Voronoi mesh. The simulation of Hirmand under quadrilateral mesh is presented in Figure 11B. It is noted that the curve for VEM is in line with the results shown by Hirmand. The contours of the normal stress is shown in Figure 11C.
FIGURE 11
Contours of displacement and pressure for the direction (Uy is the displacement in y-direction and σyy is the stress in y-direction).
4.4 Example 4: Dam With Joint
This numerical example simulates a dam problem with a cracked foundation. This example is shown in Zheng (Zheng et al., 2002) in the 2005 year. The geometric model is exhibited in Figure 12A. The model size parameters are , , , , . The Young’s modulus and Poisson’s ratio are the material parameters for this model. The is taken as volumetric weight. The displacements of the bottom left and right of the foundation are fixed. The coordinates of the joint tip are from (4,10) to (10,4). The joint end is fixed and will not propagate, and there is no friction at the crack interface.
FIGURE 12
The geometric parameters and Voronoi mesh for dam with joint. (A) Dam geometry with jointed foundation. (B) Voronoi mesh for dam with jointed foundation.
The stress situation is analyzed using two load steps. The first load step only considers the self-weight of the dam body and foundation; the second load step applies a triangularly distributed normal water pressure to the surface of the dam body to simulate the condition of the reservoir after it is full of water. The Voronoi mesh was used to discretize the model. The Voronoi mesh is presented in Figure 12B. The displacement contour along the -direction is demonstrated in Figure 13A, and the displacement contour in -direction is shown in Figure 13B. From Figures 13A,B, it is noted that the displacement contours are discontinuous at the joint. The maximum and minimum principal stress contours are presented in Figure 14. As expected, the maximum stress occurs at the joint tip. The phenomenon is consistent with Li’s research (Li et al., 2022). In their studies, the strategy derived from the meshless numerical manifold method (MNMM) is employed by Li to solve linear elastic fractures.
FIGURE 13
The contours of displacement in and directions (Ux and Uy are the displacement in displacement in x and y direction, respectively).
FIGURE 14
The contours of maximum principal stress is the maximum principal stress and minimum principal stress) is the minimum principal stress).
5 Discussion
In Example 1, when the distortion parameter is 0, the result of the FEM is more accurate than the VEM. The reason is that the function in the virtual element space satisfies the globally continuous on the element boundary. When the distortion parameter gradually increases, the downward trend of the resulting curve of the VEM is slower than that of the FEM. When the element distortion parameter exceeds a certain value, the result of the VEM is more accurate than the FEM. The result of VEM calculation is less affected by mesh quality.
In Figure 7 of example 2, the result of the VEM is better than the FEM in the contact region where is less than 1.4 because the discrete element size is small, which makes the element more twisted. And in the contact region where is more significant than 1.4, the main reason is that discrete element size is relatively large, so the element distortion is minor.
Under the 389 elements in Figure 8 of example 2, the main reason for the similar maximum normal stress of the VEM and the FEM is that the element distortion parameter is small. When the number of elements increases, the corresponding element size becomes smaller, the distortion increases, and the advantages of the VEM become more obvious. When the number of elements is 600, the difference between the VEM and the FEM results is greater than the difference between the VEM and the FEM when the amount of elements is 389. When the amount of elements is 1,155, the difference between the outcomes of the VEM and FEM is close to the difference between the outcomes of the amount of elements 600. On the whole, the better convergence and accuracy of the VEM in Hertz contact lies in the VEM is suitable for general polygons or polyhedrons, which is used flexibly for discrete complex contact surfaces (Beirão Da Veiga et al., 2013; Chen, 2015; Benedetto et al., 2016).
From example 1 and example 2, it can be known that the distortion of the element greatly influences the results. In example 3, normal contact stress in the contact interface is the same under different mesh shapes. The main reason is that the VEM test and trial space do not need to be accurately calculated, avoiding mesh dependence, and the contact interface is straight.
In example 4, a dam with cracks under Voronoi mesh was modeled. As expected, the maximum stress occurs at the joint tip in example 4, so the joint is an important cause affecting safety in engineering. So, the joint the focus of the study. Through example 4, it can be obtained that the NST contact model based on the virtual element method can solve engineering problems well under the Voronoi mesh.
6 Conclusion
A strategy to handle the contact problem is proposed in this paper, which is stemmed from the NTS model and the VEM. The effect of mesh distortion for results, the accuracy and convergence rate for Hertz contact, the impact of different mesh shapes and different elements numbers for results and the application of algorithms in engineering are implemented by several numerical examples. The results show that:
1) The VEM is insensitive to the mesh quality.
2) When the mesh on the contact interface is distorted, The VEM has high convergence and accuracy.
3) When contact problem is handled by VEM, the normal contact stress on the contact surface is slightly affected by the mesh shape.
4) The results in the fourth example show that the VEM can solve the contact problem in engineering under Voronoi mesh
Statements
Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.
Author contributions
YS performed the data analyses and wrote the manuscript. GS contributed significantly to analysis and manuscript preparation. QY and JW helped perform the analysis with constructive discussions.
Funding
Supported by the national Natural Science Foundation of China (Grant No. 11972043) and National Key R&D Program of China (2018YFE0100100).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
BéchetÉ.MoësN.WohlmuthB. (2010). A Stable Lagrange Multiplier Space for Stiff Interface Conditions within the Extended Finite Element Method. Int. J. Numer. Methods Eng.78, 931–954. 10.1002/nme.2515
2
Beirão Da VeigaL.BrezziF.CangianiA.ManziniG.MariniL. D.RussoA. (2013). Basic Principles of Virtual Element Methods. Math. Models Methods Appl. Sci.23, 199–214. 10.1142/s0218202512500492
3
Beirão da VeigaL.BrezziF.MariniL. D.RussoA.BrezziF.ManziniG. (2014). The Hitchhiker's Guide to the Virtual Element Method. Math. Models Methods Appl. Sci.24, 1541–1573. 10.1142/s021820251440003x
4
BenedettoM. F.BerroneS.ScialòS. (2016). A Globally Conforming Method for Solving Flow in Discrete Fracture Networks Using the Virtual Element Method. Finite Elem. Anal. Des.109, 23–36. 10.1016/j.finel.2015.10.003
5
ChenL. (2015). Equivalence of Weak Galerkin Methods and Virtual Element Methods for Elliptic Equations, 1–11.
6
da VeigaL. B.LipnikovK.ManziniG. (2014). The Mimetic Finite Difference Method for Elliptic Problems. Model. Simulation Appl.11, 1–389. 10.1007/978-3-319-02663-3
7
FlemischB.PusoM. A.WohlmuthB. I. (2005). A New Dual Mortar Method for Curved Interfaces: 2D Elasticity. Int. J. Numer. Meth. Engng63, 813–832. 10.1002/nme.1300
8
HautefeuilleM.AnnavarapuC.DolbowJ. E. (2012). Robust Imposition of Dirichlet Boundary Conditions on Embedded Surfaces. Int. J. Numer. Methods Eng.90, 1102–1119. 10.1002/nme.3306
9
HirmandM.VahabM.KhoeiA. R. (2015). An Augmented Lagrangian Contact Formulation for Frictional Discontinuities with the Extended Finite Element Method. Finite Elem. Anal. Des.107, 28–43. 10.1016/j.finel.2015.08.003
10
HughesT. J. R.TaylorR. L.SackmanJ. L.CurnierA.KanoknukulchaiW. (1976). A Finite Element Method for a Class of Contact-Impact Problems. Comp. Methods Appl. Mech. Eng.8, 249–276. 10.1016/0045-7825(76)90018-9
11
KhoeiA. R.ShamlooA.AzamiA. R. (2006). Extended Finite Element Method in Plasticity Forming of Powder Compaction with Contact Friction. Int. J. Sol. Structures43, 5421–5448. 10.1016/j.ijsolstr.2005.11.008
12
Krstulovic-OparaL.WriggersP.KorelcJ. (2002). A C 1 -continuous Formulation for 3D Finite Deformation Frictional Contact. Comput. Mech.29, 27–42. 10.1007/s00466-002-0317-z
13
LeeN.-S.BatheK.-J. (1993). Effects of Element Distortions on the Performance of Isoparametric Elements. Int. J. Numer. Meth. Engng.36, 3553–3576. 10.1002/nme.1620362009
14
LiW.YuX.LinS.QuX.SunX. (2022). A Numerical Integration Strategy of Meshless Numerical Manifold Method Based on Physical Cover and Applications to Linear Elastic Fractures. Eng. Anal. Boundary Elem.134, 79–95. 10.1016/j.enganabound.2021.09.028
15
LiuF.BorjaR. I. (2008). A Contact Algorithm for Frictional Crack Propagation with the Extended Finite Element Method. Int. J. Numer. Methods Eng.76, 1489–1512. 10.1002/nme.2376
16
LiuF.BorjaR. I. (2010). Finite Deformation Formulation for Embedded Frictional Crack with the Extended Finite Element Metho. Int. J. Numer. Methods Eng.82, 773–804. 10.1002/nme.2782
17
LiuF.BorjaR. I. (2010). Stabilized Low-Order Finite Elements for Frictional Contact with the Extended Finite Element Method. Comp. Methods Appl. Mech. Eng.199, 2456–2471. 10.1016/j.cma.2010.03.030
18
LiuG. R.DaiK. Y.NguyenT. T. (2007). A Smoothed Finite Element Method for Mechanics Problems. Comput. Mech.39, 859–877. 10.1007/s00466-006-0075-4
19
Nguyen-ThanhV. M.ZhuangX.Nguyen-XuanH.RabczukT.WriggersP. (2018). A Virtual Element Method for 2D Linear Elastic Fracture Analysis. Comp. Methods Appl. Mech. Eng.340, 366–395. 10.1016/j.cma.2018.05.021
20
Ortiz-BernardinA.RussoA.SukumarN. (2017). Consistent and Stable Meshfree Galerkin Methods Using the Virtual Element Decomposition. Int. J. Numer. Methods Eng.112, 655–684. 10.1002/nme.5519
21
PadmanabhanV.LaursenT. A. (2001). A Framework for Development of Surface Smoothing Procedures in Large Deformation Frictional Contact Analysis. Finite Elem. Anal. Des.37, 173–198. 10.1016/s0168-874x(00)00029-9
22
RemacleJ.-F.LambrechtsJ.SenyB.MarchandiseE.JohnenA.GeuzainetC. (2012). Blossom-Quad: A Non-uniform Quadrilateral Mesh Generator Using a Minimum-Cost Perfect-Matching Algorithm. Int. J. Numer. Meth. Engng89, 1102–1119. 10.1002/nme.3279
23
ShengZ.YuanG. (2012). An Improved Monotone Finite Volume Scheme for Diffusion Equation on Polygonal Meshes. J. Comput. Phys.231, 3739–3754. 10.1016/j.jcp.2012.01.015
24
StavroulakisG. E. (2013). Lecture Notes in Applied and Computational Mechanics: Introduction.
25
SunG.LinS.ZhengH.TanY.SuiT. (2020). The Virtual Element Method Strength Reduction Technique for the Stability Analysis of Stony Soil Slopes. Comput. Geotechnics119, 103349. 10.1016/j.compgeo.2019.103349
26
SuttonO. J. (2017). The Virtual Element Method in 50 Lines of MATLAB. Numer. Algor75, 1141–1159. 10.1007/s11075-016-0235-3
27
TalischiC.PaulinoG. H.PereiraA.MenezesI. F. M. (2012). PolyMesher: A General-Purpose Mesh Generator for Polygonal Elements Written in Matlab. Struct. Multidisc Optim45, 309–328. 10.1007/s00158-011-0706-z
28
WriggersP.Krstulovic-OparaL.KorelcJ. (2001). Smooth C1-Interpolations for Two-Dimensional Frictional Contact Problems. Int. J. Numer. Meth. Engng.51, 1469–1495. 10.1002/nme.227
29
YangY.TangX.ZhengH. (2014). A Three-Node Triangular Element with Continuous Nodal Stress. Comput. Structures141, 46–58. 10.1016/j.compstruc.2014.05.001
30
YangY.SunG.ZhengH.QiY. (2019). Investigation of the Sequential Excavation of a Soil-Rock-Mixture Slope Using the Numerical Manifold Method. Eng. Geology.256, 93–109. 10.1016/j.enggeo.2019.05.005
31
YangY.XuD.LiuF.ZhengH. (2020). Modeling the Entire Progressive Failure Process of Rock Slopes Using a Strength-Based Criterion. Comput. Geotechnics126, 103726. 10.1016/j.compgeo.2020.103726
32
YangY.SunG.ZhengH.YanC. (2020). An Improved Numerical Manifold Method with Multiple Layers of Mathematical Cover Systems for the Stability Analysis of Soil-Rock-Mixture Slopes. Eng. Geology.264, 105373. 10.1016/j.enggeo.2019.105373
33
YangY.XiaY.ZhengH.LiuZ. (2021). Investigation of Rock Slope Stability Using a 3D Nonlinear Strength-Reduction Numerical Manifold Method. Eng. Geology.292, 106285. 10.1016/j.enggeo.2021.106285
34
YangY.WuW.ZhengH. (2021). Stability Analysis of Slopes Using the Vector Sum Numerical Manifold Method. Bull. Eng. Geol. Environ.80, 345–352. 10.1007/s10064-020-01903-x
35
ZhangB. R.RajendranS. (2008). 'FE-Meshfree' QUAD4 Element for Free-Vibration Analysis. Comp. Methods Appl. Mech. Eng.197, 3595–3604. 10.1016/j.cma.2008.02.012
36
ZhengH.LiZ.GeX.YueZ. (2002). Mixed Finite Element Method for Interface Problems. Chin. J. Rock Mech. Eng.20, 1–8. 10.3321/j.issn:1000-6915.2002.01.001
37
ZhengH.YangY.ShiG. (2019). Reformulation of Dynamic Crack Propagation Using the Numerical Manifold Method. Eng. Anal. Boundary Elem.105, 279–295. 10.1016/j.enganabound.2019.04.023
Summary
Keywords
virtual element method, sensitivity of mesh, frictionless, node-to-segment contact model, voronoi mesh
Citation
Sun Y, Sun G, Yi Q and Wang J (2022) The Virtual Element Method for the Dam Foundation With Joint. Front. Earth Sci. 10:875561. doi: 10.3389/feart.2022.875561
Received
14 February 2022
Accepted
22 February 2022
Published
24 March 2022
Volume
10 - 2022
Edited by
Yongtao Yang, Institute of Rock and Soil Mechanics (CAS), China
Reviewed by
Wei Li, Linyi University, China
Shan Lin, Beijing University of Technology, China
Updates
Copyright
© 2022 Sun, Sun, Yi and Wang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Guanhua Sun, ghsun@whrsm.ac.cn
This article was submitted to Geohazards and Georisks, a section of the journal Frontiers in Earth Science
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.