Analysis of quantum-dot systems under thermal loads based on gradient elasticity

Lattice mismatch or thermal expansion difference between a burried quantum-dot (QD) and a surrounding piezoelectric matrix induces both elastic and electric fields in the structure. Since both the elastic and electric fields are equally important to the understanding of the photonic and electronic features in nanostructured semiconductors [1], a reliable analysis of these fields are crucial to the design and fabrication of such structures. Several quantum mechanical methods for modifying the material optical properties were discussed in [2]. In many previous papers, semicoupled computational approaches were applied to analyze QD systems. The initial strains are computed first based on the elasticity equations. Then the purely elastic solution is used to derive the governing equations for the electric potential and electric field. Finally, the electric field is solved under suitable purely electric boundary conditions. The semicoupled approach was applied very frequently for its simplicity [3–9], but later Pan [10, 11] made a comparison between fully coupled and semicoupled models, and demonstrated that the semicoupled model could produce very serious errors for both the elastic and electric fields. Therefore, the fully coupled approach should always be used to model QD systems where piezoelectric semiconductor inclusions are embedded in a piezoelectric matrix.

Analytical solutions for an inclusion of a polygonal shape was presented in [12, 13] using the special conformal mapping method. However, conformal mapping method is inefficient for the general polygonal shape since the conformal mapping function may involve infinite terms. Green's function method [14–16] proved to be more efficient, particularly suited for polygonal shapes. The perturbation method [17, 18] also allows considering material anisotropy and arbitrary shape of inclusion. A misfitted inclusion embedded in a finite spherical matrix was analyzed in the framework of linear piezoelectricity by Seo et al [19]. Recently, Sun et al [20] analyzed quantum wires of arbitrary polygonal cross sections with graded eigenstrain in an anisotropic piezoelectric full plane. The eigenstrain field inside the inclusion is expressed as linear functions of Cartesian coordinates. A universal and efficient analytical method was proposed by Zou and coworkers [21, 22] where the quantum wires can be of any shape and the eigenstrain can be a general polynomial function of the coordinates. Jurczak and Dluzewski [23] evaluated the effect of higher-order elastic and piezoelectric coefficients on the elastic and electric fields in III-nitride wurtzite QDs.

The shape and size of the QDs play a major role in determining the optical response of QD systems. It has been demonstrated by Liu et al [24] that the shape and the area density of the QDs can be controlled by varying the growth parameters such as the growth temperature, growth rate, annealing progress and the growth interruption. Strain and piezoelectric effects are considered as tuning parameters for optical response of semiconductor nanostructures. Coupling thermo-mechanical fields (thermal expansion coefficient) and thermo-electric fields (pyroelectric parameter) can also play a very important role in QD technology. Hence temperature can also be considered as another tuning parameter. The material composition grading factor was also used as another tuning parameter [25].

Microelectronic circuits can work at extreme thermal conditions, therefore, it is important to understand the steady state and transient effects of temperature on the optoelectronic properties. Giazotto et al [26] analyzed QDs in low temperature Peltier refrigeration, whereas Chen et al [27] considered higher temperature regime. Patil and Melnik [28] studied thermoelectromechanical effects in QDs under stationary boundary conditions and steady state thermal conditions. Thermo-piezoelectricity theory allows for modeling QDs more realistically.

The first attempt to develop a theory of thermo-piezoelectricity was done by Mindlin [29]. Later, Nowacki [30] included more physical laws into that theory. Due to the complicated coupling among the three (mechanical, electrical, thermal) fields, the theory has been applied only to very simple boundary-value problems in the beginning. Analytical methods with Fourier–Hankel transformation were successfully applied to simple problems [31, 32]. The general boundary-value problem with a complex geometry requires powerful computational tools such as the finite element method (FEM) [33–36]. In recent years, meshless formulations have also become popular due to their high adaptivity and low costs in preparation of input and output data for numerical analyses. A variety of meshless methods have been proposed and some of them were also applied to thermo-piezoelectric problems [37]. The inclusion of the thermal effects may help improve the performance characteristics of intelligent structures with distributed piezoelectric sensors and actuators [38]. Three-dimensional (3D) analysis of the coupled thermo-piezoelectro-mechanical behavior of multilayered plates was done by Vel and Batra [39] and Liew et al [40]. In many of the aforementioned studies, classical continuum mechanics proved to be accurate enough to model QD systems. This was also discussed in more details in [41].

If the dimensions of the QDs are of the same order of the material length scale, classical continuum mechanics should not be applied since it neglects the interaction of material microstructure, and the results are size-independent. The atomistic models overcome intrinsic limitations of classical elasticity. The development of multiscale approaches is preferable due to the extremely high requirements on computer memory in atomistic models. The multiscale approaches bridge atomistic and continuum subdomains. However, there are observed unphysical phenomena in multiscale approaches especially for time-dependent problems. The other alternative is to apply advanced continuum theories to account for intrinsic length scales. Advanced continuum theories can be based on the Eringen's non-local elasticity [42, 43], or Mindlin's strain-gradient elasticity [44, 45]. In the non-local theory, the stress at a reference point is a functional of strains at more points of the solid body. The higher-order strain gradients are considered in the constitutive equations of the gradient theory. The two length scales in Mindlin's theory [44] are taken equal to each other in Aifantis' theory [46–48]. This greatly simplifies further mathematical and implementational treatment. Fleck and Hutchinson [49] provided an overview of the various higher-order gradient theories.

This paper presents, for the first time, the thermo-electro-mechanical 3D analyses of QDs under transient thermal ;conditions in order to assess the thermal shock effect on the electromechanical response of the QDs in applications where cyclic thermoeletric loads are applied. The gradient theory for thermo-piezoelectricity is also applied to account for the size-dependent behavior of QDs in nano-sized piezoelectric structures. Solving this problem using the FEM requires using C¹-continuous elements to guarantee the continuity of the problem variables and their derivatives at the element boundaries. In this case, the additional complexity, as compared to using C⁰ shape functions, is considerable. Therefore, the original fourth order partial differential equation (PDE) is split into two PDEs of lower order that require only C⁰ shape functions. Numerical results obtained by the present gradient formulation are compared with those obtained by the classical coupled thermo-piezoelectric theory, and the influence of the size effect parameter is discussed.

The paper is organized as follows: section 2 presents the governing equations for the gradient thermo-piezoelectricity, section 3 details the FEM formulation, and section 4 describes the COMSOL model and numerical results. Conclusions are summarized in section 5.

The basic feature of a QD system is its inherent initial strain, which induces electric and elastic fields in a piezoelectric material. This initial strain is associated with the lattice mismatch. Usually the QD has almost vanishing thermal conduction and thermal expansion coefficients. In this analysis, a periodic distribution of QDs in the matrix is considered, and for the numerical simulations, a representative volume element (RVE) as illustrated in figure 1 (left) is used. This piezoelectric composite structure under thermal load can be described by the theory of thermo-piezoelectricity.

Zoom In Zoom Out Reset image size

Most piezoelectric materials are also temperature sensitive, i.e. an electric charge or voltage is generated when exposed to temperature variations. This effect is called the pyroelectric effect. However, for most materials, the inverse coupling, i.e., the heat generation caused by mechanical/electric fields, is very weak. Therefore, thermo-piezoelectricity with 'one-way coupling' is considered here for the QD system [28]. While the heat generation caused by mechanical and electrical fields can be neglected, the stresses and electrical displacements can be influenced by temperature variation. The governing equations of the phenomenological theory of piezoelectric materials under thermal loads consist of Maxwell's equations, the heat conduction equation and the balance of momentum. These governing equations, which are complemented by the boundary and initial conditions, should be solved for the unknown primary field variables which are the elastic displacement vector field u_i(x, τ), the electric potential ϕ(x, τ) (or, its gradient, the electric field vector E_i(x, τ)), and the temperature difference from the reference room temperature θ(x, τ), where x is the position vector and τ is the time. Taking into account the typical material coefficients, it can be found that the characteristic frequencies for thermal, elastic and electromagnetic processes are f_th = 10⁻³ Hz, f_el = 10⁴ Hz and f_elm = 10⁷ Hz, respectively. Thus, if we consider such bodies under transient loadings, with temporal changes corresponding to f_th and/or f_el, the changes of the electromagnetic fields can be assumed to be immediate, or in other words the electromagnetic fields can be considered quasi-static [50].

The strain tensor ε_ij, the electric field vector E_j and the temperature gradient vector β_j are related to the displacements u_i, the electric potential ϕ and the temperature θ, respectively, by

$$\begin{eqnarray}&&{\varepsilon }_{ij}=\displaystyle \frac{1}{2}({u}_{i,j}+{u}_{j,i});\,{E}_{j}=-{\phi }_{,j};\,{\beta }_{j}={\theta }_{,j}.\end{eqnarray} \tag{ 1 }$$

The elastic strain is given by subtraction of the eigenstrain and thermal strain from the total strain

$$\begin{eqnarray}&&{\varepsilon }_{ij}^{e}={\varepsilon }_{ij}-{\varepsilon }_{ij}^{* }-{\alpha }_{ij}\theta ,\,{\varepsilon }_{ij}^{* }=\chi \displaystyle \frac{{a}_{m}-{a}_{{\rm{Q}}{\rm{D}}}}{{a}_{{\rm{Q}}{\rm{D}}}},\\ \chi =\left\{\begin{array}{ll}1, & {\rm{in}}\,{\rm{QD}}\\ 0, & {\rm{in}}\,{\rm{matrix}}\end{array}\right.,\end{eqnarray} \tag{ 2 }$$

where α_ij is the coefficients of linear thermal expansion, ${\varepsilon }_{ij}^{* }$ is the eigenstrain tensor and a_m, a_QD are lattice constants of the matrix and QD, respectively. The reference state for displacements is specified by vanishing eigenstrains and at temperature θ = 0. As the substrate is relatively large compared to the QD, we follow common practice to neglect lattice mismatch inside the matrix and to assume homogeneous eigenstrains inside the QD [28].

The strain-gradient tensor is defined as

$$\begin{eqnarray}&&{\eta }_{ijk}={\varepsilon }_{ij,k}=\displaystyle \frac{1}{2}({u}_{i,jk}+{u}_{j,ik}).\end{eqnarray} \tag{ 3 }$$

If the size of the QDs is comparable to the size of the material microstructure, the results will be size-dependent. In this case, the classical elasticity theory should be replaced by an advanced continuum theory to account for the intrinsic length scale of the material. The one-parameter gradient theory is applied to the piezoelectric solids under thermal loads in this work. The deformation energy density in gradient theory for thermo-piezoelectricity is written as [51, 52]

$$\begin{eqnarray}&&{ \mathcal U }({\boldsymbol{\varepsilon }},{\rm{\nabla }}{\boldsymbol{\varepsilon }},{\bf{E}},\theta ,{\boldsymbol{\beta }})=\displaystyle \frac{1}{2}{c}_{ijkl}{\varepsilon }_{ij}{\varepsilon }_{kl}+\displaystyle \frac{1}{2}{g}_{ijklmn}{\eta }_{ijk}{\eta }_{lmn}\\ &&\,\,+\,\displaystyle \frac{1}{2}{h}_{ij}{E}_{i}{E}_{j}+\displaystyle \frac{1}{2}{\kappa }_{ij}{\beta }_{i}{\beta }_{j}-{e}_{ijk}{\varepsilon }_{ij}{E}_{k}-{\gamma }_{ij}{\varepsilon }_{ij}\theta \\ &&\,\,-\,{c}_{ijkl}{\varepsilon }_{ij}{\varepsilon }_{kl}^{\ast }+{p}_{i}{E}_{i}\theta ,\end{eqnarray} \tag{ 4 }$$

where c_ijkl, e_ijk, h_ij, p_i and κ_ij are the material elastic, piezoelectric, dielectric, pyroelectric and thermal conductivity tensors in the thermo-piezoelectric medium, respectively. The tensor g_ijklmn represents the purely non-local elastic effects, i.e., the elastic strain-gradient elasticity. The stress–temperature modulus γ_ij is expressed in terms of the elastic stiffness coefficients c_ijkl and the linear thermal expansion coefficients α_kl as γ_ij = c_ijklα_kl.

Under infinitesimal deformation, the Cauchy stress σ_ij, the higher-order stress τ_ijk, the electric displacement D_i can be expressed in terms of the derivatives of the deformation energy density as

$$\begin{eqnarray}&&{\sigma }_{ij}({\bf{x}},\tau )=\displaystyle \frac{\partial { \mathcal U }}{\partial {\varepsilon }_{ij}}={c}_{ijkl}({\varepsilon }_{kl}({\bf{x}},\tau )-{\varepsilon }_{kl}^{\ast }({\bf{x}},\tau ))\\ &&\,\,-\,{e}_{kij}{E}_{k}({\bf{x}},\tau )-{\gamma }_{ij}\theta ({\bf{x}},\tau );\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\tau }_{ijk}({\bf{x}},\tau )=\displaystyle \frac{\partial { \mathcal U }}{\partial {\eta }_{ijk}}={g}_{ijklmn}{\eta }_{lmn}({\bf{x}},\tau ).\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{D}_{i}({\bf{x}},\tau )=\displaystyle \frac{\partial { \mathcal U }}{\partial {E}_{i}}={e}_{ijk}{\varepsilon }_{jk}({\bf{x}},\tau )+{h}_{ij}{E}_{j}({\bf{x}},\tau )\\ &&\,\,+\,{p}_{i}\theta ({\bf{x}},\tau );\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\psi }_{i}({\bf{x}},\tau )=\displaystyle \frac{\partial { \mathcal U }}{\partial {\beta }_{i}}={\kappa }_{ij}{\beta }_{j}({\bf{x}},\tau ).\end{eqnarray} \tag{ 8 }$$

Equation (8) is the Fourier's law for thermal conduction. Note that in uncoupled thermo-elasticity, stresses are functions of temperature, but the heat flux is independent of the elastic and electrical fields. Constitutive equations (5), (7) and (8) for material with cubic symmetry can be written in matrix form, using the reduced Voigt notation, as

$$\begin{eqnarray}&&\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{23}\\ {\sigma }_{13}\\ {\sigma }_{12}\end{array}\right\}\,=\,\left[\begin{array}{cccccc}{c}_{11} & {c}_{12} & {c}_{12} & 0 & 0 & 0\\ {c}_{12} & {c}_{11} & {c}_{12} & 0 & 0 & 0\\ {c}_{12} & {c}_{12} & {c}_{11} & 0 & 0 & 0\\ 0 & 0 & 0 & {c}_{44} & 0 & 0\\ 0 & 0 & 0 & 0 & {c}_{44} & 0\\ 0 & 0 & 0 & 0 & 0 & {c}_{44}\end{array}\right]\,\left\{\begin{array}{c}{\varepsilon }_{11}^{r}\\ {\varepsilon }_{22}^{r}\\ {\varepsilon }_{33}^{r}\\ 2{\varepsilon }_{23}^{r}\\ 2{\varepsilon }_{13}^{r}\\ 2{\varepsilon }_{12}^{r}\end{array}\right\}\\ \,\,-\,\left[\begin{array}{ccc}0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0\\ {e}_{14} & 0 & 0\\ 0 & {e}_{14} & 0\\ 0 & 0 & {e}_{14}\end{array}\right]\left\{\begin{array}{c}{E}_{1}\\ {E}_{2}\\ {E}_{3}\end{array}\right\}\,-\,\left\{\begin{array}{c}{\gamma }_{11}\\ {\gamma }_{22}\\ {\gamma }_{33}\\ 0\\ 0\\ 0\end{array}\right\}\theta ,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{c}{D}_{1}\\ {D}_{2}\\ {D}_{3}\end{array}\right\}\,=\,\left[\begin{array}{cccccc}0 & 0 & 0 & {e}_{14} & 0 & 0\\ 0 & 0 & 0 & 0 & {e}_{14} & 0\\ 0 & 0 & 0 & 0 & 0 & {e}_{14}\end{array}\right]\left\{\begin{array}{c}{\varepsilon }_{11}\\ {\varepsilon }_{22}\\ {\varepsilon }_{33}\\ 2{\varepsilon }_{23}\\ 2{\varepsilon }_{13}\\ 2{\varepsilon }_{12}\end{array}\right\}\\ \,\,+\,\left[\begin{array}{ccc}{h}_{11} & 0 & 0\\ 0 & {h}_{11} & 0\\ 0 & 0 & {h}_{11}\end{array}\right]\left\{\begin{array}{c}{E}_{1}\\ {E}_{2}\\ {E}_{3}\end{array}\right\}+\left\{\begin{array}{c}{p}_{1}\\ {p}_{1}\\ {p}_{1}\end{array}\right\}\theta ,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{c}{\psi }_{1}\\ {\psi }_{2}\\ {\psi }_{3}\end{array}\right\}=\left[\begin{array}{ccc}{\kappa }_{11} & 0 & 0\\ 0 & {\kappa }_{11} & 0\\ 0 & 0 & {\kappa }_{11}\end{array}\right]\left\{\begin{array}{c}{\beta }_{1}\\ {\beta }_{2}\\ {\beta }_{3}\end{array}\right\},\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&\left\{\begin{array}{c}{\gamma }_{11}\\ {\gamma }_{22}\\ {\gamma }_{33}\end{array}\right\}=\left[\begin{array}{ccc}{c}_{11} & {c}_{12} & {c}_{12}\\ {c}_{12} & {c}_{11} & {c}_{12}\\ {c}_{12} & {c}_{12} & {c}_{11}\end{array}\right]\left\{\begin{array}{c}{\alpha }_{11}\\ {\alpha }_{11}\\ {\alpha }_{11}\end{array}\right\},\end{eqnarray} \tag{ 12 }$$

where ε_ij^r = ε_ij − ε_ij*.

Equations (9)–(11) can be written in compact form as

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}={\bf{C}}{{\boldsymbol{\varepsilon }}}^{r}-{\bf{LE}}-{\boldsymbol{\gamma }}\theta ,\,{\bf{D}}={{\bf{L}}}^{T}{\boldsymbol{\varepsilon }}+{\bf{AE}}+{\bf{p}}\theta ,\\ &&{\boldsymbol{\psi }}=\,{\boldsymbol{\kappa }}{\boldsymbol{\beta }}.\end{eqnarray} \tag{ 13 }$$

The higher-order stress tensor τ_ijk can be written in terms of the Cauchy stress tensor giving the total stress tensor T_ij in the form [48]

$$\begin{eqnarray}&&{T}_{ij}={\sigma }_{ij}-{l}^{2}{\sigma }_{ij,mm},\end{eqnarray} \tag{ 14 }$$

where l is the internal length material parameter. One can observe that in the limit, l → 0, the total stress is equal to the Cauchy stress known from the classical ('non-gradient') elasticity. In an earlier paper, the authors derived the governing equations for the gradient piezoelectricity without thermal effects [53]. The equilibrium equation for the total stress can be rewritten in terms of the Cauchy stress [54] as

$$\begin{eqnarray}&&{T}_{ij,j}={\sigma }_{ij,j}-{l}^{2}{\sigma }_{ij,jmm}=0,\end{eqnarray} \tag{ 15 }$$

where vanishing body forces are assumed. The governing equations for the electrical and thermal fields (Maxwell's and heat conduction equations), with vanishing volume density of free charges and density of heat sources, are expressed as

$$\begin{eqnarray}&&{D}_{i,i}({\bf{x}},\tau )=0,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\psi }_{i,i}({\bf{x}},\tau )+\rho c\dot{\theta }({\bf{x}},\tau )=0,\end{eqnarray} \tag{ 17 }$$

where ρ and c are the mass density and specific heat, respectively. Dots over a quantity indicate time derivatives. Considering that the size effect is restricted to only elastic strains, the constitutive equation (5) can be substituted into the equilibrium equation (15) to get

$$\begin{eqnarray}&&{c}_{ijkl}{u}_{k,lj}({\bf{x}},\tau )-{e}_{kij}{E}_{k,j}({\bf{x}},\tau )-{\gamma }_{ij}{\theta }_{,j}({\bf{x}},\tau )\\ &&\,-\,{l}^{2}{c}_{ijkl}{u}_{k,ljmm}({\bf{x}},\tau )=0.\end{eqnarray} \tag{ 18 }$$

This is a PDE of the fourth order. If traditional finite elements are used for the numerical solution of such problems, then C¹ displacement continuity is required. Askes et al [55] presented a decomposition technique, where a set of Helmholtz equations associated with the theory of gradient elasticity were solved to handle the gradient dependence. In this case, the additional complexity as compared to C⁰ shape functions is considerable. In the present paper such a technique is applied to gradient piezoelectricity under a thermal load. So, the original fourth order PDE (18) is split into two PDEs, (19) and (20), of lower order that require C⁰ shape functions in the finite element interpolation.

If we use the notation

$$\begin{eqnarray}&&{\varepsilon }_{ij}^{c}({\bf{x}},\tau )={\varepsilon }_{ij}({\bf{x}},\tau )-{l}^{2}{\varepsilon }_{ij,kk}({\bf{x}},\tau ),\end{eqnarray} \tag{ 19 }$$

the PDE (18) can be expressed as

$$\begin{eqnarray}&&{c}_{ijkl}{\varepsilon }_{kl,j}^{c}({\bf{x}},\tau )-{e}_{kij}{E}_{k,j}({\bf{x}},\tau )-{\gamma }_{ij}{\theta }_{,j}({\bf{x}},\tau )=0.\end{eqnarray} \tag{ 20 }$$

The strain ε_ij^ec can be expressed by the kinematic equation

$$\begin{eqnarray}&&{\varepsilon }_{ij}^{c}=\displaystyle \frac{1}{2}({u}_{i,j}^{c}+{u}_{j,i}^{c}).\end{eqnarray} \tag{ 21 }$$

The boundary conditions for the coupled equations are different from those valid for the original fourth order format of gradient elasticity [55]. The following essential and natural boundary conditions are considered for the mechanical field represented by equation (20):

$$\begin{eqnarray}&&{u}_{i}^{c}({\bf{x}},\tau )={\bar{u}}_{i}^{c}({\bf{x}},\tau )\,{\rm{on}}\,{{\rm{\Gamma }}}_{{u}_{i}};\\ &&{T}_{i}({\bf{x}},\tau )=({c}_{ijkl}{\varepsilon }_{kl}^{c}-{e}_{kij}{E}_{k}-{\gamma }_{ij}\theta ){n}_{j}\\ &&\,\,\,=\,{\bar{t}}_{i}({\bf{x}},\tau )\,{\rm{on}}\,{{\rm{\Gamma }}}_{{t}_{i}},\end{eqnarray} \tag{ 22 }$$

where n_i is the unit vector normal to the domain boundary, ${\bar{t}}_{i}$ are the prescribed tractions at ${{\rm{\Gamma }}}_{{t}_{i}}$ boundary, which is a part of the global domain boundary, and ${\bar{u}}_{i}^{c}$ are the prescribed displacements at Γ_u boundary. For the electrical and thermal fields, we have:

$$\begin{eqnarray}&&\phi ({\bf{x}})=\bar{\phi }({\bf{x}})\,{\rm{on}}\,{{\rm{\Gamma }}}_{p};\,{Q}_{i}({\bf{x}})={D}_{i}({\bf{x}}){n}_{i}\\ &&\,\,=\,\bar{Q}({\bf{x}}),\,{\rm{on}}\,{{\rm{\Gamma }}}_{Q},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\theta ({\bf{x}},\tau )=\bar{\theta }({\bf{x}},\tau )\,{\rm{on}}\,{{\rm{\Gamma }}}_{\theta };\\ &&q({\bf{x}},\tau )={\psi }_{i}({\bf{x}},\tau ){n}_{i}({\bf{x}})\\ &&\,\,\,=\,\bar{q}({\bf{x}},\tau )\,{\rm{on}}\,{{\rm{\Gamma }}}_{q},\end{eqnarray} \tag{ 24 }$$

where Γ_p, Γ_Q, Γ_θ, and Γ_q are the boundaries where electric potential $\bar{\phi },$ surface charge density $\bar{Q},$ temperature difference $\bar{\theta }$ and heat flux $\bar{q}$ are prescribed, respectively.

The essential and natural boundary conditions for the Helmholtz equation (19) are given by

$$\begin{eqnarray}&&{\varepsilon }_{ij}({\bf{x}},\tau )={\bar{\varepsilon }}_{ij}({\bf{x}},\tau )\,{\rm{on}}\,{{\rm{\Gamma }}}_{\varepsilon };\,{\varepsilon }_{ij,m}({\bf{x}},\tau ){n}_{m}({\bf{x}})\\ &&\,\,\,\,=\,{\bar{\varepsilon }}_{ij,m}({\bf{x}},\tau ){n}_{m}({\bf{x}})\,{\rm{on}}\,{{\rm{\Gamma }}}_{n}.\end{eqnarray} \tag{ 25 }$$

Initial conditions for the mechanical and thermal quantities are assumed as:

$$\begin{eqnarray}&&{{u}_{i}({\bf{x}},\tau )| }_{\tau =0}={u}_{i}({\bf{x}},0),\,{{\dot{u}}_{i}({\bf{x}},\tau )| }_{\tau =0}={\dot{u}}_{i}({\bf{x}},0),\\ &&{\theta ({\bf{x}},\tau )| }_{\tau =0}=\theta ({\bf{x}},0)\,{\rm{in}}\,{\rm{\Omega }}.\end{eqnarray} \tag{ 26 }$$

In the next section, a finite element formulation is developed to analyze the above stated boundary-value problem of the gradient-based thermo-piezoelectricity.

The FEM is developed here for the problem described by the system of PDEs (16), (17), (19) and (20). Since uncoupled thermo-elasticity is considered, the heat conduction problem is analyzed separately and the weak/variational forms are given by

$$\begin{eqnarray}&&\displaystyle {\int }_{{\rm{\Omega }}}({\lambda }_{i}\delta {\theta }_{,i}-\rho c\dot{\theta }\delta \theta )d{\rm{\Omega }}=\displaystyle {\int }_{{{\rm{\Gamma }}}_{q}}\bar{q}\delta \theta d{\rm{\Gamma }},\end{eqnarray} \tag{ 27 }$$

δθ is a test function such that ${\delta \theta | }_{\partial {\rm{\Omega }}-{{\rm{\Gamma }}}_{q}}=0.$

The weak forms for the rest governing equations are given by

$$\begin{eqnarray}&&\displaystyle {\int }_{{\rm{\Omega }}}{D}_{i,i}\delta \phi d{\rm{\Omega }}=\displaystyle {\int }_{{{\rm{\Gamma }}}_{Q}}\bar{Q}\delta \phi d{\rm{\Gamma }}-\displaystyle {\int }_{{\rm{\Omega }}}{D}_{i}\delta {\phi }_{,i}d{\rm{\Omega }}=0,\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\displaystyle {\int }_{{\rm{\Omega }}}[{c}_{ijkl}{\varepsilon }_{kl,j}^{c}({\bf{x}},\tau )-{e}_{kij}{E}_{k,j}({\bf{x}},\tau )-{\gamma }_{ij}{\theta }_{,j}({\bf{x}},\tau )]\delta {u}_{i}^{c}d{\rm{\Omega }}\\ &&\,=\,\displaystyle {\int }_{{{\rm{\Gamma }}}_{t}}{\bar{t}}_{i}\delta {u}_{i}^{c}d{\rm{\Gamma }}-\displaystyle {\int }_{{\rm{\Omega }}}[{c}_{ijkl}{\varepsilon }_{kl}^{c}({\bf{x}},\tau )-{e}_{kij}{E}_{k}({\bf{x}},\tau )\\ &&\,\,-\,{\gamma }_{ij}\theta ({\bf{x}},\tau )]\delta {\varepsilon }_{ij}^{c}d{\rm{\Omega }}=0,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\displaystyle {\int }_{{\rm{\Omega }}}[{\varepsilon }_{kl}({\bf{x}},\tau )-{l}^{2}{\varepsilon }_{kl,mm}({\bf{x}},\tau )-{\varepsilon }_{kl}^{c}({\bf{x}},\tau )]\delta {\varepsilon }_{kl}d{\rm{\Omega }}\\ &&\,=\,-\displaystyle {\int }_{{{\rm{\Gamma }}}_{\varepsilon }}{l}^{2}{n}_{m}{\varepsilon }_{kl,m}({\bf{x}},\tau )\delta {\varepsilon }_{kl}d{\rm{\Gamma }}\\ &&\,+\,\displaystyle {\int }_{{\rm{\Omega }}}{\varepsilon }_{kl}({\bf{x}},\tau )\delta {\varepsilon }_{kl}d{\rm{\Omega }}+\displaystyle {\int }_{{\rm{\Omega }}}{l}^{2}{\varepsilon }_{kl,m}({\bf{x}},\tau )\delta {\varepsilon }_{kl,m}d{\rm{\Omega }}\\ &&\,-\,\displaystyle {\int }_{{\rm{\Omega }}}{\varepsilon }_{kl}^{c}({\bf{x}},\tau )\delta {\varepsilon }_{kl}d{\rm{\Omega }}=0.\end{eqnarray} \tag{ 30 }$$

The used test functions are such that ${\delta \phi | }_{\partial {\rm{\Omega }}-{{\rm{\Gamma }}}_{Q}}=0,$ ${\delta {u}_{i}^{c}| }_{\partial {\rm{\Omega }}-{{\rm{\Gamma }}}_{t}}=0$ and ${\delta {\varepsilon }_{ij}| }_{\partial {\rm{\Omega }}-{{\rm{\Gamma }}}_{\varepsilon }}=0.$

The primary field variables in the above weak forms can be approximated over each standard finite element V^e in the FEM. The fields are expressed in terms of their respective nodal values and the corresponding shape function N^(a)(ξ₁, ξ₂, ξ₃) as

$$\begin{eqnarray}&&\theta \approx \theta {| }_{{V}^{e}}=\displaystyle \sum _{a=1}^{n}{\theta }^{(a)}{N}^{(a)}({\xi }_{1},{\xi }_{2},{\xi }_{3}),\\ &&\phi \approx \phi {| }_{{V}^{e}}=\displaystyle \sum _{a=1}^{n}{\phi }^{(a)}{N}^{(a)}({\xi }_{1},{\xi }_{2},{\xi }_{3}),\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{u}_{i}^{c}\approx {{u}_{i}^{c}| }_{{V}^{e}}=\displaystyle \sum _{a=1}^{n}{u}_{i}^{c(a)}{N}^{(a)}({\xi }_{1},{\xi }_{2},{\xi }_{3});\\ &&{\varepsilon }_{ij}\approx {{\varepsilon }_{ij}| }_{{V}^{e}}=\displaystyle \sum _{a=1}^{n}{\varepsilon }_{ij}^{(a)}{N}^{(a)}({\xi }_{1},{\xi }_{2},{\xi }_{3}),\end{eqnarray} \tag{ 32 }$$

where the superscripts (a) denote the global number of the ath node of V^e, (a = 1, 2, ..., n). The gradients in the Cartesian coordinate system are expressed in terms of the derivatives with respect to the intrinsic coordinates as

$$\begin{eqnarray}&&{\left.\left\{\begin{array}{c}\partial /\partial {x}_{1}\\ \partial /\partial {x}_{2}\\ \partial /\partial {x}_{3}\end{array}\right\}\right|}_{{V}^{e}}={{\bf{Y}}}^{e}\left\{\begin{array}{c}\partial /\partial {\xi }_{1}\\ \partial /\partial {\xi }_{2}\\ \partial /\partial {\xi }_{3}\end{array}\right\},\,{{\bf{Y}}}^{e}={({{\bf{J}}}^{e})}^{-1},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}{{\bf{J}}}^{e}={\left.\left[\begin{array}{ccc}\displaystyle \frac{\partial {x}_{1}}{\partial {\xi }_{1}} & \displaystyle \frac{\partial {x}_{2}}{\partial {\xi }_{1}} & \displaystyle \frac{\partial {x}_{3}}{\partial {\xi }_{1}}\\ \displaystyle \frac{\partial {x}_{1}}{\partial {\xi }_{2}} & \displaystyle \frac{\partial {x}_{2}}{\partial {\xi }_{2}} & \displaystyle \frac{\partial {x}_{3}}{\partial {\xi }_{2}}\\ \displaystyle \frac{\partial {x}_{1}}{\partial {\xi }_{3}} & \displaystyle \frac{\partial {x}_{2}}{\partial {\xi }_{3}} & \displaystyle \frac{\partial {x}_{3}}{\partial {\xi }_{3}}\end{array}\right]\right|}_{{V}^{e}}\\ =\displaystyle \sum _{a=1}^{n}\left[\begin{array}{ccc}{x}_{1}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{1}} & {x}_{2}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{1}} & {x}_{3}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{1}}\\ {x}_{1}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{2}} & {x}_{2}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{2}} & {x}_{3}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{2}}\\ {x}_{1}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{3}} & {x}_{2}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{3}} & {x}_{3}^{(a)}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{3}}\end{array}\right].\end{eqnarray} \tag{ 34 }$$

Hence,

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial f}{\partial {x}_{i}}\right|}_{{V}^{e}}=\displaystyle \sum _{a=1}^{n}{f}^{(a)}{b}_{i}^{e(a)}({\xi }_{1},{\xi }_{2},{\xi }_{3});\\ &&{b}_{i}^{e(a)}({\xi }_{1},{\xi }_{2},{\xi }_{3})={Y}_{i1}^{e}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{1}}+{Y}_{i2}^{e}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{2}}\\ &&\,\,+\,{Y}_{i3}^{e}\displaystyle \frac{\partial {N}^{(a)}}{\partial {\xi }_{3}};\end{eqnarray} \tag{ 35 }$$

and ${dV| }_{{V}^{e}}=| \det ({{\bf{J}}}^{e})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}.$

It is convenient to introduce the following column vectors:

$$\begin{eqnarray}&&{{\boldsymbol{\varepsilon }}}^{c}=\left\{\begin{array}{c}{\varepsilon }_{11}^{c}\\ {\varepsilon }_{22}^{c}\\ {\varepsilon }_{33}^{c}\\ 2{\varepsilon }_{23}^{c}\\ 2{\varepsilon }_{13}^{c}\\ 2{\varepsilon }_{12}^{c}\end{array}\right\}\,=\,\left[\begin{array}{ccc}{\hat{b}}_{1} & 0 & 0\\ 0 & {\hat{b}}_{2} & 0\\ 0 & 0 & {\hat{b}}_{3}\\ 0 & {\hat{b}}_{3} & {\hat{b}}_{2}\\ {\hat{b}}_{3} & 0 & {\hat{b}}_{1}\\ {\hat{b}}_{2} & {\hat{b}}_{1} & 0\end{array}\right]\left\{\begin{array}{c}{u}_{1}^{c}\\ {u}_{2}^{c}\\ {u}_{3}^{c}\end{array}\right\};\\ &&{\bf{E}}=\left\{\begin{array}{c}{E}_{1}\\ {E}_{2}\\ {E}_{3}\end{array}\right\}\,=\,-\left\{\begin{array}{c}{\hat{b}}_{1}\\ {\hat{b}}_{2}\\ {\hat{b}}_{3}\end{array}\right\}\phi ;\,{\boldsymbol{\beta }}=\left\{\begin{array}{c}{\beta }_{1}\\ {\beta }_{2}\\ {\beta }_{3}\end{array}\right\}\,=\,\left\{\begin{array}{c}{\hat{b}}_{1}\\ {\hat{b}}_{2}\\ {\hat{b}}_{3}\end{array}\right\}\theta ;\\ &&{\boldsymbol{\varepsilon }}=\left\{\begin{array}{c}{\varepsilon }_{11}\\ {\varepsilon }_{22}\\ {\varepsilon }_{33}\\ 2{\varepsilon }_{23}\\ 2{\varepsilon }_{13}\\ 2{\varepsilon }_{12}\end{array}\right\}.\end{eqnarray} \tag{ 36 }$$

where ${\hat{b}}_{i}=\partial /\partial {x}_{i}.$ The approximations of the secondary fields on V^e can now be written as:

$$\begin{eqnarray}&&{{\boldsymbol{\varepsilon }}}^{c}\approx {{{\boldsymbol{\varepsilon }}}^{c}| }_{{V}^{e}}=\displaystyle \sum _{a=1}^{n}{{\bf{B}}}_{\varepsilon }^{e(a)}{{\bf{u}}}^{c(a)};\\ {{\bf{B}}}_{\varepsilon }^{e(a)}={\left[\begin{array}{cccccc}{b}_{1}^{e(a)} & 0 & 0 & 0 & {b}_{3}^{e(a)} & {b}_{2}^{e(a)}\\ 0 & {b}_{2}^{e(a)} & 0 & {b}_{3}^{e(a)} & 0 & {b}_{1}^{e(a)}\\ 0 & 0 & {b}_{3}^{e(a)} & {b}_{2}^{e(a)} & {b}_{1}^{e(a)} & 0\end{array}\right]}^{T},\\ &&{{\bf{u}}}^{c(a)}=\left\{\begin{array}{c}{u}_{1}^{c(a)}\\ {u}_{2}^{c(a)}\\ {u}_{3}^{c(a)}\end{array}\right\}\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{\bf{E}}\approx {{\bf{E}}| }_{{V}^{e}}=-\displaystyle \sum _{a=1}^{n}{{\bf{b}}}^{e(a)}{\phi }^{(a)};\\ {{\bf{b}}}^{e(a)}={\left[\begin{array}{ccc}{b}_{1}^{e(a)} & {b}_{2}^{e(a)} & {b}_{3}^{e(a)}\end{array}\right]}^{T}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{\boldsymbol{\beta }}=\left\{\begin{array}{c}{\beta }_{1}\\ {\beta }_{2}\\ {\beta }_{3}\end{array}\right\}\,\approx \,{\left.\left\{\begin{array}{c}{\theta }_{,1}\\ {\theta }_{,2}\\ {\theta }_{,3}\end{array}\right\}\right|}_{{V}^{e}}=\displaystyle \sum _{a=1}^{n}{{\bf{b}}}^{e(a)}{\theta }^{(a)},\\ &&{n}_{m}{\hat{b}}_{m}{\boldsymbol{\varepsilon }}\approx {{n}_{m}{\hat{b}}_{m}{\boldsymbol{\varepsilon }}| }_{{V}^{e}}=\displaystyle \sum _{a=1}^{n}{n}_{m}{b}_{m}^{e(a)}{{\boldsymbol{\varepsilon }}}^{(a)}\\ &&\,=\,\displaystyle \sum _{a=1}^{n}({{\bf{n}}}^{T}{{\bf{b}}}^{e(a)}){{\boldsymbol{\varepsilon }}}^{(a)},\,{{\bf{n}}}^{T}=[{n}_{1}\,{n}_{2}\,{n}_{3}].\end{eqnarray} \tag{ 39 }$$

To simplify the FEM equations further, we consider vanishing pyroelectric coefficient. The variations of the primary field variables can be also expressed in terms of the variations of the nodal values. Thus, the variational statements, (27)–(29), can be rewritten as

$$\begin{eqnarray}&&\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a,b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\delta {\theta }^{(a)}[{({{\bf{b}}}^{e(a)})}^{T}{\bf{K}}{{\bf{b}}}^{e(b)}{\theta }^{(b)}\\ &&\,\,+\,\rho c{N}^{(a)}{N}^{(b)}{\dot{\theta }}^{(b)}]| \det ({{\bf{J}}}^{e})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}\\ &&\,=\,\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a=1}^{n}\delta {\theta }^{(a)}\left(\displaystyle {\int }_{{{\rm{\Gamma }}}_{q}^{e}}\bar{q}{N}^{(a)}d{\rm{\Gamma }}\right),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a,b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\delta {\phi }^{(a)}{({{\bf{b}}}^{e(a)})}^{T}\{{{\bf{L}}}^{T}{{\boldsymbol{\varepsilon }}}^{(b)}{N}^{(b)}\\ &&\,\,+\,{\bf{A}}{{\bf{b}}}^{e(b)}{\phi }^{(b)}\}| \det ({{\bf{J}}}^{e})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}\\ &&\,=\,\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a=1}^{n}\delta {\phi }^{(a)}\displaystyle {\int }_{{{\rm{\Gamma }}}_{Q}^{e}}\bar{Q}{N}^{(a)}d{\rm{\Gamma }},\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a,b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}{(\delta {{\bf{u}}}^{c(a)})}^{T}[{({{\bf{B}}}_{\varepsilon }^{e(a)})}^{T}\\ &&\,\times \,({\bf{C}}{{\bf{B}}}_{\varepsilon }^{e(b)}{{\bf{u}}}^{c(b)}+{\bf{L}}{{\bf{b}}}^{e(b)}{\phi }^{(b)})]\\ &&\,\times \,| \det ({{\bf{J}}}^{e})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}\\ &&\,=\,\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a,b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}{(\delta {{\bf{u}}}^{c(a)})}^{T}{({{\bf{B}}}_{\varepsilon }^{e(a)})}^{T}\\ &&\,\times \,{\boldsymbol{\gamma }}{N}^{(b)}{\theta }^{(b)}| \det ({{\bf{J}}}^{e})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}\\ &&\,+\,\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a=1}^{n}{(\delta {{\bf{u}}}^{c(a)})}^{T}\displaystyle {\int }_{{{\rm{\Gamma }}}_{t}^{e}}{\bf{T}}{N}^{(a)}d{\rm{\Gamma }},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a,b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}{(\delta {{\boldsymbol{\varepsilon }}}^{(a)})}^{T}[({{\boldsymbol{\varepsilon }}}^{(b)}{N}^{(b)}-{{\bf{B}}}_{\varepsilon }^{e(b)}{{\bf{u}}}^{c(b)})\\ &&\,\times \,{N}^{(a)}+{l}^{2}{({{\bf{b}}}^{e(a)})}^{T}{{\bf{b}}}^{e(b)}{{\boldsymbol{\varepsilon }}}^{(b)}]| \det ({{\bf{J}}}^{e})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}\\ &&\,=\,{l}^{2}\displaystyle \sum _{e=1}^{m}\displaystyle \sum _{a,b=1}^{n}\displaystyle {\int }_{{{\rm{\Gamma }}}_{\varepsilon }^{e}}{N}^{(a)}{(\delta {{\boldsymbol{\varepsilon }}}^{(a)})}^{T}({{\bf{n}}}^{T}{{\bf{b}}}^{e(b)}){{\boldsymbol{\varepsilon }}}^{(b)}d{\rm{\Gamma }},\end{eqnarray} \tag{ 43 }$$

where C, L and A are matrices of elastic, piezoelectric and dielectric constants, respectively. Since the variational forms (40)–(43) are, respectively, valid for any arbitrary δθ^a, δϕ^a, δu^c(a) and δε^a satisfying the Dirichlet boundary conditions, the following system of ordinary differential equations and algebraic equations is obtained:

$$\begin{eqnarray}&&\displaystyle \sum _{{e}_{k}=1}^{m}\displaystyle \sum _{b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}[{({{\bf{b}}}^{{e}_{k}({a}_{k})})}^{T}{\bf{K}}{{\bf{b}}}^{{e}_{k}(b)}{\theta }^{(b)}\\ &&\,+\,\rho c{N}^{({a}_{k})}{N}^{(b)}{\dot{\theta }}^{(b)}]| \det ({{\bf{J}}}^{{e}_{k}})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}\\ &&\,=\,\displaystyle \sum _{{e}_{k}=1}^{m}\left(\displaystyle {\int }_{{{\rm{\Gamma }}}_{q}^{{e}_{k}}}\bar{q}{N}^{({a}_{k})}d{\rm{\Gamma }}\right),\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{{e}_{k}=1}^{m}\displaystyle \sum _{b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}{({{\bf{b}}}^{{e}_{k}({a}_{k})})}^{T}\{{{\bf{L}}}^{T}{{\boldsymbol{\varepsilon }}}^{(b)}{N}^{(b)}+{\bf{A}}{{\bf{b}}}^{e(b)}{\phi }^{(b)}\}\\ &&\,\times | \det ({{\bf{J}}}^{{e}_{k}})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}=\displaystyle \sum _{{e}_{k}=1}^{m}\displaystyle {\int }_{{{\rm{\Gamma }}}_{Q}^{{e}_{k}}}\bar{Q}{N}^{({a}_{k})}d{\rm{\Gamma }},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{{e}_{k}=1}^{m}\displaystyle \sum _{b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}[{({{\bf{B}}}_{\varepsilon }^{{e}_{k}({a}_{k})})}^{T}({\bf{C}}{{\bf{B}}}_{\varepsilon }^{e(b)}{{\bf{u}}}^{c(b)}+{\bf{L}}{{\bf{b}}}^{e(b)}{\phi }^{(b)})]\\ &&\,\times \,| \det ({{\bf{J}}}^{{e}_{k}})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}=\displaystyle \sum _{{e}_{k}=1}^{m}\displaystyle \sum _{b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}{({{\bf{B}}}_{\varepsilon }^{{e}_{k}({a}_{k})})}^{T}\\ &&\,\times \,{\boldsymbol{\gamma }}{N}^{(b)}{\theta }^{(b)}| \det ({{\bf{J}}}^{{e}_{k}})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}+\displaystyle \sum _{{e}_{k}=1}^{m}\displaystyle {\int }_{{{\rm{\Gamma }}}_{t}^{{e}_{k}}}{\bf{T}}{N}^{({a}_{k})}d{\rm{\Gamma }},\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{{e}_{k}=1}^{m}\displaystyle \sum _{b=1}^{n}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}\displaystyle {\int }_{-1}^{1}[({{\boldsymbol{\varepsilon }}}^{(b)}{N}^{(b)}-{{\bf{B}}}_{\varepsilon }^{e(b)}{{\bf{u}}}^{c(b)}){N}^{({a}_{k})}\\ &&\,+\,{l}^{2}{({{\bf{b}}}^{e({a}_{k})})}^{T}{{\bf{b}}}^{e(b)}{{\boldsymbol{\varepsilon }}}^{(b)}]| \det ({{\bf{J}}}^{{e}_{k}})| d{\xi }_{1}d{\xi }_{2}d{\xi }_{3}\\ &&\,=\,{l}^{2}\displaystyle \sum _{{e}_{k}=1}^{m}\displaystyle \sum _{b=1}^{n}\displaystyle {\int }_{{{\rm{\Gamma }}}_{\varepsilon }^{{e}_{k}}}{N}^{({a}_{k})}({{\bf{n}}}^{T}{{\bf{b}}}^{e(b)}){{\boldsymbol{\varepsilon }}}^{(b)}d{\rm{\Gamma }},\\ &&(k=1,\,2,\,\mathrm{...},\,K)\end{eqnarray} \tag{ 47 }$$

where e_k denotes the finite element ${V}^{{e}_{k}}$ involving the node x^k, and a_k is the local number of the node x^k on ${V}^{{e}_{k}}.$

For a periodic distribution of QDs in a 3D piezoelectric matrix we consider a RVE with prescribed temperature on its upper surface (figure 1 (left)). The RVE is assumed to have a cubic geometry with side length of 40 nm, and a cubic QD with side length of a = 4 nm is embedded in it. The top surface of the InAs QD is 4 nm from the top surface of the matrix. Material properties of the InAs dot are [56]

$$\begin{eqnarray*}&&{c}_{11}=83.3\,{\rm{GPa}},{c}_{12}=45.3\,{\rm{GPa}},{c}_{44}=39.6\,{\rm{GPa}},\\ &&{e}_{14}=-0.045\,{\rm{C}}\,{{\rm{m}}}^{-{\rm{2}}},{h}_{11}=0.1107\times {10}^{-9}\,{{\rm{C}}}^{{\rm{2}}}\,{{\rm{N}}}^{-1}\,{{\rm{m}}}^{-2},\\ &&\alpha =6.1\times {10}^{-6}\,{\deg }^{-1},c=250\,{\rm{J}}\,{{\rm{kg}}}^{-1}\,{\deg }^{-1},\\ &&{\kappa }_{11}=46\,{\rm{W}}\,{{\rm{m}}}^{-1}\,{\deg }^{-1},\rho =5.67\times {10}^{3}\,{\rm{kg}}\,{{\rm{m}}}^{-3}.\end{eqnarray*}$$

Matrix (substrate) is made of GaAs, whose properties are

$$\begin{eqnarray*}&&{c}_{11}=118.8\,\mathrm{GPa},{c}_{12}=54\,{\rm{GPa}},{c}_{44}=59.4\,{\rm{GPa}},\\ &&{e}_{14}=-0.16\,{\rm{C}}\,{{\rm{m}}}^{-{\rm{2}}},{h}_{11}=0.1346\times {10}^{-9}\,{{\rm{C}}}^{{\rm{2}}}\,{{\rm{N}}}^{-1}\,{{\rm{m}}}^{-{\rm{2}}},\\ &&\alpha =5.1\times {10}^{-6}\,{\deg }^{-1},c=350\,{\rm{J}}\,{{\rm{k}}{\rm{g}}}^{-1}\,{{\rm{\deg }}}^{-1},\\ &&{\kappa }_{11}=0.044\,{\rm{W}}\,{{\rm{m}}}^{-1}\,{\deg }^{-1},\rho =5.32\times {10}^{3}\,{\rm{kg}}\,{{\rm{m}}}^{-3}.\end{eqnarray*}$$

It follows from the lattice constants that the eigenstrains are ${\varepsilon }_{11}^{* }={\varepsilon }_{22}^{* }={\varepsilon }_{33}^{* }=0.07.$ The pyroelectric coefficient p₁ is not considered in the numerical analyses since this quantity is small and would have only very small influence on the results.

The four side faces of the matrix are fixed along their normal direction, the bottom face is fixed along its normal direction, and the upper surface is free of tractions. Lateral sides and bottom side of the matrix cube are thermally isolated. Temperature is prescribed on the top. All surfaces have vanishing normal component of electric displacement, except for the bottom side, where a vanishing electric potential is prescribed. Then, one can write the following boundary conditions for the InAs/GaAs electronic structure (figure 1 (left)):

On surfaces BCGF and ADHE: u₁ = 0, t₂ = t₃ = 0, D₁ = 0, q = 0;

ABFE and DCGH: u₂ = 0, t₁ = t₃ = 0, D₂ = 0, q = 0;

ABCD: u₃ = 0, t₁ = t₂ = 0, ϕ = 0, q = 0;

EFGH: t₃ = 0, t₁ = t₂ = 0, D₃ = 0, θ = θ₀.

On the interface between InAs and GaAs, the continuity of displacements, electric potential and temperature has to be satisfied, namely,

$$\begin{eqnarray}&&{u}_{i}^{cm}={u}_{i}^{c\,{\rm{Q}}{\rm{D}}},\,{\phi }^{m}={\phi }^{{\rm{Q}}{\rm{D}}},\,{\theta }^{m}={\theta }^{{\rm{Q}}{\rm{D}}}\end{eqnarray} \tag{ 48 }$$

as well as the reciprocity of traction vector, electric displacement and heat flux

$$\begin{eqnarray}&&{T}_{i}^{m}+{T}_{i}^{{\rm{Q}}{\rm{D}}}=0,\,{{\rm{D}}}_{i}^{m}+{{\rm{D}}}_{i}^{{\rm{Q}}{\rm{D}}}=0,\,{q}^{m}+{q}^{{\rm{Q}}{\rm{D}}}=0.\end{eqnarray} \tag{ 49 }$$

Vanishing natural boundary conditions for the Helmholtz equation (19) are considered on all surfaces (see equation (25)).

A 3D COMSOL model was created by constructing the geometry and meshing it using 195 696 SOLID elements corresponding to 339 117 nodes as shown in figure 1 (right). COMSOL Multiphysics allows for implementing the governing equations (27)–(29) in the Weak-Form module. Computer code results were verified by ANSYS when the structural material parameter is vanishing (classical theory).

Majdoub et al [57] published material parameters for the enhanced size-dependent piezoelectricity and elasticity in nanostructures, and it follows that l = 0.5 × 10⁻⁹ m. However, in our numerical analyses we considered various values of the intrinsic material characteristic (size effect parameter), l, to illustrate its influence on strains. This material parameter represents micro-structural size, and therefore it has to be smaller than the size of the QD, a. As such, the admissible dimensionless ratio l/a should be restricted to the range from 0 to 1.

4.1. Code validation

In this section, numerical results obtained using COMSOL are compared with the ones obtained by ANSYS for classical theory, l = 0. For comparison, we have selected the applied thermal load to be 500 °C. The variation of the strain component ε₁₁ along X₁-direction (passing the QD center) and electric potential ϕ along X₃ (aligned with the edge of the QD at X₁ = 18 nm and X₂ = 22 nm) is presented in figure 2. One can observe excellent agreement between COMSOL and ANSYS results with error less than 1% for both the strain and electric potential.

Zoom In Zoom Out Reset image size

4.2. Static analysis of cubic QD in a cubic matrix

The electric potential contour plots in the X₁X₂- and X₂X₃-planes are shown in figure 3 where the prescribed temperature on the top surface of the matrix is 250 °C under stationary conditions. Electric potential on the X₁X₂ plane with X₃ = 32 nm is similar to that in figure 3 (bottom). Figure 4 shows the contour plots of the electric field intensity E_X1 at X₁X₂- , X₁X₃- and X₂X₃-planes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It can be observed from figure 3 that due to the cubic symmetry and the fact that the QD center is away from the top surface, the electric potential distributions on the X₁X₂- and X₂X₃-planes are very similar. The highest concentration of the electric potential is near the corners of the QD cube and is anti-symmetric. However, the electric field intensity, which is the derivative of the electric potential, on the corresponding planes is substantially different. This is demonstrated clearly in figure 4 for the electric field intensity E_X1 at X₁X₂- , X₂X₃- and X₁X₃-planes where one can observe field concentrations at the four corners on the X₁X₂- and X₁X₃-planes but not on the X₂X₃-plane. This strange feature is due to the fact that in the first two cases the electric field intensity E_X1 is obtained by taking the derivative of the potential with respect to X₁, which is within these planes, whilst for the last case (in the X₂X₃-plane), the derivative is out of plane.

Contour plots for strain components ε₃₃ and ε₁₂ on the X₁X₂-plane with X₃ = 34 nm are presented in figure 5. It is clearly observed that while the normal strain is concentrated within the QD region, the shear strains show field concentration at the corners.

Zoom In Zoom Out Reset image size

Next, the influence of temperature prescribed on the upper surface of the QD electronic structure EFGH in figure 1 is investigated. Figures 6 and 7 show the variations of total normal strains along the X₁- and X₃-directions, respectively, passing the QD center (X₁ = X₂ = 20 nm, X₃ = 34 nm) at various temperatures 0 °C, 250 °C, and 500 °C with classical elasticity (l = 0) (left) and gradient elasticity (l = 5 × 10⁻¹⁰) (right). It can be seen that the effect of temperature on the electromechanical quantities of GaAs/InAS QD system is not significant. A similar phenomenon was observed by Patil and Melnik [28] for cylindrical GaN/AlN QDs. The magnitude of all strain components decreases with increasing distance away from the QD. When l increases, the effect of temperature on strains increases as can be seen in figure 7. Under a given temperature, however, the magnitude of these normal strains can be significantly decreased with increasing material parameter length l, which is further demonstrated in figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 (upper left and right) shows the influence of the material characteristic length, l, on the normal strains ε₁₁ and ε₃₃ in figures 6 and 7 at a constant temperature (250 °C). As the size effect parameter, l, increases, the magnitude of the strains decreases significantly and therefore the corresponding system becomes stiffer. Figure 8 (lower) illustrates the variation of the maximum strain ε₁₁ with l/a ratio. As the QD size gets smaller with respect to the parameter l, the strains decrease. This provides a new opportunity for designing stiffer QD nanostructure by decreasing the size of the QDs to the level that experiences this size effect. Furthermore, size effect eliminates the sharp transitions of the strains at the interfaces of the boundaries, reducing the field concentration, and thus making the structure safer.

4.3. Transient analysis

To consider a transient thermal load, a 500 °C temperature shock load on the upper surface of the matrix is prescribed. The variations of strain components ε₁₁ and ε₃₃ along the X₁- and X₃-directions, respectively, are presented in figures 9 and 10. The line l = 0 corresponds to classical elasticity at time τ, under thermal shock load with τ = 7 × 10⁻¹³ s. Figure 11 shows the time evolution of maximum strains ε₁₁ and ε₃₃ at l = 5 × 10⁻¹⁰. Although steady state temperature proved to be insignificant on the strain distribution of the QD system, as was shown in section 4.2, when the system is subjected to transient or cyclic thermal loads, the induced strains vary substantially with time. On the other hand, increasing the size effect parameter is always stiffening the system and reducing the induced strains. These trends need to be carefully considered when the nanostructure is under thermal dynamic loads.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To assess the effect of a thermal shock load on a QD embedded in a piezoelectric substrate, a 3D dynamic thermo-piezoelectric finite element model has been developed. The model is a RVE that is used to simulate periodic QD array in the system and analyze the gradient elastic and piezoelectric fields induced by stationary and transient thermal loads. It was observed that the largest magnitudes of the induced field quantities, such as the induced strain fields, are within or near the QD regions, indicating that the effect of the lattice mismatch dominates. It was also found that as the size of QDs in the system gets smaller, the size effect parameter makes the system stiffer and reduces the induced strains. In addition, dynamic thermal loads proved to have a significant effect on the induced strains and the QD overall performance. Accordingly, the reported observations should be considered when designing strained QD nanostructures subjected to cyclic electrical and thermal loads.

The authors gratefully acknowledge the supports by the Slovak Science and Technology Assistance Agency APVV-14-0216 and the Slovak Grant Agency VEGA-1/0145/17.