The elastic body is composed of an open domain Ω ⊆ R 2 . The Dirichlet boundary and Neumann boundary are represented by Γ g , Γ t . The displacement field of the elastic body can be represented by u ( x ) . The Dirichlet boundary conditions are g ( x ) . The linear elastic boundary-value problem can be expressed that discovering u ( x ) satisfies the following conditions (Ortiz-Bernardin et al., 2017): { ∇ σ + b = 0 u = g g ∈ Γ g σ ⋅ n = t b ∈ Γ t (2.1) Where σ represents the stress tensor, b is the body force, n is the unit normal of boundary, the t is the external traction. The equivalent Calerkin variational formulation can be expressed that finding u ⊂ U satisfies the equation a ( u , v ) = l ( u , v ) u ⊂ U , v ⊂ V a ( u , v ) = ∫ Ω σ ( u ) : ε ( v ) d Ω l ( v ) = ∫ Ω b ⋅ v d Ω + ∫ Γ t t ⋅ v d Γ (2.2) Where V and U represent the displacement test and trial space.

The domain Ω is divided into non-overlapping polygonal elements that make up the area Γ h .

The virtual function space is defined to be V h : = { v ∈ H 1 ( E ) : v | E ∈ V h E   f o r   a l l   E ∈ Γ h } (2.3) Where V h E is the local space on the element E . V h E sustains some properties (Sutton, 2017).

The basis function of space V h E = ( V h E ) 2 can be expressed as ϕ 1 ¯ , ⋯ , ϕ N ¯ , ⋯ , ϕ 1 ¯ ⋯ ϕ N ¯ . ϕ i ¯ = [ ϕ i 0 ] , ϕ i ¯ = [ 0 ϕ i ] , i = 1 , ⋯ N

For any u = ( u 1 , u 2 ) = : ( u ¯ , u ¯ ) ∈ V h E . u = ∑ i = 1 N u i ¯ ϕ i ¯ + ∑ i = 1 N u i ¯ ϕ i ¯ , u i ¯ = χ i ( u ¯ ) , u i ¯ = χ i ( u ¯ )

Therefore d o f i ( u ) = χ i ( u 1 ) , d o f i + N ( u ) = χ i ( u 2 )

The basis function of polynomial space P E = ( P E ) 2 can be expressed as (Nguyen-Thanh et al., 2018). P E is a subspace of V h E . [ m 1 ¯ , m 2 ¯ , m 3 ¯ , m 1 ¯ , m 2 ¯ , m 3 ¯ ] (2.4) m α ¯ = [ m α 0 ] , m α ¯ = [ 0 m α ] , α = 1,2,3

We can define a projection from the virtual function space to the polynomial space Π E : V h E → P E . The defined projection needs to satisfy the following equation: a h E ( Π E v , p ) = a h E ( v , p ) , p ∈ ( P E ) 2 (2.5) ∫ ∂ E Π E v d S = ∫ ∂ E v d S (2.6)

The defined vectors form is ∫ E ∇ m T ⋅ ∇ Π E ϕ d x = ∫ E ∇ m T ⋅ ∇ ϕ d x (2.7) P 0 ( ∇ Π E ϕ ) = P 0 ( ϕ ) (2.8) Where P 0 is the constant term projection operator (Sutton, 2017).

Due to P E ⊂ V h E , The equation can be obtained m = D ϕ (2.9) Where D represents the expression of the matrix under the basis function ϕ .

Additionally, on account of Π E ϕ ∈ P E ⊂ V h E , The equation can be obtained Π E ϕ i ∑ α = 1 N E k a i , α ϕ α = ∑ j = 1 N E k s i , j m j (2.10)

From the Eq. 2.9, the equation is obtained ∑ α = 1 N E k s α = ∑ j = 1 N E k D α , j a j (2.11)

The equations of Eq. 2.11 and Eq. 2.10 are brought into Eq. 2.2, and the following equation is obtained G ¯ S = B ¯ (2.12)

The matrixes of G ¯ and B ¯ are computed as follows G α , β ¯ = { P 0 ( m ) i f   β = 1 ∫ E ∇ m α ⋅ ∇ m β d x B β , j ¯ = { P 0 ( ϕ ) i f   β = 1 ∫ E ∇ ϕ α ⋅ ∇ m j d x

Eq. 2.10 can be written as a matrix expression as G ¯ = B ¯ D (2.13)

Bring Eq. 2.13 into Eq. 2.12 S = ( B ¯ D ) − 1 B ¯ (2.14)

The element stiffness matrix is (Beirão da Veiga et al., 2014) ( K E ) i j = a ( Π E u , Π E v ) + a ( u − Π E u , v − Π E v ) = ( ∇ Π E ϕ i , ∇ Π E ϕ j ) 0 , E + ( ∇ ( 1 − ∇ Π E ϕ i , ) , ∇ ( 1 − ∇ Π E ϕ j ) ) 0 , E (2.15)

The meaning of the Eq. 2.15 can refer to the articles (Chen, 2015; Ortiz-Bernardin et al., 2017).

Figure 1 presents a two-dimensional frictionless contact model. For domain Ω 1 and   Ω 2 , the possible contact boundaries are S 1 and S 2 , g represents the contact gap. let R 1 and R 2 are the contact force vectors in the contact region.

Consider d S 1 and d S 2 are the infinitesimal region where is coming into contact. The virtual work δ W done by the contact traction is δ W = ∫ S 1 R 1 ⋅ δ u 1 d S 1 + ∫ S 2 R 2 ⋅ δ u 2 d S 2 (3.1)

In the contact area, every point should satisfy the equilibrium equation R 1 d S 1 = R 2 d S 2 (3.2)

Thus, we could consider the integral in Eq. 3.1 along the contact line S 1 or S 2 . δ W = ∫ S 1 R 1 ⋅ ( δ u 1 − δ u 2 ) d S 1 (3.3)

In this paper, the NST contact model is employed. A mapping is defined in Figure 2.

In the Figure 2, the tangent vector and normal vector of the master element establish a local coordinate system ( n , t ) . The t can be written in matrix form as t = | x m + 1 − x m y m + 1 − y m | l

The n is n = | y m − y m + 1 x m + 1 − x m | l Where l   stands for the length of node m + 1 to m and ( x m ,   y m ) is coordinates of node m .

The projection node coordinates of node s is on the surface composed of nodes m and m + 1 is ξ = ( x m + 1 − x m ) ⋅ [ x s − x m ] l (3.4)

The relative displacement from node s to the corresponding segment could be calculated as u s − u ξ = [ 10 ξ 01 − ξ 0 010 ξ 01 − ξ ] [ x s y s x m y m x m + 1 y m + 1 ] = A ( ξ ) x (3.5)

The g s indicates the normal gap, and it can be computed by the following equation g s = n T ⋅ u (3.6)

Bring the Eq. 3.5 into Eq. 3.3, The contact integral can be calculated δ W = ∫ S 1 R 1 ⋅ δ ( u s − u ξ ) d S 1 (3.7)

The δ u 1 is replaced by δ u s and δ u 2 is replaced by δ u ξ . The contact force R 1 can be split into R 1 = N + T = N ⋅ n + T ⋅ t (3.8)

In the frictionless contact problem, the tangential part disappears ( T = 0 ) . Thus, the integral of Eq. 3.7 can be written as δ W = ∫ S 1 R 1 ⋅ δ ( u s − u ξ ) d S 1 = ∫ S 1 R 1 ⋅ n δ g s d S 1 = ∫ S 1 N δ g s d S 1 (3.9)

When the PM is employed for the frictionless contact problem, contact traction N can be written as N = ω N H ( − g s ) g s   w i t h   H ( − g s ) = { 0 i f   g s > 0 1 i f   g s ≤ 0 (3.10)

Thus, the Eq. 3.9 can be written as δ W = ∫ S 1 ω N H ( − g s ) g s δ g s d S 1 (3.11)

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 δ W . This method is called linearization.

Differentiating Eq. 4–13 with respect to time t , we get d δ W d t = d δ W N d u d u d t = ω N f ( ξ 3 ) ξ 3 ˙ δ ξ 3 = ω N f ( ξ 3 ) ( v s − v m ) ⋅ n ⋅ δ x T ⋅ A T ⋅ n = ω N f ( ξ 3 ) δ x T A T n ⊗ n ⋅ A ⋅ v (3.12) Without considering the rotation of the contact body, the contact matrix is K c = D v ( δ W ) = ω N H ( − g s ) [ A ( ξ ) ] T n ⊗ n [ A ( ξ ) ] (3.13)

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 L = 40 units , b = 8 units , M = 16 units , E = 1000   units , u = 0   a n d   h = 1   uni , 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 2 d / b (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.

In this case, the value which is calculated point A is compared with the exact value which is M L 2 / ( 2 E I z ) to illustrate the sensitivity of different numerical methods to the quality of the mesh. In Figure 4, the following conclusions can be obtained:

For the twist factor 2 d / b 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 2 d / b is near 1. When the distortion parameter 2 d / b exceeds 1, the accuracy of VEM is superior to that of FEM.

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 R = 10   m is loaded by a pressure p = 100 P a . The μ = 0.3 and E = 10 3 P a is selected as material parameters. Rectangular at the bottom of the model is taken as L = 20 m , H = 10   m , E = 10 6 P a , μ = 0.3 . The quadrilateral meshes are used. Because of the model’s symmetry, half of the models are selected for research.

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 ω N 1 = 1 × 10 4 . The model is discretized by 1,155 quadrilateral elements, as exhibited in Figure 6A. The contour of normal pressure is exhibited in Figure 6B.

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:

In the contact region where x 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 x 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.

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 ( u y = − 0.1   m ) . The displacements of the bottom of the model are both fixed. The Young’s modulus is E = 10 9   P a . The Poisson’s ratio is μ = 0.3 .

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 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 ( σ y y ) is shown in Figure 11C.

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 L 1 = 25 m , L 2 = 3 m , L 3 = 2 m , H 1 = 10 m , H 2 = 5 m . The Young’s modulus E = 10 10 P a and Poisson’s ratio μ = 0.3 are the material parameters for this model. The γ = 24 k N / m 3 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.

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 x -direction is demonstrated in Figure 13A, and the displacement contour in y -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.