Assume is a homogeneous elastic domain with a uniform Young’s modulus E and Poisson’s ratio ν. The governing equations for this linear elastic problem are written as the following:with u_i represents the displacements, ε_ij denotes the strains, σ_ij represents the stresses, b_i represents the body forces, ū_i the prescribed displacement, represents the prescribed traction force, Γ_c the boundary domain for the essential boundary condition, Γ_t the boundary domain for the nature boundary condition, where Γ_t ∩ Γ_c = Ø, Γ_t ∪ Γ_c = Γ ≡ ∂Ω. λ and μ are Lamé constants withThere are numerous closed-form analytical solutions in terms of elementary functions in history. Nonetheless, when the domain Ω is complex, obtaining an analytical solution becomes challenging. The linear finite element method (LFEM) has been proven to be an effective and accurate method for obtaining an approximate solution for the equations Eq. (1). In the LFEM, the domain Ω is discretized into subdomains named “elements” where Ω = ∪_eΩ^e, Ω^e denotes the domain for eth element. Notably, in this manuscript, “element” refers to a solid or continuum element, rather than a structural element such as a beam or shell. Within one solid element, one can evaluate the so-called element stiffness matrix K^e such thatwhere subscription “I” denotes the index of the quadrature point, B_I represents the strain-displacement matrix which interpolates strain by the elemental displacement vector, J_I is the Jacobian matrix defined as J ≡ ∂x/∂ξ, and ξ represents the isoparametric coordinates. The parameter w_I represents the quadrature weight, andis the elastic stiffness tensor in Voigt notation {·}.

With given nodal connectivity of all the elements, we can assemble the global stiffness matrix K, and we may solve the following sparse linear equation systems to obtain displacements thatwhere d is the displacements to solve; ū denotes the constraint displacements; f and are the prescribed nodal forces at unknown and constraint displacement degrees-of-freedom, respectively; and r are the reaction forces at the constraint displacements. Thus, one can construct the linear equations respect to d thatOnce rigid-body motions and rotations are constrained by proper essential boundary conditions, there will be a unique solution for d.

Under the framework of the LFEM, we aim to derive the element stiffness and then the global stiffness matrix in terms of material parameters. We first rewrite the isotropic elastic stiffness tensor Eq. (4) intowhereClearly, H₀ is symmetric and positive-definite, whereas H₁ is symmetric only. We also haveIn addition, we can evaluate the spectral radius thatsince the Poisson’s ratio ν ∈ (–1, 0.5), where ρ(A) denotes the spectral radius of a matrix A. Here we define a so-called shifted Poisson’s ratio such thatto simplify the expression. Thus, we havewhere 𝕄 denotes the linear elastic compliance tensor. Then the element stiffness matrix is written aswithBoth and depend only on the nodal coordinates of this element but irrelevant to elastic parameters. In addition, both and have the same matrix dimension as K^e, and the dimensions of both and are length [L].

One can assemble the corresponding global matrix Eq. (5) with R₀ and R₁ such thatThus, Eq. (6) can be rewritten with respect to R₀ and R₁ such thatorAgain that R_0,uu, R_1,uu, R_0,uc, and R_1,uc are matrices that depends on the finite element mesh only and irrelevant to the material parameters. In other words, one can evaluate these matrices with a given finite element mesh and essential boundary Γ_c. In addition, if K_uu is invertible, we also have R_0,uu as invertible.

Similarly, we haveThe displacement solution is then given asWe can rearrange the solution into the following series form:withandThe dimensions of a_k and b_k are length [L] and force per length [F/L], respectively. In the numerical implementation, it is recommended to solve a₀, a₁, … sequentially as well as b₀, b₁, …. All the coefficient vectors a_k and b_k are solved by the same coefficient matrix R_0,uu. During the numerical implementation, one may store a proper matrix factorization of R_0,uu in memory to avoid repetitive matrix factorization.

We may continue to evaluate strain and stress tensors at one quadrature point such thatwhere d^e is the element displacement vector. Alternatively, we can rewrite Eq. (23) in matrix form such thatwhere n denotes the number of elemental degrees-of-freedom. The stress tensor is then evaluated byWe may also obtain the reaction force at constraint degrees-of-freedom byConsequently, the strains and stresses at the quadrature points as well as the reaction force at constraint degrees-of-freedom can be interpreted with series with different orders of such as Eq. (20).

The shifted Poisson’s ratio may not be straightforward enough to express the various types of solutions. Thus, we can finally rewrite Eq. (20) into a series solution on the basis of various orders of Poisson’s ratio such that:Similarly, the strains at the quadrature point arewhile the stresses are rearranged intowith . The reaction forces as follows:

In practice, we may prefer a truncated solution with a finite number of terms of different orders of the Poisson’s ratio. We defineand the truncated strains, stresses, and reaction forces can be defined in the same manner.

We start with a patch test to verify whether an exact analytical solution can be derived. We use 8-node hexahedron element with full integration scheme to mesh the plate as depicted in Fig. 1. We define the length (along the x axis), width (along the y axis), and height (along the z axis) as a = 2, b = 1, and h = 0.2, respectively. We assign u = 0 at the surface x = 0; v = 0 at the surface y = 0; w = 0 at the surface z = 0; and u = ū = 0.01 at the surface x = a = 2.

The patch test gives uniform strains and stresses over the whole domainand linearly distributed the displacement field such thatConsequently, there are only two nonzero terms ã₀ and ã₁ for the approximated analytical solution d, and match the numerical results. We may plot the values of ã₀ and ã₁ for displacement u, v, and w with respect to x, y, and z, respectively, as depicted in Fig. 2.

We continue to derive an approximate analytical solution for a purebending problem depicted in Fig. 3. The length, width, and height of the beam are L, b, and h, respectively. The beam is meshed by 8-node hexahedron elements.

We consider a beam with size L = 4, and h = 2. Thus, this beam actually is not a typical slender beam structure. Here we use cubic-shaped 8-node hexahedron element to mesh the beam with element size l_el = 1/4, 1/8, and 1/16. We set the width as the element size b = 2l_el, and keep M_y/b = 2, where M_y is the moment about the y axis. The CMAS-FE is then able to evaluate the coefficients ã_k and , k = 0, 1, … for each problem.

We first study how the truncated solution d_N converges to the reference solution via the LFEM named the d_LFEM as N increases, which can be marked as series-convergence. Because the LFEM requires values of Young’s modulus and Poisson’s ratio rather than symbols; we assign E = 2.1 × 10⁵ MPa and ν = 0.23. Accordingly, we can obtain numerical values of each d_N and compare it with the d_LFEM. We take the case l_el = 1/8 as an example. Fig. 4(b)–(d) depicts the errors between d₀ – d_LFEM, d₁ – d_LFEM, and d₂ – d_LFEM, respectively, while the reference solution d_LFEM is depicted in Fig. 4(a). Clearly the error decreases to a very small even N only reaches 2 except around the corner regions. Similar observations can be found for the truncated strain and stress fields, which are depicted in Figs. 5 and 6, respectively.

The mesh convergence can also be observed in Fig. 7, where the error between d_N, N = 0, 1, 2 and the reference solutions d_LFEM with element sizes of l_el = 1/4, 1/8, and 1/16 are depicted. We can see that the magnitude of errors is reduced by the mesh refinement or the increasing number of truncated terms. Furthermore, the relative errors of the displacement, strain, and stress fields with increasing number of truncated terms N are plotted in Fig. 8(a)–(c), respectively.

Finally, there is an analytical solution for a pure-bending problem even when the beam is not a typical slender beam thatwhere I_y = bh³/12 is the moment of inertia about the y axis. Here the reference analytical solution in Eq. (34) is organized as the same expression as Eq. (27).

We take the element size l_el = 1/16 as example. Eq. (34) actually gives an analytical solution for and thatThe comparisons of the four coefficient vectors b_I,0 for u, b_I,1 for u, b_I,1 for v, and b_I,0 for w are depicted in Fig. 9(a)–(d), respectively. In addition, all the other coefficient vectors are close to 0, which also matches the analytical solution.

In CMAS-FE, geometric parameters such as the length of the structure cannot be considered directly, i.e., these parameters will not appear in the approximate analytical solution.

Next, we consider the cantilever-beam-bending test depicted in Fig. 10. In this example, we aim to fit an approximate analytical solution for bending stiffness with respect to geometric parameters L, h, and b. These geometric parameters may not directly appear in the truncated solution d_N. However, we will repeat various tests with different L and h values, while the width b can be set as 1 all the time. For a beam structure, we assume h/L < 0.6.

We define the deflection stiffness K aswhere F is the total reaction force at the right edge of the beam. Dimensionless coefficients c_k, k = 0, 1, 2 depends on the geometric parameters h and L. We use the CMAS-FE for all the pairs of (L, h) such that L = 5, 6, …, 20 and h = 2, 3, …, 6, and obtain F in series form. Next, we use base function (h/L)³ and (h/L)⁴ to approximate the coefficients c_k, k = 0, 1, 2 and obtainorHere the base function (h/L)^p for arbitrary p ∈ can be used to enrich the interpolation function. However, the numerical optimization process shows that the terms (h/L)³ and (h/L)⁴ are dominant. From Fig. 11, we may see that the interpolated coefficients in Eq. (37) are quite accurate compared with the results given by the CMAS-FE.

Next, we consider deriving an approximate analytical expression of stress intensity factor (SIF) for the problem depicted in Fig. 12, where a crack with length 2w is embedded in a square plate under remote stress in the y direction. From linear elastic fracture mechanics we know that the SIFwhen A → ∞. Here we use the CMAS-FE to derive the SIF with finite square size A. Finite element meshes with some selected square sizes are depicted in Fig. 13. We can obtain SIFs with A as listed in Table 1. Here the CMAS-FE can determine the values for and defined as

Through the collected data with various square sizes, we can obtain the following approximated SIF with dependence on A such that:where the values of and are also listed in Table 1 for comparison with and , respectively.

Since the problem considers an embedded crack, the values of strain and stress may diverge near the crack tip. Here, we check the relative errors of the displacement, strain, and stress with various w/A ratios, as depicted in Fig. 14. The comparison shows that the approximate CAMS-FE results can converge to the classic LFEM results when the number of terms is sufficiently large.

Consider the square plate with an ellipse hole depicted in Fig. 15. The length along the x and y directions is L = 1, and the thickness is h = 0.2. The semi-major and semi-minor axes of the ellipse are a = 0.3 and b = 0.1, respectively. A rotation angle θ is introduced, as depicted in Fig. 15. We use the CMAS-FE to derive an approximated bending stiffness as a function of θ only.

To collect sufficient data, seven CMAS-FE simulations are tested with θ = {0, π/12, π/6, …, π/2}. Four selected hexahedron finite element meshes are depicted in Fig. 16. Here we only study the bending stiffness as function of θ, and the interpolation yields:In addition, we demonstrate the stress error distributions for various values of θ, as depicted in Fig. 17 with E = 2.1 × 10⁵ MPa and ν = 0.23. We can see that the maximum error occurs at the edges, where stress concentration or a high stress gradient occurs.