Semi-Infinite Structure Analysis with Bimodular Materials with Infinite Element

The modulus of elasticity of some materials changes under tensile and compressive states is simulated by constructing a typical material nonlinearity in a numerical analysis in this paper. The meshless Finite Block Method (FBM) has been developed to deal with 3D semi-infinite structures in the bimodular materials in this paper. The Lagrange polynomial interpolation is utilized to construct the meshless shape function with the mapping technique to transform the irregular finite domain or semi-infinite physical solids into a normalized domain. A shear modulus strategy is developed to present the nonlinear characteristics of bimodular material. In order to verify the efficiency and accuracy of FBM, the numerical results are compared with both analytical and numerical solutions provided by Finite Element Method (FEM) in four examples.


Introduction
It has been shown that certain materials such as composites, porous materials, rocks, cement concrete and asphalt concrete, etc., show significant differences in their strength in tension and compression states. The modulus of elasticity as well as the Poisson's ratio of the material may also change under tensile and compressive states [1][2][3]. Take the concrete material as an example; the compressive modulus is about 1.5~2 times the tensile modulus [4][5][6]. So, for an accurate numerical simulation, this characteristic of material has to be considered. It constructs a typical material nonlinear model.
In order to evaluate bearing capacity and stability, the civil structure with the soil-foundation interaction is commonly investigated numerically, including airport runways, highway pavement, stacking dock, mineral deposit, geotechnical slope and so on. The soil medium is simplified as an infinite or semi-infinite domain. The most common approach with FEM is to use massive elements to simulate an unbounded domain. The application of large-scale finite element discretization could result in an increase in computational burden [7]. Furthermore, the inaccurate results could be obtained due to the truncated boundaries in the numerical procedure. To overcome this difficulty, the Boundary Integral Equations Method (BIEM), also known as the Boundary Element Method (BEM), is coupled with the FEM [8,9]. However, it is difficult to derive the fundamental solutions in general cases, especially for non-homogeneous and nonlinearity of materials. Meanwhile, the semi-analytical finite element method was developed to reduce the time cost of 3D model simulation [10,11] and applied in pavement structural analysis [12,13], but it mainly focuses on linear analysis or problems without complicated loads. The unbounded problems can be overcome by introducing mapped infinite elements, i.e., utilizing the infinite element to extend the FEM to unbounded domain problems [14][15][16][17]. The shape function describes the far-field characteristic of the problem, which can be obtained using mapping to transform the global infinite region into a local finite domain by Bettess et al. [17][18][19][20]. As an alternative, these issues can be solved with the meshless methods coupling with an infinite-mapping technique [7].
In engineering analysis, the linear elasticity of material is not valid for general issue. The material mechanical properties are closely related to their micro structure. The scanning images of the building materials are shown in Figure 1 and present similar mottled patterns at different scales. The heterogeneity is manifested in the micro-scale for the metal materials, and its mechanical properties accord with the linear elastic hypothesis. For the rock or concrete materials, their heterogeneity is displayed in the mesoscale and the assumption of linear elasticity sometimes produces computational errors which cannot be ignored. problems can be overcome by introducing mapped infinite elements, i.e., utilizing the infinite element to extend the FEM to unbounded domain problems [14][15][16][17]. The shape function describes the far-field characteristic of the problem, which can be obtained using mapping to transform the global infinite region into a local finite domain by Bettess et al. [17][18][19][20]. As an alternative, these issues can be solved with the meshless methods coupling with an infinite-mapping technique [7].
In engineering analysis, the linear elasticity of material is not valid for general issue. The material mechanical properties are closely related to their micro structure. The scanning images of the building materials are shown in Figure 1 and present similar mottled patterns at different scales. The heterogeneity is manifested in the micro-scale for the metal materials, and its mechanical properties accord with the linear elastic hypothesis. For the rock or concrete materials, their heterogeneity is displayed in the mesoscale and the assumption of linear elasticity sometimes produces computational errors which cannot be ignored. Commercial numerical software in engineering, including ABAQUS, are widely used in engineering and manufacturing. However, it is still a challenging task to solve bimodular problems efficiently [21][22][23][24][25][26]. Nevertheless, the development of new numerical methods is always attractive to solve difficult and complicated engineering problems. Unlike the traditional numerical method, the computational framework of the meshless method was based on the scattered nodes. In the 1990s, the meshless method was developed based upon the Galerkin method. In 1992, the diffuse element method (DEM) was proposed by Nayroles et al. [27]. The Moving-Least Square (MLS) method was introduced to construct the meshless shape functions with Galerkin method in numerical discretization. In 1994, Belytschkoet al. presented the Element-Free Galerkin method (EFGM) [28], in which Lagrange was employed to ensure the boundary conditions were being satisfied. Since then, the EFGM has been widely used to simulate the fracture failure of materials and to show its superiority over the traditional FEM [29,30]. In 1996, Belytschko et al. published a comprehensive review [31] which attracted exclusive attention in computational mechanics. This can be regarded as the beginning of the meshless method in numerical engineering. Another important development was the introduction of the local weak form methods. In 1998, Atluri et al. proposed the Meshless Local Petrov-Galerkin(MLPG) method [32]. The discrete system equation is based on a nodal assembly with more conciseness in numerical implementation. In 1995, Liu et al. proposed a Reproducing Kernel Particle Method (RKPM) approximation [33][34][35]. Thereafter, several meshless methods were developed such as the Method of Fundamental Solution (MFS) [36][37][38], the local Radial Point Interpolation Method (RPIM) [39][40][41], the local Radial Basis Function (RBF) collocation method [42][43][44] and the Meshless Intervention-Point (MIP) method [45], etc. In 2014, Wen et al. proposed the meshless FBM [46]. In the finite block method, the mapping technique is implemented numerically with the infinite elements for the infinite domain problems [7]. Afterwards, the FBM is successfully applied to nonlinear elasticity problems, contact Commercial numerical software in engineering, including ABAQUS, are widely used in engineering and manufacturing. However, it is still a challenging task to solve bimodular problems efficiently [21][22][23][24][25][26]. Nevertheless, the development of new numerical methods is always attractive to solve difficult and complicated engineering problems. Unlike the traditional numerical method, the computational framework of the meshless method was based on the scattered nodes. In the 1990s, the meshless method was developed based upon the Galerkin method. In 1992, the diffuse element method (DEM) was proposed by Nayroles et al. [27]. The Moving-Least Square (MLS) method was introduced to construct the meshless shape functions with Galerkin method in numerical discretization. In 1994, Belytschkoet al. presented the Element-Free Galerkin method (EFGM) [28], in which Lagrange was employed to ensure the boundary conditions were being satisfied. Since then, the EFGM has been widely used to simulate the fracture failure of materials and to show its superiority over the traditional FEM [29,30]. In 1996, Belytschko et al. published a comprehensive review [31] which attracted exclusive attention in computational mechanics. This can be regarded as the beginning of the meshless method in numerical engineering. Another important development was the introduction of the local weak form methods. In 1998, Atluri et al. proposed the Meshless Local Petrov-Galerkin (MLPG) method [32]. The discrete system equation is based on a nodal assembly with more conciseness in numerical implementation. In 1995, Liu et al. proposed a Reproducing Kernel Particle Method (RKPM) approximation [33][34][35]. Thereafter, several meshless methods were developed such as the Method of Fundamental Solution (MFS) [36][37][38], the local Radial Point Interpolation Method (RPIM) [39][40][41], the local Radial Basis Function (RBF) collocation method [42][43][44] and the Meshless Intervention-Point (MIP) method [45], etc. In 2014, Wen et al. proposed the meshless FBM [46]. In the finite block method, the mapping technique is implemented numerically with the infinite elements for the infinite domain problems [7]. Afterwards, the FBM is successfully applied to nonlinear elasticity problems, contact problems and heat conduction problems [47][48][49]. It has been demonstrated in the analysis of bimodular problems for two-dimensional problems [50].
In this paper, the FBM is extended to three-dimensional semi-infinite structures in bimodular materials. The infinite block mapping technique is introduced to present the semi-infinite structure and implemented with the meshless finite block method to construct the intrinsic constitutive equations in iterative analysis. The meshless finite block method with the infinite block mapping technique is formulated for 3D bimodular problems. The FEM solution is considered as a benchmark for numerical analysis, and the accuracy of the proposed method is observed by ABAQUS with subroutine UMAT developed for bimodular materials.

Bimodular Material Constitutive Equations
Suppose σ α , σ β and σ γ are principal stresses, as shown in Figure 2. The generalized Hooke's law, in matrix form, asε = Aσ · · · or · · ·σ = Q Iε (1) where A is the flexibility matrix, Q I is the elasticity matrix,ε is the nodal strain vector in the principal directions andσ is the nodal stress vector in the principal directions, which are defined as There are three cases to obtain  G ,  G and  G , (1) If all three principal stresses are equal, i.e. With the analytical theory proposed by Ambartsumyan and complemented with shear moduli [1,21,22], it is assumed that a ij = −v − /E − = −v + /E + , a jj = 1/E + or 1/E − , (I = 1, j = 1, 2, 3), where E + and E − present as the tensile and compressive moduli respectively, v + and v − are the tensile and compressive Poisson's ratio, respectively; a 44 = 1/G βγ , a 55 = 1/G αγ , a 66 = 1/G αβ , in which, G βγ , G αγ and G αβ are the shear moduli. The shear stresses or strains in the principal directions are zero. According to the shear moduli algorithm [13], it is assumed that the axes x, y and z tend to axes α, β and γ, respectively. Then, we have There are three cases to obtain G βγ , G αγ and G αβ , (1) If all three principal stresses are equal, i.e., σ α = σ β = σ γ , we have a.

Lagrange Polynomial Interpolation
Consider a 3D square in normalized domain mapping to the physical domain, as shown in Figure 3. The Lagrange polynomials in the coordinate system (ξ, η, ζ) are used to interpolate function u where u p indicates the nodal value, subscript p denotes the number of node at P(ξ i , η i , ζ i ) in the global system and functions where N ξ , N η and N ζ denote the numbers of node distributed along the axes ξ, η and ζ, respectively. The shape function is obtained simply as Materials 2022, 14, x FOR PEER REVIEW 6 of 28

Partial Differential Matrix
The partial derivative of function u in Equation (15) can be arranged in a vector. For example, the nodal first order partial derivative of function u can be written, in the vector form, as where p is the number of node P(i, j, k) in the global system; ) ( the number of nodes in the local coordinate system, and The partial differential with respect to axis ξ can be obtained directly

Partial Differential Matrix
The partial derivative of function u in Equation (15) can be arranged in a vector. For example, the nodal first order partial derivative of function u can be written, in the vector form, as where p is the number of node P(i, j, k) in the global system; M(= N ξ × N η × N ζ ) indicates the number of nodes in the local coordinate system, and In addition, the L-th order partial derivative with respect to the coordinates ξ, η and ζ can be approximated as Therefore, the higher-order partial differentials in Equation (22) can be obtained, in terms of the first-order partial derivative matrices D ξ , D η and D ζ , as

Mapping Differential Matrix
For three-dimensional problems, a hexahedron block with 20 seeds is selected in order to transform the coordination (x, y, z) to (ξ, η, ζ) as shown in Figure 3. The mapping function is expressed as The partial differentials of function u(x, y, z) with subject to axis ξ, η or ζ can be written as  (25) Then the partial differentials of the function u(x, y, z) with respect to x, y and z are given by  (26) in which β ij express the terms in the cofactor of Jacobi matrix J, and Therefore, the first order partial differential in the physical domain can be written as in which where β (1) ij / J (1) can be determined from Equation (27) at each node in the normalized domain, and the first order differentials matrix is determined by the Lagrange interpolation functions in normalized domain (|ξ| ≤ 1, |η| ≤ 1, |ζ| ≤ 1).

Mapping Technology with 3D Blocks
For the semi-infinite structure shown in Figure 4a, the semi-infinite domain is divided into several subdomains with two 20-seed-finite blocks, two 12-seed-one-infinite-edge blocks, two 7-seed-two-infinite-edge blocks and two 8-seed-three-infinite-edge blocks as shown in figures from Figure 4b-e. The infinite blocks in different directions can be obtained by rotating the initial mapping function. The mapping function for the finite block and infinite blocks in a general form is written as where q is the seed number shown in Figure 4. The details of the mapping function and their partial differentials can be presented in Appendix A in different categories.
can be determined from Equation (27) at each node in the normalized domain, and the first order differentials matrix is determined by the Lagrange interpolation functions in normalized domain

Mapping Technology with 3D Blocks
For the semi-infinite structure shown in Figure 4a, the semi-infinite domain is divided into several subdomains with two 20-seed-finite blocks, two 12-seed-one-infiniteedge blocks, two 7-seed-two-infinite-edge blocks and two 8-seed-three-infinite-edge blocks as shown in figures from Figure 4b-e. The infinite blocks in different directions can be obtained by rotating the initial mapping function. The mapping function for the finite block and infinite blocks in a general form is written as where q is the seed number shown in Figure 4. The details of the mapping function and their partial differentials can be presented in appendix A in different categories.

Formulations for Bimodular Material with Meshless FBM
The equilibrium equation, in the domain, gives

Formulations for Bimodular Material with Meshless FBM
The equilibrium equation, in the domain, gives in which σ αβ , (α, β = x, y, z) denotes stress; f α are body force. Substituting the constitutive equation, Equation (1), into the kinematic equation in Equation (33) without body forces yields C 11 u x + C 12 u y + C 13 u z = 0, where u x , u y , u z are vectors of nodal displacements, and C ij , (i, j = 1, 2, 3) are coefficients by the constitutive and equilibrium equations, and given by where Q ij , (i, j = 1, 2, · · · , 6, Q ij = Q ji ) are the terms in elasticity matrix Q and given by Consider the following boundary conditions defined as where t(x) and u(x) are given traction and displacement on the boundary, t(x) = t x , t y , t z T , x is the collocation point on the boundary. Traction t(x) can be rewritten as where matrix B ij , (i, j = 1, 2, 3) is associated with the boundary collocation point B 11 = D x (Q 11 n x + Q 16 n y + Q 15 n z ) + D y (Q 16 n x + Q 66 n y + Q 56 n z ) + D z (Q 15 n x + Q 56 n y + Q 55 n z ), B 22 = D x (Q 66 n x + Q 26 n y + Q 46 n z ) + D y (Q 26 n x + Q 22 n y + Q 24 n z ) + D z (Q 46 n x + Q 24 n y + Q 44 n z ), B 33 = D x (Q 55 n x + Q 45 n y + Q 35 n z ) + D y (Q 45 n x + Q 44 n y + Q 34 n z ) + D z (Q 35 n x + Q 34 n y + Q 33 n z ), B 12 = D x (Q 16 n x + Q 66 n y + Q 56 n z ) + D y (Q 12 n x + Q 26 n y + Q 25 n z ) + D z (Q 14 n x + Q 46 n y + Q 45 n z ), B 13 = D x (Q 15 n x + Q 56 n y + Q 55 n z ) + D y (Q 14 n x + Q 46 n y + Q 45 n z ) + D z (Q 13 n x + Q 36 n y + Q 35 n z ), B 21 = D x (Q 16 n x + Q 12 n y + Q 14 n z ) + D y (Q 66 n x + Q 26 n y + Q 46 n z ) + D z (Q 56 n x + Q 25 n y + Q 45 n z ), B 23 = D x (Q 56 n x + Q 25 n y + Q 45 n z ) + D y (Q 46 n x + Q 24 n y + Q 44 n z ) + D z (Q 36 n x + Q 23 n y + Q 34 n z ), B 31 = D x (Q 15 n x + Q 14 n y + Q 13 n z ) + D y (Q 56 n x + Q 46 n y + Q 36 n z ) + D z (Q 55 n x + Q 45 n y + Q 35 n z ), B 32 = D x (Q 56 n x + Q 46 n y + Q 36 n z ) + D y (Q 25 n x + Q 24 n y + Q 23 n z ) + D z (Q 35 n x + Q 34 n y + Q 33 n z ) where n α , (α = x, y, z) is the boundary outwards normal. Therefore, 3 × M linear algebraic equations are obtained in total from Equations (33) and (38). In addition, on the interfaces between blocks, the following continued conditions should be taken into account where u (i) α and t (i) α represent the displacement and traction on the interface between block i and block j. Finally, a set of linear algebraic equations is established the in global system as follows Where K is the stiffness matrix, U is the vector of displacements and F is the vector consisting of the boundary value of the displacement, tractions and domain body forces. The following nonlinear iterative algorithm is adopted in this paper.
Step 1: m = 0, take either tensile or compressive modulus at all collocation points. Solve the global stiffness matrix to obtain the initial displacements, stresses and strains.
Step 3: Modify the stiffness matrix K and vector F based on the current step. Solve the equations again to obtain the displacements, stresses and strains at each node.
Step 4: Calculate the average error from all collocation points where U i (m) presents the displacement at step m. if κ < 10 −6 , terminate the iteration and print out the result. Otherwise, let m = m + 1; go to step 2.

Numerical Examples
In this section, four examples are presented to demonstrate the accuracy of the meshless FBM with bimodular materials. A 3D tensile column with gravity is investigated in the first example. Then, FBM is applied to an arch bridge model, a single-layer semi-infinite model and a multi-layer pavement foundation under different loadings. All codes were written with Matlab (R2021b, The MathWorks, Inc., Natick, MA, USA).and Fortran in subroutine UMAT using ABAQUS (2019, Dassault Systèmes Simulia Corp., Vélizy-Villacoublay, France).

Tensile Column with Gravity
Consider a gravitational column of the length l = 2; the dimension of the cross-section is normalized as 1 × 1, and the mass density γ = 2 as shown in Figure 5a. It is fixed on the bottom and a tensile force P of 2 units is applied to the top. It is assumed that a compression modulus is 5000 units, and the Poisson's ratios in tension and compression is zero. The exact solution of the displacement [1] along the z-axis is given as where c = l − P/γ. The numbers of node in x-axis and y-axis are 9, and in the z-axis is 14. The locations of node along different axes in the normalized domain are chosen  9.0 × 10 −3 8.9 × 10 −3 2 2  The total number of nodes for the FBM is 1134 (= 9 × 9 × 14), and 396 C3D20R elements are used in FEM. The node distribution of FBM is shown in Figure 5b. Comparison between the exact solution and FBM solution at point z = 1.96 and the number of iterations for convergence between FEM and FBM are presented in Table 1. With different ratios of tensile and compression modulus, the vertical displacement changes along the z-axis and exact solution are shown in Figure 6. Obviously, the FBM can give an accurate solution for the problem and shows a similar convergence rate when compared with the FEM method. To investigate the accuracy for different node density, the average relative errors are defined as 9.0 × 10 −3 8.9 × 10 −3 2 2  The numerical results presented in Table 2 demonstrate the average errors with iteration numbers of convergence over all collocation points when E − /E + = 10. Observing the results in Table 2, it is evident that increasing the node density improves the degrees of accuracy, and convergency is easily approached in iterations when the node number N ξ is more than 3.

Arch Bridge in Bimodular Materials
Consider a simplified arch bridge as shown in Figure 7. Due to the symmetry of the structure, half of the model is taken for analysis. The radius of the arc is a = 1 unit. There is a vertical pressure load p 0 of 1 unit applied on the top, and the lengths in both y-axis and xaxis are w(=2a). The displacement u y is fixed on the bottom face (y = 0), and u x is zero on the surface x = 0. The ratios of Young's moduli are selected as E − /E + = 1, 2, 5, compression modulus E − = 1 unit and Poisson's ratio v − = 0.4 in the computation procedure.
structure, half of the model is taken for analysis. The radius of the arc is a = 1 unit. There is a vertical pressure load  The bridge is divided into three blocks using FBM shown in Figure 7, where blocks I and II are finite blocks with 20-seed and block III is one semi-infinite block with a 12-seedone-infinite-edge. In the discretization of each block, there are 12 and 14 collocation nodes along finite and infinite directions, respectively. The distribution of nodes along each axis is the same as Equation (45)   The bridge is divided into three blocks using FBM shown in Figure 7, where blocks I and II are finite blocks with 20-seed and block III is one semi-infinite block with a 12seed-one-infinite-edge. In the discretization of each block, there are 12 and 14 collocation nodes along finite and infinite directions, respectively. The distribution of nodes along each axis is the same as Equation (45) in Section 5.1, as shown in Figure 8a. Stresses along two segments, AB and CD, shown in Figure 7 are plotted to illustrate the degree of accuracy. Simulation with FEM is complemented with 90,912 C3D10 elements as shown in Figure 8b. The length in the x-axis is w = 40 unit. The normalized stress σ x along AB, CD and AC by FBM and FEM is plotted in Figure 9 to show the difference between these two methods with bimodular materials. Reasonable agreements can clearly be observed. It is also noticed that there are several kinks in Figure 9 for stress distributions by the FEM due to the discontinuity of the Young's modulus.

A Semi-Infinite Solid with Bimodular Materials
The semi-infinite structures are introduced to simulate soil foundations. Consider a semi-infinite body as shown in Figure 10a with the linear distributed vertical load in a square area of width 1 unit on the surface. The linear distributed load is plotted in Figure 10b with a unit maximum absolute value of q in compression and tension. Bimodular materials are selected with three different ratios of tensile and compressive moduli, as shown in Table 3. Due to the symmetry of the structure and loading, only a half model is analyzed as shown in Figure 10a. To accurately capture the stress near the loading area, the structure is subdivided into two layers. The first layer includes one 20-seed finite block III, three 12-seed-one-infinite-edge blocks I, IV, V and two 7-seed-two-infinite-edge block II and VI. However, in the second layer, one 12-seed-one-infinite-edge block, three 7-seed-two-infinite-edge blocks and two 8-seed-three-infinite-edges blocks are used. For each block, 9 collocation points are used on the finite edge and 12 points for the infinite edge. Normalized stress σ x along AB and AC are presented to demonstrate the accuracy of the FBM shown in Figure 11a,b versus the different ratios of tensile and compressive moduli, and Poisson ratios. In this example, FEM simulation is complemented by use of 362,484 C3D10 elements with dimensions of 20 units in length and height and 10 units in the width. A reasonable agreement was clearly achieved.
AC by FBM and FEM is plotted in Figure 9 to show the difference between these two methods with bimodular materials. Reasonable agreements can clearly be observed. It is also noticed that there are several kinks in Figure 9 for stress distributions by the FEM due to the discontinuity of the Young's modulus.

Multi-Layered Infinite Model with Bimodular Materials
Consider a multi-layered infinite structure, as shown in Figure 12, to simulate a highway pavement structure under two symmetric circular pressure loads. The pressure is assumed to be 0.7 MPa with a radius of 0.1065 m. The distance between two centers of loads is 0.3195 m. The model contains four layers; the first and second layers are bimodulus materials and the third and fourth layers are isotropic materials. The details of material parameters and dimensions of each layer are listed in Table 4. Again, due to the symmetry of the structure and load condition, quarter of the structure is analyzed as shown in Figure 12. During the numerical process, each layer is divided into four blocks. For the first layer, the top layer contains one 20-seed finite block, two 12-seed-infinite-edge blocks II and III and one 7-node-two-infinite-edge block IV. In the second and third layers, the same block distribution is applied as in the first layer. The bottom layer contains one 12-seed-one-infinite-edge block I, two 7-seed-two-infinite-edge blocks II and III and one 8-seed-three-infinite-edge block IV.

A Semi-Infinite Solid with Bimodular Materials
The semi-infinite structures are introduced to simulate soil foundations. Consider a semi-infinite body as shown in Figure 10a with the linear distributed vertical load in a square area of width 1 unit on the surface. The linear distributed load is plotted in Figure  10b with a unit maximum absolute value of q in compression and tension. Bimodular materials are selected with three different ratios of tensile and compressive moduli, as shown in Table 3. Due to the symmetry of the structure and loading, only a half model is analyzed as shown in Figure 10a. To accurately capture the stress near the loading area, the structure is subdivided into two layers. The first layer includes one 20-seed finite block III, three 12-seed-one-infinite-edge blocks I, IV, V and two 7-seed-two-infinite-edge block II and VI. However, in the second layer, one 12-seed-one-infinite-edge block, three 7-seedtwo-infinite-edge blocks and two 8-seed-three-infinite-edges blocks are used. For each block, 9 collocation points are used on the finite edge and 12 points for the infinite edge. Normalized stress x  along AB and AC are presented to demonstrate the accuracy of the FBM shown in Figure 11a,b versus the different ratios of tensile and compressive moduli, and Poisson ratios. In this example, FEM simulation is complemented by use of 362,484 C3D10 elements with dimensions of 20 units in length and height and 10 units in the width. A reasonable agreement was clearly achieved.

Multi-Layered Infinite Model with Bimodular Materials
Consider a multi-layered infinite structure, as shown in Figure 12, to simulate a highway pavement structure under two symmetric circular pressure loads. The pressure is assumed to be 0.7MPa with a radius of 0.1065m. The distance between two centers of loads is 0.3195m. The model contains four layers; the first and second layers are bimodulus materials and the third and fourth layers are isotropic materials. The details of material parameters and dimensions of each layer are listed in Table 4. Again, due to the symmetry of the structure and load condition, quarter of the structure is analyzed as shown in Figure  12. During the numerical process, each layer is divided into four blocks. For the first layer, the top layer contains one 20-seed finite block, two 12-seed-infinite-edge blocks II and III and one 7-node-two-infinite-edge block IV. In the second and third layers, the same block distribution is applied as in the first layer. The bottom layer contains one 12-seed-oneinfinite-edge block I, two 7-seed-two-infinite-edge blocks II and III and one 8-seed-threeinfinite-edge block IV.
Like the node distributions in Section 5.3, the 8 seeds are used on the finite edges and 14 seeds on the infinite edges. The total number of collocation nodes by FBM is 12288. To   Like the node distributions in Section 5.3, the 8 seeds are used on the finite edges and 14 seeds on the infinite edges. The total number of collocation nodes by FBM is 12288. To validate the computational accuracy, the results of stresses σ z by FBM and FEM along segment AB and segment CD are compared in Figure 13. The contours of von Mises stress with bimodular materials on y = 0 are presented by using FBM in Figure 14. FEM is also used for analysis with no dimension of 20 × 20 × 20 and 127,760 C3D20R elements used in this example. It can be seen that the position of the maximum von Mises stress with these two methods is the almost the same, and the values are also very close to each other. In addition, the FBM results are smoother.

Conclusions
A meshless finite block method with infinite block analyzing three-dimensional solids of bimodular materials was demonstrated in this paper. A mapping technique was applied to determine the first order of derivatives. The 20-node finite block, 12-seed-one-edgeinfinite block, 7-seed-two-edge-infinite block and 8-seed-three-edge-infinite block were introduced to simulate all semi-infinite domains. An iterative process for the meshless finite block method with a shear-modulus-complemented algorithm to solve bimodular problems was proposed. The numerical algorithm was validated with four examples. The finite element software ABAQUS was used for comparison. The following conclusions can be summarized: (1) FBM easily tackles nonlinear problems with semi-infinite boundaries; (2) a shear modulus algorithm efficiently and accurately describes the bimodular mechanical behavior of materials; (3) the proposed method shows efficiency and accuracy for semiinfinite problems with bimodular materials. Compared to FEM, FBM is more accurate with the same computational effort; (4) FBM can be applied to more complicated problems, such as 3D elastoplasticity, thermoelasticity and elastodynamics.
FE methods are rather general and efficient numerical tools to deal with complicated problems in engineering. However, as an alternative, the meshless finite block method with an infinite-mapping technique provides a new approach in solving unbounded bimodular material problems, with many advantages including efficiency and simplicity. As ABAQUS is a commercialized package, the CPU times used by different approaches are not comparable in this work. At present, dividing blocks is still a manual process in FBM; the versatility needs to be further improved with complex regional models. In future work, the FBM is expected to be extended to apply to more complicated problems, such as 3D elastoplasticity, thermoelasticity and elastodynamics.

Conflicts of Interest:
The authors declare no conflict of interest.

20-node finite block
For this type of finite element, the physical domain is mapped to a cube with 20 seeds in coordination system (ξ, η, ζ) in the region |ξ| ≤ 1, |η| ≤ 1 and |ζ| ≤ 1, as shown in Figure 4b. The mapping function can be written as follows [51]:  Their partial differentials of the mapping function are listed below: 2.

12-seed-one-edge-infinite block
In the normalized domain, the face of the upper side (ζ = 1) is mapped to infinite area as shown in Figure 4c. The mapping functions [17] are The first-order partial differentials of Equations (A9)-(A12) are