A Comprehensive $p$ -VEM Framework for Advanced
Variable Stiffness Plates with Arbitrary Shapes

Paola Pia Foligno Affiliation: Dipartimento di Scienze e Tecnologie Aerospaziali, Politecnico di Milano, 20156 Milano, Italy Daniele Boffi Affiliation: CEMSE Division, King Abdullah University of Science and Technology, 23955 Thuwal, Saudi Arabia Affiliation: Dipartimento di Matematica “F. Casorati”, Università degli Studi di Pavia, 27100 Pavia, Italy Affiliation: IMATI “E. Magenes”, CNR, 27100 Pavia, Italy Fabio Credali Note: Corresponding author.
Email addresses:paolapia.foligno@polimi.it (Paola Pia Foligno), daniele.boffi@kaust.edu.sa (daniele Boffi), fabio.credali@kaust.edu.sa (fabio Credali), riccardo.vescovini@polimi.it (Riccardo Vescovini) Affiliation: CEMSE Division, King Abdullah University of Science and Technology, 23955 Thuwal, Saudi Arabia Riccardo Vescovini Affiliation: Dipartimento di Scienze e Tecnologie Aerospaziali, Politecnico di Milano, 20156 Milano, Italy

Abstract

This paper presents a comprehensive, high-order ( $p$ -version) Virtual Element Method (VEM) framework for the structural analysis of innovative variable stiffness plates. VEM is particularly suited for complex configurations due to its ability to handle arbitrary polygonal meshes, including curved edges. However, its mathematical formulation may hinder its spread in the engineering community. This work illustrates a formulation with an accessible implementation using well-known FEM notation and integrating at the same time a set of new advanced capabilities. Specifically, both standard stabilized and advanced self-stabilized strategies are adopted. To further improve the robustness of VEM in the presence of variable coefficients, polynomial projections taking into account the coefficients are employed. This approach is referred to as Variable Coefficients-VEM approach ( $\mathrm{VC}$ -VEM). This unified framework is applied to linear static, free-vibration, and buckling analyses, and validated against analytical solutions and numerical benchmarks. In particular, plates with cutouts and problems featuring high-gradient solutions are investigated, demonstrating that the proposed comprehensive approach provides a flexible and ready-to-implement tool for advanced structural design.

1 Introduction

Innovative structural configurations feature plates characterized by a variable stiffness skin, where the fibers follow curvilinear paths. These variable stiffness (VS) plates offer significant potential for weight minimization and improved efficiency compared to classical designs. Thus, increasing attention has been devoted to their study.

The literature on these innovative configurations is largely limited to the use of classical numerical methods, namely the Finite Element Method and the Ritz method. Recent applications in finite element frameworks include [41], where shear buckling and postbuckling responses were compared to classical straight fiber configurations, and [33], which addressed thermal buckling. On the other hand, applications of the Ritz method can be found for the mechanical [53, 18] and thermal [47] buckling, where the use of simple geometries, such as rectangular domains, is due to the global nature of the trial functions. Interesting exceptions where the Ritz method is constructed on more complex geometries include [26, 27], in which plates with arbitrary cutouts have been investigated, and [28], which accounted for arbitrary geometries. Innovative Ritz-based formulations for the linear and geometrically nonlinear responses of arbitrary domains have been proposed in [50, 48].

It is therefore clear that the study of VS plates requires advanced and effective numerical methods capable of tackling the challenges arising in this field. In this regard, the Virtual Element Method (VEM) represents a promising alternative, as it employs general polygonal elements (even with curved edges [12]), thereby significantly simplifying the mesh generation for complex geometries.

Pioneering works [5, 6] introduced the innovative formulation of the method and its implementation for the Laplace problem, while its accuracy for general second-order elliptic problems has been assessed in [7].

The mathematical foundations of VEM have been extensively investigated over the last decade. However, most of the existing literature is formulated within a rigorous mathematical framework that may hinder its spread in the engineering community, this consideration being exacerbated by a relatively steep learning curve. Investigations of elasticity problems within the VEM framework can be found in [8, 9, 2, 35, 21], while Kirchhoff-Love and Reissner-Mindlin plate models have been addressed in [17, 37] and [11, 22], respectively. Applications to buckling problems can be found in [38]. A key limitation of the VEM is the need for a stabilization term to handle the non-polynomial residual component of the trial functions and ensure the well-posedness of the discrete problem. The stability of the method has been investigated in [10, 34], where several stabilization formulas have also been introduced. Generally, in the standard VEM, the stabilization term does not conform with the physics of the problem under consideration and it is scaled by a user-defined tuning parameter. This may reduce the generality of the method and lead to over-stabilization, resulting in degraded solution quality. This is particularly relevant for anisotropic problems [14] and eigenvalue problems [16]. Alternative approaches have therefore been proposed to avoid arbitrary stabilization terms, such as strategies based on higher-order polynomial projections [15], divergence-free polynomial spaces [13] and enriched VEM spaces with additional internal degrees of freedom (DOFs) [29]. The effect of the stabilization parameter has been studied in [25], showing that choosing it based on the bending strain energy density yields more accurate results compared to classical approaches. Self-stabilized formulations for mixed linear elasticity have been proposed in [31].

The treatment of variable coefficients in the bilinear form is another relevant aspect of the VEM that requires further investigation, particularly in advanced applications. In [7], the standard $L^{2}$ projection was adopted and the variable coefficients were introduced only in the construction of the final discrete bilinear form. This approach was shown to perform better in the presence of variable coefficients compared to other alternatives. A VEM application to composites with spatially varying fiber directions was investigated in [42]. The authors demonstrated that approximating the fiber direction as a constant equal to the average at the nodal values yields more accurate results than using the element’s centroid value. Spatially varying material properties have been studied within a Hellinger–Reissner VEM framework in 2D [3] and in 3D [32]. A recent work [23] introduces a novel formulation, denoted as $\mathrm{VC}$ -VEM, in which the coefficient is fully incorporated into the definition of the projection operators, both in stabilized and self-stabilized settings.

This paper aims to bridge the gap between theoretical VEM formulations and engineering practice by presenting a robust, high-order $p$ -VEM computational framework with focus on the analysis of variable stiffness plates with complex geometries. For this purpose, the main objective of this work is to provide a self-contained reference, illustrating VEM formulation with an implementation-oriented matrix notation, familiar to the engineering computational mechanics community. In addition to the above mentioned objectives, the proposed advanced formulation exploits the geometric flexibility of VEM to consider polygonal elements with curved edges and hanging nodes, hence simplifying the mesh generation for plates with complex boundaries and cutouts. To address the critical choice of the stabilization term, both stabilized and self-stabilized strategies are included and discussed. Furthermore, to enhance numerical robustness in the presence of non-uniform elasticity properties, a novel Variable Coefficient-VEM ( $\mathrm{VC}$ -VEM) approach is adopted to integrate spatial variability within the projection operators.

The manuscript is organized as follows. Section 2 presents the structural modeling of the variable stiffness plate. In Section 3, the VEM formulation is introduced, including the stabilized and self-stabilized strategies, as well as the $\mathrm{VC}$ -VEM approach. Section 4 provides the numerical results, in which the developed formulation is validated against analytical solutions and numerical results from the literature. Lastly, the conclusions are drawn in Section 5.

2 Structural modeling

The structures under investigation are variable stiffness (VS) plates, where the fiber orientation varies across the domain, rendering the stiffness a function of the planar position. A two-dimensional model is employed, with a Cartesian reference system $xyz$ , whose origin is located on the plate midsurface. The $x$ and $y$ axes are directed in the longitudinal and transverse directions, respectively, and the $z$ axis is in the thickness direction. The geometry is arbitrary, with maximum dimensions equal to $a$ and $b$ and thickness $t$ . The plate domain is denoted by $\Omega$ and its Lipschitz boundary $\Gamma$ consists of a finite number of smooth curves $\left\{\Gamma_{i}\right\}_{i=1,\dots,N_{e}}$ , with $N_{e}$ being the number of edges. Each curve $\Gamma_{i}$ is of class $\mathcal{C}^{m+1}$ for $m\geq 0$ and it is parametrized by an invertible $\mathcal{C}^{m+1}$ map $\gamma_{i}\mathrel{\mathop{\mathchar 58\relax}}I_{i}=\left[a_{i},b_{i}\right]\rightarrow\Gamma_{i}$ [12].

2.1 Variational statement

The formulation is developed within a displacement-based variational framework.

Linear static, buckling, and free-vibration analyses are of concern. Hence, the variational statement can be expressed in a unified form as [30]:

\left(\beta_{1}+\beta_{2}+\beta_{3}\right)\delta W_{i}+\beta_{2}\delta W_{b}+\beta_{3}\delta W_{k}=\beta_{1}\delta W_{e},

(1)

where $W_{\alpha}$ , with $\alpha=i,b,k,e$ , are used to denote the internal virtual work, the pre-buckling energy contribution, the first variation of the kinetic energy contribution and the external virtual work, respectively. The Boolean flags $\beta_{i}$ , with $i=1,2,3$ , are chosen dependently on the analysis of interest, following the summary reported in Table 1.

Analysis type	$\beta_{1}$	$\beta_{2}$	$\beta_{3}$
Linear static	1	0	0
Buckling	0	1	0
Free-vibration	0	0	1

Table 1: Variational statement coefficients for the different analysis types.

2.2 Plate model

Hereafter, the two-dimensional plate model is presented. First, the kinematics, the strain measure, and the constitutive law, which accounts for the spatially varying fiber orientations, are detailed. Lastly, the energy terms are derived within the variational framework.

2.2.1 Kinematics

The plate kinematics is modeled according to the First-order Shear Deformation Theory (FSDT), which allows thin and relatively thick panels to be considered. The displacement of a generic point on the plate is expressed as [43]:

	$\displaystyle\bm{d}\left(x,y,z\right)=\begin{Bmatrix}d_{x}\left(x,y,z\right)\\ d_{y}\left(x,y,z\right)\\ d_{z}\left(x,y,z\right)\end{Bmatrix}$	$\displaystyle=\begin{Bmatrix}u\left(x,y\right)\\ v\left(x,y\right)\\ w\left(x,y\right)\end{Bmatrix}+z\begin{bmatrix}1&0\\ 0&1\\ 0&0\\ \end{bmatrix}\begin{Bmatrix}\theta_{x}\left(x,y\right)\\ \theta_{y}\left(x,y\right)\\ 0\end{Bmatrix}=\bm{u}^{0}\left(x,y\right)+z\bm{L}^{\bm{\theta}}\,\bm{\theta}\left(x,y\right)$		(2)
		$\displaystyle=\begin{bmatrix}\bm{I}&z\,\bm{L}^{\bm{\theta}}\end{bmatrix}\begin{Bmatrix}\bm{u}^{0}\left(x,y\right)\\ \bm{\theta}\left(x,y\right)\end{Bmatrix}=\begin{bmatrix}\bm{I}&z\,\bm{L}^{\bm{\theta}}\end{bmatrix}\bm{u},$		(2)

where $\bm{u}^{0}$ and $\bm{\theta}$ are the generalized displacements and rotations of the midsurface, and $\bm{u}$ is the vector collecting them. The conventions of the plate kinematics are shown in Figure 1(a).

The strains $\hat{\bm{\varepsilon}}=\begin{Bmatrix}\varepsilon_{xx},\;\varepsilon_{yy},\;\gamma_{xy}\end{Bmatrix}^{T}$ and $\bm{\gamma}=\begin{Bmatrix}\gamma_{yz},\;\gamma_{xy}\end{Bmatrix}^{T}$ are expressed in terms of the displacements as:

\displaystyle\hat{\bm{\varepsilon}}=\bm{\varepsilon}^{0}\left(\bm{u}\right)+z\,\bm{k}\left(\bm{u}\right),\qquad\bm{\gamma}=\bm{\gamma}\left(\bm{u}\right),

(3)

where $\bm{\varepsilon}^{0}$ denotes the membrane strains, $\bm{k}$ the curvatures, and $\bm{\gamma}$ the transverse shear strains, and their expression reads:

\begin{split}\bm{\varepsilon}^{0}\left(\bm{u}\right)=\begin{Bmatrix}u_{,x}\\ v_{,y}\\ u_{,y}+v_{,x}\end{Bmatrix},\qquad\bm{k}\left(\bm{u}\right)=\begin{Bmatrix}\theta_{x,x}\\ \theta_{y,y}\\ \theta_{x,y}+\theta_{y,x}\end{Bmatrix},\qquad\bm{\gamma}\left(\bm{u}\right)=\begin{Bmatrix}\theta_{y}+w_{,y}\\ \theta_{x}+w_{,x}\\ \end{Bmatrix},\end{split}

(4)

where $(\cdot)_{,x}$ and $(\cdot)_{,y}$ are the derivatives with respect to the in-plane coordinates.

In view of future developments, the generalized strains can be conveniently collected into a single vector:

\bm{\varepsilon}\left(\bm{u}\right)=\begin{Bmatrix}\bm{\varepsilon}^{0}\left(\bm{u}\right)\\ \bm{k}\left(\bm{u}\right)\\ \bm{\gamma}\left(\bm{u}\right)\end{Bmatrix}.

(5)

2.2.2 Constitutive law

In variable stiffness plates, the fiber orientation angles are a function of the position. Different representations have been proposed in the literature. Linear variations [40], Lobatto distributions [1] and NURBS [39] are examples. In this work, the angles $\theta_{mn}$ are specified on a grid of M $\times$ N points over one quarter of the plate domain, and the angles at a generic point of the domain are retrieved via Lagrange polynomials interpolation [53, 51, 56, 49]:

\theta\left(x,y\right)=\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}\theta_{mn}\prod_{m\neq i}\left(\frac{|x|-x_{i}}{x_{m}-x_{i}}\right)\prod_{n\neq j}\left(\frac{|y|-y_{j}}{y_{n}-y_{j}}\right).

(6)

The proposed distribution is illustrative of a specific choice, but other choices could be easily accommodated within the present framework. From Eq. (6), the angle $\theta\left(x,y\right)$ is nonlinear with respect to the coordinates $x$ and $y$ . The linear variation is retrieved as a special case by considering only two points, e.g. the center and the edge of the plate. The fiber orientation distribution is illustrated in Figure 1(b).

Refer to caption — (a) Conventions for generalized displacements.

For each ply $k$ , the constitutive law is expressed in global coordinates as:

\displaystyle\bm{\sigma}_{k}=\bar{\bm{Q}}_{k}\left(x,y\right)\,\left(\hat{\bm{\varepsilon}}_{k}-\bm{\varepsilon}_{k}^{\mathrm{th}}\right),

\displaystyle\bm{\tau}_{k}=\bar{\bm{Q}}_{nk}\left(x,y\right)\,\bm{\gamma},

(7)

where $\bar{\bm{Q}}_{k}\left(x,y\right)$ and $\bar{\bm{Q}}_{nk}\left(x,y\right)$ are the constitutive matrices expressed in laminate axes. The vector $\hat{\bm{\varepsilon}}_{k}$ is the total deformation, and $\bm{\varepsilon}_{k}^{\mathrm{th}}$ is the thermal contribution, accounted only for the in-plane behavior, as out-of-plane thermal stresses are typically negligible for the plates under consideration [49]. The thermal contribution $\bm{\varepsilon}_{k}^{\mathrm{th}}$ is expressed as in [49]:

\bm{\varepsilon}_{k}^{\mathrm{th}}=\bar{\bm{\alpha}}_{k}\left(x,y\right)\Delta T,

(8)

where $\bar{\bm{\alpha}}_{k}$ is the ply thermal expansion coefficient in global coordinates and $\Delta T$ is the temperature gradient. The coefficients of $\bar{\bm{\alpha}}_{k}$ are assumed to be temperature-independent.

The thermo-elastic constitutive law is then obtained as:

\bm{S}=\begin{Bmatrix}\bm{N}\\ \bm{M}\\ \bm{Q}\end{Bmatrix}=\begin{bmatrix}\bm{A}\left(x,y\right)&\bm{B}\left(x,y\right)&\bm{0}\\ \bm{B}\left(x,y\right)&\bm{D}\left(x,y\right)&\bm{0}\\ \bm{0}&\bm{0}&\bm{A}^{\mathrm{s}}\left(x,y\right)\\ \end{bmatrix}\bm{\varepsilon}-\begin{Bmatrix}\hat{\bm{N}}\\ \hat{\bm{M}}\\ \bm{0}\end{Bmatrix}\Delta T=\mathbb{C}\left(x,y\right)\bm{\varepsilon}-\hat{\bm{R}}\Delta T,

(9)

where $\bm{N}$ , $\bm{M}$ and $\bm{Q}$ are the forces and moments per unit length and $\bm{S}$ is the vector collecting them. The resulting matrices $\bm{A}$ , $\bm{B}$ , $\bm{D}$ and $\bm{A}^{\mathrm{s}}$ are the laminate stiffness matrices, and $\hat{\bm{N}}$ and $\hat{\bm{M}}$ are the thermal forces and moments per unit length, defined as [49]:

\displaystyle\hat{\bm{N}}\left(x,y\right)=\sum_{k=1}^{N}\int_{t_{k}}^{t_{k+1}}\bar{\bm{Q}}_{k}\,\bar{\bm{\alpha}}_{k}\left(x,y\right)\,\mathrm{d}z,

\displaystyle\hat{\bm{M}}\left(x,y\right)=\sum_{k=1}^{N}\int_{t_{k}}^{t_{k+1}}z\,\bar{\bm{Q}}_{k}\,\bar{\bm{\alpha}}_{k}\left(x,y\right)\,\mathrm{d}z.

(10)

2.2.3 Energy terms

Having defined the kinematics, the strain measure, and the constitutive relation, it is now convenient to derive the energy terms in the Principle of Virtual Work framework. In particular, the internal virtual work expression reads:

\delta W_{i}=\int_{\Omega}\bm{\varepsilon}\left(\delta\bm{u}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{u}\right)\,\mathrm{d}\Omega.

(11)

The buckling contribution is expressed as:

\delta W_{b}=\int_{\Omega}\bm{\varepsilon}_{\mathrm{b}}\left(\delta\bm{u}\right)^{T}\begin{bmatrix}N_{xx}&N_{xy}\\ N_{xy}&N_{yy}\end{bmatrix}\bm{\varepsilon}_{\mathrm{b}}\left(\bm{u}\right)\,\mathrm{d}\Omega=\int_{\Omega}\bm{\varepsilon}_{\mathrm{b}}\left(\delta\bm{u}\right)^{T}\mathbb{N}\left(x,y\right)\bm{\varepsilon}_{\mathrm{b}}\left(\bm{u}\right)\,\mathrm{d}\Omega,

(12)

with:

\bm{\varepsilon}_{\mathrm{b}}\left(\cdot\right)=\begin{bmatrix}0&0&,x&0&0\\ 0&0&,y&0&0\end{bmatrix}.

(13)

For the free-vibration analysis, the contribution due to inertial forces is:

\displaystyle\delta W_{k}

\displaystyle=\int_{\Omega}\delta\bm{u}^{T}\int_{t}\rho_{0}\begin{bmatrix}\bm{I}&z\bm{L}\\ z\bm{L}&z^{2}\bm{L}\end{bmatrix}\,\mathrm{d}z\bm{\ddot{u}}\,\mathrm{d}\Omega=\int_{\Omega}\delta\bm{u}^{T}\begin{bmatrix}\bm{I}_{0}&\bm{I}_{1}\\ \bm{I}_{1}&\bm{I}_{2}\end{bmatrix}\bm{\ddot{u}}\,\mathrm{d}\Omega=\int_{\Omega}\delta\bm{u}^{T}\mathbb{M}\bm{\ddot{u}}\,\mathrm{d}\Omega,

(14)

where $\bm{\ddot{u}}$ denotes the second time derivative of $\bm{u}$ , $\bm{I}_{0},\bm{I}_{1},\bm{I}_{2}$ are the inertial moments and $\mathbb{M}$ is the mass matrix of the plate.

The external virtual work accounts for body forces, boundary traction forces, and thermal loads, in the form:

\delta W_{e}=\int_{\Omega}\delta\bm{u}^{T}\bar{\bm{b}}\,\mathrm{d}\Omega+\sum_{j}\int_{\Gamma_{j}}\delta\bm{u}^{T}\bar{\bm{t}}_{j}\,\mathrm{d}\Gamma_{j}+\int_{\Omega}\bm{\varepsilon}\left(\delta\bm{u}\right)^{T}\hat{\bm{R}}\,\mathrm{d}\Omega,

(15)

where $\bar{\bm{b}}$ and $\bar{\bm{t}}_{j}$ are the prescribed body and the traction forces on the generic edge $\Gamma_{j}$ , respectively. Prescribed displacements are enforced as in standard finite element strategies.

3 Virtual Element Method

The numerical approximation is carried out via the Virtual Element Method (VEM). The method allows the use of arbitrary polygonal elements, which is a desirable feature to simplify the mesh generation on complex domains. Moreover, the VEM naturally handles hanging nodes, allowing the treatment of nonconforming discretizations as conforming ones.

By following the approach in [12], the standard VEM framework is extended here to account for curved edges, effectively enabling the treatment of arbitrary curved boundaries. Moreover, the $p$ -version of VEM is considered, where $p\geq 1$ is the order of accuracy. The combined effect of curved edges representation capability and higher-order approximations is particularly effective to achieve simple yet accurate models, eliminating the need to operate geometry-induced mesh refinement strategies.

In the following, the VEM discretization, space, and associated degrees of freedom (DOFs), are recalled. Then, the discrete forms and corresponding projector operators of the energy contributions are presented, for both standard and $\mathrm{VC}$ -VEM [23]. The reader is referred to the Supporting material Section for more details on the steps to build the relevant matrices, as well as for the definitions of the different operators and quantities not otherwise specified.

3.1 Discretization

The domain $\Omega$ is decomposed into a finite set $\mathcal{T}_{h}$ of non-overlapping star-shaped polygons with possibly curved edges $E\in\mathcal{T}_{h}$ [5], with $h=\max_{E\in\mathcal{T}_{h}}h_{E}$ and $h_{E}$ the diameter of $E$ . The boundary of $E$ is denoted by $\Gamma^{E}=\partial E$ , which can be partitioned into straight segments $\Gamma^{E}_{e}$ with $e=1,\dots,n_{e}$ and curved edges $\tilde{\Gamma}^{E}_{e}$ with $e=1,\dots,\tilde{n}_{e}$ , the latter satisfying the regularity assumptions introduced earlier.

The generic continuous bilinear form $a\left(\cdot,\cdot\right)$ is now introduced. It is obtained as the sum of the contributions of the elements $E$ in $\mathcal{T}_{h}$ as:

a\left(\delta\bm{u},\bm{u}\right)=\sum_{E\in\mathcal{T}_{h}}a^{E}\left(\delta\bm{u},\bm{u}\right)\quad\forall\delta\bm{u},\bm{u}\in\bm{\mathcal{V}},

(16)

where $\bm{\mathcal{V}}\equiv\left[\mathcal{V}^{u}\times\mathcal{V}^{v}\times\mathcal{V}^{w}\times\mathcal{V}^{\theta_{x}}\times\mathcal{V}^{\theta_{y}}\right]$ is the continuous space in which the generalized displacement components lie. The discretization uses the space $\bm{\mathcal{V}}_{h}\equiv\left[\mathcal{V}_{h}^{u}\times\mathcal{V}_{h}^{v}\times\mathcal{V}_{h}^{w}\times\mathcal{V}_{h}^{\theta_{x}}\times\mathcal{V}_{h}^{\theta_{y}}\right]$ , with $\bm{\mathcal{V}}_{h}\subset\bm{\mathcal{V}}$ . Consequently, the discrete version of Eq. (16) translates into:

a_{h}\left(\delta\bm{u}_{h},\bm{u}_{h}\right)=\sum_{E\in\mathcal{T}_{h}}a_{h}^{E}\left(\delta\bm{u}_{h},\bm{u}_{h}\right)\quad\forall\delta\bm{u}_{h},\bm{u}_{h}\in\bm{\mathcal{V}}_{h},

(17)

where $\bm{u}_{h}$ is the discrete counterpart of $\bm{u}$ .

3.2 Virtual Element space and degrees of freedom

The local virtual element spaces, on a generic polygon $E$ , with order of accuracy $p\geq 1$ and $m\geq p$ , are hereafter defined. For a generic $\phi$ representing a generalized displacement, the space is defined as:

$\displaystyle\mathcal{V}_{h}^{r}\left(E\right)=\left\{\right.$	$\displaystyle{\displaystyle\phi}_{h}\in H^{1}\left(E\right)\cap C^{0}\left(E\right)\mathrel{\mathop{\mathchar 58\relax}}\bm{L}\left[\mathbb{C}\bm{\varepsilon}\left({\phi}_{h}\right)\right]\big\|_{E}\in\left[\mathcal{P}_{p}\left(E\right)\right]^{5},$	(18)
	$\displaystyle\left.{\phi}_{h}\big\|_{\Gamma^{E}_{e}}\in\mathcal{P}_{p}\left(\Gamma^{E}_{e}\right)\quad\forall e=1,\dots,n_{e},\;{\phi}_{h}\big\|_{\tilde{\Gamma}^{E}_{e}}\in\mathcal{\tilde{P}}_{p}\left(\tilde{\Gamma}^{E}_{e}\right)\quad\forall e=1,\dots,\tilde{n}_{e},\right.$
	$\displaystyle\left.\int_{E}{\phi}_{h}\;q\;\,\mathrm{d}E=\int_{E}{\Pi^{\mathcal{K}}_{p}}^{\phi}{\phi}_{h}\;q\;\,\mathrm{d}E\quad\forall q\in\mathcal{P}_{p}\left(E\right)\setminus\mathcal{P}_{p-r}\left(E\right)\right\},$

where $\bm{\varepsilon}\left({\phi}_{h}\right)$ is the strain operator obtained by considering only the corresponding displacement component. $\mathcal{P}_{p}\left(E\right)$ is the polynomial space of degree less than or equal to $p$ and $\mathcal{\tilde{P}}_{p}\left(\tilde{\Gamma}^{E}_{e}\right)$ is the polynomial space on a curved edge, which, following [12], can be expressed as:

\mathcal{\tilde{P}}_{p}\left(\tilde{\Gamma}^{E}_{e}\right)=\left\{\tilde{q}=q\circ\gamma_{e}^{-1}\mathrel{\mathop{\mathchar 58\relax}}q\in\mathcal{P}_{p}\left(I_{e}\right)\right\}.

(19)

Hence, the functions on the curved edges are polynomials with respect to the parametrization $\gamma_{e}$ . The operator ${\Pi^{\mathcal{K}}_{p}}^{\phi}$ is the scalar counterpart, referred to the generic generalized displacement $\phi$ , of the operator $\bm{\Pi}^{\mathcal{K}}_{p}$ that projects the trial functions from the VEM to the polynomial space, which will be defined later. The spaces for the different displacement components are then obtained as:

\mathcal{V}_{h}^{u}\left(E\right)=\mathcal{V}_{h}^{v}\left(E\right)=\mathcal{V}_{h}^{2}\left(E\right),\qquad\mathcal{V}_{h}^{w}\left(E\right)=\mathcal{V}_{h}^{1}\left(E\right),\qquad\mathcal{V}_{h}^{\theta_{x}}\left(E\right)=\mathcal{V}_{h}^{\theta_{y}}\left(E\right)=\mathcal{V}_{h}^{0}\left(E\right).

(20)

The total local virtual element space for a generic element $E$ of dimension $N_{d}\big|_{E}$ reads:

\bm{\mathcal{V}}_{h}\left(E\right)=\left[\mathcal{V}_{h}^{u}\times\mathcal{V}_{h}^{v}\times\mathcal{V}_{h}^{w}\times\mathcal{V}_{h}^{\theta_{x}}\times\mathcal{V}_{h}^{\theta_{y}}\right].

(21)

The unknown $\bm{u}_{h}\in\bm{\mathcal{V}}_{h}\left(E\right)$ is uniquely identified by the following set of degrees of freedom ( $\textbf{DOF}_{i}$ ):

•

$\textbf{DOF}_{1}$ : the values of $\bm{u}_{h}$ at the $n_{v}$ vertices of $E$ ,
•

$\textbf{DOF}_{2}$ : for $p>1$ , the values of $\bm{u}_{h}$ at the $p-1$ internal Gauss-Lobatto quadrature points on each straight edge $\Gamma^{E}_{e}$ , and the values of $\bm{u}_{h}$ at the $p-1$ points on each curved edge $\tilde{\Gamma}^{E}_{e}$ , which are images through $\gamma_{e}$ of the $p-1$ internal Gauss-Lobatto quadrature points on $I_{e}$ [12],

•

$\textbf{DOF}_{3,4}$ : for $p>1$ , the internal moments of $\mathcal{V}_{h}^{u}\left(E\right)$ and $\mathcal{V}_{h}^{v}\left(E\right)$ up to order $p-2$ :

\frac{1}{\left|E\right|}\int_{E}q\left(x,y\right)u_{h}\left(x,y\right)\,\mathrm{d}E\quad\forall q\in\mathcal{Q}_{p-2}\left(E\right),\quad\frac{1}{\left|E\right|}\int_{E}q\left(x,y\right)v_{h}\left(x,y\right)\,\mathrm{d}E\quad\forall q\in\mathcal{Q}_{p-2}\left(E\right).

(22)

•

$\textbf{DOF}_{5}$ : the internal moments of $\mathcal{V}_{h}^{w}\left(E\right)$ up to order $p-1$ :

\frac{1}{\left|E\right|}\int_{E}q\left(x,y\right)w_{h}\left(x,y\right)\,\mathrm{d}E\quad\forall q\in\mathcal{Q}_{p-1}\left(E\right).

(23)

•

$\textbf{DOF}_{6,7}$ : the internal moments of $\mathcal{V}_{h}^{\theta_{x}}\left(E\right)$ and $\mathcal{V}_{h}^{\theta_{y}}\left(E\right)$ up to order $p$ :

\frac{1}{\left|E\right|}\int_{E}q\left(x,y\right){\theta_{x}}_{h}\left(x,y\right)\,\mathrm{d}E\quad\forall q\in\mathcal{Q}_{p}\left(E\right),\quad\frac{1}{\left|E\right|}\int_{E}q\left(x,y\right){\theta_{y}}_{h}\left(x,y\right)\,\mathrm{d}E\quad\forall q\in\mathcal{Q}_{p}\left(E\right),

(24)

where $\left|E\right|$ is the area of the element $E$ and $\mathcal{Q}_{p}\left(E\right)$ is a basis for the polynomial space $\mathcal{P}_{p}\left(E\right)$ . The displacements are polynomials of order $p$ on the edges. Conversely, in the element interior they differ for the FSDT formulation [22], and are associated to polynomials of degree $p-2$ , $p-1$ and $p$ for the in-plane displacements, out-of-plane displacement and in-plane rotations, respectively.

A schematic representation of the degrees of freedom (DOFs) for $p=2$ is shown in Figure 2, where circles and squares represent the vertex and edge DOFs, respectively, while triangles correspond to the internal moments.

The total global virtual element space is then obtained as:

\bm{\mathcal{V}}_{h}=\left\{\bm{u}\in\left[H_{0}^{1}\left(\Omega\right)\right]^{5}\mathrel{\mathop{\mathchar 58\relax}}\bm{u}\in\bm{\mathcal{V}}_{h}\left(E\right)\quad\forall E\in\mathcal{T}_{h}\right\}=\left[\mathcal{V}_{h}^{u}\times\mathcal{V}_{h}^{v}\times\mathcal{V}_{h}^{w}\times\mathcal{V}_{h}^{\theta_{x}}\times\mathcal{V}_{h}^{\theta_{y}}\right].

(25)

The dependence of the constitutive law $\mathbb{C}$ on the spatial coordinates $x,y$ is here neglected. Otherwise, the definition of the discrete spaces is inevitably more involved. The extension to this case is addressed at the end of this section.

3.3 Discrete bilinear forms and projection operators

As usual with VEM, the trial functions are solutions of a PDE inside each element and they are not explicitly computed. Consequently, only a projection from the VEM onto the polynomial space of order $p$ is computable through the degrees of freedom:

\bm{\Pi}_{E}\mathrel{\mathop{\mathchar 58\relax}}\bm{\mathcal{V}}_{h}\left(E\right)\rightarrow\bm{\mathcal{P}}_{p}\left(E\right).

(26)

For the generic bilinear form, the projection is the solution of:

a^{E}\left(\bm{u}_{h},\bm{q}\right)=a^{E}\left(\bm{\Pi}_{E}\bm{u}_{h},\bm{q}\right)\quad\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right).

(27)

The displacement field can be written using the Lagrangian-type interpolation identity as:

\bm{u}_{h}=\sum_{i=1}^{N_{d}\big|_{E}}\text{dof}_{i}\left(\bm{u}_{h}\right)\bm{\psi}_{i}\quad\forall\bm{u}_{h}\in\bm{\mathcal{V}}_{h}\left(E\right),

(28)

with $\text{dof}_{j}\left(\bm{\psi}_{i}\right)=\bm{\delta}_{ij}$ for all $i,j=1,\dots,N_{d}\big|_{E}$ , where $\bm{\delta}_{ij}$ is the $5\times 1$ Kronecker- $\delta$ vector having unit value on the corresponding displacement component when $i=j$ and zeros elsewhere. Substitution of Eq. (28) into Eq. (27) yields:

a^{E}\left(\bm{\psi}_{i},\bm{q}\right)=a^{E}\left(\bm{\Pi}_{E}\bm{\psi}_{i},\bm{q}\right)\quad\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right).

(29)

In the following, the projections are defined for the stiffness matrix, the geometric stiffness matrix, the mass matrix, and the body and thermal force vectors.

3.3.1 Stiffness matrix

The evaluation of the stiffness matrix is first presented by assuming that the constitutive tensor is constant within each element, with a value equal to its average at the integration points. This approximation facilitates the computation of the projection using the degrees of freedom. The generalization to spatially varying constitutive tensors $\mathbb{C}\left(x,y\right)$ is outlined at the end of the section.

Accordingly, the bilinear form associated with the stiffness matrix is of elliptic type and is defined as:

\displaystyle\mathcal{K}^{E}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)=\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{\psi}_{j}\right)\,\mathrm{d}E,

(30)

and the projection operator $\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i}$ is the solution of:

\begin{cases}\begin{aligned} &\mathcal{K}^{E}\left(\bm{\psi}_{i},\bm{q}\right)=\mathcal{K}^{E}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)\\ &P_{0}\left(\bm{\psi}_{i},\bm{q}\right)=P_{0}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{0}\left(E\right)\oplus\hat{\bm{P}}\left(E\right),\end{aligned}\end{cases}

(31)

for all $i=1,\dots,N_{d}\big|_{E}$ ; $\hat{\bm{P}}\left(E\right)$ corresponds to $\begin{Bmatrix}x\;0\;0\;0\;0\end{Bmatrix}^{T}$ and $P_{0}\left(\bm{\psi}_{i},\bm{q}\right)$ , introduced to take care of the rigid body motions, is defined as:

\displaystyle P_{0}\left(\bm{\psi}_{i},\bm{q}\right)=\frac{1}{n_{v}}\sum_{j=1}^{5n_{v}}\text{dof}_{j}\left(\bm{\psi}_{i}\right)\text{dof}_{j}\left(\bm{q}\right).

(32)

The evaluation of the right-hand side of Eq. (31) is straightforward as it involves the product of known polynomials. Regarding the left-hand side, integration by parts is applied. Subsequently, line integrals are easily computed as the trial functions are known on the boundary. The surface integrals are evaluated from the internal degrees of freedom.

Once the projection operator is available, following [6], a generic virtual function can be written as:

\bm{\psi}_{i}=\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i}+\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{i}.

(33)

Substituting Eq. (33) into the stiffness bilinear form, and omitting the mixed terms, yields:

{\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)={\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i},\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{j}\right)+{\mathcal{K}_{\mathrm{VC}}}^{E}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{i},\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{j}\right),

(34)

where:

\displaystyle{\displaystyle\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)=\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{\psi}_{j}\right)\,\mathrm{d}E.

(35)

In Eq. (34), if $\mathbb{C}$ is constant within the elements, the two mixed terms vanish by the definition of $\bm{\Pi}^{\mathcal{K}}_{p}$ ; otherwise, they are omitted.

The first term on the right-hand side of Eq. (34) is the consistency term that ensures the accuracy of the solution. The second term is a contribution required to guarantee the stability of the solution; however, since the trial functions are not explicitly known inside the element, it cannot be computed explicitly.

Stabilized VEM

Following [5], the stabilization term is a symmetric and semi-positive definite bilinear form defined in such a way that there exist two positive constants $\alpha_{*}\left(p\right)$ and $\alpha^{*}\left(p\right)$ independent of $h$ :

\alpha_{*}a^{E}\left(\delta\bm{u}_{h},\delta\bm{u}_{h}\right)\leq S^{E}\left(\delta\bm{u}_{h},\delta\bm{u}_{h}\right)\leq\alpha^{*}a^{E}\left(\delta\bm{u}_{h},\delta\bm{u}_{h}\right)\quad\forall\delta\bm{u}_{h}\in\bm{\mathcal{V}}_{h}\left(E\right)\quad\text{with}\quad\bm{\Pi}^{\mathcal{K}}_{p}\delta\bm{u}_{h}=0.

(36)

Consequently, the discrete version of Eq. (34) is defined as:

{\mathcal{K}_{\mathrm{VC}}}^{E}_{h}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)={\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i},\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{j}\right)+S^{E}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{i},\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{j}\right).

(37)

In [23], several stabilization recipes from the literature have been compared for different PDEs and approximation orders to identify the most robust one. Relying upon the outcomes of the mentioned study, the stabilization term is selected as in [34]. In particular, the stabilization is defined block-wise according to the different components:

\left\{\begin{aligned} &S^{E}=\tau\sum_{r\in\mathcal{I}^{u,v}}\text{dof}_{r}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{i}\right)\cdot\left(s^{E}\right)_{r}\cdot\text{dof}_{r}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{j}\right),&&\text{if}\quad i,j\in\mathcal{I}^{u,v}\\ &S^{E}=\tau\sum_{r\in\mathcal{I}^{w}}\text{dof}_{r}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{i}\right)\cdot\left(s^{E}\right)_{r}\cdot\text{dof}_{r}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{j}\right),&&\text{if}\quad i,j\in\mathcal{I}^{w}\\ &S^{E}=\tau\sum_{r\in\mathcal{I}^{\theta_{x},\theta_{y}}}\text{dof}_{r}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{i}\right)\cdot\left(s^{E}\right)_{r}\cdot\text{dof}_{r}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{j}\right),&&\text{if}\quad i,j\in\mathcal{I}^{\theta_{x},\theta_{y}}.\end{aligned}\right.

(38)

where $\mathcal{I}^{u,v}$ , $\mathcal{I}^{w}$ and $\mathcal{I}^{\theta_{x},\theta_{y}}$ refer to the subset of degrees of freedom associated to the in-plane displacements, out-of-plane displacements and in-plane rotations. Therefore, the stabilization term is block diagonal.

Moreover, $\tau$ is a user-defined parameter, taken as $0.5$ in linear elasticity [35], and $\left(s^{E}\right)_{r}=\max\left(1,{\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\psi}_{r},\bm{\psi}_{r}\right)\right)$ , which scales the stabilization term as the diagonal of the consistency term.

Consequently, in the standard stabilized VEM, the second term in Eq. (34) is substituted with Eq. (38). Since this heuristic choice may spoil the accuracy of the method, self-stabilized formulations have been introduced in the literature. These formulations, by adopting higher-order polynomial projections, do not rely on the arbitrary stabilization term. Owing to the results in [23, 24], in which several stabilized and self-stabilized formulations from the literature have been compared in terms of accuracy and conditioning for different classes of PDEs, one stabilized and one self-stabilized strategy are adopted in the present work, corresponding to the most robust and suitable for elasticity problems.

Self-stabilized VEM

Self-stabilized VEM formulations do not rely on ad hoc stabilization terms. In contrast, they employ higher-order polynomial projections to ensure the stability of the linear system.

In [7, 23], it was demonstrated that strategies based on $L^{2}$ projections are generally more robust in the presence of variable coefficients. Owing to these outcomes, just one self-stabilized formulation based on the $L^{2}$ projection is used in the present work. Specifically, the one introduced in [15] for the lowest-order VEM in the context of the Laplace problem is here adapted to linear elasticity and generic order $p$ .

In order to compute higher-order polynomial projections, the spaces in Eq. (20) must be enlarged. For a generic $\phi$ , it holds that:

$\displaystyle\mathcal{\tilde{V}}_{h}^{r}\left(E\right)=\left\{\right.$	$\displaystyle{\displaystyle\phi}_{h}\in H^{1}\left(E\right)\cap C^{0}\left(E\right)\mathrel{\mathop{\mathchar 58\relax}}\bm{L}\left[\mathbb{C}\bm{\varepsilon}\left({\phi}_{h}\right)\right]\big\|_{E}\in\left[\mathcal{P}_{p+\ell_{E}}\left(E\right)\right]^{5},$	(39)
	$\displaystyle\left.{\phi}_{h}\big\|_{\Gamma^{E}_{e}}\in\mathcal{P}_{p}\left(\Gamma^{E}_{e}\right)\quad\forall e=1,\dots,n_{e},\;{\phi}_{h}\big\|_{\tilde{\Gamma}^{E}_{e}}\in\mathcal{\tilde{P}}_{p}\left(\tilde{\Gamma}^{E}_{e}\right)\quad\forall e=1,\dots,\tilde{n}_{e}\right.,$
	$\displaystyle\int_{E}{\phi}_{h}\;q\;\,\mathrm{d}E=\int_{E}{\Pi^{\mathcal{K}}_{p}}^{\phi}{\phi}_{h}\;q\;\,\mathrm{d}E\quad\forall q\in\mathcal{P}_{p}\left(E\right)\setminus\mathcal{P}_{p-r}\left(E\right),$
	$\displaystyle\left.\int_{E}{\phi}_{h}\;q\;\,\mathrm{d}E=\int_{E}{\Pi^{0}_{p}}^{\phi}{\phi}_{h}\;q\;\,\mathrm{d}E\quad\forall q\in\mathcal{P}_{p+\ell_{E}}\left(E\right)\setminus\mathcal{P}_{p}\left(E\right)\right\},$

and the spaces for the different displacement components are obtained as:

\mathcal{\tilde{V}}_{h}^{u}\left(E\right)=\mathcal{\tilde{V}}_{h}^{v}\left(E\right)=\mathcal{\tilde{V}}_{h}^{2}\left(E\right),\qquad\mathcal{\tilde{V}}_{h}^{w}\left(E\right)=\mathcal{\tilde{V}}_{h}^{1}\left(E\right),\qquad\mathcal{\tilde{V}}_{h}^{\theta_{x}}\left(E\right)=\mathcal{\tilde{V}}_{h}^{\theta_{y}}\left(E\right)=\mathcal{\tilde{V}}_{h}^{0}\left(E\right).

(40)

The subscript $\ell_{E}$ denotes the required increase in polynomial order to obtain a non-singular system. The operators ${\Pi^{\mathcal{K}}_{p}}^{\phi}$ and ${\Pi^{0}_{p}}^{\phi}$ are the scalar versions of $\bm{\Pi}^{\mathcal{K}}_{p}$ and $\bm{\Pi}_{p}^{0}$ associated with the generic displacement component.

The $L^{2}$ projection for the self-stabilized VEM reads:

\displaystyle\left(\bm{\varepsilon}\left(\bm{\psi}_{i}\right),\mathbb{C}\bm{q}^{\star}\right)_{E}=\left(\bm{\Pi}^{\mathcal{K}}_{p+\ell_{E}-1}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T},\mathbb{C}\bm{q}^{\star}\right)_{E}\quad\forall\bm{q}^{\star}\in\bm{\mathcal{P}}^{\star}_{p+\ell_{E}-1}\left(E\right),

(41)

for all $i=1,\dots,N_{d}\big|_{E}$ , where $\bm{q}^{\star}$ and $\bm{\mathcal{P}}^{\star}_{p+\ell_{E}-1}\left(E\right)$ are the polynomial vector and polynomial space referred to the self-stabilized formulation. The projector operator $\bm{\Pi}^{\mathcal{K}}_{p+\ell_{E}-1}$ maps the strains and not the functions themselves, as in the stabilized VEM, in the polynomial space.

The procedure to obtain the projection is the same as in standard stabilized VEM; hence, integration by parts is applied to the left-hand side. For what concerns line integrals, they are evaluated as in stabilized VEM. Conversely, surface integrals cannot be evaluated solely from the internal degrees of freedom, as in stabilized formulations. Therefore, the enlarged spaces in Eq. (40) are used to evaluate higher-order polynomials. Once the projector is found, the bilinear form is evaluated as in Eq. (34), by considering only the consistency term and the self-stabilized projection operator.

To determine the augmented polynomial order $\ell_{E}$ , two conditions are verified. First, the number of generalized strain modes must be greater than or equal to the number of degrees of freedom minus the number of rigid body motions. Second, a full-rank condition is enforced.

With this purpose, an algorithm is implemented that iteratively increases $\ell_{E}$ until this condition is satisfied. More details can be found in [23, 24].

3.3.2 Remaining forms

The projection procedure and the resulting discrete forms for the geometric stiffness and mass matrices, as well as for the body and thermal load vectors, follow the same construction adopted for the stiffness matrix. The stabilization term is required only for the stiffness matrix to ensure invertibility of the linear system. For all the other contributions, only the consistency term is retained. A summary of the discrete forms is provided in Table 2.

Continuous form	Projection	Discrete form
Stiffness matrix (stabilized VEM)
$\begin{aligned} \mathcal{K}^{E}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)=\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{\psi}_{j}\right)\,\mathrm{d}E\end{aligned}$ $\begin{aligned} {\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)=\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{\psi}_{j}\right)\,\mathrm{d}E\end{aligned}$	for all $i=1,\dots,N_{d}\big\|_{E}$ : $\left\{\begin{aligned} &\mathcal{K}^{E}\left(\bm{\psi}_{i},\bm{q}\right)=\mathcal{K}^{E}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)\\ &P_{0}\left(\bm{\psi}_{i},\bm{q}\right)=P_{0}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{0}\left(E\right)\oplus\hat{\bm{P}}\left(E\right)\end{aligned}\right.$	$\begin{aligned} {\mathcal{K}_{\mathrm{VC}}}^{E}_{h}&\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)={\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i},\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{j}\right)\\ &+S^{E}\left(\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{i},\left(\bm{I}-\bm{\Pi}^{\mathcal{K}}_{p}\right)\bm{\psi}_{j}\right)\end{aligned}$
Stiffness matrix (self-stabilized VEM)
$\begin{aligned} {\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)=\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{\psi}_{j}\right)\,\mathrm{d}E\end{aligned}$	for all $i=1,\dots,N_{d}\big\|_{E}$ : $\begin{aligned} \left(\bm{\varepsilon}\left(\bm{\psi}_{i}\right),\mathbb{C}\bm{q}^{\star}\right)_{E}=\left(\bm{\Pi}^{\mathcal{K}}_{p+\ell_{E}-1}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T},\mathbb{C}\bm{q}^{\star}\right)_{E}\end{aligned}$ $\begin{aligned} \forall\bm{q}^{\star}\in\bm{\mathcal{P}}^{\star}_{p+\ell_{E}-1}\left(E\right)\end{aligned}$	${\mathcal{K}_{\mathrm{VC}}}^{E}_{h}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)=\left(\bm{\Pi}^{\mathcal{K}}_{p+\ell_{E}-1}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T},\mathbb{C}\left(x,y\right)\bm{\Pi}^{\mathcal{K}}_{p+\ell_{E}-1}\bm{\varepsilon}\left(\bm{\psi}_{j}\right)\right)_{E}$
Geometric stiffness matrix
$\begin{aligned} \mathcal{G}^{E}\left(\psi_{i},\psi_{j}\right)=\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(\psi_{i}\right)^{T}\mathbb{N}\bm{\varepsilon}_{\mathrm{b}}\left(\psi_{j}\right)\,\mathrm{d}E\end{aligned}$ $\begin{aligned} {\mathcal{G}_{\mathrm{VC}}}^{E}\left(\psi_{i},\psi_{j}\right)=\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(\psi_{i}\right)^{T}\mathbb{N}\left(x,y\right)\bm{\varepsilon}_{\mathrm{b}}\left(\psi_{j}\right)\,\mathrm{d}E\end{aligned}$	for all $i=1,\dots,N_{d}\big\|_{E}^{w}$ : $\left\{\begin{aligned} &\mathcal{G}^{E}\left(\psi_{i},q\right)=\mathcal{G}^{E}\left(\bm{\Pi}^{\mathcal{G}}_{p}\psi_{i},q\right)&&\forall q\in\mathcal{P}_{p}\left(E\right)\\ &P_{0}\left(\psi_{i},q\right)=P_{0}\left(\bm{\Pi}^{\mathcal{G}}_{p}\psi_{i},q\right)&&\forall q\in\mathcal{P}_{0}\left(E\right)\end{aligned}\right.$	${\mathcal{G}_{\mathrm{VC}}}^{E}_{h}\left(\psi_{i},\psi_{j}\right)={\mathcal{G}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{\mathcal{G}}_{p}\psi_{i},\bm{\Pi}^{\mathcal{G}}_{p}\psi_{j}\right)$
Mass matrix
$\begin{aligned} \mathcal{M}^{E}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)=\int_{E}\bm{\psi}_{i}^{T}\mathbb{M}\bm{\psi}_{j}\,\mathrm{d}E\end{aligned}$	for all $i=1,\dots,N_{d}\big\|_{E}$ : $\begin{aligned} &\mathcal{M}^{E}\left(\bm{\psi}_{i},\bm{q}\right)=\mathcal{M}^{E}\left(\bm{\Pi}^{\mathcal{M}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)\end{aligned}$	$\begin{aligned} \mathcal{M}^{E}_{h}\left(\bm{\psi}_{i},\bm{\psi}_{j}\right)=\mathcal{M}^{E}\left(\bm{\Pi}^{\mathcal{M}}_{p}\bm{\psi}_{i},\bm{\Pi}^{\mathcal{M}}_{p}\bm{\psi}_{j}\right)\end{aligned}$
Body forces vector
$\mathcal{F}^{E}\left(\bm{\psi}_{i}\right)=\int_{E}\bm{\psi}_{i}^{T}\bm{\bar{b}}\,\mathrm{d}E$	for all $i=1,\dots,N_{d}\big\|_{E}$ : $\left(\bm{\psi}_{i}^{T},\bm{q}\right)_{E}=\left(\bm{\Pi}_{p}^{0}\bm{\psi}_{i}^{T},\bm{q}\right)_{E}\hskip 8.50012pt\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)$	$\mathcal{F}^{E}_{h}\left(\bm{\psi}_{i}\right)=\mathcal{F}^{E}\left(\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\right)$
Thermal forces vector
$\begin{aligned} \mathcal{T}^{E}\left(\bm{\psi}^{}_{i},\bm{\psi}^{}_{j}\right)=\int_{E}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\psi}^{}_{i}\right)^{T}\check{\bm{R}}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\psi}^{}_{j}\right)\,\mathrm{d}E\end{aligned}$ $\begin{aligned} {\mathcal{T}_{\mathrm{VC}}}^{E}\left(\bm{\psi}^{}_{i}\right)=\int_{E}\check{\bm{R}}\left(x,y\right)\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\psi}^{}_{i}\right)\,\mathrm{d}E\end{aligned}$	for all $i=1,\dots,N_{d}\big\|_{E}^{}$ : $\left\{\begin{aligned} &\mathcal{T}^{E}\left(\bm{\psi}^{}_{i},\bm{q}^{}\right)=\mathcal{T}^{E}\left(\bm{\Pi}^{\mathcal{T}}_{p}\bm{\psi}^{}_{i},\bm{q}^{}\right)&&\forall\bm{q}^{}\in\bm{\mathcal{P}}^{}_{p}\left(E\right)\\ &P_{0}\left(\bm{\psi}^{}_{i},\bm{q}^{}\right)=P_{0}\left(\bm{\Pi}^{\mathcal{T}}_{p}\bm{\psi}^{}_{i},\bm{q}^{}\right)&&\forall\bm{q}^{}\in\bm{\mathcal{P}}^{}_{0}\left(E\right)\oplus\hat{\bm{P}}^{}\left(E\right)\end{aligned}\right.$	${\mathcal{T}_{\mathrm{VC}}}^{E}_{h}\left(\bm{\psi}^{}_{i}\right)={\mathcal{T}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{\mathcal{T}}_{p}\bm{\psi}^{}_{i}\right)$

Table 2: Summary of remaining forms.

The projection is of elliptic type for the stabilized stiffness matrix, the geometric stiffness matrix and the thermal load vector, and of $L^{2}$ type for the self-stabilized stiffness matrix, the mass matrix and the body force vector. For the geometric stiffness matrix, only the contributions associated with the out-of-plane displacement $w$ are retained, as the remaining terms are neglected in the buckling formulation and would otherwise lead to a singular system.

For the mass matrix, the density is assumed constant. So, the inertial coefficient matrix $\mathbb{M}$ is independent of $x$ and $y$ . For the thermal load vector, the out-of-plane displacement $w$ is excluded, as it is assumed not to contribute to thermal effects.

The line load vector is constructed following the standard finite element procedure, since the trial functions are polynomials, or polynomial images, along the element edges. Prescribed displacements are enforced at the assembled system level, as in standard finite elements.

3.4 Discrete forms and projection operators with variable coefficients

In variable stiffness laminates, the elastic properties are not constant over the domain, so proper handling is required for the projection operators. In the previous section, the elastic coefficients were approximated as a constant within each element, which is appropriate for low-order approaches based on $h$ -refinement, but is, in general, not suitable within a $p$ -refinement framework.

To overcome this issue, the Variable Coefficients-VEM approach ( $\mathrm{VC}$ -VEM) proposed by the authors in the recent work [23] is here employed. As opposed to standard VEM practice, the coefficient variability is directly included in the projection, and non-computable terms are approximated using the $L^{2}$ VEM projector $\bm{\Pi}_{p}^{0}$ . The $\mathrm{VC}$ -VEM applies to the forms with a non-constant coefficient, i.e. stiffness matrix, geometric stiffness matrix, and thermal forces vector. A summary of the corresponding projections is provided in Table 3.

Term	Projection
Stiffness matrix (stabilized)	for all $i=1,\dots,N_{d}\big\|_{E}$ : $\left\{\begin{aligned} &{\mathcal{K}_{\mathrm{VC}}}_{*}^{E}\left(\bm{\psi}_{i},\bm{q}\right)={\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)\\ &P_{0}\left(\bm{\psi}_{i},\bm{q}\right)=P_{0}\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{0}\left(E\right)\oplus\hat{\bm{P}}\left(E\right)\end{aligned}\right.$
Stiffness matrix (self-stabilized)	for all $i=1,\dots,N_{d}\big\|_{E}$ : $\begin{aligned} \left(\bm{\varepsilon}^{}\left(\bm{\psi}_{i}\right),\mathbb{C}\left(x,y\right)\bm{q}^{\star}\right)_{,E}=\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p+\ell_{E}-1}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T},\mathbb{C}\left(x,y\right)\bm{q}^{\star}\right)_{E}\hskip 9.24994pt\forall\bm{q}^{\star}\in\bm{\mathcal{P}}^{\star}_{p+\ell_{E}-1}\left(E\right)\end{aligned}$
Geometric stiffness matrix	for all $i=1,\dots,N_{d}\big\|_{E}^{w}$ : $\left\{\begin{aligned} &{\mathcal{G}_{\mathrm{VC}}}_{*}^{E}\left(\psi_{i},q\right)={\mathcal{G}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{{\mathcal{G}_{\mathrm{VC}}}}_{p}\psi_{i},q\right)&&\forall q\in\mathcal{P}_{p}\left(E\right)\\ &P_{0}\left(\psi_{i},q\right)=P_{0}\left(\bm{\Pi}^{{\mathcal{G}_{\mathrm{VC}}}}_{p}\psi_{i},q\right)&&\forall q\in\mathcal{P}_{0}\left(E\right)\end{aligned}\right.$
Thermal forces vector	for all $i=1,\dots,N_{d}\big\|_{E}^{}$ : $\left\{\begin{aligned} &{\mathcal{T}_{\mathrm{VC}}}_{}^{E}\left(\bm{\psi}^{}_{i},\bm{q}^{}\right)={\mathcal{T}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{{\mathcal{T}_{\mathrm{VC}}}}_{p}\bm{\psi}^{}_{i},\bm{q}^{}\right)&&\forall\bm{q}^{}\in\bm{\mathcal{P}}^{}_{p}\left(E\right)\\ &P_{0}\left(\bm{\psi}^{}_{i},\bm{q}^{}\right)=P_{0}\left(\bm{\Pi}^{{\mathcal{T}_{\mathrm{VC}}}}_{p}\bm{\psi}^{}_{i},\bm{q}^{}\right)&&\forall\bm{q}^{}\in\bm{\mathcal{P}}^{}_{0}\left(E\right)\oplus\hat{\bm{P}}^{*}\left(E\right)\end{aligned}\right.$

Table 3: Projection for VC-VEM.

For the stiffness matrix, the procedure described above is extended to the variable stiffness case, i.e., spatially varying constitutive laws $\mathbb{C}\left(x,y\right)$ . The idea behind $\mathrm{VC}$ -VEM is to compute the polynomial projection $\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}$ as the solution of the following problem, which directly involves the bilinear form:

\left\{\begin{aligned} &{\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\psi}_{i},\bm{q}\right)={\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)\\ &P_{0}\left(\bm{\psi}_{i},\bm{q}\right)=P_{0}\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{0}\left(E\right)\oplus\hat{\bm{P}}\left(E\right).\end{aligned}\right.

(42)

However, while the right-hand side is fully computable as it involves only polynomials (and the coefficients), the left-hand side is not. Indeed, integration by parts yields the following identity:

	$\displaystyle{\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\psi}_{i},\bm{q}\right)=$	$\displaystyle-\int_{E}\bm{\psi}_{i}\cdot\bm{L}_{d}\left[\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}\right)\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E}+$		(43)
		$\displaystyle+\int_{E}\bm{\varepsilon}^{\mathrm{\theta}}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}\right)\,\mathrm{d}E,$		(43)

where it is clear that the surface integrals are not computable through the degrees of freedom. Eq. (42) is then modified as follows:

\left\{\begin{aligned} &{\mathcal{K}_{\mathrm{VC}}}_{*}^{E}\left(\bm{\psi}_{i},\bm{q}\right)={\mathcal{K}_{\mathrm{VC}}}^{E}\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)\\ &P_{0}\left(\bm{\psi}_{i},\bm{q}\right)=P_{0}\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\bm{\psi}_{i},\bm{q}\right)&&\forall\bm{q}\in\bm{\mathcal{P}}_{0}\left(E\right)\oplus\hat{\bm{P}}\left(E\right),\end{aligned}\right.

(44)

where the left-hand side is the modified form defined as:

	$\displaystyle{\mathcal{K}_{\mathrm{VC}}}_{*}^{E}\left(\bm{\psi}_{i},\bm{q}\right)\mathrel{\mathop{\mathchar 58\relax}}=$	$\displaystyle-\int_{E}\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\cdot\bm{L}_{d}\left[\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}\right)\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E}+$		(45)
		$\displaystyle+\int_{E}\bm{\varepsilon}^{\mathrm{\theta}}\left(\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}\right)\,\mathrm{d}E.$		(45)

${\mathcal{K}_{\mathrm{VC}}}_{*}^{E}$ is a computable approximation of the original bilinear form through the $L^{2}$ projection operator $\bm{\Pi}_{p}^{0}$ , and $\hat{\bm{S}}$ contains the spatial variation of $\mathbb{C}\left(x,y\right)$ . The evaluation of the line integral does not require special care, although the non-polynomial nature of $\mathbb{C}\left(x,y\right)$ leads to a non-exact integration.

As implied by Eq. (45), the coefficient $\mathbb{C}\left(x,y\right)$ should be at least of class $\mathcal{C}^{1}$ within each element. This is a suitable assumption for the problems under consideration, as the constitutive law $\mathbb{C}\left(x,y\right)$ is obtained through standard laminate stiffness matrices assembly and the fiber orientations are interpolated via Lagrange polynomials.

The same strategy applies to the self-stabilized VEM formulation, as well as to the geometric stiffness matrix and the thermal load vector. The related $\mathrm{VC}$ -VEM projections are collected in Table 3. The subscript “ $*$ ” denotes the modified forms defined by replacing the trial functions with their polynomial projection in the non-computable integrals, after having applied integration by parts. In addition, $\mathbb{N}\left(x,y\right)$ and $\check{\bm{R}}\left(x,y\right)$ are required to be of class $\mathcal{C}^{1}$ within each element, a condition easily satisfied for the class of problems considered in this investigation.

The spatial variation of $\mathbb{N}\left(x,y\right)$ depends not only on the constitutive law $\mathbb{C}\left(x,y\right)$ , but also on the strain field. Therefore, even for constant fiber orientation, $\mathbb{N}\left(x,y\right)$ remains spatially dependent.

The final discrete bilinear forms are constructed following the standard VEM formulation, with the stabilization term of the stiffness matrix defined using the standard projector $\bm{\Pi}^{\mathcal{K}}_{p}$ , which does not include the spatial variation of the constitutive tensor in the projection.

4 Results

In this section, the accuracy and robustness of the proposed VEM formulation are assessed by comparison with analytic solutions and benchmark problems from the literature, as well as with commercial finite element simulations conducted using Abaqus. First, the convergence of the method is assessed through three representative test cases in a $p$ -refinement setting. Specifically, a first test case serves to investigate the performance of standard and $\mathrm{VC}$ -VEM, in their stabilized and self-stabilized versions, in the presence of different layups, ranging from constant to curvilinear fiber orientations. Moreover, the static and free-vibration responses of problems featuring high-gradient solutions are investigated to demonstrate the effectiveness of performing local mesh refinements. Lastly, the method is validated against results from the literature and Abaqus for both free-vibration and buckling analyses. Geometries of varying complexity are considered to illustrate the ability of the method to handle general geometries. Different meshes and approximation orders are considered throughout the section. In all the cases, the number of integration points is selected to ensure sufficient accuracy of the results. For each element, all curve types are interpolated as Bézier curves, which are selected in this work for their robustness and simplicity of implementation.

4.1 Test case 1

The first test case investigates the convergence properties of the standard VEM and its $\mathrm{VC}$ -version in a $p$ -refinement framework. In particular, this analysis extends the results of [23] to a more complex domain with a cutout, curvilinear fiber orientations, and both membrane and bending behavior.

The plate is square with side $a=1000$ mm, and a circular cutout of radius $R=150$ mm is centered in the middle of the plate. A schematic representation is available in Figure 3, where the essential boundary conditions, prescribed to all the displacement components $\bm{u}=\{u\,v\,w\,\theta_{x}\,\theta_{y}\}^{\mathrm{T}}$ , are reported.

An orthotropic material with properties $E_{11}=181000$ MPa, $E_{22}=10273$ MPa, $G_{12}=G_{13}=G_{23}=7170.5$ MPa and $\nu_{12}=0.28$ is considered. The laminate has sixteen plies with thickness $h=0.1272$ mm each. Two different layups are investigated:

\displaystyle\text{layup 1}\mathrel{\mathop{\mathchar 58\relax}}[\pm 90\langle 45|45\rangle]_{4s},\qquad\text{layup 2}\mathrel{\mathop{\mathchar 58\relax}}[\langle 90|0\rangle_{8},\langle-90|0\rangle_{8}].

(46)

These layups feature an increasing level of complexity. The first one displays uniform properties over the domain, while the second accounts for variable properties through a linear variation of the fiber angle and adds membrane-bending coupling due to the asymmetry of the stacking sequence.

The errors of the VEM solution are evaluated with respect to the exact solution of the problem, formulated via an in inverse approach: the exact solution is postulated in advance and the corresponding loading conditions are retrieved. In particular, the solution is imposed to be:

\bm{u}\left(x,y\right)=\sin\left(\pi x\right)\sin\left(\pi y\right)\begin{Bmatrix}1&1&1&1&1\end{Bmatrix}^{T},

(47)

where fixed boundary conditions are assumed along the outer edges.

Error estimates are conducted using the energy norm:

\displaystyle e_{\bm{\varepsilon}}^{\bm{u}}=\frac{\left(\sum_{E\in\mathcal{T}_{h}}\mathinner{\!\left\lVert\sqrt{\mathbb{C}\left(x,y\right)}\bm{\varepsilon}\left(\bm{u}-\bm{\Pi}_{p}^{0}\bm{u}_{h}\right)\right\rVert}_{0,E}^{2}\right)^{\frac{1}{2}}}{\mathinner{\!\left\lVert\sqrt{\mathbb{C}\left(x,y\right)}\bm{\varepsilon}\left(\bm{u}\right)\right\rVert}_{0,E}}.

(48)

A $p$ -refinement strategy is adopted, with $p=1,\dots,10$ . Both standard and $\mathrm{VC}$ -VEM, as well as stabilized and self-stabilized VEM, are employed, thereby allowing a comprehensive comparison of the different strategies in the presence of curved edges, high approximation orders, and complex layup configurations.

The error estimates are reported in Figure 4 for the two layups at hand. The mesh, which features curved edges, is plotted in the same figures next to the legend.

For the first layup, Figure 4(a), the results demonstrate equal accuracy across the different strategies. In contrast, when the fiber orientation is no longer constant (layup 2 in Figure 4(b)), the standard stabilized VEM loses its accuracy as the order $p$ increases. Instead, both the stabilized and self-stabilized $\mathrm{VC}$ -VEM, and the standard self-stabilized VEM, which adopts a $L^{2}$ projection, maintain a good level of accuracy, despite the increased complexity given by the non-uniform stiffness distribution and membrane-bending coupling. Across all layups, some oscillations are present at lower orders; while, as the order increases, a smoother trend is observed. For layup 2, the error is slightly higher than for layup 1. This behavior is ascribed to the complexity introduced by the presence of the curvilinear fiber orientations.

These findings are in agreement with those presented in [7], in which it is shown that $L^{2}$ projections typically perform better in the presence of variable coefficients, and with those in [23], in which it is demonstrated that $\mathrm{VC}$ -VEM maintains optimal accuracy even with elliptic projections. The present test case broadens the ones analyzed in [23] to a more complex scenario. Indeed, while [23] investigates the performance of $\mathrm{VC}$ -VEM by analyzing the membrane behavior of a square domain with a polynomial variation of the elastic properties, the present work assesses the effectiveness of the method by considering a plate with a cutout, curvilinear fiber orientations, and membrane-bending coupling.

4.2 Test case 2

The second test case regards a cracked panel under tension load, with a singular stress state at the crack tip. This is a classical benchmark to assess the convergence properties of numerical methods, as the presence of the singularity makes the convergence particularly challenging [55]. In this work, this test case is of interest to demonstrate the potential of VEM to perform mesh refinements where desired, thereby enhancing the convergence properties in the presence of singularities. Hanging nodes, i.e. nodes located along the edges of an element, are used to preserve mesh conformity in the refined areas.

The configuration is reported in Figure 5. By exploiting the symmetry of the problem, only half of the panel is considered and the crack is simulated via Dirichlet boundary conditions. The panel has a half-length $a=2$ mm and is made of isotropic material with $E=10$ MPa and $\nu=0.3$ . Plane strain conditions are assumed.

The analytical stress field in polar coordinates $\left(\theta,r\right)$ is given by [55]:

		$\displaystyle\sigma_{xx}=\frac{K_{1}}{\sqrt{2\pi r}}\cos{\frac{\theta}{2}}\left(1-\sin{\frac{\theta}{2}}\sin{\frac{3\theta}{2}}\right),$		$\displaystyle\sigma_{yy}=\frac{K_{1}}{\sqrt{2\pi r}}\cos{\frac{\theta}{2}}\left(1+\sin{\frac{\theta}{2}}\sin{\frac{3\theta}{2}}\right),$		(49)
		$\displaystyle\sigma_{xy}=\frac{K_{1}}{\sqrt{2\pi r}}\sin{\frac{\theta}{2}}\cos{\frac{\theta}{2}}\cos{\frac{3\theta}{2}},$				(49)

where $K_{1}=1$ is the stress intensity factor. Neumann boundary conditions are applied along the free edges [55]:

		$\displaystyle\bm{t}_{1}=\bm{\sigma}\cdot\bm{n}=-\begin{bmatrix}\sigma_{xx}\\ \sigma_{xy}\end{bmatrix}\quad\forall\bm{x}\in\Gamma_{1},$		$\displaystyle\bm{t}_{2}=\bm{\sigma}\cdot\bm{n}=+\begin{bmatrix}\sigma_{xy}\\ \sigma_{yy}\end{bmatrix}\quad\forall\bm{x}\in\Gamma_{2},$		(50)
		$\displaystyle\bm{t}_{3}=\bm{\sigma}\cdot\bm{n}=+\begin{bmatrix}\sigma_{xx}\\ \sigma_{xy}\end{bmatrix}\quad\forall\bm{x}\in\Gamma_{3}.$				(50)

The convergence properties are evaluated using the energy norm error:

e_{\bm{\varepsilon}}^{\bm{u}}=\frac{\sum_{E\in\mathcal{T}_{h}}\mathinner{\!\left\lVert\sqrt{\bm{A}}\bm{\varepsilon}^{0}\left(\bm{u}-\bm{\Pi}_{p}^{0}\bm{u}_{h}\right)\right\rVert}_{0,E}^{2}}{\mathinner{\!\left\lVert\sqrt{\bm{A}}\bm{\varepsilon}^{0}\left(\bm{u}\right)\right\rVert}_{0}}.

(51)

In the following, different refinement strategies are employed to assess the convergence properties of the method. Firstly, uniform $h$ -refinement and $p$ -refinement are adopted. Then, local $h$ -refinement is performed at the crack tip. The effect of mesh distortion is investigated, too. With this purpose, both structured quadrilateral and Voronoi meshes are considered.

All the simulations reported below are based on standard stabilized VEM, as preliminary analyses showed no significant differences compared to the self-stabilized variant.

Uniform $h$ -refinement and $p$ -refinement

The first investigation deals with the $h$ - and $p$ -refinement strategies. For the latter, the order $p$ is progressively increased from $1$ to $4$ . The meshes used in the simulations are presented in Figure 6. They consist of rectangular elements of increasing density.

A summary of the results is available in Figure 7, where the errors in the energy norm are plotted against the total number of degrees of freedom.

As shown in Figure 7(a), the convergence rate observed by $h$ -refinement is essentially independent of the order $p$ . This behavior is due to the presence of the singularity in the stress field, hence any increase in the order $p$ does not provide any benefit, and the convergence remains limited by this fixed rate.

Similar conclusions are drawn by inspection of Figure 7(b), where the results of the $p$ -refinement strategy are reported. In particular, an algebraic convergence rate with respect to $p$ can be observed. The convergence is not exponential due to the presence of the singularity at the crack tip. Furthermore, a similar convergence rate is achieved for the three different meshes considered here.

As a further assessment, the stress profile at the coordinate $y=0$ is reported for the different refinement strategies in Figure 8. This comparison aims at illustrating the quality of the predictions in the proximity of the crack tip. Three different models are considered for this purpose. The first one corresponds to the $h$ -refined mesh with $12\times 6$ elements, the second one is obtained via $p$ -refinement of the $4\times 2$ mesh, and the third one by combining these two strategies. The VEM stress is recovered from the polynomial projection via $\bm{\Pi}_{p}^{0}$ of the numerical solution.

As shown in Figure 8(a), the first strategy provides a relatively inaccurate piecewise constant description in the proximity of the stress peak. An improved stress profile prediction is observed with the second one, as revealed by Figure 8(b). The $hp$ -refined strategy, whose results are presented in Figure 8(c), further improves the quality of the predictions. However, noticeable discrepancies are still present.

Local $h$ -refinement

Moving from the results obtained in the previous section, the potential of VEM is now exploited to locally refine the grid in the proximity of the stress singularity. In particular, the size of the elements is progressively halved and mesh conformity is preserved by employing hanging nodes at the center of the edge of the interface elements with non-matching dimensions. This strategy is defined as local $h$ -refinement.

By defining as level 1 the uniform mesh, six more refinements are considered, ranging from level 2 up to level 7. In particular, level 2 corresponds to a mesh where two different sizes of $h$ coexist: one for the broad field and another in correspondence of the crack tip. Similarly, the other refinement levels correspond to an increasing number of progressively refined meshes. For clarity, different levels in the case of $h=12\times 6$ are presented in Figure 9.

The error curves are shown in Figure 10. Within the same local $h$ -refinement paradigm, two different studies are conducted. In the first case, three different mesh densities are considered, and refinement is operated while keeping fixed the order at $p=1$ . In the second case, the initial mesh density is fixed at $12\times 6$ elements. Different orders $p=1,\ldots,4$ are considered and each model is progressively refined locally.

Compared to the uniform refinement investigated earlier, the local $h$ -refinement of Figure 10(a) illustrates an improved convergence rate.

While the results above refer to the global energy response, it is interesting to address the local behavior in terms of stress distribution in the most critical region, i.e. the crack tip. For this reason, the stress profile and the contour of the stress component $\sigma_{xx}$ are reported in Figure 11 for the mesh with $h=12\times 6$ , $p=4$ , and $7$ refinement levels.

The quality of the stress prediction is excellent, as seen from Figure 11(b). A slight overshoot of the solution can be noted at the crack tip. This is a common effect of high-order approximations, but it does not affect the overall accuracy of the solution.

To conclude, the local $h$ -refinement strategy guarantees improved convergence and local stress prediction capability when compared to uniform $h$ - and $p$ -refinements approaches presented earlier. Owing to the VEM inherent ability to handle hanging nodes, ease of modeling and excellent accuracy-to-degrees of freedom ratio are guaranteed.

Non-structured mesh

The investigation is further extended to verify how the element regularity may affect the conclusions drawn in the previous section. So, a distorted Voronoi mesh is used. Two refinements are considered: $p$ -refinement performed on the level $1$ mesh, and local $h$ -refinement up to seven levels. Refinement levels $1$ , $4$ , and $7$ are shown in Figure 12. On the contrary, uniform $h$ -refinement is excluded in this analysis owing to inherent restrictions associated with the type of mesh at hand, since non-structured meshes may lose accuracy under uniform refinements due to increasing mesh distortion. The errors curves for $p$ and local $h$ -refinement are reported in Figures 13(a) and 13(b), respectively.

The curves in Figure 13 illustrate that diffused distortion of the elements has a detrimental effect on the convergence of the method. On the one hand, the results demonstrate that convergence is achieved despite the high degree of distortion. This robustness is a key feature of VEM. On the other hand, no benefits are now experienced by the local $h$ -refinement strategy against the pure $p$ -refinement. Indeed, the rate of convergence is similar in both cases. Moreover, some oscillations are present, especially for lower orders.

Remarks

The results obtained in this test case demonstrate the difficulties of the standard $h$ - and $p$ -refinement techniques to handle scenarios characterized by drastic stress gradients. In contrast, local $h$ -refinements are more suitable for these cases, as shown by noticeable benefits on the convergence plots. Furthermore, this example demonstrates the effectiveness of employing general polygonal elements: local mesh refinements are easily introduced where needed, still preserving the mesh conformity. Moreover, the effects of mesh distortion are investigated by adopting a non-structured Voronoi mesh, illustrating that the presence of distorted elements affects the convergence of the method to some extent, leading to oscillations while preserving the overall trend.

4.3 Test case 3

In this third test case, the convergence of the method is assessed for the free-vibration response of a highly anisotropic plate. The benchmark has been studied in [52, 46] and proves to be useful to investigate the VEM ability to capture localization induced by drastic material anisotropy.

A simply supported square plate of dimension $a=100$ mm and thickness $h=0.01$ mm is considered. The domain and the corresponding boundary conditions at the edges are specified in Figure 14.

The material is a highly anisotropic pre-preg P100/AS3501 carbon/epoxy, with elastic properties: $E_{11}=369000$ MPa, $E_{22}=5030$ MPa, $G_{12}=G_{13}=G_{23}=5240$ MPa, $\nu_{12}=\nu_{13}=\nu_{23}=0.31$ and $\rho=1.5\times 10^{-6}$ kg/mm³, exhibiting an artificially high orthotropic ratio. The laminate is composed of a single ply oriented at $45^{\circ}$ . So, membrane and flexural anisotropy are maximized. The combination of these material and layup features exacerbates the convergence challenges, making this benchmark of special interest.

The free vibration response of the plate is investigated in terms of the first nondimensional circular frequency, defined as:

\bar{\omega}=\omega\frac{a^{2}}{h}\sqrt{\frac{\rho}{E_{22}}}.

(52)

The error is measured as:

e_{\bar{\omega}}=\frac{\left|\bar{\omega}_{\text{VEM}}-\bar{\omega}_{\text{REF}}\right|}{\bar{\omega}_{\text{REF}}}.

(53)

The subscripts VEM and REF refer to the present VEM and the reference solution reported in [54], $\bar{\omega}=21.9263$ , obtained by application of a refined $ps$ -FEM based model.

As done previously, convergence is studied in terms of uniform $h$ -, $p$ - and local $h$ -refinement. Moreover, both the standard stabilized and self-stabilized VEM are employed to assess the influence of the stabilization term on the accuracy of the solution. This aspect is particularly relevant, as eigenvalue problems are sensitive to the choice of the stabilization term.

Uniform $h$ -refinement and $p$ -refinement

The meshes used for the $h$ -refinement are shown in Figure 15. The $p$ -refinement is conducted for orders $p=2,\dots,5$ .

For both the uniform $h$ and $p$ -refinement, the error plots are reported in Figure 16. The number of degrees of freedom refers solely to the flexural ones.

By adopting a uniform $h$ -refinement strategy, both the stabilized and the self-stabilized formulation exhibit a relatively flat trend, as revealed by Figures 16(a) and 16(b). This trend is ascribed to the presence of strong anisotropy-induced shear gradients in correspondence of the corners. Compared to the self-stabilized variant, the stabilized VEM yields higher errors. This behavior is motivated by the sensitivity of eigenvalue problems to the choice of the stabilization term. Therefore, self-stabilized formulations are preferable in these cases.

Regarding the $p$ -refinement in Figures 16(c) and 16(d), the stabilized VEM exhibits a faster error decay, but larger errors. Conversely, the self-stabilized formulation features some oscillations, but lower overall errors.

Local $h$ -refinement

Improved results can be obtained by application of the local $h$ -refinement strategy already presented in the previous test cases. For the problem at hand, local refinement is beneficial due to internal shear gradients in the proximity of the corners. For this reason, ten refinement levels are considered with increasing mesh density toward the corners. The different levels of refinement are shown in Figure 17, where the mesh $4\times 4$ is used as the initial reference discretization.

The error curves are reported in Figure 18, where a fixed mesh $h=4\times 4$ and different polynomial orders up to $p=5$ are considered.

The stabilized VEM features a flattening trend, indicating a slow convergence rate, as shown in Figure 18(a). Conversely, self-stabilized VEM achieves lower errors overall, but with oscillations and occasional error increases, as seen by inspection of Figure 18(b). The convergence rate improves significantly as the order is increased up to $p=5$ . This is consistent with the trend observed for progressively refined models, which tend to yield lower frequency estimates, suggesting convergence toward a more accurate solution.

Remarks

Consistently with the previous test case, this example further demonstrates the effectiveness of the virtual element method in enabling local $h$ -refinement to capture localized behavior. Since eigenvalue problems are particularly sensitive to the choice of the stabilization term, some differences can be observed between the stabilized and self-stabilized formulation. In particular, the stabilized VEM shows smoother convergence but with higher errors, whereas the self-stabilized VEM achieves lower errors although it exhibits some oscillations.

4.4 Test case 4

Having established the convergence properties of the method, its numerical performance is further assessed by investigating the response of plates with complex geometries. This benchmark, taken from [20, 19], considers the free-vibration and thermal buckling analysis of a square plate with a heart-shaped cutout.

The plate is square with side length $a=10$ m. The heart-shaped cutout is defined as a NURBS curve in the reference, whereas in the present work it is constructed by combining three simple geometries: a square of side $a_{1}=4$ m, centered at $\left(x_{s},y_{s}\right)=\left(4,6\right)$ m, and two circles of radius $R_{1}=R_{2}=a_{1}/2$ , centered at $\left(x_{c1},y_{c1}\right)=\left(4,4\right)$ m and $\left(x_{c2},y_{c2}\right)=\left(6,6\right)$ m, respectively.

Two configurations are considered with equal geometry, but different boundary conditions depending on the analysis type. In particular, the plate is simply supported or fully clamped, dependently on whether free-vibration or buckling analysis are considered. The dimensions and boundary conditions are summarized in Figure 19.

The material properties and layups are specified in the respective subsections.

Two structured meshes are employed, as shown in Figure 20. The first mesh is relatively coarse and is employed in combination with high approximation orders, whereas the second is finer and coupled with lower approximation orders. This setup enables an assessment of the influence of both mesh density and approximation order $p$ on the accuracy and computational efficiency of the method. Particular attention should be paid to the presence of elements with curved edges. Their treatment is particularly straightforward within the VEM framework: in the coarse mesh, such elements are relatively irregular, yet accurately captured thanks to the proposed curved-edge formulation; in the finer mesh, a different scenario is observed, where hanging nodes prove useful during the meshing process.

Free-vibrations

For the free-vibration analysis, a composite material with properties $E_{11}/E_{22}=2.45$ , $G_{12}/E_{22}=0.48$ , $G_{13}/E_{22}=G_{23}/E_{22}=0.2$ , $\nu_{12}=0.23$ , $\rho=8000$ kg/m³, and thickness $h=0.06$ m, is considered. Three different layups are analyzed:

\displaystyle\text{layup 1}\mathrel{\mathop{\mathchar 58\relax}}[15^{\circ}/-15^{\circ}/15^{\circ}],

\displaystyle\text{layup 2}\mathrel{\mathop{\mathchar 58\relax}}[30^{\circ}/-30^{\circ}/30^{\circ}],

\displaystyle\text{layup 3}\mathrel{\mathop{\mathchar 58\relax}}[45^{\circ}/-45^{\circ}/45^{\circ}].

(54)

The frequencies are normalized as [20]:

\bar{\omega}_{i}=\left(\frac{\rho\,h\,\omega_{i}\,a^{4}}{\overline{D}}\right)^{\frac{1}{2}},\quad\text{with}\quad\overline{D}=\frac{E_{11}\,h^{3}}{12\left(1-\nu_{12}\,\frac{E_{22}}{E_{11}}\,\nu_{12}\right)}.

(55)

In order to simplify the comparison across the three layups, only Mesh 2 with approximation order $p=5$ is employed. The standard stabilized VEM is employed.

The normalized natural frequencies for the three layups are reported in Table 4 and compared with those obtained via isogeometric analysis (IGA) in [20].

Mode	layup 1			layup 2			layup 3
Mode	VEM	IGA [20]	Error [%]	VEM	IGA [20]	Error [%]	VEM	IGA [20]	Error [%]
1	19.04	18.91	0.69	20.53	20.40	0.64	21.24	21.10	0.66
2	31.70	31.83	0.41	33.51	33.66	0.45	34.30	34.45	0.44
3	35.69	36.09	1.11	36.82	37.23	1.10	37.49	37.91	1.11
4	56.98	57.00	0.04	59.19	59.20	0.02	60.29	60.29	0.00
5	61.93	62.73	1.28	64.09	65.03	1.45	65.20	66.22	1.54
6	83.78	83.93	0.18	87.84	87.92	0.09	91.02	90.68	0.37

Table 4: Test case 4: nondimensional natural frequencies

\bar{\omega}

and normalized percentage errors with respect to the reference values.

In all cases, the results obtained with the present VEM formulation show good agreement with the reference results from [20], with a relative error below $1\%$ in most cases. In particular, for all layups, the highest error is observed for Mode 5, with a value slightly higher than $1\%$ .

Thermal buckling

For the thermal buckling benchmark [19], the material properties are $E_{11}/E_{22}=15$ , $G_{12}/E_{22}=0.5$ , $G_{13}/E_{22}=G_{23}/E_{22}=0.3356$ , $\nu_{12}=0.3$ , with $E_{22}=1$ GPa, and thermal expansion coefficients $\alpha_{11}/\alpha_{0}=0.015$ and $\alpha_{22}/\alpha_{0}=1$ , with $\alpha_{0}=0.025$ $1/^{\circ}$ C. The thickness is $h=0.1$ m, and the laminate layup is $[0^{\circ}/90^{\circ}/90^{\circ}/0^{\circ}]$ . The plate is loaded with a temperature gradient $\Delta T$ .

Both meshes shown in Figure 20 are employed in combination with the stabilized VEM. The approximation orders are selected to balance accuracy and computational efficiency, and the two following models are considered: Mesh 1 with a higher-order approximation ( $p=8$ ), and Mesh 2 with a lower-order approximation ( $p=5$ ).

The comparison is performed with an Abaqus model consisting of $2032$ S4R shell elements.

The nondimensional critical temperatures corresponding to the first five buckling modes are summarized in Table 5. The results are obtained with the standard stabilized VEM approach. So, the spatial dependency of the pre-stress distribution is neglected in the computation of the geometric stiffness matrix projection operator.

Mode	Mesh 1 and $p=8$		Mesh 2 and $p=5$		Abaqus
Mode	VEM	Error [%]	VEM	Error [%]	Abaqus
1	0.01197	0.50	0.01206	0.25	0.01203
2	0.01242	0.64	0.01254	0.32	0.01250
3	0.01308	0.15	0.01315	0.38	0.01310
4	0.01320	0.00	0.01324	0.30	0.01320
5	0.01438	0.69	0.01451	0.21	0.01448

Table 5: Test case 4: nondimensional critical temperatures

\alpha_{0}\Delta T^{cr}

and normalized percentage errors with respect to the reference values.

The results obtained with the present VEM formulation closely match the Abaqus ones. In all cases, the relative error in the nondimensional critical temperature is below $1\%$ for every mode. Comparing the two meshes, the combination of the coarser Mesh 1 with the higher approximation order $p=8$ yields slightly lower nondimensional critical temperatures than those obtained with the finer Mesh 2 with $p=5$ . This highlights that coarse discretization with high-order approximations can provide accurate results with a reduced number of degrees of freedom. Indeed, for Mesh 2 with $p=5$ , the number of degrees of freedom is $9684$ , whereas for Mesh 1 with $p=8$ it is $5877$ .

The VEM-Abaqus comparison in terms of buckling modes is provided in Figure 21. Mesh 1 with $p=8$ is considered, although similar results are obtained with Mesh 1, but they are omitted here for the sake of brevity.

The contours illustrate the close agreement between the two computational strategies, further demonstrating the correctness of the proposed implementation.

4.5 Test case 5

To further assess the effectiveness of the developed tool, the validation is now extended to variable stiffness plates with complex geometries. The test case is taken from [36], where the Ritz method is employed to investigate the buckling of a square plate with a circular cutout. The plate is square with dimensions $a=254$ mm, with a circular cutout of radius $R=50.8$ mm centered in the middle of the plate. The plate is simply supported: the out-of-plane deflections are prevented, while the in-plane conditions are free apart from the normal component on the loaded edges, as shown in Figure 22.

The material properties are $E_{11}=181000$ MPa, $E_{22}=10273$ MPa, $G_{12}=G_{13}=G_{23}=7170.5$ MPa and $\nu_{12}=0.28$ . The laminate consists of sixteen plies, each with thickness $h=0.1272$ mm. Four different layups are considered. In particular:

\displaystyle\text{layup 1}\mathrel{\mathop{\mathchar 58\relax}}[\pm 90\langle 0|75\rangle]_{4s},\;\;\text{layup 2}\mathrel{\mathop{\mathchar 58\relax}}[\pm 0\langle 60|30\rangle]_{4s},\;\;\text{layup 3}\mathrel{\mathop{\mathchar 58\relax}}[\pm 45\langle 15|0\rangle]_{4s},\;\;\text{layup 4}\mathrel{\mathop{\mathchar 58\relax}}[\pm 90\langle 45|45\rangle]_{4s}.

(56)

The first three layups are characterized by stiffness variability, while the last one is a constant stiffness configuration. The presence of both variable and constant stiffness plates is useful to compare standard and $\mathrm{VC}$ -VEM.

Four different meshes are chosen, as shown in Figure 23.

In particular, both structured and non-structured meshes are employed. The first two meshes are regular, while the third and fourth are characterized by the presence of elements with a relatively high degree of distortion. In all the cases, the curved edge feature is exploited to accurately represent the cutout.

The four meshes are associated with different orders $p$ to account for the different refinement levels. Mesh 1 is used with $p=7$ , Mesh 2 with $p=4$ , Mesh 3 with $p=6$ , and Mesh 4 with $p=4$ . Both standard and $\mathrm{VC}$ -VEM are employed, in their stabilized and self-stabilized versions.

A summary of the results is provided in Table 6 for the different meshes and stabilization techniques, and the comparison is presented against Ref. [36]. In particular, the results are presented in terms of the nondimensional ratio $\bar{N}_{xx}^{cr}/\bar{N}_{iso,xx}^{cr}$ , the former being the buckling resultant of the plate under investigation, the latter the buckling resultant of a corresponding isotropic plate without cutout and elastic properties $E=69668$ MPa and $\nu=0.296$ .

			layup 1	layup 2	layup 3	layup 4
			${\bar{N}_{xx}^{cr}/\bar{N}_{xx,iso}^{cr}}$ [-]
Mesh 1 & ${p=7}$	Stabilized	Standard	2.03	1.06	1.07	1.07
	Stabilized	VC	2.06	1.03	1.08	1.06
	Self-Stabilized	Standard	2.02	1.02	1.02	1.04
	Self-Stabilized	VC	2.01	1.01	1.04	1.04
Mesh 2 & ${p=4}$	Stabilized	Standard	2.01	1.04	1.04	1.06
	Stabilized	VC	2.05	1.04	1.05	1.07
	Self-Stabilized	Standard	2.03	1.03	1.04	1.06
	Self-Stabilized	VC	2.02	1.02	1.04	1.06
Mesh 3 & ${p=6}$	Stabilized	Standard	1.60	1.33	0.46	1.51
	Stabilized	VC	2.05	1.05	1.12	1.05
	Self-Stabilized	Standard	2.04	1.10	1.07	1.08
	Self-Stabilized	VC	2.04	1.09^a	1.13^a	1.08
Mesh 4 & ${p=4}$	Stabilized	Standard	1.81^a	1.18^a	1.04	0.96^a
	Stabilized	VC	2.03	0.99	1.04	1.06
	Self-Stabilized	Standard	2.05^a	1.04	1.05	1.07^a
	Self-Stabilized	VC	2.04^a	1.03	1.05	1.07
Ref. [36]			2.06	1.02	1.05	1.08

a

spurious modes were removed

Table 6: Test case 5: normalized buckling resultant

\bar{N}_{xx}^{cr}/\bar{N}_{xx,iso}^{cr}

for different meshes and stabilization techniques.

A first consideration regards the substantial agreement between the results obtained with the two structured meshes, i.e. Mesh 1 and Mesh 2, and the reference results. No significant discrepancies are observed in these cases, irrespective of the layup and adopted VEM strategy.

On the other hand, some noticeable deviations can be seen for the distorted Meshes 3 and 4. In particular, the standard stabilized VEM exhibits incorrect results for some configurations, both in terms of buckling load and mode shape. A first observation regards the modes predicted by Mesh 3 for the standard stabilized VEM: for all the layups, they deviate from the expected shape. Instead, regarding Mesh 4, several buckling loads are slightly underestimated and the corresponding mode shapes feature a milder amplitude, hence an increase in the approximation order $p$ may be necessary. Nonetheless, the order is not further augmented since $p=4$ provides sufficiently accurate results for the $\mathrm{VC}$ -VEM, as well as for both self-stabilized variants. Therefore, the same order $p=4$ is retained to ensure a proper comparison of the different formulations. In contrast, both stabilized and self-stabilized $\mathrm{VC}$ -VEM approaches, as well as the standard self-stabilized VEM, are able to correctly predict the buckling load even in the presence of mesh distortion. These results are in agreement with those presented in the first test case, further highlighting the superior robustness of the $L^{2}$ projection and $\mathrm{VC}$ -VEM in the presence of more complex scenarios.

Lastly, the buckling modes of the different layups are shown in Figure 24. The VEM results are obtained using Mesh 1 with stabilized $\mathrm{VC}$ -VEM.

As shown in Figure 24, the VEM predictions are in good agreement with the reference ones for all the layups considered. Overall, this test case further proves the effectiveness of the developed tool, and in particular of the $\mathrm{VC}$ -VEM, in predicting the buckling behavior of complex plate domains with variable stiffness properties.

5 Conclusions

This work presented a comprehensive high-order Virtual Element Method ( $p$ -VEM) framework for the static, free-vibration, and buckling analysis of variable stiffness plates with arbitrary shapes. The proposed formulation aimed at bridging the gap between existing VEM mathematical formulations and engineering applications to plate problems in structural mechanics. For this purpose, an implementation-oriented formulation based on standard finite element notation was developed. Furthermore, a number of advanced VEM capabilities have been integrated within a unified computational framework: arbitrary polygonal elements with curved edges, hanging nodes, stabilized and self-stabilized formulations, and a newly proposed Variable Coefficients VEM (VC-VEM) approach.

Five test cases were presented for problems involving both standard and innovative variable stiffness configurations. The numerical investigations demonstrate the accuracy and robustness of the method for a wide range of structural applications. In particular, the combination of high-order approximation with local mesh refinement is effective for problems characterized by localized phenomena, such as stress concentration; the use of arbitrary polygonal elements, curved edges, and hanging nodes offers an excellent potential to simplify the discretization of complex geometries. Furthermore, owing to the $p$ -VEM capabilities, accurate solutions are obtained on relatively coarse and distorted meshes, highlighting the robustness of the formulation and the effectiveness of $p$ -refinement. For variable stiffness laminates, the proposed $\mathrm{VC}$ -VEM formulation and the standard self-stabilized VEM, which employs a $L^{2}$ projection, provide improved treatment of the spatially varying constitutive properties, leading to higher accuracy than the standard stabilized VEM formulation, particularly for distorted meshes and higher approximation orders.

Overall, the test cases highlight both the general capabilities of the framework and its specific advantages for the analysis of variable stiffness structures. The proposed framework successfully combines the geometric flexibility of VEM with an accurate treatment of spatially varying constitutive properties, while maintaining an implementation-oriented formulation suitable for practical engineering applications. The resulting methodology therefore represents a versatile computational tool for the analysis of advanced composite structures characterized by complex geometries and non-uniform stiffness distributions.

References

[1] A. Alhajahmad, M.M. Abdalla, and Z. Gürdal (2008) Design tailoring for pressure pillowing using tow-placed steered fibers. Journal of Aircraft 45 (2), pp. 630–640. Cited by: §2.2.2.
[2] E. Artioli, L. Beirão da Veiga, C. Lovadina, and E. Sacco (2017) Arbitrary order 2D virtual elements for polygonal meshes: part I, elastic problem. Computational Mechanics 60, pp. 355–377. External Links: Document Cited by: §1.
[3] E. Artioli, S. de Miranda, C. Lovadina, and L. Patruno (2018) A family of virtual element methods for plane elasticity problems based on the Hellinger–Reissner principle. Computer Methods in Applied Mechanics and Engineering 340, pp. 978–999. External Links: Document Cited by: §1.
[4] F. Bassi, L. Botti, A. Colombo, D. A. Di Pietro, and P. Tesini (2012) On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations. Journal of Computational Physics 231 (1), pp. 45–65. External Links: Document Cited by: Polynomial space.
[5] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo (2013) BASIC principles of virtual element methods. Mathematical Models and Methods in Applied Sciences 23 (01), pp. 199–214. External Links: Document Cited by: §1, §3.1, §3.3.
[6] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo (2014) The hitchhiker’s guide to the virtual element method. Mathematical Models and Methods in Applied Sciences 24 (08), pp. 1541–1573. External Links: Document Cited by: §1, §3.3.1, Stiffness matrix.
[7] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo (2016) Virtual element method for general second-order elliptic problems on polygonal meshes. Mathematical Models and Methods in Applied Sciences 26 (04), pp. 729–750. External Links: Document Cited by: §1, §1, §3.3, §4.1.
[8] L. Beirão da Veiga, F. Brezzi, and L. D. Marini (2013) Virtual elements for linear elasticity problems. SIAM Journal on Numerical Analysis 51 (2), pp. 794–812. External Links: Document Cited by: §1.
[9] L. Beirão da Veiga, C. Lovadina, and D. Mora (2015) A virtual element method for elastic and inelastic problems on polytope meshes. Computer Methods in Applied Mechanics and Engineering 295, pp. 327–346. External Links: Document Cited by: §1.
[10] L. Beirão da Veiga, C. Lovadina, and A. Russo (2017) Stability analysis for the virtual element method. Mathematical Models and Methods in Applied Sciences 27 (13), pp. 2557–2594. External Links: Document Cited by: §1.
[11] L. Beirão da Veiga, D. Mora, and G. Rivera (2019) Virtual elements for a shear-deflection formulation of Reissner–Mindlin plates. Mathematics of Computation 88 (315), pp. 149–178. External Links: Document Cited by: §1.
[12] L. Beirão da Veiga, A. Russo, and G. Vacca (2019) The virtual element method with curved edges. ESAIM: Mathematical Modelling and Numerical Analysis 53 (2), pp. 375–404. External Links: Document Cited by: §1, §2, 2nd item, §3.2, §3, Stiffness matrix.
[13] S. Berrone, A. Borio, D. Fassino, and F. Marcon (2025) Stabilization-free Virtual Element Method for 2D second order elliptic equations. Computer Methods in Applied Mechanics and Engineering 438, pp. 117839. External Links: Document Cited by: §1.
[14] S. Berrone, A. Borio, and F. Marcon (2022) Comparison of standard and stabilization free virtual elements on anisotropic elliptic problems. Applied Mathematics Letters 129, pp. 107971. External Links: Document Cited by: §1.
[15] S. Berrone, A. Borio, and F. Marcon (2025) Lowest order stabilization free virtual element method for the 2D Poisson equation. Computers & Mathematics with Applications 177, pp. 78–99. External Links: Document Cited by: §1, §3.3.
[16] D. Boffi, F. Gardini, and L. Gastaldi (2020) Approximation of PDE eigenvalue problems involving parameter dependent matrices. Calcolo 57 (4), pp. 41. External Links: Document Cited by: §1.
[17] F. Brezzi and L. D. Marini (2013) Virtual element methods for plate bending problems. Computer Methods in Applied Mechanics and Engineering 253, pp. 455–462. External Links: Document Cited by: §1.
[18] B. H. Coburn and P. M. Weaver (2016) Buckling analysis, design and optimisation of variable-stiffness sandwich panels. International Journal of Solids and Structures 96, pp. 217–228. External Links: Document Cited by: §1.
[19] B. Devarajan and R. K. Kapania (2020) Thermal buckling of curvilinearly stiffened laminated composite plates with cutouts using isogeometric analysis. Composite Structures 238, pp. 111881. External Links: Document Cited by: §4.4, §4.4.
[20] B. Devarajan (2021) Free vibration analysis of curvilinearly stiffened composite plates with an arbitrarily shaped cutout using isogeometric analysis. arXiv preprint arXiv:2104.12856. Cited by: §4.4, §4.4, §4.4, §4.4, Table 4, Table 4, Table 4.
[21] A. M. D’Altri, S. de Miranda, L. Patruno, E. Artioli, and C. Lovadina (2020) Error estimation and mesh adaptivity for the virtual element method based on recovery by compatibility in patches. International Journal for Numerical Methods in Engineering 121 (19), pp. 4374–4405. External Links: Document Cited by: §1.
[22] A. M. D’Altri, L. Patruno, S. de Miranda, and E. Sacco (2022) First-order VEM for Reissner–Mindlin plates. Computational Mechanics 69, pp. 315–333. External Links: Document Cited by: §1, §3.2.
[23] P. P. Foligno, D. Boffi, F. Credali, and R. Vescovini (2026) Benchmarking stabilized and self-stabilized $p$ -virtual element methods with variable coefficients. Computer Methods in Applied Mechanics and Engineering 455, pp. 118863. Cited by: §1, §3.3, §3.3, §3.3, §3.3, §3.4, §3, §4.1, §4.1.
[24] P. P. Foligno (2026) Novel p-virtual element method for the nonlinear analysis of curvilinearly stiffened panels. Ph.D. Thesis, Politecnico di Milano. Cited by: §3.3, §3.3, Supporting material.
[25] R. Fujimoto and I. Saiki (2024) Study of the stabilization parameter in the virtual element method. Computer Methods in Applied Mechanics and Engineering 428, pp. 117106. External Links: Document Cited by: §1.
[26] T. A. Janssens and S. G. Castro (2021) Semi-analytical modelling of variable stiffness laminates with cut-outs. In AIAA Scitech 2021 Forum, pp. 0440. Cited by: §1.
[27] Z. Jing and L. Duan (2023) Discrete Ritz method for buckling analysis of arbitrarily shaped plates with arbitrary cutouts. Thin-Walled Structures 193, pp. 111294. External Links: Document Cited by: §1.
[28] Z. Jing and L. Duan (2024) Free vibration analysis of three-dimensional solids with arbitrary geometries using discrete Ritz method. Journal of Sound and Vibration 571, pp. 118132. External Links: Document Cited by: §1.
[29] A. Lamperti, M. Cremonesi, U. Perego, A. Russo, and C. Lovadina (2023) A Hu–Washizu variational approach to self-stabilized virtual elements: 2D linear elastostatics. Computational Mechanics 71 (5), pp. 935–955. External Links: Document Cited by: §1.
[30] K. M. Liew and C. M. Wang (1993) Pb-2 Rayleigh-Ritz method for general plate analysis. Engineering Structures 15 (1), pp. 55–60. External Links: Document Cited by: §2.1.
[31] F. S. Liguori, A. Madeo, S. Marfia, G. Garcea, and E. Sacco (2024) A stabilization-free hybrid virtual element formulation for the accurate analysis of 2D elasto-plastic problems. Computer Methods in Applied Mechanics and Engineering 431, pp. 117281. External Links: Document Cited by: §1.
[32] M.Visinoni (2024) A family of three-dimensional virtual elements for Hellinger-Reissner elasticity problems. Computers & Mathematics with Applications 155, pp. 97–109. External Links: Document Cited by: §1.
[33] G. Manickam, A. Bharath, A. N. Das, A. Chandra, and P. Barua (2018) Thermal buckling behaviour of variable stiffness laminated composite plates. Materials Today Communications 16, pp. 142–151. External Links: Document Cited by: §1.
[34] L. Mascotto (2018) Ill-conditioning in the virtual element method: Stabilizations and bases. Numerical Methods for Partial Differential Equations 34 (4), pp. 1258–1281. External Links: Document Cited by: §1, §3.3, Polynomial space, Polynomial space, Stiffness matrix.
[35] M. Mengolini, M. F. Benedetto, and A. M. Aragón (2019) An engineering perspective to the virtual element method and its interplay with the standard finite element method. Computer Methods in Applied Mechanics and Engineering 350, pp. 995–1023. External Links: Document Cited by: §1, §3.3, Stiffness matrix, Stiffness matrix, Stiffness matrix, Stiffness matrix.
[36] A. Milazzo, G. Guarino, and V. Gulizzi (2023) Buckling and post-buckling of variable stiffness plates with cutouts by a single-domain Ritz method. Thin-Walled Structures 182, pp. 110282. External Links: Document Cited by: 24(e), 24(e), 24(f), 24(f), 24(g), 24(g), 24(h), 24(h), §4.5, §4.5, Table 6.
[37] D. Mora, G. Rivera, and I. Velásquez (2018) A virtual element method for the vibration problem of Kirchhoff plates. ESAIM: Mathematical Modelling and Numerical Analysis 52 (4), pp. 1437–1456. External Links: Document Cited by: §1.
[38] D. Mora and I. Velásquez (2020) Virtual element for the buckling problem of Kirchhoff–Love plates. Computer Methods in Applied Mechanics and Engineering 360, pp. 112687. External Links: Document Cited by: §1.
[39] S. Nagendra, S. Kodiyalam, J.E. Davis, and V. Parthasaraty (1995) Optimization of tow fiber paths for composite design. In $36^{th}$ AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Material Conference, AIAA-1995-1275-CP, New Orleans, LA. Cited by: §2.2.2.
[40] R. Olmedo and Z. Gürdal (1993) Buckling response of laminates with spatially varying fiber orientations. In $34^{th}$ AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Material Conference, AIAA-93-1567-CP, La Jolla, CA, pp. 2261–2269. Cited by: §2.2.2.
[41] G. Raju, Z. Wu, and P. M. Weaver (2015) Buckling and postbuckling of variable angle tow composite plates under in-plane shear loading. International Journal of Solids and Structures 58, pp. 270–287. External Links: Document Cited by: §1.
[42] B. D. Reddy and D. van Huyssteen (2019) A virtual element method for transversely isotropic elasticity. Computational Mechanics 64 (4), pp. 971–988. External Links: Document Cited by: §1.
[43] J. N. Reddy (2004) Mechanics of Laminated Composite Plates and Shells. CRC Press LCC. Cited by: §2.2.1.
[44] A. Sommariva and M. Vianello (2007) Product Gauss cubature over polygons based on Green’s integration formula. BIT Numerical Mathematics 47, pp. 441–453. External Links: Document Cited by: Stiffness matrix.
[45] A. Sommariva and M. Vianello (2023) Low cardinality Positive Interior cubature on NURBS-shaped domains. BIT Numerical Mathematics 63 (2), pp. 22. External Links: Document Cited by: Stiffness matrix.
[46] R. Vescovini, L. Dozio, M. d’Ottavio, and O. Polit (2018) On the application of the Ritz method to free vibration and buckling analysis of highly anisotropic plates. Composite Structures 192, pp. 460–474. External Links: Document Cited by: §4.3.
[47] R. Vescovini and L. Dozio (2018) Thermal buckling behaviour of thin and thick variable-stiffness panels. Journal of Composites Science 2 (4), pp. 58. External Links: Document Cited by: §1.
[48] R. Vescovini and P. P. Foligno (2025) Geometrically nonlinear analysis of variable-stiffness plates using the R-functions combined with the Ritz method. Acta Mechanica, pp. 1–29. External Links: Document Cited by: §1.
[49] R. Vescovini, V. Oliveri, D. Pizzi, L. Dozio, and P. M. Weaver (2020) A semi-analytical approach for the analysis of variable-stiffness panels with curvilinear stiffeners. International Journal of Solids and Structures 188-189, pp. 244–260. External Links: Document Cited by: §2.2.2, §2.2.2, §2.2.2.
[50] R. Vescovini (2023) Ritz R-function method for the analysis of variable-stiffness plates. AIAA Journal 61 (6), pp. 2689–2701. External Links: Document Cited by: §1.
[51] Z. Wu, G. Raju, and P. M. Weaver (2018) Optimization of postbuckling behaviour of variable thickness composite panels with variable angle tows: Towards “Buckle-Free” design concept. International Journal of Solids and Structures 132, pp. 66–79. External Links: Document Cited by: §2.2.2.
[52] Z. Wu, G. Raju, and P.M. Weaver (2012) Comparison of variational, differential quadrature, and approximate closed-form solution methods for buckling of highly flexurally anisotropic laminates. Journal of Engineering Mechanics 139 (8), pp. 1073–1083. Cited by: §4.3.
[53] Z. Wu, P. M. Weaver, G. Raju, and B. C. Kim (2012) Buckling analysis and optimisation of variable angle tow composite plates. Thin-Walled Structures 60, pp. 163–172. External Links: Document Cited by: §1, §2.2.2.
[54] C. A. Yan and R. Vescovini (2023) Application of the $ps$ - version of the finite element method to the analysis of laminated shells. Materials 16 (4), pp. 1395. External Links: Document Cited by: §4.3.
[55] N. Zander, M. Ruess, T. Bog, S. Kollmannsberger, and E. Rank (2017) Multi-level $hp$ -adaptivity for cohesive fracture modeling. International Journal for Numerical Methods in Engineering 109 (13), pp. 1723–1755. External Links: Document Cited by: §4.2, §4.2, §4.2.
[56] W. Zhao and R. K. Kapania (2019) Prestressed vibration of stiffened variable-angle two laminated plates. AIAA Journal 57 (6), pp. 2575–2593. External Links: Document Cited by: §2.2.2.

Supporting material

Hereafter, the construction of the relevant quantities of the VEM formulation is presented in detail. A more in-depth derivation is available in [24].

Supporting material for Subsection 3.2: Spaces definition

Virtual Element space

The differential operator $\bm{L}\left[\cdot\right]$ defines the problem in its strong form, and its expression reads:

\bm{L}\left[\cdot\right]=\begin{bmatrix},x&0&,y&0&0&0&0&0\\ 0&,y&,x&0&0&0&0&0\\ 0&0&0&0&0&0&,x&,y\\ 0&0&0&,x&0&,y&-1&0\\ 0&0&0&0&,y&,x&0&-1\end{bmatrix}.

(57)

The dimensions of the local spaces defined in Eq. (20) are:

		$\displaystyle N_{d}\big\|_{E}^{s}=\text{dim}\left(\mathcal{V}_{h}^{s}\left(E\right)\right)=pn_{v}+\frac{\left(p-1\right)p}{2}\ \text{for }s=u,v,$		(58)
		$\displaystyle N_{d}\big\|_{E}^{w}=\text{dim}\left(\mathcal{V}_{h}^{w}\left(E\right)\right)=pn_{v}+\frac{p\left(p+1\right)}{2},$
		$\displaystyle N_{d}\big\|_{E}^{s}=\text{dim}\left(\mathcal{V}_{h}^{s}\left(E\right)\right)=pn_{v}+\frac{\left(p+1\right)\left(p+2\right)}{2}\ \text{for }s=\theta_{x},\theta_{y}.$

The dimension of the total local space of Eq. (21) is:

N_{d}\big|_{E}=\text{dim}\left(\bm{\mathcal{V}}_{h}\left(E\right)\right)=5pn_{v}+2\frac{\left(p-1\right)p}{2}+\frac{p\left(p+1\right)}{2}+2\frac{\left(p+1\right)\left(p+2\right)}{2}.

(59)

Polynomial space

The most common choice for $\mathcal{Q}_{p}\left(E\right)$ , which serves as a basis for $\mathcal{P}_{p}\left(E\right)$ , is the set of scaled monomials, here denoted as $\bar{q}_{\bm{n}}\left(x,y\right)$ . By denoting with $\left(x_{C},y_{C}\right)$ the coordinates of the centroid of $E$ , the basis is defined as [34]:

\bar{q}_{\bm{n}}\left(x,y\right)\big|_{E}=\left(\frac{x-x_{C}}{h_{E}}\right)^{n_{1}}\left(\frac{y-y_{C}}{h_{E}}\right)^{n_{2}}\quad\forall\bm{n}=\left[n_{1},n_{2}\right]=\left[\left(0,\dots,p\right),\left(0,\dots,p\right)\right].

(60)

However, scaled monomials suffer from numerical instability for higher values of $p$ . To improve stability, $L^{2}\left(E\right)$ orthonormal polynomials $q_{\bm{n}}\left(x,y\right)$ are employed in this work, obtained via Modified Gram-Schmidt (MGS) orthonormalization [34]:

q_{\bm{n}}\left(x,y\right)\big|_{E}=\sum_{\bm{\beta}=1}^{\bm{n}}\bm{GS}_{\bm{n},\bm{\beta}}\big|_{E}\bar{q}_{\bm{n}}\left(x,y\right)\big|_{E}\quad\forall\bm{n}=\left[n_{1},n_{2}\right]=\left[\left(0,\dots,p\right),\left(0,\dots,p\right)\right],

(61)

where $\bm{GS}\big|_{E}$ is the matrix of the orthonormalization coefficients, which are obtained for each element $E$ [4].

To build a complete polynomial basis of order $p$ , $\frac{\left(p+1\right)\left(p+2\right)}{2}$ independent polynomials are required. Hence, for the five displacement components, the dimension of the polynomial space is:

N_{p}\big|_{E}=\text{dim}\left(\bm{\mathcal{P}}_{p}\left(E\right)\right)=5\frac{\left(p+1\right)\left(p+2\right)}{2}.

(62)

Supporting material for Subsection 3.3: Matrices and vectors construction

Stiffness matrix

By using Eq. (28), the unknowns $\bm{u}_{h}$ can be compactly written as a function of the VEM trial functions as:

\bm{u}_{h}=\begin{Bmatrix}u_{h}\quad v_{h}\quad w_{h}\quad{\theta_{x}}_{h}\quad{\theta_{y}}_{h}\end{Bmatrix}^{T}=\begin{Bmatrix}{\bm{\Psi}^{u}}\quad{\bm{\Psi}^{v}}\quad{\bm{\Psi}^{w}}\quad{\bm{\Psi}^{\theta_{x}}}\quad{\bm{\Psi}^{\theta_{y}}}\end{Bmatrix}^{T}{\bm{c}_{h}}=\sum_{i=1}^{N_{d}\big|_{E}}\bm{\psi}_{i}{{\bm{c}_{h}}}_{i}=\bm{\Psi}{\bm{c}_{h}}.

(63)

where $\bm{\Psi}^{k}$ are $N_{d}\big|_{E}\times 1$ column vectors associated to the $k$ -th displacement component, and ${\bm{c}_{h}}$ is the vector of dimension $N_{d}\big|_{E}\times 1$ , collecting the degrees of freedom. The projection of the generic VEM trial functions $\bm{\psi}_{i}\big|_{i=1,\dots,N_{d}\big|_{E}}$ onto the polynomial space is defined as [35]:

\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i}=\sum_{\beta=1}^{N_{p}\big|_{E}}s_{i,\beta}\bm{q}_{\beta}\quad\forall i=1,\dots,N_{d}\big|_{E}.

(64)

By assembling the VEM trial functions and the polynomial basis functions column-wise, Eq. (64) becomes:

\bm{\Pi}^{\mathcal{K}}_{p}\bm{\Psi}=\bm{Q}\tilde{\bm{\Pi}}^{\mathcal{K}}_{p},

(65)

where $\bm{\Pi}^{\mathcal{K}}_{p}\bm{\Psi}$ is the matrix whose columns contain the projections of the VEM trial functions, and $\bm{Q}$ contains the polynomial basis functions. Matrix $\bm{\Pi}^{\mathcal{K}}_{p}\bm{\Psi}$ has dimension $5\times N_{d}\big|_{E}$ , while matrix $\bm{Q}$ has dimension $5\times N_{p}\big|_{E}$ . Accordingly, matrix $\tilde{\bm{\Pi}}^{\mathcal{K}}_{p}$ has dimension $N_{p}\big|_{E}\times N_{d}\big|_{E}$ and its $i$ -th column contains the polynomial coefficients of the projection of the $i$ -th VEM trial function. The first six polynomials, for which the strain $\bm{\varepsilon}\left(\bm{q}\right)=0$ , correspond to the rigid body motions. Therefore, the following invertible augmented system is considered [35]:

\left\{\begin{aligned} &\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}\right)\,\mathrm{d}E=\int_{E}\bm{\varepsilon}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}\right)\,\mathrm{d}E&&\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)\\ &\frac{1}{n_{v}}\sum_{j=1}^{5n_{v}}\text{dof}_{j}\left(\bm{\psi}_{i}\right)\text{dof}_{j}\left(\bm{q}_{\alpha}\right)=\frac{1}{n_{v}}\sum_{j=1}^{5n_{v}}\text{dof}_{j}\left(\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i}\right)\text{dof}_{j}\left(\bm{q}_{\alpha}\right)&&\forall\alpha=1,\dots,6,\end{aligned}\right.

(66)

for all $i=1,\dots,N_{d}\big|_{E}$ . Here, $\text{dof}_{j}\left(\bm{q}_{\alpha}\right)$ represents the value of the polynomial at the vertex associated with the degree of freedom $j$ . The system in Eq. (66) can be expressed in matrix form as [35]:

\bm{B}^{\mathcal{K}}=\tilde{\bm{G}}^{\mathcal{K}}\tilde{\bm{\Pi}}^{\mathcal{K}}_{p}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{\mathcal{K}}_{p}=\left(\tilde{\bm{G}}^{\mathcal{K}}\right)^{-1}\bm{B}^{\mathcal{K}}.

(67)

Matrices $\bm{B}^{\mathcal{K}}$ and $\tilde{\bm{G}}^{\mathcal{K}}$ have dimensions $N_{p}\big|_{E}\times N_{d}\big|_{E}$ and $N_{p}\big|_{E}\times N_{p}\big|_{E}$ , respectively, and are defined as:

\displaystyle\bm{B}^{\mathcal{K}}_{\alpha,i}=\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E,\qquad\tilde{\bm{G}}^{\mathcal{K}}_{\alpha,\beta}=\int_{E}\bm{\varepsilon}\left(\bm{q}_{\beta}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E.

(68)

Notice that, according to the second line of Eq. (66), additional rows are included in $\tilde{\bm{G}}^{\mathcal{K}}$ and $\bm{B}^{\mathcal{K}}$ to account for the rigid body motions. The same applies to the subsequent matrices and vectors that require the augmented rows. To construct these matrices, it is convenient first to define polynomial vectors. Let $q_{n}\in\mathcal{P}_{p}\left(E\right)$ with $n=1,\dots,\frac{N_{p}\big|_{E}}{5}$ denote the scalar polynomial basis. The corresponding vector is defined as:

\bm{q}_{\alpha=5\left(n-1\right)+l}=\bm{e}_{l}q_{n},\quad\forall l=1,\dots,5,\quad\forall n=1,\dots,\frac{N_{p}\big|_{E}}{5},\quad\forall\alpha=1,\dots,N_{p}\big|_{E},

(69)

where $\bm{e}_{l}$ is the $l$ -th versor in $\mathbb{R}^{5}$ . The construction of matrix $\tilde{\bm{G}}^{\mathcal{K}}_{\alpha,\beta}$ in Eq. (68) is straightforward, as it involves only polynomial terms. For $\alpha=1,\dots,6$ , the rigid body motions defined in Eq. (66) must also be included. The construction of $\bm{B}^{\mathcal{K}}$ is more involved because the VEM trial functions are not explicitly known inside the element. First, the strain operator $\bm{\varepsilon}\left(\cdot\right)$ can be decomposed as:

\displaystyle\bm{\varepsilon}\left(\cdot\right)

\displaystyle=\begin{bmatrix},x&0&0&0&0\\ 0&,y&0&0&0\\ ,y&,x&0&0&0\\ 0&0&0&,x&0\\ 0&0&0&0&,y\\ 0&0&0&,y&,x\\ 0&0&,x&1&0\\ 0&0&,y&0&1\end{bmatrix}=\begin{bmatrix},x&0&0&0&0\\ 0&,y&0&0&0\\ ,y&,x&0&0&0\\ 0&0&0&,x&0\\ 0&0&0&0&,y\\ 0&0&0&,y&,x\\ 0&0&,x&0&0\\ 0&0&,y&0&0\end{bmatrix}+\begin{bmatrix}0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&1&0\\ 0&0&0&0&1\end{bmatrix}=\bm{\varepsilon}^{\mathrm{\partial}}\left(\cdot\right)+\bm{\varepsilon}^{\mathrm{\theta}}\left(\cdot\right).

(70)

Using this decomposition, the matrix $\bm{B}^{\mathcal{K}}$ can be written as:

	$\displaystyle\bm{B}^{\mathcal{K}}_{\alpha,i}=\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E$	$\displaystyle=\int_{E}\left[\bm{\varepsilon}^{\mathrm{\partial}}\left(\bm{\psi}_{i}\right)+\bm{\varepsilon}^{\mathrm{\theta}}\left(\bm{\psi}_{i}\right)\right]^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E$		(71)
		$\displaystyle=\int_{E}\bm{\varepsilon}^{\mathrm{\partial}}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E+\int_{E}\bm{\varepsilon}^{\mathrm{\theta}}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E.$		(71)

Integration by parts can now be applied to the first term, obtaining:

\displaystyle\bm{B}^{\mathcal{K}}_{\alpha,i}

\displaystyle=-\int_{E}\bm{\psi}_{i}\cdot\bm{L}_{d}\left[\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E}+\int_{E}\bm{\varepsilon}^{\mathrm{\theta}}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E,

(72)

where $\bm{\hat{n}}_{\Gamma^{E}}$ is the outward unit normal on the element boundary, and $\bm{L}_{d}\left[\cdot\right]$ and $\hat{\bm{S}}$ are defined as:

\displaystyle\bm{L}_{d}\left[\cdot\right]=\begin{bmatrix},x&0&,y&0&0&0&0&0\\ 0&,y&,x&0&0&0&0&0\\ 0&0&0&0&0&0&,x&,y\\ 0&0&0&,x&0&,y&0&0\\ 0&0&0&0&,y&,x&0&0\end{bmatrix},\qquad\hat{\bm{S}}=\begin{bmatrix}N_{xx}&N_{xy}\\ N_{xy}&N_{yy}\\ Q_{xx}&Q_{yy}\\ M_{xx}&M_{xy}\\ M_{xy}&M_{yy}\end{bmatrix}.

(73)

The three terms defined in Eq. (72) can be computed entirely from the degrees of freedom. Starting from the first term and approximating $\mathbb{C}$ as constant within the element, with a value equal to its average at the integration points, it follows that:

\left\{\begin{aligned} &\bm{L}_{d}\left[\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\right]\in\bm{\mathcal{P}}_{p-1}\left(E\right)&&\forall\alpha=3+5r,\;\;\text{with}\;\;r=0,\dots,\frac{N_{p}\big|_{E}-5}{5}\\ &\bm{L}_{d}\left[\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\right]\in\bm{\mathcal{P}}_{p-2}\left(E\right)&&\text{otherwise}.\end{aligned}\right.

(74)

This reduces to an integral of a VEM trial function multiplied by a polynomial of degree $p-1$ for $w$ , and of degree $p-2$ for the other displacement components. Thus, the first term depends only on the internal degrees of freedom of $\bm{\mathcal{V}}_{h}$ . Following [35], the coefficients can be expressed in terms of the polynomial basis as:

\left\{\begin{aligned} &\bm{L}_{d}\left[\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\right]=\sum_{\beta=1}^{N_{p}\big|_{E}-1}d_{\alpha,\beta}\bm{q}_{\beta}&&\forall\alpha=3+5r,\;\;\text{with}\;\;r=0,\dots,\frac{N_{p}\big|_{E}-5}{5}\\ &\bm{L}_{d}\left[\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\right]=\sum_{\beta=1}^{N_{p}\big|_{E}-2}d_{\alpha,\beta}\bm{q}_{\beta}&&\text{otherwise}.\end{aligned}\right.

(75)

For a scaled monomial basis, these coefficients are trivial to obtain. In the case of the orthonormal basis adopted in this work, the coefficients are computed using the orthonormality property [34]:

d_{\alpha,\beta}=\int_{E}\bm{L}_{d}\left[\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\right]\cdot\bm{q}_{\beta}\,\mathrm{d}E.

(76)

Regarding the third term, it holds that:

\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\in\bm{\mathcal{P}}_{p}\left(E\right)\quad\forall\alpha=4+s+5r,\;\;\text{with}\;\;r=0,\dots,\frac{N_{p}\big|_{E}-5}{5},\;\;s=0,1.

(77)

Therefore, this matrix also depends solely on the internal degrees of freedom, as it corresponds to the integral of a VEM trial function times a polynomial of degree $p$ for $\theta_{x}$ and $\theta_{y}$ . The coefficients can be written as:

\mathbb{C}\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)=\sum_{\beta=1}^{N_{p}\big|_{E}}d_{\alpha,\beta}\bm{q}_{\beta}\quad\forall\alpha=4+s+5r,\;\;\text{with}\;\;r=0,\dots,\frac{N_{p}\big|_{E}-5}{5},\;\;s=0,1.

(78)

This step is straightforward, and the coefficients are obtained in the same way as in Eq. (75). Lastly, the boundary integrals can be easily evaluated because the trial functions are known along the element boundary:

\displaystyle\int_{\Gamma^{E}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E}

\displaystyle=\sum_{j=1}^{n_{e}}\int_{\Gamma^{E}_{e}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\bm{\hat{n}}_{\Gamma_{e}}\,\mathrm{d}\Gamma^{E}_{e}+\sum_{j=1}^{\tilde{n}_{e}}\int_{\tilde{\Gamma}^{E}_{e}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\bm{\hat{n}}_{\tilde{\Gamma}_{e}}\,\mathrm{d}\tilde{\Gamma}^{E}_{e}.

(79)

For straight edges, the integrals are computed using the Gauss-Lobatto quadrature:

\sum_{j=1}^{n_{e}}\int_{\Gamma^{E}_{e}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\bm{\hat{n}}_{\Gamma_{e}}\,\mathrm{d}\Gamma^{E}_{e}=\sum_{j=1}^{n_{e}}\big|e_{j}\big|\sum_{r=1}^{n_{int}}w_{r}\bm{\psi}_{i}\big|_{x_{r},y_{r}}\cdot\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\big|_{x_{r},y_{r}}\bm{\hat{n}}_{\Gamma_{e}},

(80)

where $\big|e_{j}\big|$ is the edge length, $w_{r}$ and $x_{r},y_{r}$ are the quadrature weights and points, and $n_{int}$ is the number of integration points required for exact integration. For curved edges, following [12], the mapping is exploited:

$\displaystyle\sum_{j=1}^{\tilde{n}_{e}}\int_{\tilde{\Gamma}^{E}_{e}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\bm{\hat{n}}_{\tilde{\Gamma}_{e}}\,\mathrm{d}\tilde{\Gamma}^{E}_{e}$	$\displaystyle=\sum_{j=1}^{\tilde{n}_{e}}\int_{\tilde{\Gamma}^{E}_{e}}\left[\tilde{\bm{\psi}}_{i}\cdot\bm{\tilde{\hat{S}}\left(\bm{q}_{\alpha}\right)}\tilde{\bm{\hat{n}}}_{\tilde{\Gamma}_{e}}\right]\circ\gamma_{e}^{-1}\,\mathrm{d}\tilde{\Gamma}^{E}_{e}$	(81)
	$\displaystyle=\sum_{j=1}^{\tilde{n}_{e}}\int_{I_{e}}\left[\tilde{\bm{\psi}}_{i}\cdot\bm{\tilde{\hat{S}}\left(\bm{q}_{\alpha}\right)}\tilde{\bm{\hat{n}}}_{\tilde{\Gamma}_{e}}\right]\mathinner{\!\left\lVert\gamma_{e}^{{}^{\prime}}\right\rVert}\,\mathrm{d}I_{e}$
	$\displaystyle=\sum_{j=1}^{\tilde{n}_{e}}\sum_{r=1}^{n_{int}}w_{r}\left[\tilde{\bm{\psi}}_{i}\cdot\bm{\tilde{\hat{S}}\left(\bm{q}_{\alpha}\right)}\tilde{\bm{\hat{n}}}_{\tilde{\Gamma}_{e}}\right]\big\|_{x_{r},y_{r}}\mathinner{\!\left\lVert\gamma_{e}^{{}^{\prime}}\right\rVert}\big\|_{x_{r},y_{r}}$
	$\displaystyle=\sum_{j=1}^{\tilde{n}_{e}}\sum_{r=1}^{n_{int}}w_{r}\tilde{\bm{\psi}}_{i}\big\|_{x_{r},y_{r}}\cdot\left[\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\bm{\hat{n}}_{\tilde{\Gamma}_{e}}\right]\big\|_{\gamma_{e}\left(x_{r},y_{r}\right)}\mathinner{\!\left\lVert\gamma_{e}^{{}^{\prime}}\right\rVert}\big\|_{x_{r},y_{r}}.$

Here, $\tilde{\bm{\psi}}_{i}$ is the image of $\bm{\psi}_{i}$ through $\gamma_{e}$ and is a polynomial, see Eq. (19); specifically, Lagrange polynomials are employed. Therefore, matrix $\bm{B}^{\mathcal{K}}_{\alpha,i}$ is fully computable. For $\alpha=1,\dots,6$ , the rigid body motions defined in Eq. (66) must also be included. It is convenient to also define the matrix:

\bm{D}_{i,\alpha}=\text{dof}_{i}\left(\bm{q}_{\alpha}\right).

(82)

For the boundary degrees of freedom, this corresponds to evaluating the polynomial at the given node, while for internal moments it reduces to integrating the product of known polynomials. Numerical integration over a generic polygon with curved Bézier edges is carried out using the open-access code from [45, 44].

By defining the same operator as in Eq. (64), this time expressed in terms of the trial functions themselves [6], the following expression is obtained:

\displaystyle\bm{\Pi}^{\mathcal{K}}_{p}\bm{\psi}_{i}=\sum_{j=1}^{N_{d}\big|_{E}}\pi_{i,j}\bm{\psi}_{j}\quad\forall i=1,\dots,N_{d}\big|_{E},

\displaystyle\text{with}\quad\pi_{i,j}=\sum_{\beta=1}^{N_{p}\big|_{E}}s_{i,\beta}\text{dof}_{j}\left(\bm{q}_{\beta}\right).

(83)

In matrix form:

\bm{\Pi}=\bm{D}\tilde{\bm{\Pi}}^{\mathcal{K}}_{p}.

(84)

The consistency part of the stiffness matrix is constructed as:

\bm{K}_{c}^{E}=\left(\tilde{\bm{\Pi}}^{\mathcal{K}}_{p}\right)^{T}\bm{G}^{{\mathcal{K}_{\mathrm{VC}}}}\tilde{\bm{\Pi}}^{\mathcal{K}}_{p}.

(85)

Notice that, to construct the stiffness matrix, $\bm{G}^{{\mathcal{K}_{\mathrm{VC}}}}$ accounts for the variable stiffness:

\bm{G}^{{\mathcal{K}_{\mathrm{VC}}}}_{\alpha,\beta}=\int_{E}\bm{\varepsilon}\left(\bm{q}_{\beta}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E.

(86)

To evaluate Eq. (86), for a spatially varying $\mathbb{C}\left(x,y\right)$ , a higher-order quadrature rule is required compared to the constant case.

Stabilized VEM

In matrix form, Eq. (38) can be written as:

\bm{S}^{E}=\begin{bmatrix}{\bm{S}^{E}}^{u,v}&\bm{0}&\bm{0}\\ \bm{0}&{\bm{S}^{E}}^{w}&\bm{0}\\ \bm{0}&\bm{0}&{\bm{S}^{E}}^{\theta_{x},\theta_{y}}\end{bmatrix},\qquad\text{with}\quad{\bm{S}^{E}}^{\alpha}=\tau\left(\bm{I}^{\alpha}-\bm{\Pi}^{\alpha}\right)^{T}{\bm{s}^{E}}^{\alpha}\left(\bm{I}^{\alpha}-\bm{\Pi}^{\alpha}\right),

(87)

where $\bm{\Pi}^{\alpha}$ and ${\bm{s}^{E}}^{\alpha}$ are the sub-matrices associated to the degrees of freedom of the sub-block $\alpha$ , see Eq. 38.

So that, in the case of stabilized VEM, the stiffness matrix can be written as:

\bm{K}^{E}=\bm{K}_{c}^{E}+\bm{S}^{E}.

(88)

Self-stabilized VEM

The VEM space has the same dimension as the classical stabilized formulation, i.e. $N_{d+\ell_{E}}\big|_{E}=N_{d}\big|_{E}$ . Conversely, the dimension of the enlarged polynomial space is given by:

N_{p+\ell_{E}}\big|_{E}=\text{dim}\left(\bm{\mathcal{P}}^{\star}_{p+\ell_{E}-1}\left(E\right)\right)=8\frac{\left(p+\ell_{E}\right)\left(p+\ell_{E}+1\right)}{2}.

(89)

The projection in the $L^{2}$ norm is obtained by solving the system:

\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{q}^{\star}\,\mathrm{d}E=\int_{E}\bm{\Pi}^{\mathcal{K}}_{p+\ell_{E}-1}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{q}^{\star}\,\mathrm{d}E\quad\forall\bm{q}^{\star}\in\bm{\mathcal{P}}^{\star}_{p+\ell_{E}-1}\left(E\right),

(90)

for all $i=1,\dots,N_{d}\big|_{E}$ , where $\bm{q}^{\star}$ is the polynomial vector defined as:

\bm{q}^{\star}_{\alpha=8\left(n-1\right)+l}=\bm{e}_{l}q_{n},\quad\forall l=1,\dots,8,\quad\forall n=1,\dots,\frac{N_{p}\big|_{E}}{5},\quad\alpha=1,\dots,N_{p+\ell_{E}}\big|_{E},

(91)

with $\bm{e}_{l}$ being the $l-th$ versor in $\mathbb{R}^{8}$ and $\bm{\mathcal{P}}^{\star}_{p+\ell_{E}-1}\left(E\right)=\left[\mathcal{P}_{p+\ell_{E}-1}\left(E\right)\right]^{8}$ . The projector $\bm{\Pi}^{\mathcal{K}}_{p+\ell_{E}-1}$ maps the strain $\bm{\varepsilon}\left(\bm{\psi}_{i}\right)$ , rather than the trial function $\bm{\psi}_{i}$ itself, into this new polynomial space.

Eq. (90) in matrix form becomes:

\bm{B}^{\mathcal{K}\star}=\bm{G}^{\mathcal{K}\star}\tilde{\bm{\Pi}}^{\mathcal{K}}_{p+\ell_{E}-1}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{\mathcal{K}}_{p+\ell_{E}-1}=\left(\bm{G}^{\mathcal{K}\star}\right)^{-1}\bm{B}^{\mathcal{K}\star}.

(92)

The matrix $\bm{B}^{\mathcal{K}\star}$ has dimension $N_{p+\ell_{E}}\big|_{E}\times N_{d}\big|_{E}$ and the matrix $\bm{G}^{\mathcal{K}\star}$ has dimension $N_{p+\ell_{E}}\big|_{E}\times N_{p+\ell_{E}}\big|_{E}$ . The projection operator $\tilde{\bm{\Pi}}^{\mathcal{K}}_{p+\ell_{E}-1}$ has dimension $N_{p+\ell_{E}}\big|_{E}\times N_{d}\big|_{E}$ . The matrices $\bm{B}^{\mathcal{K}\star}$ and $\bm{G}^{\mathcal{K}\star}$ are defined as:

\displaystyle\bm{B}^{\mathcal{K}\star}_{\alpha,i}=\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{q}^{\star}_{\alpha}\,\mathrm{d}E,\qquad\bm{G}^{\mathcal{K}\star}_{\alpha,\beta}=\int_{E}{\bm{q}^{\star}_{\beta}}^{T}\mathbb{C}\bm{q}^{\star}_{\alpha}\,\mathrm{d}E.

(93)

These matrices are constructed following the same procedure as in the stabilized VEM. After applying integration by parts, to construct matrix $\bm{B}^{\mathcal{K}\star}_{\alpha,i}$ , the enlarged enhanced space in Eq. (40) is used, giving:

\displaystyle\bm{B}^{\mathcal{K}\star}_{\alpha,i}

\displaystyle=-\int_{E}\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\cdot\bm{L}_{d}\left[\mathbb{C}\bm{q}^{\star}_{\alpha}\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}^{\star}_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E}+\int_{E}\bm{\varepsilon}^{\mathrm{\theta}}\left(\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\right)^{T}\mathbb{C}\bm{q}^{\star}_{\alpha}\,\mathrm{d}E,

(94)

where $\hat{\bm{S}}\left(\bm{q}^{\star}_{\alpha}\right)$ is constructed starting from $\mathbb{C}\bm{q}^{\star}_{\alpha}$ . The operator $\bm{\Pi}_{p}^{0}$ is the $L^{2}$ projector and is defined later in Eq. (107). The stiffness matrix is then obtained in the form:

\bm{K}^{E}=\left(\tilde{\bm{\Pi}}^{\mathcal{K}}_{p+\ell_{E}-1}\right)^{T}\bm{G}^{{\mathcal{K}_{\mathrm{VC}}\star}}\tilde{\bm{\Pi}}^{\mathcal{K}}_{p+\ell_{E}-1},

(95)

where $\bm{G}^{{\mathcal{K}_{\mathrm{VC}}\star}}$ is the variable coefficient counterpart of $\bm{G}^{\mathcal{K}\star}$ .

Geometric stiffness matrix

Following the same logical flow adopted for the stiffness matrix, the following system is obtained:

\left\{\begin{aligned} &\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(\psi_{i}\right)^{T}\mathbb{N}\bm{\varepsilon}_{\mathrm{b}}\left(q\right)\,\mathrm{d}E=\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(\bm{\Pi}^{\mathcal{G}}_{p}\psi_{i}\right)^{T}\mathbb{N}\bm{\varepsilon}_{\mathrm{b}}\left(q\right)\,\mathrm{d}E&&\forall q\in\mathcal{P}_{p}\left(E\right)\\ &\frac{1}{n_{v}}\sum_{j=1}^{n_{v}}\text{dof}_{j}\left(\psi_{i}\right)\text{dof}_{j}\left(q_{\alpha}\right)=\frac{1}{n_{v}}\sum_{j=1}^{n_{v}}\text{dof}_{j}\left(\bm{\Pi}^{\mathcal{G}}_{p}\psi_{i}\right)\text{dof}_{j}\left(q_{\alpha}\right)&&\alpha=1,\end{aligned}\right.

(96)

for all $i=1,\dots,N_{d}\big|_{E}^{w}$ . Since only the terms related to $w$ are retained, the rigid body motion corresponds to $\alpha=1$ and $\psi_{i}$ and $q$ are both scalars. The strain operator is defined as $\bm{\varepsilon}_{\mathrm{b}}=\begin{bmatrix},x,\quad,y\end{bmatrix}^{T}$ , corresponding to the scalar version of Eq. (13). The system in Eq. (96) can be written in matrix form as:

\bm{B}^{\mathcal{G}}=\tilde{\bm{G}}^{\mathcal{G}}\tilde{\bm{\Pi}}^{\mathcal{G}}_{p}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{\mathcal{G}}_{p}=\left(\tilde{\bm{G}}^{\mathcal{G}}\right)^{-1}\bm{B}^{\mathcal{G}}.

(97)

Matrix $\tilde{\bm{G}}^{\mathcal{G}}$ has dimension $\frac{N_{p}\big|_{E}}{5}\times\frac{N_{p}\big|_{E}}{5}$ , while matrix $\bm{B}^{\mathcal{G}}$ has dimension $\frac{N_{p}\big|_{E}}{5}\times N_{d}\big|_{E}^{w}$ , and they read:

\displaystyle\bm{B}^{\mathcal{G}}_{\alpha,i}=\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(\psi_{i}\right)^{T}\mathbb{N}\bm{\varepsilon}_{\mathrm{b}}\left(q_{\alpha}\right)\,\mathrm{d}E,

\displaystyle\tilde{\bm{G}}^{\mathcal{G}}_{\alpha,\beta}=\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(q_{\beta}\right)^{T}\mathbb{N}\bm{\varepsilon}_{\mathrm{b}}\left(q_{\alpha}\right)\,\mathrm{d}E.

(98)

In this context, $q_{\alpha}\in\mathcal{P}_{p}\left(E\right)$ for $\alpha=1,\dots,\frac{N_{p}\big|_{E}}{5}$ . Regarding matrix $\bm{B}^{\mathcal{G}}$ , integration by parts yields:

\bm{B}^{\mathcal{G}}_{\alpha,i}=-\int_{E}\psi_{i}\cdot\bm{L}_{b}\left[\mathbb{N}\bm{\varepsilon}_{\mathrm{b}}\left(q_{\alpha}\right)\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\psi_{i}\cdot\mathbb{N}\bm{\varepsilon}_{\mathrm{b}}\left(q_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E},

(99)

where $\bm{L}_{b}\left[\cdot\right]=\begin{bmatrix},x&,y\end{bmatrix}$ . The steps to compute Eq. (99) follow the same procedure detailed for the stiffness matrix, requiring only the degrees of freedom. The construction of matrix $\tilde{\bm{G}}^{\mathcal{G}}$ is straightforward, allowing the projection in Eq. (97) to be computed. Therefore, the matrix is constructed as:

\bm{K}_{g}^{E}=\left(\tilde{\bm{\Pi}}^{\mathcal{G}}_{p}\right)^{T}\bm{G}^{{\mathcal{G}_{\mathrm{VC}}}}\tilde{\bm{\Pi}}^{\mathcal{G}}_{p},

(100)

with:

\displaystyle\bm{G}^{{\mathcal{G}_{\mathrm{VC}}}}_{\alpha,\beta}=\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(q_{\beta}\right)^{T}\mathbb{N}\left(x,y\right)\bm{\varepsilon}_{\mathrm{b}}\left(q_{\alpha}\right)\,\mathrm{d}E.

(101)

Mass matrix

The steps for the construction of the mass matrix are detailed hereafter. The projection is of $L^{2}$ type and its extended form reads:

\int_{E}\bm{\psi}_{i}^{T}\mathbb{M}\bm{q}\,\mathrm{d}E=\int_{E}\bm{\Pi}^{\mathcal{M}}_{p}\bm{\psi}_{i}^{T}\mathbb{M}\bm{q}\,\mathrm{d}E\quad\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right),

(102)

for all $i=1,\dots,N_{d}\big|_{E}$ . In matrix form:

\bm{Q}^{\mathcal{M}}=\bm{H}^{\mathcal{M}}\tilde{\bm{\Pi}}^{\mathcal{M}}_{p}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{\mathcal{M}}_{p}=\left(\bm{H}^{\mathcal{M}}\right)^{-1}\bm{Q}^{\mathcal{M}},

(103)

where matrix $\bm{H}^{\mathcal{M}}$ has dimension $N_{p}\big|_{E}\times N_{p}\big|_{E}$ and matrix $\bm{Q}^{\mathcal{M}}$ has dimension $N_{p}\big|_{E}\times N_{d}\big|_{E}$ , with entries:

\displaystyle\bm{Q}^{\mathcal{M}}_{\alpha,i}=\int_{E}\bm{\psi}_{i}^{T}\mathbb{M}\bm{q}_{\alpha}\,\mathrm{d}E,\qquad\bm{H}^{\mathcal{M}}_{\alpha,\beta}=\int_{E}\bm{q}_{\beta}^{T}\mathbb{M}\bm{q}_{\alpha}\,\mathrm{d}E.

(104)

By using the enhanced space, matrix $\bm{Q}^{\mathcal{M}}$ is constructed as follows:

\bm{Q}^{\mathcal{M}}_{\alpha,i}=\begin{cases}\int_{E}\bm{\psi}_{i}^{T}\mathbb{M}\bm{q}_{\alpha}\,\mathrm{d}E\qquad\;\,\forall\alpha=1+s+5r,\;\text{with}\;r=0,\dots,\frac{N_{p-2}\big|_{E}-5}{5},\;s=0,1\\ \int_{E}\bm{\psi}_{i}^{T}\mathbb{M}\bm{q}_{\alpha}\,\mathrm{d}E\qquad\;\,\forall\alpha=3+5r,\;\text{with}\;r=0,\dots,\frac{N_{p-1}\big|_{E}-5}{5}\\ \int_{E}\bm{\psi}_{i}^{T}\mathbb{M}\bm{q}_{\alpha}\,\mathrm{d}E\qquad\;\,\forall\alpha=4+s+5r,\;\text{with}\;r=0,\dots,\frac{N_{p}\big|_{E}-5}{5},\;s=0,1\\ \int_{E}\tilde{\bm{\Pi}}^{\mathcal{K}}_{p}\bm{\psi}_{i}^{T}\mathbb{M}\bm{q}_{\alpha}\,\mathrm{d}E\quad\text{otherwise}.\end{cases}

(105)

The first three rows correspond to the internal degrees of freedom and are computed exactly. The remaining terms are evaluated using the enhanced condition defined in Eqs. (18) and (20). The construction of the matrix $\bm{H}^{\mathcal{M}}$ is straightforward. The mass matrix is then obtained as:

\bm{M}^{E}=\left(\tilde{\bm{\Pi}}^{\mathcal{M}}_{p}\right)^{T}\bm{H}^{\mathcal{M}}\tilde{\bm{\Pi}}^{\mathcal{M}}_{p}.

(106)

Body forces vector

The $L^{2}$ projection for the body forces vector reads:

\int_{E}\bm{\psi}_{i}^{T}\bm{q}\,\mathrm{d}E=\int_{E}\bm{\Pi}_{p}^{0}\bm{\psi}_{i}^{T}\bm{q}\,\mathrm{d}E\quad\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right),

(107)

for all $i=1,\dots,N_{d}\big|_{E}$ . In matrix form:

\bm{Q}_{p}^{0}=\bm{H}_{p}^{0}\tilde{\bm{\Pi}}_{p}^{0}\quad\longrightarrow\quad\tilde{\bm{\Pi}}_{p}^{0}=\left(\bm{H}_{p}^{0}\right)^{-1}\bm{Q}_{p}^{0},

(108)

where:

\displaystyle{\bm{Q}_{p}^{0}}_{\alpha,i}=\int_{E}\bm{\psi}_{i}^{T}\bm{q}_{\alpha}\,\mathrm{d}E,\qquad{\bm{H}_{p}^{0}}_{\alpha,\beta}=\int_{E}\bm{q}_{\beta}^{T}\bm{q}_{\alpha}\,\mathrm{d}E.

(109)

Matrix $\bm{Q}_{p}^{0}$ is straightforward as it is the product of known polynomials, whereas matrix $\bm{H}_{p}^{0}$ is built by using the enhanced condition in Eq. (20). Using this projection, the body forces vector can be computed as:

\left(\bm{f}_{b}^{E}\right)_{i}=\int_{E}\bm{\Pi}_{p}^{0}\bm{\psi}_{i}^{T}\bm{\bar{b}}\,\mathrm{d}E,

(110)

where $\bm{\bar{b}}$ are the body forces.

Thermal forces vector

The projection for the thermal forces vector can be written as:

\left\{\begin{aligned} &\int_{E}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\psi}^{*}_{i}\right)^{T}\check{\bm{R}}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}\right)\,\mathrm{d}E=\int_{E}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\Pi}^{\mathcal{T}}_{p}\bm{\psi}^{*}_{i}\right)^{T}\check{\bm{R}}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}\right)\,\mathrm{d}E&&\;\bm{q}^{*}\in\bm{\mathcal{P}}^{*}_{p}\left(E\right)\\ &\frac{1}{n_{v}}\sum_{j=1}^{n_{v}}\text{dof}_{j}\left(\bm{\psi}^{*}_{i}\right)\text{dof}_{j}\left(\bm{q}^{*}_{\alpha}\right)=\frac{1}{n_{v}}\sum_{j=1}^{n_{v}}\text{dof}_{j}\left(\bm{\Pi}^{\mathcal{T}}_{p}\bm{\psi}^{*}_{i}\right)\text{dof}_{j}\left(\bm{q}^{*}_{\alpha}\right)&&\;\alpha=1,\dots,5,\end{aligned}\right.

(111)

for all $i=1,\dots,N_{d}\big|_{E}^{*}$ , with $N_{d}\big|_{E}^{*}=N_{d}\big|_{E}^{u}+N_{d}\big|_{E}^{v}+N_{d}\big|_{E}^{\theta_{x}}+N_{d}\big|_{E}^{\theta_{y}}$ , and $\bm{\mathcal{P}}^{*}_{p}=\left[\mathcal{P}_{p}\right]^{4}$ . In $\bm{\psi}^{*}$ , the out-of-plane displacement $w$ is excluded, as the shear contribution is not accounted for in the thermal forces. Notice that in Table 2, $\hat{\bm{P}}^{*}\left(E\right)$ corresponds to $\begin{Bmatrix}x\;0\;0\;0\end{Bmatrix}^{T}$ . The strain operator $\bm{\varepsilon}_{\mathrm{th}}\left(\cdot\right)$ and matrix $\check{\bm{R}}$ are defined as:

\displaystyle\bm{\varepsilon}_{\mathrm{th}}\left(\cdot\right)=\begin{bmatrix},x&0&0&0\\ 0&,y&0&0\\ ,y&,x&0&0\\ 0&0&,x&0\\ 0&0&0&,y\\ 0&0&,y&,x\\ 0&0&0&0\\ 0&0&0&0\end{bmatrix},\qquad\check{\bm{R}}\left(x,y\right)=\begin{bmatrix}\check{\bm{N}}\left(x,y\right)&\bm{0}&\bm{0}\\ \bm{0}&\check{\bm{M}}\left(x,y\right)&\bm{0}\\ \bm{0}&\bm{0}&\check{\bm{Q}}\left(x,y\right)\end{bmatrix}.

(112)

The operator $\bm{\varepsilon}_{\mathrm{th}}\left(\cdot\right)$ corresponds to a reduced version of $\bm{\varepsilon}\left(\cdot\right)$ , in which shear contributions and the out-of-plane displacement $w$ are omitted. In $\check{\bm{R}}\left(x,y\right)$ :

\check{\bm{N}}\left(x,y\right)=\begin{bmatrix}\hat{N}_{xx}&\hat{N}_{xy}\\ \hat{N}_{xy}&\hat{N}_{yy}\end{bmatrix},\qquad\check{\bm{M}}\left(x,y\right)=\begin{bmatrix}\hat{M}_{xx}&\hat{M}_{xy}\\ \hat{M}_{xy}&\hat{M}_{yy}\end{bmatrix},\qquad\check{\bm{Q}}=\begin{bmatrix}0&0\\ 0&0\end{bmatrix}.

(113)

In matrix form, the projection is expressed as:

\bm{B}^{\mathcal{T}}=\tilde{\bm{G}}^{\mathcal{T}}\tilde{\bm{\Pi}}^{\mathcal{T}}_{p}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{\mathcal{T}}_{p}=\left(\tilde{\bm{G}}^{\mathcal{T}}\right)^{-1}\bm{B}^{\mathcal{T}}.

(114)

Matrices $\bm{B}^{\mathcal{T}}$ and $\tilde{\bm{G}}^{\mathcal{T}}$ have dimensions $\frac{4N_{p}\big|_{E}}{5}\times N_{d}\big|_{E}^{*}$ and $\frac{4N_{p}\big|_{E}}{5}\times\frac{4N_{p}\big|_{E}}{5}$ , respectively, and are defined as:

\displaystyle\bm{B}^{\mathcal{T}}_{\alpha,i}=\int_{E}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\psi}^{*}_{i}\right)^{T}\check{\bm{R}}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}_{\alpha}\right)\,\mathrm{d}E,\qquad\tilde{\bm{G}}^{\mathcal{T}}_{\alpha,\beta}=\int_{E}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}_{\beta}\right)^{T}\check{\bm{R}}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}_{\alpha}\right)\,\mathrm{d}E.

(115)

In this case, the polynomial space is defined as:

\bm{q}^{*}_{\alpha=4\left(n-1\right)+l}=\bm{e}_{l}q_{n},\quad\forall l=1,\dots,4,\quad\forall n=1,\dots,\frac{N_{p}\big|_{E}}{5},\quad\forall\alpha=1,\dots,\frac{4N_{p}\big|_{E}}{5},

(116)

where $\bm{e}_{l}$ is the $l-th$ versor in $\mathbb{R}^{4}$ . The matrix $\bm{B}^{\mathcal{T}}_{\alpha,i}$ is constructed following the same steps detailed for the stiffness matrix. Specifically, by applying integration by parts, it can be expressed as:

\bm{B}^{\mathcal{T}}_{\alpha,i}=-\int_{E}\bm{\psi}^{*}_{i}\cdot\bm{L}_{\mathrm{th}}\left[\check{\bm{R}}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}_{\alpha}\right)\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\bm{\psi}^{*}_{i}\cdot\tilde{\bm{R}}\left(\bm{q}^{*}_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E},

(117)

where the operator $\bm{L}_{\mathrm{th}}\left[\cdot\right]$ and the matrix $\tilde{\bm{R}}$ are defined as:

\displaystyle\bm{L}_{\mathrm{th}}\left[\cdot\right]=\begin{bmatrix},x&0&,y&0&0&0&0&0\\ 0&,y&,x&0&0&0&0&0\\ 0&0&0&,x&0&,y&0&0\\ 0&0&0&0&,y&,x&0&0\end{bmatrix},\qquad\tilde{\bm{R}}=\begin{bmatrix}\hat{N}_{xx}&\hat{N}_{xy}\\ \hat{N}_{xy}&\hat{N}_{yy}\\ \hat{M}_{xx}&\hat{M}_{xy}\\ \hat{M}_{xy}&\hat{M}_{yy}\end{bmatrix}.

(118)

Therefore, the matrix $\bm{B}^{\mathcal{T}}_{\alpha,i}$ is fully computable from the degrees of freedom. The matrix $\tilde{\bm{G}}^{\mathcal{T}}_{\alpha,\beta}$ is easily evaluated since it involves only polynomial components. The thermal load vector is then obtained as:

\left(\bm{f}_{\mathrm{th}}^{E}\right)_{i}=\int_{E}\check{\bm{R}}\left(x,y\right)\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\Pi}^{\mathcal{T}}_{p}{\bm{\psi}^{*}_{i}}^{T}\right)\,\mathrm{d}E.

(119)

Line loads and prescribed displacements

The construction of the line load vector follows the standard FEM procedure, as trial functions are polynomials along the element edges. The line load vector is expressed as:

\left(\bm{f}_{\mathrm{l}}^{E}\right)_{i}=\sum_{j}\int_{\Gamma^{E}_{j}}\bm{\psi}_{i}\bm{\bar{t}}_{j}\,\mathrm{d}\Gamma^{E}_{j},

(120)

where $\bm{\bar{t}}_{j}$ is the prescribed traction on the generic edge $\Gamma^{E}_{j}$ . The prescribed displacements are also imposed as in standard FEM.

Supporting material for Subsection 3.4: Matrix and vectors construction

Stiffness matrix

To obtain a projection that accounts for the spatial variation of the fiber orientation, the following strategy is proposed. The objective is to construct a projection operator $\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}$ such that:

\left\{\begin{aligned} &\left[\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}\right)\,\mathrm{d}E\right]^{\mathrm{VC}}=\int_{E}\bm{\varepsilon}\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}\right)\,\mathrm{d}E&&\forall\bm{q}\in\bm{\mathcal{P}}_{p}\left(E\right)\\ &\frac{1}{n_{v}}\sum_{j=1}^{5n_{v}}\text{dof}_{j}\left(\bm{\psi}_{i}\right)\text{dof}_{j}\left(\bm{q}_{\alpha}\right)=\frac{1}{n_{v}}\sum_{j=1}^{5n_{v}}\text{dof}_{j}\left(\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\bm{\psi}_{i}\right)\text{dof}_{j}\left(\bm{q}_{\alpha}\right)&&\forall\alpha=1,\dots,6,\end{aligned}\right.

(121)

for all $i=1,\dots,N_{d}\big|_{E}$ . This leads to the standard matrix form:

\bm{B}^{{\mathcal{K}_{\mathrm{VC}}}}=\tilde{\bm{G}}^{{\mathcal{K}_{\mathrm{VC}}}}\tilde{\bm{\Pi}}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{{\mathcal{K}_{\mathrm{VC}}}}_{p}=\left(\tilde{\bm{G}}^{{\mathcal{K}_{\mathrm{VC}}}}\right)^{-1}\bm{B}^{{\mathcal{K}_{\mathrm{VC}}}},

(122)

where matrix $\tilde{\bm{G}}^{{\mathcal{K}_{\mathrm{VC}}}}$ reads:

\displaystyle\tilde{\bm{G}}^{{\mathcal{K}_{\mathrm{VC}}}}_{\alpha,\beta}=\int_{E}\bm{\varepsilon}\left(\bm{q}_{\beta}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E,

(123)

and its computation only requires a sufficient number of integration points to accurately evaluate the non-constant constitutive law. The notation $[\cdot]^{\mathrm{VC}}$ used at the left-hand side of Eq. (121) refers to the $\mathrm{VC}$ -VEM technique and hides the manipulations described in Section 3.4. The left-hand side of Eq. (121) is thus defined as:

	$\displaystyle\bm{B}^{{\mathcal{K}_{\mathrm{VC}}}}_{\alpha,i}=$	$\displaystyle-\int_{E}\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\cdot\bm{L}_{d}\left[\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E}$		(124)
		$\displaystyle+\int_{E}\bm{\varepsilon}^{\mathrm{\theta}}\left(\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{\varepsilon}\left(\bm{q}_{\alpha}\right)\,\mathrm{d}E,$		(124)

where $\hat{\bm{S}}$ also contains the spatial variation of $\mathbb{C}\left(x,y\right)$ .

The subsequent steps for constructing the final stiffness matrix are identical to those of the standard VEM, with the stabilization term built using the standard projector.

Self-stabilized VEM

In the case of self-stabilized VEM, the corresponding projection is obtained:

\left[\int_{E}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{q}^{\star}\,\mathrm{d}E\right]^{\mathrm{VC}}=\int_{E}\bm{\Pi}^{{\mathcal{K}_{\mathrm{VC}}}}_{p+\ell_{E}-1}\bm{\varepsilon}\left(\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{q}^{\star}\,\mathrm{d}E\quad\forall\bm{q}^{\star}\in\bm{\mathcal{P}}^{\star}_{p+\ell_{E}-1}\left(E\right),

(125)

for all $i=1,\dots,N_{d}\big|_{E}$ . In matrix form:

\bm{B}^{{\mathcal{K}_{\mathrm{VC}}\star}}=\bm{G}^{{\mathcal{K}_{\mathrm{VC}}\star}}\tilde{\bm{\Pi}}^{{\mathcal{K}_{\mathrm{VC}}}}_{p+\ell_{E}-1}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{{\mathcal{K}_{\mathrm{VC}}}}_{p+\ell_{E}-1}=\left(\bm{G}^{{\mathcal{K}_{\mathrm{VC}}\star}}\right)^{-1}\bm{B}^{{\mathcal{K}_{\mathrm{VC}}\star}}.

(126)

The matrix $\bm{B}^{{\mathcal{K}_{\mathrm{VC}}\star}}$ is computed as:

	$\displaystyle\bm{B}^{{\mathcal{K}_{\mathrm{VC}}\star}}_{\alpha,i}$	$\displaystyle=-\int_{E}\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\cdot\bm{L}_{d}\left[\mathbb{C}\left(x,y\right)\bm{q}^{\star}_{\alpha}\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\bm{\psi}_{i}\cdot\hat{\bm{S}}\left(\bm{q}^{\star}_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E}+$		(127)
		$\displaystyle+\int_{E}\bm{\varepsilon}^{\mathrm{\theta}}\left(\bm{\Pi}_{p}^{0}\bm{\psi}_{i}\right)^{T}\mathbb{C}\left(x,y\right)\bm{q}^{\star}_{\alpha}\,\mathrm{d}E,$		(127)

where $\hat{\bm{S}}$ also contains the spatial variation of $\mathbb{C}\left(x,y\right)$ . The matrix $\bm{G}^{{\mathcal{K}_{\mathrm{VC}}\star}}$ is defined as:

\bm{G}^{{\mathcal{K}_{\mathrm{VC}}\star}}_{\alpha,\beta}=\int_{E}{\bm{q}^{\star}_{\beta}}^{T}\mathbb{C}\left(x,y\right)\bm{q}^{\star}_{\alpha}\,\mathrm{d}E.

(128)

Geometric stiffness matrix

The projection that accounts for the spatial variability of $\mathbb{N}\left(x,y\right)$ is defined such that:

\left\{\begin{aligned} &\left[\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(\psi_{i}\right)^{T}\mathbb{N}\left(x,y\right)\bm{\varepsilon}_{\mathrm{b}}\left(q\right)\,\mathrm{d}E\right]^{\mathrm{VC}}=\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(\bm{\Pi}^{{\mathcal{G}_{\mathrm{VC}}}}_{p}\psi_{i}\right)^{T}\mathbb{N}\left(x,y\right)\bm{\varepsilon}_{\mathrm{b}}\left(q\right)\,\mathrm{d}E&&\forall q\in\mathcal{P}_{p}\left(E\right)\\ &\frac{1}{n_{v}}\sum_{j=1}^{n_{v}}\text{dof}_{j}\left(\psi_{i}\right)\text{dof}_{j}\left(q_{\alpha}\right)=\frac{1}{n_{v}}\sum_{j=1}^{n_{v}}\text{dof}_{j}\left(\bm{\Pi}^{{\mathcal{G}_{\mathrm{VC}}}}_{p}\psi_{i}\right)\text{dof}_{j}\left(q_{\alpha}\right)&&\alpha=1,\end{aligned}\right.

(129)

for all $i=1,\dots,N_{d}\big|_{E}^{w}$ . In matrix form, this can be written as:

\bm{B}^{{\mathcal{G}_{\mathrm{VC}}}}=\tilde{\bm{G}}^{{\mathcal{G}_{\mathrm{VC}}}}\tilde{\bm{\Pi}}^{{\mathcal{G}_{\mathrm{VC}}}}_{p}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{{\mathcal{G}_{\mathrm{VC}}}}_{p}=\left(\tilde{\bm{G}}^{{\mathcal{G}_{\mathrm{VC}}}}\right)^{-1}\bm{B}^{{\mathcal{G}_{\mathrm{VC}}}},

(130)

with:

\displaystyle\tilde{\bm{G}}^{{\mathcal{G}_{\mathrm{VC}}}}_{\alpha,\beta}=\int_{E}\bm{\varepsilon}_{\mathrm{b}}\left(q_{\beta}\right)^{T}\mathbb{N}\left(x,y\right)\bm{\varepsilon}_{\mathrm{b}}\left(q_{\alpha}\right)\,\mathrm{d}E.

(131)

To handle the surface integrals involving unknown terms, the matrix $\bm{B}^{{\mathcal{G}_{\mathrm{VC}}}}_{\alpha,i}$ is written as:

\bm{B}^{{\mathcal{G}_{\mathrm{VC}}}}_{\alpha,i}=-\int_{E}\Pi^{0}_{p}\psi_{i}\cdot\bm{L}_{b}\left[\mathbb{N}\left(x,y\right)\bm{\varepsilon}_{\mathrm{b}}\left(q_{\alpha}\right)\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\psi_{i}\cdot\mathbb{N}\left(x,y\right)\bm{\varepsilon}_{\mathrm{b}}\left(q_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E}.

(132)

Thermal forces vector

For the thermal forces vector, the projection operator is defined such that:

\left\{\begin{aligned} &\left[\int_{E}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\psi}^{*}_{i}\right)^{T}\check{\bm{R}}\left(x,y\right)\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}\right)\,\mathrm{d}E\right]^{\mathrm{VC}}=\int_{E}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{\Pi}^{{\mathcal{T}_{\mathrm{VC}}}}_{p}\bm{\psi}^{*}_{i}\right)^{T}\check{\bm{R}}\left(x,y\right)\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}\right)\,\mathrm{d}E&&\forall\bm{q}^{*}\in\bm{\mathcal{P}}^{*}_{p}\left(E\right)\\ &\frac{1}{n_{v}}\sum_{j=1}^{n_{v}}\text{dof}_{j}\left(\bm{\psi}^{*}_{i}\right)\text{dof}_{j}\left(\bm{q}^{*}_{\alpha}\right)=\frac{1}{n_{v}}\sum_{j=1}^{n_{v}}\text{dof}_{j}\left(\bm{\Pi}^{{\mathcal{T}_{\mathrm{VC}}}}_{p}\bm{\psi}^{*}_{i}\right)\text{dof}_{j}\left(\bm{q}^{*}_{\alpha}\right)&&\forall\alpha=1,\dots,5,\end{aligned}\right.

(133)

for all $i=1,\dots,N_{d}\big|_{E}^{*}$ . This can be compactly rewritten as:

\bm{B}^{{\mathcal{T}_{\mathrm{VC}}}}=\tilde{\bm{G}}^{{\mathcal{T}_{\mathrm{VC}}}}\tilde{\bm{\Pi}}^{{\mathcal{T}_{\mathrm{VC}}}}_{p}\quad\longrightarrow\quad\tilde{\bm{\Pi}}^{{\mathcal{T}_{\mathrm{VC}}}}_{p}=\left(\tilde{\bm{G}}^{{\mathcal{T}_{\mathrm{VC}}}}\right)^{-1}\bm{B}^{{\mathcal{T}_{\mathrm{VC}}}}.

(134)

In order to compute $\tilde{\bm{\Pi}}^{{\mathcal{T}_{\mathrm{VC}}}}_{p}$ , matrix $\tilde{\bm{G}}^{{\mathcal{G}_{\mathrm{VC}}}}_{\alpha,\beta}$ is constructed as:

\displaystyle\tilde{\bm{G}}^{{\mathcal{G}_{\mathrm{VC}}}}_{\alpha,\beta}=\int_{E}\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}_{\beta}\right)^{T}\check{\bm{R}}\left(x,y\right)\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}_{\alpha}\right)\,\mathrm{d}E.

(135)

To make the expression of $\bm{B}^{{\mathcal{T}_{\mathrm{VC}}}}_{\alpha,i}$ computable, the matrix is rewritten as:

\bm{B}^{{\mathcal{T}_{\mathrm{VC}}}}_{\alpha,i}=-\int_{E}\bm{\Pi}_{p}^{0,*}\bm{\psi}^{*}_{i}\cdot\bm{L}_{\mathrm{th}}\left[\check{\bm{R}}\left(x,y\right)\bm{\varepsilon}_{\mathrm{th}}\left(\bm{q}^{*}_{\alpha}\right)\right]\,\mathrm{d}E+\int_{\Gamma^{E}}\bm{\psi}^{*}_{i}\cdot\tilde{\bm{R}}\left(\bm{q}^{*}_{\alpha}\right)\bm{\hat{n}}_{\Gamma^{E}}\,\mathrm{d}\Gamma^{E},

(136)

where $\tilde{\bm{R}}$ accounts for the dependency on $x,y$ and $\bm{\Pi}_{p}^{0,*}$ is $\bm{\Pi}_{p}^{0}$ referred to $\bm{\psi}^{*}$ .

		$\displaystyle N_{d}\big\|_{E}^{s}=\text{dim}\left(\mathcal{V}_{h}^{s}\left(E\right)\right)=pn_{v}+\frac{\left(p-1\right)p}{2}\ \text{for }s=u,v,$		(58)
		$\displaystyle N_{d}\big\|_{E}^{w}=\text{dim}\left(\mathcal{V}_{h}^{w}\left(E\right)\right)=pn_{v}+\frac{p\left(p+1\right)}{2},$
		$\displaystyle N_{d}\big\|_{E}^{s}=\text{dim}\left(\mathcal{V}_{h}^{s}\left(E\right)\right)=pn_{v}+\frac{\left(p+1\right)\left(p+2\right)}{2}\ \text{for }s=\theta_{x},\theta_{y}.$

A Comprehensive pp-VEM Framework for Advanced Variable Stiffness Plates with Arbitrary Shapes

Abstract

1 Introduction

2 Structural modeling

2.1 Variational statement

2.2 Plate model

2.2.1 Kinematics

2.2.2 Constitutive law

2.2.3 Energy terms

3 Virtual Element Method

3.1 Discretization

3.2 Virtual Element space and degrees of freedom

3.3 Discrete bilinear forms and projection operators

3.3.1 Stiffness matrix

Stabilized VEM

Self-stabilized VEM

3.3.2 Remaining forms

3.4 Discrete forms and projection operators with variable coefficients

4 Results

4.1 Test case 1

4.2 Test case 2

Uniform hh-refinement and pp-refinement

Local hh-refinement

Non-structured mesh

Remarks

4.3 Test case 3

Uniform hh-refinement and pp-refinement

Local hh-refinement

Remarks

4.4 Test case 4

Free-vibrations

Thermal buckling

4.5 Test case 5

5 Conclusions

References

Supporting material

Supporting material for Subsection 3.2: Spaces definition

Virtual Element space

Polynomial space

Supporting material for Subsection 3.3: Matrices and vectors construction

Stiffness matrix

Stabilized VEM

Self-stabilized VEM

Geometric stiffness matrix

Mass matrix

Body forces vector

Thermal forces vector

Line loads and prescribed displacements

Supporting material for Subsection 3.4: Matrix and vectors construction

Stiffness matrix

Self-stabilized VEM

Geometric stiffness matrix

Thermal forces vector

A Comprehensive $p$ -VEM Framework for Advanced
Variable Stiffness Plates with Arbitrary Shapes

Uniform $h$ -refinement and $p$ -refinement

Local $h$ -refinement

Uniform $h$ -refinement and $p$ -refinement

Local $h$ -refinement