Enhanced Multiphysics Simulation via a Tailored PINN Architecture for Piezoelectric Cantilever Beam Energy Harvesters

1. Introduction

The Finite Element Method (FEM) has long been the standard approach for modeling and simulating complex physical systems [1]. FEM is known for its accuracy and reliability, making it indispensable for engineering and applied physics, among other application fields [2]. However, it is computationally expensive, and as computational methods continue to evolve, researchers have begun exploring alternative approaches that can address emerging challenges in scientific applications and new problem domains, driving new scientific discoveries and engineering innovations [3,4].

In recent years, the accessibility and advancement of Machine Learning (ML) technologies have opened new doors for simulation methodologies [5]. Among them, Physics-Informed Neural Networks (PINNs) have been developed as a promising alternative [6]. Unlike traditional Artificial Neural Networks (ANN), PINNs do not require any measurement data for training [7]. Instead, they leverage the governing physical equations and the universal approximation theorem to simulate the behavior of specific domains or systems [8]. PINNs can approximate general solutions to the governing equations of physical systems, making them suited for nonlinearities and complex dynamic simulations [9].

The current state of research in ML-based simulations has seen significant advances, with key studies demonstrating the potential of PINNs in various applications, from fluid dynamics to structural analysis [10,11]. However, the complexities and nonlinearities of piezoelectric multiphysics simulations, particularly those involving the direct piezoelectric effect, present significant challenges [12].

PINNs have been applied to various problems involving the piezoelectric effect modeling semiconductor nanowires [13] and piezoelectric Energy Harvesters (EH) [14]. Such studies demonstrate that PINNs can accurately simulate physical fields and successfully estimate key material and system parameters related to the device. Their validity as a tool opens the door to their use in solving practical problems efficiently – e.g. by combining PINNs with experimental voltage measurements for efficient inverse modeling [14]. Elsewhere, PINNs have been successfully applied to: estimate piezoelectric EH in cantilevered metamaterial concrete structures subjected to wind excitation, incorporating complex interactions such as gravity, foundation effects, and aerodynamic forces [15]; thin-walled structural problems involving piezoelectric thin films, overcoming traditional meshing limitations and achieving accurate results even for extremely small thickness-to-length ratios [16]; and the accurate solution of multiphysics piezoelectric 2D problems via a modification of the PINN framework called the Deep Energy Method (DEM) [17].

Instead of simulating the more general piezoelectric governing equations, many of the aforementioned methods rely on Partial Differential Equation (PDE) formulations unique to specific piezoelectric devices, such as beams or thin films [14,15,16]. Moreover, the approach that addressed the fully general coupled PDE system employs DEM [17]. This paper specifically addresses the problem of effectively modeling the fully coupled piezoelectric PDE system within a PINN framework. The model is based on general coupled piezoelectric PDEs rather than beam-specific simplifications, even though the physical device under study is a 2D bimorph piezoelectric cantilever beam.

2. Materials and Methods

This section presents the main baseline used to develop the proposed Physics-Informed Neural Network (PINN) model. It begins with a brief description of the piezoelectric constitutive equations, which enable modeling a piezoelectric device. The characteristics of the simulated system are detailed, specifically, a 2D piezoelectric bimorph polyvinylidene fluoride (PVDF) beam. The proposed PINN model architecture is then presented, including the formulation of the neural network’s loss function, the training strategy employed, and the error metrics used to validate the model’s performance against a Finite Element Method (FEM) simulation. Finally, the implementation details and computational setup are presented to allow reproducibility of the results.

2.1. Piezoelectric Problem Statement

The piezoelectric effect refers to the coupling between mechanical and electrical behavior of certain crystalline materials. This phenomenon manifests in two effects [18]. First, the direct piezoelectric effect is the generation of an electrical displacement or voltage when mechanical stress is applied to the crystal. Second, the indirect piezoelectric effect describes the opposite behavior, which is a physical deformation of the crystal when an external electric field is applied.

Some crystals have a specific physical structure which causes the piezoelectric effect [19]. When they are under no mechanical strain, their internal electrical charges are balanced, which means that they do not produce any electric field. However, since these crystals belong to a non-centrosymmetric space group (i.e. their crystal unit structure lacks a center of symmetry), a reorganization of the electrical charge distribution occurs when a mechanical force is applied to the crystal, resulting in the generation of an electric field. Similarly, when an external electric field is applied, the same asymmetry causes rearrangements in the physical structure of the material, generating a mechanical deformation.

2.1.1. Piezoelectricity Constitutive Equations and Material Properties

The piezoelectric effect is described by a system of coupled constitutive equations which capture the relation between the electric and mechanical behavior [20]. These equations, expressed in their linear form, are as follows:

$$\overrightarrow{\sigma }=c^E\overrightarrow{S}-e^T\overrightarrow{E},$$

(1)

$$\overrightarrow{D}=e\overrightarrow{S}+{ϵ}^S\overrightarrow{E}.$$

(2)

where Equation (1) describes the indirect piezoelectric effect, in which an applied electric field produces mechanical stress in the piezoelectric material. The mechanical stress $\overrightarrow{\sigma }$ depends on a combination of the interaction between the material stiffness, represented by the tensor $c^E$, and the mechanical strain $\overrightarrow{S}$ ; and the interaction of the piezoelectric properties of the material, represented by the tensor $e^T,$and the electric field $\overrightarrow{E}$. Both tensors $e^T$ and $c^E$ contain constants that express the specific electromechanical properties of the crystal.

Equation (2) describes the direct piezoelectric effect, in which a dielectric displacement $\overrightarrow{D}$ results from mechanical deformation. The dielectric displacement $\overrightarrow{D}$ depends on the mechanical strain of the material $\overrightarrow{S}$ caused by the applied forces, mediated by the piezoelectric tensor $e$ and the interaction between the dielectric permittivity tensor ${ϵ}^S$ and the electric field $\overrightarrow{E}$.

In order to model piezoelectric devices, the constitutive equations are coupled with the governing equations of mechanical and electrostatic equilibrium in their time-independent form since, for EH applications, we are only interested in the steady-state solution. The resulting set of Partial Differential Equations (PDEs) is as follows [21]:

$$\nabla \cdot \sigma +\overrightarrow{f}=0,$$

(3)

$$\nabla \cdot \overrightarrow{D}=0.$$

(4)

Equation (3) is the mechanical equilibrium equation, in the static case. Here $\sigma$ represents the stress tensor and $f$ is the vector of applied body forces for the direct piezoelectric effect. In the case of the indirect effect, no forces are applied, therefore the values of the vector are zero. Equation (4) corresponds to Gauss’ law for the electric displacement $D$, with zero free charge density due to the insulation nature of piezoelectric materials.

Depending on the characteristics of each material, the coupled piezoelectric constitutive equations for the 2D case can be written in matrix form [22]. The material selected for this study is PVDF, a piezoelectric semi-crystalline polymer, used for several applications [23,24,25]. The electromechanical relations for this material in matrix form are as follows [26]:

$$\left(\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{12}\\ D_1\\ D_2\end{array}\right)=\left(\begin{array}{ccccc}{c}_{11}& c_{12}& 0& e_{11}& e_{31}\\ c_{12}& c_{22}& 0& e_{13}& e_{33}\\ 0& 0& c_{33}& e_{14}& e_{34}\\ e_{11}& e_{13}& e_{14}& -{ϵ}_1& 0\\ e_{31}& e_{33}& e_{34}& 0& -{ϵ}_2\end{array}\right)\left(\begin{array}{c}{S}_1\\ S_2\\ S_3\\ \begin{array}{c}{E}_1\\ E_2\end{array}\end{array}\right).$$

(5)

To complement and simplify the system of equations, kinematic relations linking strain to mechanical displacement and electric field to electric potential are required. The strain $\overrightarrow{S}$ can be expressed in terms of the derivatives of the Cartesian mechanical displacement $\overrightarrow{d}$, with components $u$ and $v$. The electric field can be expressed in terms of the derivatives of the electric potential $\phi$, as shown in Equations (6) and (7). These kinematic relations can be substituted into Equation (5), yielding a system of PDEs in terms of only the mechanical displacements and the electric potential:

$$\left(\begin{array}{c}{S}_1\\ S_2\\ S_3\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{\partial u}{\partial x}}\\ {\displaystyle \frac{\partial v}{\partial y}}\\ {\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\end{array}\right),$$

(6)

$$\left(\begin{array}{c}{E}_1\\ E_2\end{array}\right)=-\left(\begin{array}{c}{\displaystyle \frac{\partial \phi }{\partial x}}\\ {\displaystyle \frac{\partial \phi }{\partial y}}\end{array}\right).$$

(7)

2.1.2. Modeled Piezoelectric Device

The device to be modeled in the present study is a parallel piezoelectric bimorph cantilever beam in 2D. Figure 1 presents the characteristics and the general configuration for the simulations. The length of the beam is 0.1 meters, and its height is 0.001 meters. The lower end of the beam is connected to a grounded electrode while the upper layer is connected to a voltage electrode used to measure or induce electrical potential. For the direct piezoelectric effect simulation, an upward force of 1 N is applied at the free end of the beam.

In addition to the physical characteristics of the beam, specific boundary conditions need to be specified for a unique solution to exist. To complement the previously presented PDEs and to model the piezoelectric beam, the following set of boundary conditions [27] is applied:

$$\overrightarrow{d}\left(\overrightarrow{x}\right)=\overrightarrow{\overline{d}} , \forall \overrightarrow{x}\in {\mathsf{\Omega }}_u,$$

(8)

$$\phi \left(\overrightarrow{x}\right)=\overline{\phi }, \forall \overrightarrow{x}\in {\mathsf{\Omega }}_{\phi },$$

(9)

$$\sigma \left(\overrightarrow{x}\right)\cdot \overrightarrow{n}=0, \forall \overrightarrow{x}\in {\mathsf{\Omega }}_t,$$

(10)

$$\overrightarrow{D}\left(\overrightarrow{x}\right)\cdot \overrightarrow{n}=0, \forall \overrightarrow{x}\in {\mathsf{\Omega }}_p.$$

(11)

Equation (8) describes a Dirichlet boundary condition where the mechanical displacement $\overrightarrow{u}$ of the modeled device in segment ${\mathsf{\Omega }}_u$ has a prescribed value $\overrightarrow{\overline{u}}$. Similarly, Equation (9) describes a Dirichlet boundary condition that prescribes a specific value of the electric potential $\phi$ along the segment ${\mathsf{\Omega }}_{\phi }$. Equation (10) represents a Neumann boundary condition that enforces zero traction (zero applied stress) along segment ${\mathsf{\Omega }}_t$. Equation (11) is also a Neumann boundary condition that enforces zero normal dielectric displacement on boundary segment ${\mathsf{\Omega }}_p$.

The specific segments where the boundary conditions are applied and the prescribed values of the mechanical displacements and electric potential for the indirect and direct piezoelectric effect simulations are illustrated in Figure 2. For the simulation of both piezoelectric effects, the beam is clamped at its left end and is free at its right end. Additionally, both right and left ends are under traction-free conditions. The right end maintains zero normal dielectric displacement.

Figure 2a shows the configuration of the beam when modelling the indirect piezoelectric effect. The prescribed voltage of the upper edge of the beam is 100 V and of the lower edge is 0 V. Figure 2b illustrates the direct piezoelectric effect, where only the lower edge is grounded, while the voltage of the upper edge is measured, as shown in Figure 1.

Finally, in order to give a complete description of the piezoelectric beam, the values of the elastic, piezoelectric, and dielectric permittivity coefficients, appearing in Equations (5) and (6), are required. The values of the coefficients common to both layers can be found in Table 1. Table 2 lists the specific layer values of certain piezoelectric coefficients that differ because of the polarization orientations of the layers [17].

2.2. Physics-Informed Neural Networks

Physics-Informed Neural Networks (PINNs) provide a framework for solving PDEs by embedding physical laws into the loss function of a feed-forward Deep Neural Network (DNN) [6]. To this end, we consider the general form of a time-independent PDE and its boundary conditions:

$$\mathcal{N}\left[u\left(x\right)\right]=0, \hspace{1em}x \in \mathsf{\Omega },$$

(12)

$$B\left[u\left(x\right)\right]=0,\hspace{1em}x\in \partial \mathsf{\Omega }.$$

(13)

where $u\left(x\right)$ is the solution to the PDE, $\mathcal{N}$ is a differential operator corresponding to the governing equations of the problem, $\mathsf{\Omega }$ is the domain of interest of the problem, $B$ represents the boundary condition operator, $\partial \mathsf{\Omega }$ is the boundary where boundary conditions are enforced, and $x$ denotes the spatial coordinates.

The solution of the PDE $u\left(x\right)$ is approximated by a DNN $u_{\theta }\left(x\right)$, where $\theta$ denotes the set of trainable weights and biases of the DNN. To train the DNN in the PINN framework, a four-step procedure is followed [28]:

• Sampling of collocation points and boundary condition points

A training dataset is constructed by sampling N collocation points ${\left\{x_i\right\}}_{i=1}^N$ from the interior domain $\mathsf{\Omega }$ and $M$ boundary points ${\{x_j, y_j\}}_{j=1}^M$ from $\partial \mathsf{\Omega }$. The boundary condition points include target values $y_j$ representing the prescribed boundary values for the Dirichlet or Neumann boundary conditions.

2.

Residual computation

The core of the learning process of a PINN is the minimization of the norm of the PDE residual $r(x;\theta )$. The residual is defined as the output of the differential operator $\mathcal{N}$ to the neural network approximation $u_{\theta }\left(x\right)$:

$$r\left(x;\theta \right):=\mathcal{N}\left[u_{\theta }\left(x\right)\right].$$

(14)

The residual represents the extent to which the neural network fails to satisfy the governing PDEs. If $\left[u_{\theta }\left(x\right)\right]$ were the exact solution, the residual would be zero throughout $\mathsf{\Omega }$. To compute the residual, the neural network is evaluated at the collocation points and its partial derivatives with respect to the independent variables $x$ are calculated using automatic differentiation. These derivatives are then substituted into the differential operator $\mathcal{N}$. To ensure that $u_{\theta }\left(x\right)$ approximates the solution of the PDE, the Mean Squared Error (MSE) of the residual evaluated at the collocation points is minimized:

$${\mathcal{L}}_{PDE}(\theta )={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\Vert r\left(x_i;\theta \right)\Vert }^2.$$

(15)

3.

Boundary conditions enforcement

To obtain a unique solution, the boundary condition operator $B$ from Equation (13) must be satisfied. To this end, a collective boundary loss term ${\mathcal{L}}_{BC}\left(\theta \right)$ that penalizes deviations from Dirichlet and/or Neumann conditions prescribed values is defined:

$${\mathcal{L}}_{BC}(\theta )={\displaystyle \frac{1}{M}}\sum _{j=1}^M{\Vert B\left[u_{\theta }\left(x_j\right)\right]\Vert }^2.$$

(16)

4.

Total loss minimization

Finally, the total loss function combines the error components in the form:

$$\mathcal{L}\left(\theta \right)= {\mathcal{L}}_{PDE}\left(\theta \right) + {\mathcal{L}}_{BC}\left(\theta \right).$$

(17)

The total loss in Equation (17) is minimized and the network parameters $\theta$ are optimized using gradient-based methods.

2.3. PINN Architecture

In our case, the input layer of the neural network consists of two neurons, one for each of the spatial 2D coordinates of the points of the beam domain. The neural network uses two hidden layers with 100 and 250 nodes respectively. Both hidden layers use hyperbolic tangent (tanh) as an activation function. Finally, the output layer consists of eight nodes which represent the physical variables relevant to the piezoelectric problem modeled: mechanical displacements, electric potential, stress, and dielectric displacements. Table 3 presents a summary of the architecture of the proposed PINN in terms of its layers.

To better illustrate the important design features and decisions of the proposed architecture, Figure 3 shows a more detailed diagram of the PINN. Unlike traditional PINN architectures, which typically use a uniform number of nodes per layer [29,30], our model uses a non-uniform number of neurons per layer. Additionally, we use two hidden layers. These design choices were inspired by previous studies suggesting that PINNs with different numbers of neurons per hidden layer could be more effective, and that deeper neural networks do not necessarily lead to improved accuracy [31,32].

Another crucial feature is the use of not only mechanical displacements and electric potential as outputs of the neural network, but also the stress tensor and the dielectric displacements, an approach inspired by a previous application of PINNs to elastodynamic problems [33]. From Equations (5)–(7), it follows that the latter can be expressed in terms of the derivatives of the mechanical displacements and electric potential. If only the mechanical displacements and electric potential had been modeled as the neural network outputs, enforcing the PDEs in Equations (3) and (4) would have required computing second-order derivatives. This would have significantly increased the computational cost of the training [34]. However, by modeling stress and dielectric displacement explicitly, this computational burden is avoided. In this case, only the first derivatives of the stress and dielectric displacement are needed to satisfy Equations (3) and (4). Similarly, only the first derivatives of the mechanical displacement and electric potential are needed to satisfy the constitutive relations in Equations (1) and (2). This design choice simplified the enforcement of the governing equations and contributed to the capacity of the PINN to accurately model the coupled multiphysics problem.

2.4. Loss Formulation Based on Governing Equations and Boundary Conditions

The first contribution to the PDE loss ${\mathcal{L}}_{PDE}$ is based on Equations (1), (2), and (5), and the design of the PINN, presented in Section 2.3. First, stress and dielectric displacement are expressed in terms of the derivatives of mechanical displacements and the electric potential, modeled by the PINN following Equation (5). Then, these values are compared with the explicit outputs of the stress and dielectric displacement of the neural network. This imposes a constraint which ensures the PINN models the constitutive piezoelectric equations. This is done by the following set of equations, where the mechanical displacements, electric potential, stress and dielectric displacements correspond to the outputs of the neural network, and the $x$ and $y$ independent variables correspond to the inputs of the network:

$$r_{{\mathsf{\sigma }}_{11}}={\mathsf{\sigma }}_{11,\theta }-\left(c_{11}{\displaystyle \frac{\partial u_{\theta }}{\partial x}}+c_{12}{\displaystyle \frac{\partial v_{\theta }}{\partial y}}+e_{11}{\displaystyle \frac{\partial {\varphi }_{\theta }}{\partial x}}+e_{31}{\displaystyle \frac{\partial {\varphi }_{\theta }}{\partial y}}\right),$$

(18)

$$r_{{\mathsf{\sigma }}_{22}}={\mathsf{\sigma }}_{22,\theta }-\left(c_{12}{\displaystyle \frac{\partial u_{\theta }}{\partial x}}+c_{22}{\displaystyle \frac{\partial v_{\theta }}{\partial y}}+e_{13}{\displaystyle \frac{\partial {\varphi }_{\theta }}{\partial x}}+e_{33}{\displaystyle \frac{\partial {\varphi }_{\theta }}{\partial y}}\right),$$

(19)

$$r_{{\mathsf{\sigma }}_{12}}={\mathsf{\sigma }}_{12,\theta }-\left(c_{33}{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial u_{\theta }}{\partial y}}+{\displaystyle \frac{\partial v_{\theta }}{\partial x}}\right)+e_{14}{\displaystyle \frac{\partial {\varphi }_{\theta }}{\partial x}}+e_{34}{\displaystyle \frac{\partial {\varphi }_{\theta }}{\partial y}}\right),$$

(20)

$$r_{D_x}=D_{x,\theta }-\left(e_{11}{\displaystyle \frac{\partial u_{\theta }}{\partial x}}+e_{13}{\displaystyle \frac{\partial v_{\theta }}{\partial y}}+e_{14}\left({\displaystyle \frac{\partial u_{\theta }}{\partial y}}+{\displaystyle \frac{\partial v_{\theta }}{\partial x}}\right)-{\mathsf{ϵ}}_1{\displaystyle \frac{\partial {\varphi }_{\theta }}{\partial x}}\right),$$

(21)

$$r_{D_y}=D_{y,\theta }-\left(e_{31}{\displaystyle \frac{\partial u_{\theta }}{\partial x}}+e_{33}{\displaystyle \frac{\partial v_{\theta }}{\partial y}}+e_{34}\left({\displaystyle \frac{\partial u_{\theta }}{\partial y}}+{\displaystyle \frac{\partial v_{\theta }}{\partial x}}\right)-{\mathsf{ϵ}}_2{\displaystyle \frac{\partial {\varphi }_{\theta }}{\partial y}}\right).$$

(22)

The second contribution to ${\mathcal{L}}_{PDE}$ ensures that the PINN complies with Equations (3) and (4), by computing the divergence of the stress and dielectric displacement outputs of the neural network. Depending on the modelling scenario, the exact equation used will vary. For the indirect piezoelectric effect, both force components are zero. For the direct piezoelectric effect, the force component in $x,F_X$ remains zero while the component in $y,F_y$ is $1 N$, corresponding to an upward applied force at the tip of the beam:

$${\mathrm{r}}_{\sigma ,1}={\displaystyle \frac{\partial {\mathsf{\sigma }}_{11,\theta }}{\partial x}}+{\displaystyle \frac{\partial {\mathsf{\sigma }}_{12,\theta }}{\partial y}}+ F_x,$$

(23)

$${\mathrm{r}}_{\sigma ,2}={\displaystyle \frac{\partial {\mathsf{\sigma }}_{12,\theta }}{\partial x}}+{\displaystyle \frac{\partial {\mathsf{\sigma }}_{22,\theta }}{\partial y}}+F_y,$$

(24)

$${\mathrm{r}}_D={\displaystyle \frac{\partial D_{x,\theta }}{\partial x}}+{\displaystyle \frac{\partial D_{y,\theta }}{\partial y}}.$$

(25)

The complete PDE loss of the PINN combines these terms, computing an error function by applying the Mean Squared Error (MSE) to each term, using $N_r$ collocation points:

$$\begin{gathered} {\mathcal{L}}_c\left(\theta \right)={\displaystyle \frac{1}{N_r}}\sum _{i=1}^{N_r}{\Vert r_{{\sigma }_{11}}\left(x_i,y_i\right) +r_{{\sigma }_{22}}\left(x_i,y_i\right) + r_{{\sigma }_{12}}\left(x_i,y_i\right) + r_{D_x}\left(x_i,y_i\right) \\ + r_{D_y}\left(x_i,y_i\right) \Vert }^2, \end{gathered}$$

(26)

$${\mathcal{L}}_d\left(\theta \right)={\displaystyle \frac{1}{N_r}}\sum _{i=1}^{N_r}{\Vert r_{\sigma ,1}\left(x_i,y_i\right)+r_{\sigma ,2}\left(x_i,y_i\right)+r_D\left(x_i,y_i\right) \Vert }^2,$$

(27)

$${\mathcal{L}}_{PDE}\left(\theta \right)={\mathcal{L}}_c\left(\theta \right)+{\mathcal{L}}_d\left(\theta \right).$$

(28)

Additionally, to comply with the boundary conditions of Equations (10) and (11), as mentioned in the procedure of Section 2.2, a set of boundary conditions points specific to the boundary subdomains presented in Figure 2 is evaluated in the neural network. The output is used to enforce the traction-free and zero dielectric displacement conditions as follows:

$${\mathcal{L}}_T\left(\theta \right)= {\displaystyle \frac{1}{N_T}}\sum _{i=1}^{N_T}{\Vert {\mathsf{\sigma }}_{11,\theta }\left(x_i,y_i\right)+ {\mathsf{\sigma }}_{12,\theta }\left(x_i,y_i\right)\Vert }^2 ,$$

(29)

$${\mathcal{L}}_{D_{right}}\left(\theta \right)={\displaystyle \frac{1}{N_{right}}}\sum _{i=1}^{N_{right}}{\Vert D_{x,\theta }\left(x_i,y_i\right) \Vert }^2,$$

(30)

$${\mathcal{L}}_{D_{left}}\left(\theta \right)={\displaystyle \frac{1}{N_{left}}}\sum _{i=1}^{N_{left}}{\Vert D_{x,\theta }\left(x_i,y_i\right) \Vert }^2,$$

(31)

$${\mathcal{L}}_{BC}\left(\theta \right)={\mathcal{L}}_T\left(\theta \right)+{\mathcal{L}}_{D_{right}}\left(\theta \right)+{\mathcal{L}}_{D_{left}}\left(\theta \right),$$

(32)

where $N_{right}$ is the number of boundary condition points at the right end of the beam and $N_{left}$ is the number of boundary condition points at the left end of the beam.

The combination of the PDE loss and the boundary conditions loss results in the total loss:

$$\mathcal{L}\left(\theta \right) = {\mathcal{L}}_{PDE}\left(\theta \right) + {\mathcal{L}}_{BC}\left(\theta \right).$$

(33)

Finally, boundary conditions given by Equations (8) and (9) are directly enforced in the neural network outputs as hard constraints in order to automatically satisfy them [35], as follows:

$$u_{modified} =x·u,$$

(34)

$$v_{modified}=x·v,$$

(35)

$${\phi }_{modified}=y ·\left(y-1\right)· \phi +\left({\displaystyle \frac{y}{height}} ·voltage\right).$$

(36)

$${\phi }_{modified}=y · \phi .$$

(37)

In these equations, $x$ and $y$ represent the spatial coordinates where the neural network is evaluated. Equations (34) and (35) guarantee zero mechanical displacements at the left end of the beam, representing clamped boundary conditions. Equation (36) is used to satisfy the electric potential constraint for the indirect piezoelectric effect, where 0 V are applied at the bottom of the beam, and 100 volts are applied at the top. Equation (37) is used for the direct piezoelectric effect, where the bottom of the beam is grounded (0 V).

2.5. PINN Training

This section describes the training procedure employed for the proposed PINN, including the generation of collocation and boundary condition points, the evaluation strategy, and the optimization process used to solve the piezoelectric problem.

2.5.1. Training and evaluation: Data Generation

Data sampling and generation were performed using Latin Hypercube Sampling (LHS), for both boundary conditions and collocation points. Figure 4 illustrates the points generated to enforce the traction-free and zero dielectric displacements conditions at the right and left ends of the beam. 150 points were sampled for each of the right and left ends, generating 300 boundary points in total.

Figure 5 shows the collocation points used to train the PINN on the domain of the beam. 22,500 points were sampled in order to evaluate the PDE residuals according to the loss functions previously defined. The high density of the points helped to represent accurately the piezoelectric equations.

In order to visually analyze the predictions of the PINN, 4000 testing points were generated, depicted in Figure 6. This amount and distribution of points enabled a complete evaluation of the predictions generated by the PINN for the entire beam domain.

In order to perform quantitative testing against an “exact” solution obtained via the classical FEM technique, we used the testing points shown in Figure 7, at each of which we calculated a comparison value using FEM to compare against values generated using the trained neural network. The distribution of the points was selected based on the FEM tool employed, described later.

2.5.2. Training Procedure and Optimizer Configuration

The training of the neural network was done using a two-phase strategy to minimize the loss function Equation (30) [28]. In the first phase, we employed the ADAM optimizer with a learning rate of 0.001 for 1,000 epochs. During each epoch, the full set of collocation and boundary points was passed through the network as a single batch. The second phase of the training employed the L-BFGS optimizer, with a learning rate of 0.01, iterating over 200 epochs.

2.6. Validation of PINN Results

The indirect piezoelectric predictions of the PINN were validated against data generated using a Finite Element Method (FEM) solution. The FEM solution used a triangular 2-D mesh with approximately 3605 nodes [26]. The spatial coordinates generated to compare both methods are shown in Figure 7. These points were fed into the PINN and the predictions for mechanical displacements and electric potential were compared with the FEM results. Let $\mathsf{\Psi }\left(x\right)$ denote the ground truth values (e.g. the FEM solution) and ${\mathsf{\Psi }}_{\theta }\left(x\right)$ the approximation generated by the neural network. First, the comparison between PINN and FEM results was done, by computing and plotting the pointwise relative error, computed as:

$$\mathrm{P}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e} \mathrm{R}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e} \mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}={\displaystyle \frac{\left|{\mathsf{\Psi }}_{\theta }\left(x\right)- \mathsf{\Psi }\left(x\right)\right|}{\left|\mathsf{\Psi }\left(x\right)+{10}^{-25}\right|}},$$

(38)

where a small constant term ${10}^{-25}$ was added to the denominator to avoid numerical instability and division by zero in boundary regions where $\mathsf{\Psi }\left(x\right)$ approaches zero.

Second, to obtain a global comparison, the relative $L^2$ error metric was computed over the full set of points compared as:

$${Relative L}^{2}error={\displaystyle \frac{{\Vert {\mathsf{\Psi }}_{\theta }\left(x\right)- \mathsf{\Psi }\left(x\right)\Vert }_2}{{\Vert \mathsf{\Psi }\left(x\right)\Vert }_2}}.$$

(39)

2.7. Implementation Details and Reproducibility

The simulations and training were performed using Python. PyTorch (v2.6.0+cu124) was the main library employed for automatic differentiation, neural network training, and numerical computations. The experiments were carried out on the Google Colab platform. An Nvidia Tesla T4 GPU with 15 GB of RAM was used, along with a CPU featuring 51 GB of system memory. The environment was configured with CUDA 12.4 and cuDNN version 90100.

3. Results

This section presents the results obtained with the proposed Physics-Informed Neural Network (PINN) model for predicting piezoelectric behaviour in a bimorph cantilever beam. A brief description of the model’s performance when modeling the indirect effect, as well as qualitative and quantitative comparisons between the PINN predictions and Finite Element Method (FEM) simulations, are presented. This includes visualizations of the predicted fields and the relative $L^2$ error. For the direct effect, the evaluation focuses on the predicted distributions of displacement and electric potential.

3.1. Results: Indirect Piezoelectric Effect

This subsection presents the results obtained with the PINN when modeling the indirect piezoelectric effect. A visual comparison, supported by a relative error computation, is provided between the predicted horizontal and vertical displacements, as well as the electric potential, and the corresponding FEM results. Finally, the relative $L^2$ error norms are reported to assess the accuracy of the model.

3.1.1. Comparison Between PINN and FEM

Figure 8 illustrates the bimorph piezoelectric cantilever beam vertical cross-section with the results obtained for the deflection in x under the indirect effect. Panel (a) shows the prediction obtained by the proposed PINN, where the red color represents a positive deflection in x, and the blue color a negative deflection in x. PINN predicted positive deflection in x in the upper layer and negative deflection in x in the lower layer. As indicated by the PINN prediction, the displacement is on the order of ${10}^{-7}$ m. Panel (b) presents the expected deflection in x computed via FEM, which shows similar qualitative behavior to the PINN prediction: positive deflection in x in the upper layer and negative deflection in x in the lower layer. Furthermore, the deflection in x given by the FEM simulation is also on the order of ${10}^{-7}$ m. Panel (c) displays the pointwise relative error between both methods. The main disparity occurs at the left end of the beam, where a boundary condition with value 0 is enforced, and at certain regions of the midline. The greatest relative error reaches a value close to 2.4.

Figure 9 presents the deflection in $y$ of the beam. Panel (a) depicts the predictions of the PINN, while panel (b) shows the expected behavior of the deflection in y. For both results, the left part of the beam has almost no deflection in y, indicated by the red color, as expected by the clamped left end of the beam. Towards the right end, the deflection in y gradually increases negatively, with the maximum negative deflection in y occurring at the left tip of the beam. Panel (c) shows the pointwise relative error between the PINN and FEM results. The relative error on the left side of the beam is low. However, the error increases progressively towards the right side of the beam. The largest difference can be found at the left end of the beam, where a Dirichlet boundary condition of value 0 is enforced. It can be observed that the relative error decreases as the deflection in $y$ increases negatively. The highest relative error reaches a value close to 1.0.

Figure 10 shows the electric potential ϕ distribution applied in the piezoelectric beam under the indirect effect. Panel (a) presents the prediction obtained with the proposed PINN, while panel (b) displays the FEM reference solution. In both results, the electric potential increases along the vertical direction, with the minimum value at the bottom surface and the maximum at the top, with a value of approximately 100 V. Both results show smooth gradients across the height of the beam. Panel (c) displays the pointwise relative error between the two fields. The main source of discrepancy is visible at the bottom of the beam, where a boundary condition with value 0 is enforced. It can be observed that the relative error in the rest of the domain is 0. The maximum value of error is 1.0.

3.1.2. Relative $L^2$

Error Evaluation

Table 5 lists the relative $L^2$ error measures for the PINN model evaluated, concerning the FEM reference solutions. The errors for the horizontal and vertical displacements are on the order of ${10}^{-1}$, and the error for the electric potential is in the order of ${10}^{-3}$.

3.2. PINN Results for the Direct Piezoelectric Effect

In this subsection, the results obtained for modeling the direct piezoelectric effect, using the proposed PINN, are presented. No comparison is made with FEM results due to the absence of available experimental or FEM reference data.

PINN Predictions of the Direct Piezoelectric Effect

Figure 11 presents the predictions from the PINN model for the behavior of the beam under the direct piezoelectric effect. Panel (a) depicts the deflection in $x$ of the beam. The negative values of the deflection in $x$ of the upper layer, which are associated with the blue region, indicate the compression of the beam in the upper layer. The positive values of the lower layer, represented by the red region, indicate the beam under tension.

Panel (b) shows the deflection in $y$ of the beam. Due to the clamped boundary condition on the left end, the deflection in $y$ is zero in that region. Moving towards the right, the deflection increases gradually, with the right tip of the beam experiencing the greatest deflection in $y$ This deflection in $y$ is positive, indicating that the beam bends upward.

Panel (c) displays the electric potential, which is the primary quantity of interest when considering the direct piezoelectric effect. The electric potential is zero at the lower surface of the beam, as imposed by the boundary condition. The PINN predicts a clear voltage gradient between the lower and upper surfaces, with the potential decreasing from bottom to top. At the top the electric potential is on the order of ${-10}^{-3}V$.

Finally, in Figure 12, the beam deformation is illustrated. The original undeformed beam is shown in blue, while the deformed shape is shown in orange. The deformation is scaled by a factor of 50 for visualization purposes. As observed, the beam is bent upwards.

4. Discussion

The obtained results demonstrate that the proposed Physics-Informed Neural Network (PINN) design, explicitly including stress and dielectric displacement as outputs, is able to model the behavior of both the indirect and direct piezoelectric effect on a polyvinylidene fluoride (PVDF) bimorph cantilever beam. For the indirect effect, the model generates accurate predictions when compared with FEM reference solutions in both visual analysis and low metric errors. For the direct piezoelectric effect, the model predictions exhibit physically consistent behavior, even without FEM validation data.

For the indirect piezoelectric effect, the pointwise comparison of the PINN predictions for mechanical displacements with the reference FEM results (Figure 9c and Figure 10c) shows good correspondence between both methods. This indicates that the PINN was able to accurately model both variables. This is quantitatively supported by the $L^2$ errors of Table 5 which are low. However, the figures also evidence that the PINN accuracy decreases in the regions where a Dirichlet boundary condition with value 0 is enforced, in contrast with regions with larger displacement values. This behavior is expected, as the relative error becomes more sensitive in areas where the reference solution approaches zero, amplifying discrepancies.

Regarding the direct piezoelectric effect, even if there was no FEM data available for comparison, the predictions of the PINN show physically consistent behavior. The beam bends upward under the applied 1 N force in the positive y-direction, as expected. Additionally, the electric potential shows a smooth gradient from 0 V at the lower end of the beam to approximately −0.001 V at the upper end. This indicates a difference in electric potential generated by the applied load, as typical for the direct piezoelectric effect. The consistent behavior of the predictions for the direct effect suggests that the proposed PINN design can generalize to different configurations of the piezoelectric problem.

The analysis of the results evidenced that the proposed PINN model was able to accurately model the behavior of the bimorph cantilever beam under the indirect piezoelectric effect. The PINN also showed physically consistent behavior for the direct piezoelectric effect. Additionally, the decision to include stress and dielectric displacement explicitly in the neural network outputs proved key to the efficiency of the approach and the modeling of the PDEs.

5. Conclusions

This study proposes a Physics Informed Neural Network (PINN) model to simulate both the indirect and direct piezoelectric effects in a 2D polyvinylidene fluoride (PVDF) bimorph cantilever beam under static conditions. A key contribution to the approach is the inclusion of stress and dielectric displacement as explicit outputs of the Neural Network (NN), in addition to the mechanical displacements and electric potential. This proposal avoids the need for second-order derivatives and simplifies the formulation of the Partial Differential Equations (PDEs) involved in the NN loss function. The results show that PINNs accurately modeled the indirect piezoelectric effect, by achieving low pointwise absolute and $L^2$ errors. Furthermore, the model generated physically consistent predictions for the direct effect, even in the absence of Finite Element Method (FEM) validation data. This work demonstrates the potential of PINNs as an alternative to traditional FEM approaches for multiphysics piezoelectric simulations.

6. Future Work

First, experimental data should be obtained to validate the reliability of the PINN predictions for the direct piezoelectric effect. Second, the proposed architecture could be tested on more complex piezoelectric configurations and devices. In addition, emerging neural network architectures such as Physics-Informed Kolmogorov-Arnold Networks (PIKANs) could be explored to potentially improve convergence, interpretability, and approximation capabilities. Finally, the approach could be extended to time-dependent and frequency-domain formulations of the piezoelectric equations, enabling applications in different domains such as MEMS devices and energy harvesting systems.