Deformations of Single-Crystal Silicon Circular Plate: Theory and Experiment

1. Introduction

Nowadays, silicon membrane technology is of a great interest for modern electronics, since it opens the way for the creation of sensitive elements and supporting structures of a number of promising micro/nanoelectromechanical system (MEMS/NEMS)-based devices, such as accelerometers [1,2], pressure sensors [3,4], gas flow meters [5], microbolometers [6,7], digital micromirrors [8,9], resonators [10], microactuators [11], etc. The main task of this technology is to fabricate high-quality thin-film membranes with specified thermo-mechanical properties that ensure appropriate operating characteristics of MEMS/NEMS devices, which, depending on its purpose, is achieved by using various methods: 1) the formation of stress-compensated (or stress-free) multilayer films for thermal sensors by combining materials with compressive (SiO₂, TEOS) and tensile (Si₃N₄, Fe₆₅Co₃₅) stress [5,12,13], 2) selection of layer materials in optical and infrared (IR) membrane sensors with effective absorption of radiation in a given wavelength range [14,15], 3) fabrication of ultrathin extreme ultraviolet (EUV) transparent pellicles to protect the photomask from various contaminants [16], 4) synthesis of new compositions of highly reflective multilayer X-ray mirrors with a roughness below 1 nm [17], 5) use of silicon-on-insulator (SOI) structures with a functional layer based on single-crystalline highly-doped silicon (Si) as a field-emission electrode material in nanoscale vacuum channel transistors [18] or as a thermocouple material with an increased Seebeck coefficient in microbolometer IR detector arrays [19], etc.

To implement the above task in all the described methods, it is important to know the laws of evolution of mechanical stresses in the membrane structure depending on the conditions of its formation - technological regimes (operation temperature, growth rate, elemental composition, etc.), history of the fabrication and mechanical properties of the underlying layers, as well as the physical parameters of the crystalline substrate [20,21,22]. In particular, one of these parameters that influences the residual mechanical stresses and the resulting deformation of the membrane under the action of an applied external (mechanical or electromagnetic) force is the crystallographic orientation of the substrate or its crystalline anisotropy. As applied to SOI structures, this parameter not only specifies the mechanical behaviour of single-crystal Si membranes [23], but also determines the functional (thermal and electrical) properties of the electronic components integrated in the membrane device layer [24]. At the same time, there is no analytical model that allows one to reliably predict the effect of crystalline anisotropy of the substrate on the nature of the distribution of internal mechanical fields in the presence of an applied load.

In the present work, we compare the bending profiles of a single-crystal Si membrane of the SOI structure with crystallographic orientation <100>, which are measured in the experiment at various pressure (in the range from 0 to about 300 kPa) and calculated within the framework of non-linear elasticity theory based on iterative solving three-dimensional Lamé equations with non-linear RHS part. Based on this comparison, the great influence of the initial tension in the single-crystal Si membrane on the deformation of the SOI substrate were established.

2. Materials and methods

2.1. Design and technology of the fabrication of SOI-based thin-film Si membranes

A description of the process flow for the fabrication of a single-crystalline silicon (Si) membrane formed on the basis of a silicon-on-insulator (SOI) structure is presented in Figure 1. The initial structure of the SOI substrate (Figure 1a) consists of the following sequence of layers: top (thin) Si layer with a thickness of 2.86 µm, from which a single-crystal membrane is formed, an intermediate dielectric layer of SiO₂ with a thickness of 1 µm, and the bottom (thick) Si layer, the thickness of which is in accordance with the total SOI thickness (608 µm) is 604.14 µm.

In the first step of the fabrication process, a mask is formed on the SOI substrate for etching silicon (Si) to a depth comparable to its thickness (608 µm). A photoresist layer with a thickness of several microns is etched off before etching of a SOI substrate with a thickness of 608 µm occurs, and therefore the standard photoresist layer is replaced by an aluminum (Al) layer, which, compared to the standard photoresist, has a greater selectivity of etching to silicon (Si). Thus, a layer of aluminum (Al) with a thickness of 1.1 µm is formed on the reverse side of the SOI substrate by magnetron sputtering technique, after which it is covered with a polymethyl methacrylate (PMMA) photoresist layer using the spin coating method for subsequent contact photolithography (Figure 1b). The photoresist is then exposed through a mask and developed, after which the aluminum (Al) layer in the exposed area is removed by wet etching in a specially selected etchant consisting of 80% H₃PO₄, 5% HNO₃, 5% HAc and 10% H₂O (Figure 1c). After removing the photoresist layer, the next step is plasma-chemical etching of Si from the back side of the SOI substrate (Bosch process), where a layer of silicon oxide (SiO₂) serves as a etching-stop layer. The last step is to remove the remaining SiO2 layer in hydrofluoric acid (HF) vapor, thereby freeing the single-crystalline silicon (Si) membrane (Figure 1d). To fabricate a single-crystalline Si membrane with a round shape with a thickness of 2.86 µm and a diameter of 0.936 mm, a BSOI 600-1.5-7 substrate with a diameter of 150 mm and crystallographic orientation (100) was used. To experimentally demonstrate the change in the profile of a thin-film Si membrane under the influence of applied external pressure, several of the most characteristic representatives were selected from the original SOI substrate with membrane structures manufactured according to the fabrication process flow (see Fig. 1).

2.2. Experimental setup for measuring the profile of fabricated Si membranes at non-zero pressure

The experimental setup for measuring the mechanical properties of thin-film structures, which is illustrated in Figure 2, consists of the following components: a pipe line with the ability to supply air or any other gas with a pressure in the range of up to 7 atm., a receiver cylinder with an analog pressure gauge, a filter regulator MC104- D00, pressure regulator ER104-5PAP, Arduino Mega 2560 board, personal computer (PC), Veeco Wyko NT9300 profilometer, vibration isolation platform.

The receiver cylinder acts as a buffer tank to smooth out gas pressure differences and ensure the overall stability of the measuring system. The supplied air enters the filter regulator MC104-D00 with a filter element made of high-density polyethylene (HDPE) membrane with a filtration fineness of 25 µm according to the ISO 8573-1:2010 standard. The pressure regulator is responsible for maintaining a constant and controlled level of excess pressure. Several pressure gauges are located at various points on the bench to provide real-time monitoring of pressure levels. There are also built-in safety valves to automatically relieve excess pressure in the event of a diaphragm rupture. Software for the Arduino Mega 2560 board provides control of the pressure applied to the membrane. The accuracy of overpressure supply is 0.01 atm. The accuracy of measuring the relief (deflection) of the membrane along the vertical axis is no worse than 5 nm. A schematic diagram of the experimental setup is presented in Figure 3.

Using the experimental setup described above, the deflection profile of the fabricated silicon (Si) membranes with a diameter of 0.936 mm was measured. Figure 4 illustrates the results of measurements of the shape of two samples of fabricated single-crystal Si membrane under applied pressure, obtained using the experimental setup shown in Figure 2. As can be seen from the figure, the vertical coordinate of the membrane profile in the center changes nonlinearly with pressure, while the profile remains symmetrical relative to the center.

3. Theoretical statement of the problem

In this part we construct a closed solution of boundary value problem for linear elastic anisotropic circular plate. We will seek the solution in the class of square-integrable vector-functions, defined over the domain, occupied by the plate. It is known that under certain boundary conditions the separation of variables can be used. We’ll take advantage of this opportunity. To this end the specific boundary conditions of roller contact type must be specified on the lateral surface of the plate. This limitation does not seem burdensome, since it corresponds to the conditions under which the experiment is provided (as will be shown below, they can be taken into account in the model very accurately by choosing a specific distribution of volumetric force). These considerations allow the solution to be represented in the form of an expansion over a system of vector-valued eigenfunctions, generated by the auxiliary problem of Bessel type.

First let’s make some notes about the representation of coordinates and fields. We will represent all kinematic fields and stresses in cylindrical coordinates $(r,\phi ,z)$:

$$r=\sqrt{x^2+y^2},\phi =arctan\left(y,x\right),z=z,$$

and corresponding local physical bases $({\mathit{e}}_r,{\mathit{e}}_{\phi },{\mathit{e}}_z)$:

$${\mathit{e}}_r=\mathit{i}cos\phi +\mathit{j}sin\phi ,{\mathit{e}}_{\phi }=-\mathit{i}sin\phi +\mathit{j}cos\phi ,{\mathit{e}}_z=\mathit{k}.$$

Here $(x,y,z)$ denote standard rectangular coordinates, while $(\mathit{i},\mathit{j},\mathit{k})$ stand for the elements of Cartesian orthonormal basis. Since the problem at the first stage is considered in a linear formulation (geometric non-linearity will be taken into account in section 6), we will not distinguish the reference and actual shapes of the body, assuming that both coincide with the region

$$\mathcal{D}=\left\{\left(r,\phi ,z\right)\in {\mathbb{R}}^3:r\in \left(0,R\right),\phi \in \left(0,2\pi \right),z\in \left(0,h\right)\right\}.$$

Due to the peculiarities of the experimental data, we will need to study only axisymmetric deformation. In this regard, we define the general representation of the desired displacement field as

$$\mathit{u}=u{\mathit{e}}_r+w{\mathit{e}}_z,$$

where $u,w$ are the functions to be defined:

$$u:\left(0,R\right)\times \left(0,h\right)\ni (r,z)\mapsto u(r,z)\in \mathbb{R},w:\left(0,R\right)\times \left(0,h\right)\ni (r,z)\mapsto w(r,z)\in \mathbb{R}.$$

The design diagram of the plate is shown in Figure 5.

Let’s move on to the discussion of boundary conditions that meet the method of sample fixing in an experiment. Suppose that the lateral boundary $\Gamma$ of the cylindrical plate is free from tangential stresses and fixed in radial direction. From the viewpoint of structural mechanics these conditions correspond to ”roller contact” of cylindrical part of the plate boundary with absolutely rigid tube without friction. Some restrictions, that has to be taken into account in conditions on the cylindrical part of the boundary, give complete freedom to choose conditions on the flat part, i.e., on the bases of the cylinder (plate face surfaces). In order to reproduce the experimental conditions, we assume that the bases ${\Omega }_1,{\Omega }_2$ are not fixed and believe that a volumetric distributed force $\mathit{K}$ is applied on the interior of the cylinder. Note that, under given above boundary conditions, the plate can move as a rigid body along its axis. In order to exclude such motion, we assume that, besides axial symmetry, the external loads (both surface and volumetric) are self-balanced, that is, their net force and net moment vanish, i.e.

$$\underset{0}{\overset{h}{\int }}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{R}{\int }}\mathit{K}rdrd\phi dz=0.$$

(1)

We especially note one of the particular distribution of volumetric force $\mathit{K}$

$$\begin{array}{cc}\hfill & \mathit{K}=\left[p-\frac{p{R}^2}{R^2-a^2}\theta \left(r-a\right)\right]{\mathit{e}}_z.\hfill \end{array}$$

Here $\theta (·)$ is Heaviside step function. The parameter a is chosen to be sufficiently close to the plate radius. These distributions correspond to an uniformly distributed volume force in the main area of the plate and distributed reaction in a narrow annular area. The latter integrally characterises the reaction of the support in experimental setup. It should be noted that, on the one hand, this formalization allows to remain within the framework of the boundary conditions discussed above. On the other hand, because it’s actually very difficult to experimentally identify in detail the specific distributions of reaction at the cylindrical part of sample boundary, they are still considered in integral form.

Now we turn to the formulation of the constitutive equation for an transversely isotropic material. Firstly note, that in the case of axial symmetry infinitesimal deformations are determined by the formulae

$${\epsilon }_{rr}=\frac{\partial u}{\partial r},{\epsilon }_{\phi \phi }=\frac{u}{r},{\epsilon }_{zz}=\frac{\partial w}{\partial z},{\epsilon }_{rz}={\epsilon }_{zr}=\frac{1}{2}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right).$$

(2)

These deformations cause stresses, defined by the fourth-rank elasticity tensor $\mathit{C}$:

$$\sigma =\mathit{C}:\epsilon ,$$

(3)

which in the case of a transversely isotropic material can be written in the Voigt form as the following matrix $6\times 6$:

$$\mathit{C}=\left(\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& 0\\ & c_{11}& c_{13}& 0& 0& 0\\ & *& c_{33}& 0& 0& 0\\ & *& *& c_{44}& 0& 0\\ & *& *& *& c_{44}& 0\\ & *& *& *& *& \frac{1}{2}\left(c_{11}-c_{12}\right)\end{array}\right).$$

Thus one can obtain relation (3) in component form:

$$\begin{array}{cc}\hfill {\sigma }_{rr}& =c_{11}{\epsilon }_{rr}+c_{12}{\epsilon }_{\phi \phi }+c_{13}{\epsilon }_{zz},{\sigma }_{\phi \phi }=c_{12}{\epsilon }_{rr}+c_{11}{\epsilon }_{\phi \phi }+c_{13}{\epsilon }_{zz},\hfill \\ \hfill {\sigma }_{zz}& =c_{13}{\epsilon }_{rr}+c_{13}{\epsilon }_{\phi \phi }+c_{33}{\epsilon }_{zz},{\sigma }_{rz}={\sigma }_{zr}=2c_{44}{\epsilon }_{rz},\hfill \end{array}$$

where $c_{11},c_{12},c_{13},c_{33},c_{44}$ are material constants, which characterise transverse isotropy in most general form. The only axis of such anisotropy is supposed to be oriented along the axis of the cylindrical plate. For the convenience of subsequent computational analysis, we present the anisotropic elastic moduli as the result of a continuous morphism of isotropic elastic one:

$$\begin{array}{c}{c}_{11}\left(t\right)=\left(2\mu +\lambda \right)\left(1-t\right)+c_{11}t,c_{12}\left(t\right)=\lambda \left(1-t\right)+c_{12}t,c_{13}\left(t\right)=\lambda \left(1-t\right)+c_{13}t,\\ c_{33}\left(t\right)=\left(2\mu +\lambda \right)\left(1-t\right)+c_{33}t,c_{44}\left(t\right)=\mu \left(1-t\right)+c_{44}t,\end{array}$$

(4)

where $\mu ,\lambda$ are Lamé material constants of polycrystalline material, which correspond to those of crystalline material by Voigt formulae [25]:

$$\mu =\frac{1}{30}\left(7c_{11}-5c_{12}-4c_{13}+2c_{33}+12c_{44}\right),\lambda =\frac{1}{15}\left(c_{11}+5c_{12}+8c_{13}+c_{33}-4c_{44}\right).$$

(5)

Substituting (2), (3) into the equilibrium equations

$$\begin{array}{c}\hfill \frac{\partial {\sigma }_{rr}}{\partial r}+\frac{{\sigma }_{rr}-{\sigma }_{\phi \phi }}{r}+\frac{\partial {\sigma }_{rz}}{\partial z}=\mathit{K}·{\mathit{e}}_r,\frac{\partial {\sigma }_{zr}}{\partial r}+\frac{{\sigma }_{zr}}{r}+\frac{\partial {\sigma }_{zz}}{\partial z}=\mathit{K}·{\mathit{e}}_z.\end{array}$$

we arrive at the system of partial differential equations:

$$\left\{\begin{array}{c}{\displaystyle c_{11}\left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}-\frac{u}{r^2}\right)+c_{13}\frac{{\partial }^{2}w}{\partial z\partial r}+c_{44}\left(\frac{{\partial }^{2}u}{\partial z^2}+\frac{{\partial }^{2}w}{\partial r\partial z}\right)=\mathit{K}·{\mathit{e}}_r,}\hfill \\ {\displaystyle c_{13}\left(\frac{{\partial }^{2}u}{\partial r\partial z}+\frac{1}{r}\frac{\partial u}{\partial z}\right)+c_{33}\frac{{\partial }^{2}w}{\partial z^2}+c_{44}\left(\frac{{\partial }^{2}u}{\partial z\partial r}+\frac{1}{r}\frac{\partial u}{\partial z}+\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{r}\frac{\partial w}{\partial r}\right)=\mathit{K}·{\mathit{e}}_z.}\hfill \end{array}\right.$$

(6)

In light of the above the boundary conditions can be stated as follows:

$$\begin{array}{c}{\left.{\sigma }_{zz}\right|}_{z=0}={\left.c_{13}\left(\frac{\partial u}{\partial r}+\frac{u}{r}\right)+c_{33}\frac{\partial w}{\partial z}\right|}_{z=0}=0,{\left.{\sigma }_{zz}\right|}_{z=h}={\left.c_{13}\left(\frac{\partial u}{\partial r}+\frac{u}{r}\right)+c_{33}\frac{\partial w}{\partial z}\right|}_{z=h}=0,\\ {\left.{\sigma }_{rz}\right|}_{z=0}={\left.\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right|}_{z=0}=0,{\left.{\sigma }_{rz}\right|}_{z=h}={\left.\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right|}_{z=h}=0,\\ {\left.{\sigma }_{rz}\right|}_{r=1}={\left.\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right|}_{r=1}=0,{\left.u\right|}_{r=1}=0.\end{array}$$

(7)

The equations (6) and boundary conditions (7) together with requirement of self-balance (1) gives linear statement for the boundary-valued problem under consideration.

It is convenient to carry out further calculations in dimensionless variables (in the next line they are denoted with twiddle), which are related to their dimensional counterparts by formulae:

$$\tilde{R}=1,\tilde{h}=\frac{h}{R},\tilde{r}=\frac{r}{R},\tilde{z}=\frac{z}{R},\tilde{u}=\frac{z}{R},\tilde{w}=\frac{z}{R},{\tilde{c}}_{ij}=\frac{c_{ij}}{\mu },\tilde{\mathit{K}}=\frac{R}{\mu }\mathit{K}.$$

In what follows, we will omit the tilde sign, since the question of which formulation is used, dimensional or dimensionless, is easy to answer based on the context.

4. Solution of the linear boundary-value problem

Before proceeding to the solution, we note that so formulated problem allows for complete separation of spatial variables. Indeed, if one represent sought variables u, w in the form of multiplicative decomposition

$$u=f_1\left(r\right)g_1\left(z\right),w=f_2\left(r\right)g_2\left(z\right),$$

that meet the conditions

$$f_2^{\prime }=-\alpha f_1,f_1^{\prime }+\frac{f_1}{r}=\alpha f_2,$$

(8)

where $\alpha$ is an arbitrary constant, then the system of homogeneous equation, defined by LHS of (6) can be written in the form

$$\left\{\begin{array}{c}{\displaystyle -\left(f_1^{\prime \prime }+\frac{f_1^{\prime }}{r}-\frac{f_1}{r^2}\right)\frac{1}{f_1}=\frac{c_{44}}{c_{11}}\frac{g_1^{\prime \prime }}{g_1}-\alpha \frac{c_{13}+c_{44}}{c_{11}}\frac{g_2^{\prime }}{g_1}={\lambda }_1,}\hfill \\ {\displaystyle -\left(f_2^{\prime \prime }+\frac{f_2^{\prime }}{r}\right)\frac{1}{f_2}=\frac{c_{33}}{c_{44}}\frac{g_2^{\prime \prime }}{g_2}+\alpha \frac{c_{13}+c_{44}}{c_{44}}\frac{g_1^{\prime }}{g_2}={\lambda }_2.}\hfill \end{array}\right.$$

(9)

Here ${\lambda }_1,{\lambda }_2$ are constants of separation [26], and the prime means differentiation with respect to the variable on which a specific function depends. Actually, these constants turn out to be equal, which is easy to see from the relations, that can be obtained by differentiation of equations (8):

$$f_1^{\prime \prime }+\frac{f_1^{\prime }}{r}-\frac{f_1}{r^2}=-{\alpha }^2{f}_1,f_2^{\prime \prime }+\frac{f_2^{\prime }}{r}=-{\alpha }^2{f}_2.$$

Taking them into account, from (9) one gets

$${\lambda }_1={\lambda }_2={\alpha }^2.$$

Therefore, a homogeneous system of partial differential equations splits into two independent systems of ordinary differential equations

$$\begin{array}{cc}\hfill & \left\{{\displaystyle \begin{array}{c}{r}^2{f}_1^{\prime \prime }+r{f}_1^{\prime }+\left({\alpha }^2{r}^2-1\right)f_1=0,\hfill \\ r^2{f}_2^{\prime \prime }+r{f}_2^{\prime }+{\alpha }^2{r}^2{f}_2=0,\hfill \end{array}}\right.\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & \left\{{\displaystyle \begin{array}{c}{c}_{44}{g}_1^{\prime \prime }-{\alpha }^2{c}_{11}{g}_1-\alpha \left(c_{13}+c_{44}\right)g_2^{\prime }=0,\hfill \\ c_{33}{g}_2^{\prime \prime }-{\alpha }^2{c}_{44}{g}_2+\alpha \left(c_{13}+c_{44}\right)g_1^{\prime }=0.\hfill \end{array}}\right.\hfill \end{array}$$

(11)

Separating the variables in equations is only half a matter. To complete the matter, it is also necessary to separate the variables in the boundary conditions. But this is where the considerations, discussed above for the cylindrical part of the boundary, come in handy. Substituting (4) into LHS of relations (7) one obtains:

$${\left.c_{33}\frac{g_2^{\prime }}{g_1}+\frac{c_{13}}{f_2}\left(f_1^{\prime }+\frac{f_1}{r}\right)\right|}_{z=0,h}=0,{\left.\frac{g_1^{\prime }}{g_2}+\frac{f_2^{\prime }}{f_1}\right|}_{z=0,h}=0,{\left.f_1\right|}_{r=1}=0,{\left.f_2^{\prime }\right|}_{r=1}=0.$$

With account of (8) these equalities can be rewritten as:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\displaystyle {\left.f_1\right|}_{r=1}=0,}\\ {\left.f_2^{\prime }\right|}_{r=1}=0,\end{array}\right.\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\displaystyle {\left.c_{33}{g}_2^{\prime }+\alpha c_{13}{g}_1\right|}_{z=0,h}=0,}\\ {\left.g_1^{\prime }-\alpha g_2\right|}_{z=0,h}=0.\end{array}\right.\hfill \end{array}$$

(13)

Thus, we arrive at two independent boundary-value problems containing only ordinary differential equations (ODE). It is this fact, that allows us to further construct representations of the solution in terms of known elementary and special functions, since all of them are solutions of corresponding ODE’s.

Let’s now start constructing the solution of (10) – (13) itself. First, consider the auxiliary problem (10), (12). Note that, because the corresponding differential operator is self-conjugate, it can be classified as Sturm–Liouville problem on the interval $(0,1)$ [27]. It follows that all it`s non-trivial solutions constitute a basis in Hilbert space of vector-valued functions over this interval. The solution of the equations (10) (for $\alpha \ne 0$) can be given in the form of combinations of Bessel functions:

$$f_1=C_1{J}_1\left(\alpha r\right)+C_2{Y}_1\left(\alpha r\right),f_2=C_3{J}_0\left(\alpha r\right)+C_4{Y}_0\left(\alpha r\right),$$

where $C_1,\dots C_4$ and $\alpha$ are arbitrary constants, $J_n(·),Y_n(·)$ are Bessel functions of the first and second kind correspondingly [28]. Keeping in mind that the sought displacement functions has to be bounded at the pole ($r=0$), we set constants $C_2$, $C_4$ equal to zero. The two remaining boundary conditions (12) will be satisfied if any root of the transcendental equation

$$J_1\left(\alpha \right)=0$$

is taken as $\alpha$. In fact, each such a value is two-fold eigenvalue of the Sturm-Liouville problem, which corresponds to the two-dimensional eigensubspace. The bases in each subspace can be chosen by setting the constants $(C_2,C_4)$ equal to $(0,1)$ and $(1,0)$. It remains only to note that in the special case $\alpha =0$ the solution to the boundary value problem has the form $f_1=0,f_2=1$. Thus the eigenfunctions can be defined as countable sequence

$$\left(\begin{array}{c}0\\ 1\end{array}\right),\left(\begin{array}{c}{J}_1\left({\alpha }_{1}r\right)\\ 0\end{array}\right),\left(\begin{array}{c}0\\ J_0\left({\alpha }_{1}r\right)\end{array}\right),\left(\begin{array}{c}{J}_1\left({\alpha }_{2}r\right)\\ 0\end{array}\right),\left(\begin{array}{c}0\\ J_0\left({\alpha }_{2}r\right)\end{array}\right),\dots$$

where

$${\alpha }_n\in \left\{\alpha \in \mathbb{R}:J_1\left(\alpha \right)=0\right\},n\in \mathbb{N},$$

and any square-integrable vector-function can be represent (in weak form) as decomposition

$$\mathit{u}=\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}0\\ y_0\end{array}\right)+\sum _{n=1}^{\infty }\left(\begin{array}{c}{y}_{n1}{J}_1\left({\alpha }_{n}r\right)\\ y_{n2}{J}_0\left({\alpha }_{n}r\right)\end{array}\right)$$

with suitable sequence of Fourier coefficients ${\left(y_n\right)}_{n=1}^{\infty }$. To simplify calculations, it is advisable to use normalized eigenfunctions

$$\sqrt{2}\left(\begin{array}{c}0\\ 1\end{array}\right),\frac{\sqrt{2}}{J_0\left({\alpha }_1\right)}\left(\begin{array}{c}{J}_1\left({\alpha }_{1}r\right)\\ 0\end{array}\right),\frac{\sqrt{2}}{J_0\left({\alpha }_1\right)}\left(\begin{array}{c}0\\ J_0\left({\alpha }_{1}r\right)\end{array}\right),\frac{\sqrt{2}}{J_0\left({\alpha }_2\right)}\left(\begin{array}{c}{J}_1\left({\alpha }_{2}r\right)\\ 0\end{array}\right),\frac{\sqrt{2}}{J_0\left({\alpha }_2\right)}\left(\begin{array}{c}0\\ J_0\left({\alpha }_{2}r\right)\end{array}\right),\dots$$

The decomposition procedure can be formalized with orthogonal projectors onto eigensubspaces, represented in the following form:

$$P_0\left(\mathit{u}\right)=\sqrt{2}\underset{0}{\overset{1}{\int }}{u}_{2}rdr,P_{n1}\left(\mathit{u}\right)=\frac{\sqrt{2}}{J_0\left({\alpha }_n\right)}\underset{0}{\overset{1}{\int }}{u}_1{J}_1\left({\alpha }_{n}r\right)rdr,P_{n2}\left(\mathit{u}\right)=\frac{\sqrt{2}}{J_0\left({\alpha }_n\right)}\underset{0}{\overset{1}{\int }}{u}_2{J}_0\left({\alpha }_{n}r\right)rdr.$$

(14)

Now any square-integrable function can be represented by a sequence of Fourier coefficients, and vice versa, from the sequence of Fourier coefficients a vector function can be reconstructed (in the weak sense)

$$\begin{array}{cc}\hfill & \mathit{u}\sim \left(y_0,y_{11},y_{12},y_{21},y_{22},\dots \right)\in l^2,y_0=P_0\left(\mathit{u}\right),y_{n1}=P_{n1}\left(\mathit{u}\right),y_{n2}=P_{n2}\left(\mathit{u}\right);\hfill \\ \hfill & \left(y_0,y_{11},y_{12},y_{21},y_{22},\dots \right)\sim \widehat{\mathit{u}}=\left(\begin{array}{c}0\\ \sqrt{2}{y}_0\end{array}\right)+\sum _{n=1}^{\infty }\frac{\sqrt{2}}{J_0\left({\alpha }_{n}r\right)}\left(\begin{array}{c}{y}_{n1}{J}_1\left({\alpha }_{n}r\right)\\ y_{n2}{J}_0\left({\alpha }_{n}r\right)\end{array}\right)\in L^2.\hfill \end{array}$$

Here $l^2$ denotes the space of square summable sequences, while $L^2$ is the space of quadratically integrable (with weight r) two-component vector-functions, and $\widehat{\mathit{u}}$ stands for the result of successive mutual mappings between these spaces, which may differ from the original $\mathit{u}$ on a set of measure zero [27].

Note that direct numerical computation for a long sequences of transcendental equation roots ${\alpha }_n$ gives rise to a well-known computational problem (in a view of possible missing roots or incorrect determination of their multiplicity). In application to the problem under consideration it can be easily solved if one applies the asymptotic approximation for ${\alpha }_n$:

$${\alpha }_n\approx \frac{{(\pi +4n\pi )}^2{J}_2\left(\left(1/4+n\right)\pi \right)}{16J_1\left(\left(1/4+n\right)\pi \right)-4(\pi +4n\pi )J_0\left(\left(1/4+n\right)\pi \right)}.$$

The result of the action of projection operators (14) on the right-hand sides of the equations (6) will be denoted by the symbols

$${\mathcal{B}}_0=P_0\left(\mathit{K}\right),{\mathcal{A}}_n=P_{n1}\left(\mathit{K}\right),{\mathcal{A}}_n=P_{n1}\left(\mathit{K}\right).$$

For further analysis, we will need projections of two special cases of distributions of volumetric forces relative to the radial coordinate, namely

• stepwise self-balanced distribution taking a unit value on the interval $(0,a)$ and a constant value of the opposite sign on the interval $(a,1)$, i.e.

  $${\mathit{K}}_1=\left(\begin{array}{c}0\\ \left\{\begin{array}{cc}1\hfill & 0\le r\le a\hfill \\ \frac{a^2}{a^2-1}\hfill & a<r<1\hfill \end{array}\right.\end{array}\right);$$
• self-balanced distribution (see (1)), taking a unit value on the interval $(c,d)\subset (0,h)$ and a constant value of the opposite sign on the interval $(a,1)$, i.e.

  $${\mathit{K}}_2=\left(\begin{array}{c}0\\ \left\{\begin{array}{cc}1\hfill & c\le r\le d\hfill \\ \frac{d^2-c^2}{a^2-1}\hfill & a<r<1\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}\right).$$

The formulations of such distributions and the result of the action of projectors on them are given below:

$$P_0\left({\mathit{K}}_1\right)=P_{n1}\left({\mathit{K}}_1\right)=0,P_{n2}\left({\mathit{K}}_1\right)=\frac{\sqrt{2}a{J}_1\left(a\alpha \right)}{\alpha \left(1-a^2\right)J_0\left(\alpha \right)},$$

$$P_0\left({\mathit{K}}_2\right)=P_{n1}\left({\mathit{K}}_2\right)=0,P_{n2}\left({\mathit{K}}_2\right)=\frac{\sqrt{2}\left(a(d^2-c^2)J_1\left(a\alpha \right)+(1-a^2)\left(c{J}_1\left(c\alpha \right)-d{J}_1\left(d\alpha \right)\right)\right)}{\alpha (1-a^2)J_0\left(\alpha \right)}.$$

(15)

Passing to the limit $a\to 1$ in these relations, and taking into account the equality

$$\underset{a\to 1}{lim}\frac{J_1\left(a{j}_n^1\right)}{1-a^2}=\frac{j_n^1}{2}{J}_2\left(j_n^1\right),$$

where $j_n^1$ is a zero of $J_1(·)$, one can obtain the following expressions for projections, valid for the case $\alpha =j_n^1$:

$$P_{n2}\left({\mathit{K}}_1\right)=\frac{J_2\left(j_n^1\right)}{\sqrt{2}{J}_0\left(j_n^1\right)},P_{n2}\left({\mathit{K}}_2\right)=\frac{\left((d^2-c^2)J_2\left(j_n^1\right)+c{J}_1\left(c{j}_n^1\right)-d{J}_1\left(d{j}_n^1\right)\right)}{\sqrt{2}{J}_0\left(j_n^1\right)}.$$

This corresponds to the case when the reactive body forces are distributed in an infinitely thin boundary ring (and can be formally represented in terms of Dirac delta functions).

Using constructed above orthonormal basis, one can represent the sought displacement functions in the form of expansions

$$u=\sum _{n=1}^{\infty }{y}_{1n}{J}_1\left({\alpha }_{n}r\right),w=y_{20}+\sum _{n=1}^{\infty }{y}_{2n}{J}_0\left({\alpha }_{n}r\right),$$

(16)

where $y_{20}=y_{20}\left(z\right)$, $y_{1n}=y_{1n}\left(z\right)$, $y_{2n}=y_{2n}\left(z\right)$ are Fourier coefficients which themselves are functions, but of another spatial variable z. To find these coefficients, one has to substitute the formal expansion (16) into equations (6), boundary conditions (7) and then apply projectors (14) to the left and right sides of the resulting functional equations. This results in countable sequence of independent ordinary differential equations (with respect to z),

$$\left\{\begin{array}{c}{c}_{44}{y}_{1n}^{\prime \prime }-{\alpha }^2{c}_{11}{y}_{1n}-\alpha \left(c_{13}+c_{44}\right)y_{2n}^{\prime }={\mathcal{A}}_n,\hfill \\ c_{33}{y}_{2n}^{\prime \prime }-{\alpha }^2{c}_{44}{y}_{2n}+\alpha \left(c_{13}+c_{44}\right)y_{1n}^{\prime }={\mathcal{B}}_n,\hfill \end{array}\right.$$

(17)

and two-point boundary conditions

$${\left.g_{1n}^{\prime }-\alpha y_{2n}\right|}_{z=0,h}=0,{\left.c_{33}{y}_{2n}^{\prime }+\alpha c_{13}{y}_{1n}\right|}_{z=0,h}=0.$$

(18)

As above, we will separately consider the case $\alpha =0$

$$c_{33}{y}_{20}^{\prime \prime }={\mathcal{B}}_0,$$

(19)

$${\left.y_{20}^{\prime }\right|}_{z=0}=0,{\left.y_{20}^{\prime }\right|}_{z=h}=0.$$

One can see that something is wrong with this boundary-value problem. Indeed, the solution of differential equation (19) can be obtained in the form

$$y_{20}=\frac{1}{c_{33}}\underset{0}{\overset{z}{\int }}(z-\zeta ){\mathcal{B}}_0\left(\zeta \right)d\zeta +C_1+C_{2}z,$$

but from the boundary conditions it is impossible to determine the constant $C_1$. On the other hand, differentiating the resulting solution, we find

$$y_{20}^{\prime }=\frac{1}{c_{33}}\underset{0}{\overset{z}{\int }}{\mathcal{B}}_0\left(\zeta \right)d\zeta +C_2.$$

Therefore, to satisfy first boundary condition, one should set $C_2=0$ and the second boundary condition will be satisfied only if

$$\underset{0}{\overset{h}{\int }}{\mathcal{B}}_0\left(\zeta \right)d\zeta =0.$$

It easy to see that the condition for the existence of a solution is equivalent to the requirement for external fields to be self-balanced.

All we have to do is to obtain solutions for two-point boundary value problems with respect to expansion coefficients. These boundary-value problems are of the same type and are determined by ordinary differential equations with constant coefficients. To this end we firstly get the general solution of differential equation (17) in the form (to shorten the notation, we will omit the index n):

$$\mathit{y}={\mathit{y}}^{*}-\mathit{Y}.\left({\mathit{m}}_1{C}_1+{\mathit{m}}_2{C}_2\right).$$

(20)

Here $\mathit{Y}$ stands for the matrix containing fundamental system of solutions for homogeneous counterpart of equation (17):

$$\mathit{Y}=\left(\begin{array}{cccc}As{e}^{A\alpha z}& As{e}^{-A\alpha z}& Bs{e}^{B\alpha z}& Bs{e}^{-B\alpha z}\\ -g{e}^{A\alpha z}& g{e}^{-A\alpha z}& -j{e}^{B\alpha z}& j{e}^{-B\alpha z}\end{array}\right),$$

(21)

where A, B are the roots of corresponding characteristic equation,

$$A=\sqrt{\frac{r-2c_{13}{c}_{44}-\sqrt{r}\sqrt{r-4c_{44}s}}{2c_{44}{c}_{33}}},B=\sqrt{\frac{r-2c_{13}{c}_{44}+\sqrt{r}\sqrt{r-4c_{44}s}}{2c_{44}{c}_{33}}},$$

(22)

and ${\mathit{y}}^{*}$ denotes particular solution, obtained with Lagrange method:

$${\mathit{y}}^{*}=\frac{1}{2c_{11}{c}_{44}{c}_{33}\alpha sAB\left(A^2-B^2\right)}\mathit{Y}.\underset{0}{\overset{z}{\int }}\left(\begin{array}{c}B{e}^{-B\alpha \zeta }\left(Bj{c}_{33}\mathcal{A}\left(\zeta \right)+s{c}_{11}\mathcal{B}\left(\zeta \right)\right)\\ B{e}^{A\alpha \zeta }\left(s{c}_{11}\mathcal{B}\left(\zeta \right)-Aj{c}_{33}\mathcal{A}\left(\zeta \right)\right)\\ -A{e}^{-B\alpha \zeta }\left(Bg{c}_{33}\mathcal{A}\left(\zeta \right)+s{c}_{11}\mathcal{B}\left(\zeta \right)\right)\\ -A{e}^{B\alpha \zeta }\left(s{c}_{11}\mathcal{B}\left(\zeta \right)-Bg{c}_{33}\mathcal{A}\left(\zeta \right)\right)\end{array}\right)d\zeta .$$

(23)

Constants of integration, defined by the expression in parentheses (20), are chosen in such a way that the boundary conditions (18) are satisfied. They are given by the following expressions

$${\mathit{m}}_1=\frac{{\alpha }^3}{\Delta }\left(\begin{array}{c}-2bp{e}^{-\alpha Ah}+b(p+q)e^{-\alpha Bh}+b(p-q)e^{\alpha Bh}\\ -2bp{e}^{\alpha Ah}+b(p+q)e^{\alpha Bh}+b(p-q)e^{-\alpha Bh}\\ a(q-p)e^{\alpha Ah}-2aq{e}^{-\alpha Bh}+at{e}^{-\alpha Ah}\\ a(q-p)e^{-\alpha Ah}+a(p+q)e^{\alpha Ah}-2aq{e}^{\alpha Bh}\end{array}\right),$$

$${\mathit{m}}_2=\frac{{\alpha }^3}{\Delta }\left(\begin{array}{c}B\left(cp\left(e^{-\alpha Bh}-e^{\alpha Bh}\right)+cq\left(-2e^{-\alpha Ah}+e^{-\alpha Bh}+e^{\alpha Bh}\right)\right)\\ B\left(cp\left(e^{-\alpha Bh}-e^{\alpha Bh}\right)+cq\left(2e^{\alpha Ah}-e^{\alpha Bh}-e^{-\alpha Bh}\right)\right)\\ A\left(dp\left(e^{-\alpha Ah}+e^{\alpha Ah}-2e^{-\alpha Bh}\right)+dq\left(e^{-\alpha Ah}-e^{\alpha Ah}\right)\right)\\ Ad\left(p\left(-e^{-\alpha Ah}-e^{\alpha Ah}+2e^{\alpha Bh}\right)+q\left(e^{-\alpha Ah}-e^{\alpha Ah}\right)\right)\end{array}\right),$$

$$\Delta =\frac{4c_{11}{\alpha }^4{r}^2{s}^2}{c_{33}}\left(2AB(cosh\left(\alpha Ah\right)cosh\left(\alpha Bh\right)-1)-\left(A^2+B^2\right)sinh\left(\alpha Ah\right)sinh\left(\alpha Bh\right)\right),$$

$$C_1=\alpha c_{13}{y}_1^{*}+c_{33}{\left(y_2^{*}\right)}^{\prime }{|}_{z=h},C_2={\left(y_1^{*}\right)}^{\prime }-\alpha {\left(y_2^{*}\right)}^{\prime }{|}_{z=h}.$$

In the above expressions, the following notation was used for abbreviation:

$$\begin{array}{cccccc}\hfill r& =c_{11}{c}_{33}-c_{13}^2,\hfill & \hfill s& =c_{13}+c_{44},\hfill & \hfill g& =c_{11}-A^2{c}_{44},\hfill \\ \hfill j& =c_{11}-B^2{c}_{44},\hfill & \hfill b& =c_{11}+B^2{c}_{13},\hfill & \hfill a& =c_{11}+A^2{c}_{13},\hfill \\ \hfill c& =c_{13}s-c_{33}j,\hfill & \hfill d& =c_{13}s-c_{33}g,\hfill & \hfill p& =Bac,q=Abd,\hfill \end{array}$$

and

$$\begin{array}{c}t=(A+B)(-c_{11}^2{c}_{33}+c_{11}\left(c_{44}{c}_{33}\left(A^2-AB+B^2\right)+c_{13}\left(c_{44}-AB{c}_{33}\right)+c_{13}^2\right)+\hfill \\ \hfill +AB{c}_{13}\left(AB{c}_{44}{c}_{33}+c_{13}^2+c_{13}{c}_{44}\right)).\end{array}$$

In the case when the density of volumetric force does not depend on z (and equal to $(0,1)$), the particular solution (23) takes the form:

$${\mathit{y}}^{*}=\frac{1}{{\alpha }^2{c}_{44}{c}_{33}}\left(\begin{array}{c}\frac{\left(c_{13}+c_{44}\right)}{A\left(A^2-B^2\right)}sinh\left(\alpha Az\right)-\frac{\left(c_{13}+c_{44}\right)}{B\left(A^2-B^2\right)}sinh\left(\alpha Bz\right)\\ \frac{\left(c_{11}-B^2{c}_{44}\right)}{B^2\left(A^2-B^2\right)}cosh\left(\alpha Bz\right)-\frac{\left(c_{11}-A^2{c}_{44}\right)}{A^2\left(A^2-B^2\right)}cosh\left(\alpha Az\right)-\frac{c_{11}}{A^2{B}^2}\end{array}\right).$$

Corresponding values for constants of integration are

$$C_1=\frac{Bdsinh\left(\alpha Ah\right)-Acsinh\left(\alpha Bh\right)}{\alpha AB{c}_{44}{c}_{33}\left(A^2-B^2\right)},C_2=\frac{1}{\alpha c_{44}{c}_{33}}\left(\frac{c_{11}}{A^2{B}^2}+\frac{acosh\left(\alpha Ah\right)}{A^2\left(A^2-B^2\right)}-\frac{bcosh\left(\alpha Bh\right)}{B^2\left(A^2-B^2\right)}\right).$$

Another important case, which will be used in section 6, corresponds to the unit distribution on interval $(a,b)\subset (0,h)$, i.e.

$$\mathcal{B}\left(z\right)=\left\{\begin{array}{cc}1,\hfill & z\in (a,b),\hfill \\ 0,\hfill & z\notin (a,b).\hfill \end{array}\right.$$

In this case, the integral (23) is calculated in the form

$${\mathit{y}}^{*}=\left\{\begin{array}{cc}0,\hfill & 0\le z<a,\hfill \\ \stackrel{1}{\mathit{y}}*,\hfill & a\le z\le b,\hfill \\ \stackrel{2}{\mathit{y}}*,\hfill & b<z\le h.\hfill \end{array}\right.$$

(24)

where

$$\stackrel{1}{\mathit{y}}*=\left(\begin{array}{c}\frac{(c_{13}+c_{44})}{2{\alpha }^{2}AB{c}_{44}{c}_{33}\left(A^2-B^2\right)}\left(A{e}^{\alpha B(a-z)}-A{e}^{\alpha B(z-a)}-B{e}^{\alpha A(a-z)}+B{e}^{\alpha A(z-a)}\right)\\ \frac{1}{2{\alpha }^2{A}^2{B}^2{c}_{44}{c}_{33}\left(A^2-B^2\right)}\left({\mathbb{A}}^{2}j\left(e^{\alpha B(z-a)}+e^{\alpha B(a-z)}\right)-2c_{11}\left(A^2-B^2\right)-B^{2}g\left(e^{\alpha A(z-a)}+e^{\alpha A(a-z)}\right)\right)\end{array}\right),$$

$$\stackrel{2}{\mathit{y}}*=\left(\begin{array}{c}\frac{(c_{13}+c_{44})\left(A\left(e^{\alpha B(a-z)}-e^{\alpha B(z-a)}-e^{\alpha B(b-z)}+e^{\alpha B(z-b)}\right)+B\left(e^{\alpha A(z-a)}-e^{\alpha A(a-z)}+e^{\alpha A(b-z)}-e^{\alpha A(z-b)}\right)\right)}{2{\alpha }^{2}AB{c}_{44}{c}_{33}\left(A^2-B^2\right)}\\ \frac{B^{2}g\left(-e^{\alpha A(z-a)}-e^{\alpha A(a-z)}+e^{\alpha A(z-b)}+e^{\alpha A(b-z)}\right)+A^{2}j\left(e^{\alpha B(a-z)}+e^{\alpha B(z-a)}-e^{\alpha B(b-z)}-e^{\alpha B(z-b)}\right)}{2{\alpha }^2{A}^2{B}^2{c}_{44}{c}_{33}\left({\mathbb{A}}^2-{\mathbb{B}}^2\right)}\end{array}\right).$$

It should be noted that in the case of an isotropic material ($t=0$), the values of the roots (22) of the characteristic equation coincide, $A=B$, and to construct a solution, instead of the fundamental matrix (21), one should use the matrix, which is obtained from the latter by passing to the limit and contains the ODE associated solutions:

$$\mathit{Y}=\left(\begin{array}{cccc}{e}^{\alpha z}& e^{-\alpha z}& -k(\alpha z+1)e^{\alpha z}& (\alpha kz-2)e^{\alpha z}\\ -e^{\alpha z}& e^{-\alpha z}& (\alpha kz-2)e^{\alpha z}& (\alpha kz+2)e^{-\alpha z}\end{array}\right),$$

where $k=1+\lambda /\mu$ ($\lambda$, $\mu$ are the Voigt average moduli (5)). All other relations are obtained similarly by passing to the limit.

5. Algorithmic solution implementation

From the above constructions it follows that the desired solution for an arbitrary right-hand side can be represented by the series

$$\begin{array}{cc}\hfill u& =\sum _{n=1}^{\infty }{J}_1\left({\alpha }_{n}r\right)\left({\mathcal{E}}_{1n}{e}^{{\alpha }_{n}Az}+{\mathcal{E}}_{2n}{e}^{-{\alpha }_{n}Az}+{\mathcal{E}}_{3n}{e}^{{\alpha }_{n}Bz}+{\mathcal{E}}_{4n}{e}^{-{\alpha }_{n}Bz}+{\mathcal{E}}_{5n}\right),\hfill \\ \hfill w& =\sum _{n=1}^{\infty }{J}_0\left({\alpha }_{n}r\right)\left({\mathcal{D}}_{1n}{e}^{{\alpha }_{n}Az}+{\mathcal{D}}_{2n}{e}^{-{\alpha }_{n}Az}+{\mathcal{D}}_{3n}{e}^{{\alpha }_{n}Bz}+{\mathcal{D}}_{4n}{e}^{-{\alpha }_{n}Bz}+{\mathcal{D}}_{5n}\right),\hfill \end{array}$$

(25)

where $E_{1n},\dots ,E_{4n}$ and $D_{1n},\dots ,D_{4n}$ are constants, while $E_{5n}$, $D_{5n}$ are functions on z (they become constant in the case of homogeneous volumetric force distribution with respect to z). Explicit expressions for these quantities can be obtained from the above formulas. Their form depends on the result of integration (23), (14) and can be very cumbersome (due to known difficulties of integration in a closed form). But to implement a computational algorithm, the use of closed analytical expressions at each stage seems to be superfluous. In order to use the simple solution obtained above for a constant density of volumetric forces on a rectangular support, we represent arbitrary RHS in (6) as piecewise constant vector-function

$$\mathbf{f}=\sum _{i,j}{\theta }_{ij}{\mathbf{F}}_{ij}$$

where ${\theta }_{ij}$ denotes the characteristic function

$${\theta }_{ij}=\left\{\begin{array}{cc}0,\hfill & (r,z)\notin {\mathcal{R}}_{ij},\hfill \\ 1,\hfill & (r,z)\in {\mathcal{R}}_{ij},\hfill \end{array}\right.$$

of rectangular domain

$${\mathcal{R}}_{ij}=(c_i,d_i)\times (a_i,b_i)\subset (0,1)\times (0,h),$$

and $F_{ij}$ are numerical two-component vectors, which values determine the levels of the RHS components in corresponding domain. Using the solutions (24), (15) one can find the values for each individual partial load, identified with pair numbers $(i,j)$, as follows:

$$\left({\mathcal{R}}_{ij},F_{ij}\right)\mapsto \left({\mathcal{E}}_{1n}^{ij},\dots ,{\mathcal{E}}_{5n}^{ij};{\mathcal{D}}_{1n}^{ij},\dots ,{\mathcal{D}}_{5n}^{ij}\right)$$

Since at this stage the problem is solved in a linear formulation, it is possible to use the superposition of solutions to determine the response on the full field of body forces. Thus the desired solution of the problem (6) can be presented in the same form (25) with following parameters

$$\mathbf{f}\mapsto \left(\sum _{ij}{\mathcal{E}}_{1n}^{ij},\dots ,\sum _{ij}{\mathcal{E}}_{5n}^{ij};\sum _{ij}{\mathcal{D}}_{1n}^{ij},\dots ,\sum _{ij}{\mathcal{D}}_{5n}^{ij}.\right)$$

To implement the computational algorithm, it is convenient to pre-calculate solutions for partial unit volumetric forces, thereby determining the discrete analogue of the Green’s function. This algorithm is a computational implementation of the linear operator inverse to the operator, generated by differential equations (6) and boundary conditions (7). This linear operator will be further denoted by the symbol ${\mathcal{L}}^{-1}$.

6. Geometric non-linearity

To take into account the geometric non-linearity, one has to treat stresses, introduced previously in the framework of linear elasticity, as second Piola – Kirchhoff stresses $\mathit{S}$. This allows formal expression for constitutive law, similar to (3)

$$\mathit{S}=\mathit{C}:\mathit{E}=\mathit{C}:\epsilon +2\mathit{C}:\left({\mathit{D}}^{\mathrm{T}}·\mathit{D}\right),$$

(26)

where $\mathit{E}$ is Green – Saint-Venant strain tensor

$$\mathit{E}=\frac{1}{2}\left(\mathit{D}+{\mathit{D}}^{\mathrm{T}}+{\mathit{D}}^{\mathrm{T}}·\mathit{D}\right),$$

$\epsilon$ is infinitesimal strain tensor, defined in (2), $\mathit{D}={\left(\nabla \otimes \mathit{u}\right)}^{\mathrm{T}}$ denotes displacement gradient and $\mathit{C}$ stands for the forth rank tensor of linear elastic moduli, used in (3). Understanding that this is only the one of the possible, and probably not the best, type of constitutive law in non-linear elasticity, we use it to take advantage of the previously constructed linear solution.

Let us now proceed to the construction of an iterative algorithm for determining displacements taking into account geometric non-linearity. The starting point will be the equilibrium equation formulated in reference coordinates with respect to second Piola-Kirchhoff tensor $\mathit{S}$:

$$div\left(\left(\mathit{I}+\mathit{D}\right)·\mathit{S}\right)+\mathit{K}=\mathbf{0}.$$

Bearing in mind the linearity of the divergence operator and taking into account the additive expansion (26), we rewrite the equation in the form

$$div\left(\mathit{C}:\epsilon \right)+div\left(\mathit{D}·\mathit{S}+2\mathit{C}:\left({\mathit{D}}^{\mathrm{T}}·\mathit{D}\right)\right)+\mathit{K}=\mathbf{0}.$$

Now we can implement a fixed-point iteration process:

$$div\left(\mathit{C}:{\epsilon }_k\right)+{\mathit{K}}_{k-1}^{*}=\mathbf{0}.$$

Here

$${\mathit{K}}_0^{*}=\mathit{K},\forall k>0{\mathit{K}}_{k-1}^{*}=div\left({\mathit{D}}_{k-1}·\mathit{C}:{\epsilon }_{k-1}+2(\mathit{I}+{\mathit{D}}_{k-1})·\mathit{C}:\left({\mathit{D}}_{k-1}^{\mathrm{T}}·{\mathit{D}}_{k-1}\right)\right)+\mathit{K}.$$

(27)

At each step of the iterative process, it is necessary to solve a linear problem (6), (7) for the distribution of volumetric forces specified by the relation (27). For this purpose, we use the analytical solution, obtained above, which determines the linear (inverse) operator ${\mathcal{L}}^{-1}$. Using it, we get the following recurrent sequence:

$${\mathit{u}}_k={\mathcal{L}}^{-1}\left[\mathit{K}\left[{\mathit{u}}_{k-1}\right]\right],\mathit{K}\left[{\mathit{u}}_{k-1}\right]={\mathit{K}}_{k-1}^{*}.$$

(28)

This sequence converges quickly enough, if at each step the right-hand side of the equation (28) sligtly differs from $\mathit{K}$. At the same time, it is the cases when the nonlinear term manifests itself significantly that are most interesting. To overcome the divergence problem, it is proposed to use regularization of the iterative process, the idea of which is similar to the Fejér summation method (with arithmetic means). In this case, at each step we calculate

$${\mathit{u}}_k={\mathcal{L}}^{-1}\left[\mathcal{R}\left[{\left(\mathit{K}\left[{\mathit{u}}_{n-1}\right]\right)}_{n=1}^k\right]\right],$$

where $\mathcal{R}\left[{\left(\mathit{K}\left[{\mathit{u}}_{n-1}\right]\right)}_{n=1}^k\right]$ is the regularizer of the iterative process, using the ”memory” of all previous iterations and, thereby, creating ”computational viscosity” in the algorithm:

$$\mathcal{R}\left[{\left(\mathit{K}\left[{\mathit{u}}_{n-1}\right]\right)}_{n=1}^k\right]=\sum _{n=1}^k{\zeta }_n\mathit{K}\left[{\mathit{u}}_{n-1}\right].$$

Sequence of weight factors ${\left({\zeta }_n\right)}_{n=1}^k$ specifies the ”memory” depth. If

$${\zeta }_n=\left\{\begin{array}{cc}1,\hfill & n=k,\hfill \\ 0,\hfill & n<k,\hfill \end{array}\right.$$

then the regularizer does not make any corrections to the iterations, but if $\forall n{\zeta }_n=1/k$, then the regularizer takes into account the entire history of the iterative process and performs the best smoothing of its variability. Of course, in this case more iterations are required. Choice of a specific type of sequence of weighting factors ${\left({\zeta }_n\right)}_{n=1}^k$ depends on the contribution of nonlinear terms to the solution and, ultimately, on the ratio of the expected deflection to the characteristic size of the plate.

In order to take into account the initial tension, one should replace displacements $\mathit{u}$ with $\mathit{u}+{\mathit{u}}_0$ in the above formulas. Here ${\mathit{u}}_0=\xi {\mathit{e}}_r$ and $\xi$ defines predetermined homogeneous initial stretching in radial direction.

7. Sensitivity of proposed theoretical model on the mechanical parameters

Before identification of mechanical characteristics upon the experimental data, we analyse the sensitivity of proposed theoretical model on the parameters, which characterise the elastic properties of the material, and on the loading intensity.

Above solutions were constructed for both isotropic and anisotropic cases. Moreover, anisotropic elastic moduli can be continuously transformed into isotropic counterparts, when the anisotropy factor t changes from 1 to 0 (see (4)). Now, to estimate the impact of anisotropy factor, we find the dependence on t of plate displacements under a unit distributed volumetric force. Recall that the case $t=1$ corresponds to a full account of anisotropy, at $t=0.5$ we have a partial account, and the case $t=0$ corresponds to a completely isotropic material which elastic moduli are obtained by Voigt averaging (5). For test calculations, we deliberately significantly increased the thickness of the plate to 20% of the radius. This allows, firstly, to illustrate the impact of the anisotropy factor more clearly, and secondly, to present the calculation results on a natural two-dimensional picture. Keeping in mind the linearity of the model, we perform calculations for a unit intensity of body forces, and present results in dimensionless form.

The following values of elastic moduli for crystalline silicon were used in the calculations:

Table 1. Elastic moduli for crystalline silicon, GPa

${\mathit{c}}_{11}$	${\mathit{c}}_{12}$	${\mathit{c}}_{13}$	${\mathit{c}}_{33}$	${\mathit{c}}_{44}$
186	64	64	80	186

Table 1. Elastic moduli for crystalline silicon, GPa

${\mathit{c}}_{11}$	${\mathit{c}}_{12}$	${\mathit{c}}_{13}$	${\mathit{c}}_{33}$	${\mathit{c}}_{44}$
186	64	64	80	186

They correspond to the following Voigt-averaged moduli for polycrystallite

$$\mu =68.6\mathrm{GPa},\lambda =58.9\mathrm{GPa}.$$

In Figure 6 the graphs that illustrate such dependencies are given. In Figure 6 one can find the difference between distorted mesh with and without accounting of anisotropy, while Figure 6 demonstrates the distinction in displacements of the bottom of the plate. It can be seen that the account of anisotropy does not have a significant effect on the theoretical flexural rigidity of the plate: the largest difference in displacements does not exceed 8.7%.

To assess the impact of the anisotropy factor on the stress state, the stress values were found at the lower central point (point A in Figure 5) and at the corner point on the boundary (point B in Figure 5). The results are shown in the Table 2. The tensions found with account of anisotropy and without it differ by no more than 8%.

Much more significant is the influence of geometric nonlinearity caused by an increase in tensile stresses on the cylindrical part of the boundary when the plate is deformed. To illustrate this, consider the dependence of the displacements of the center of the plate bottom with increasing intensity of a uniformly distributed pressure on the plate top. We will change the load intensity in the range $[0,100]$ kPa. The dependence over the entire range is shown in Figure 7. It clearly has non-linear nature. At the same time, in certain intervals it can be interpreted with sufficient accuracy as linear. The first, and most natural, such interval is rather short, it is $[0,0.1]$ kPa. It characterizes deformation range, in which all non-linear effects are negligibly small. Corresponding pressure-displacement curve is shown on Figure 7 in purple shading. Already at the next interval, $[0.1,1]$ kPa, nonlinear effects manifest themselves significantly. This can be seen on the pressure-displacement curve on Figure 7 (in green shading), where the linear dependence that continues the dependence from previous interval is shown in black. For comparison the tangent to the pressure-displacement curve (at the middle of the interval) is shown in red. Keeping in mind the concept of tangent modulus used in elasticity theory, we can interpret this tangent as some effective stiffness. It must only be borne in mind that such effective stiffness, due to the fact that it depends on the tension on boundary, reflects a property of the structure and not of the material alone. The dependence in the next interval $[1,10]$ kPa is most nonlinear and is shown in blue shading on Figure 7 Just as before, the black line represents the continuation of the linear relationship, and the tangent in the middle of the interval is shown in red. It should be noted that, unlike the previous intervals, the displacement computed within the linear theory differs from the displacements calculated with account of geometric non-linearity by several times. The slope of the black and red lines are also significantly different, which means that the effective stiffness varies considerably from what can be calculated in the linear approximation. The last interval is the longest, $[10,100]$ kPa (Figure 7). Displacements calculated within the framework of linear theory differ by an order of magnitude from those calculated taking into account geometric non-linearity. At the same time, the very dependence of displacements on pressure in this interval is close to linear, as can be seen from its comparison with the tangent (red straight line). One just needs to keep in mind that such a linear dependence does not characterize the linearity of the elastic properties of the material, but the linearity of the rigidity of the structure in this load range.

The increase in the effective stiffness of the plate due to nonlinear effects can be quantitatively assessed by comparing the slope of tangents at midpoints of different intervals of the pressure-displacement curve. Their values are given in Table 3. Note that the minimal and maximal slope differ by a factor of 32.

Now consider the effect of the initial tension in the radial direction on the effective stiffness of the plate. For different values of the initial tension in the radial direction, we will apply the same load 100 kPa to the upper surface of the plate. Displacements of the bottom surface are shown on Figure 8. Note that by changing the initial tension one can change the maximal deflections by an order of magnitude.

Therefore, the preliminary analysis of the mathematical model, applied to the plate tested in the experiment, allows us to state the following:

• Displacements weakly depend on whether or not the anisotropy factor is taken into account.
• Geometric nonlinearity can be neglected at a load of no more than 100 Pa. At pressures exceeding this value, the influence of geometric nonlinearity increases dramatically.
• Displacements can vary by more than an order of magnitude depending on the initial in-plane tension.

8. Identification of experimental data

Now let’s move on to the analysis of experimental data. Two series of experiments were carried out, during which the profiles of the plate deformed by uniform pressure were determined. The intensity of the pressure varies from 20 kPa up to 300 kPa with increment 10 kPa. The difference between the profile measured at the n-th step and the profile measured before loading determines the displacement along the z axis of the top surface of the plate. The experimentally determined displacements in the first series of measurements are shown in Figure 9, and in the second in Figure 10, where the distribution of measured values are shown as polygons of alternating colours.

The displacements of the central point are determined by averaging over a centered strip of $0.1$ mm wide (one tenth of the diameter). It is shown in yellow in the pictures. The results of averaging are shown with brown lines and pointed on lateral plane with brown dots. These dots together represent the pressure-displacement curve for the central point of plate bottom (shown in black). Confidence intervals (with 0.95 level) are also shown in brown.

As expected, the dependencies for the displacements of the central points on pressure are close to linear. This, however, does not mean at all that a linear model of an elastic body can be used to interpret the experiment. Based on the preliminary calculations carried out above, we can expect that the pressure-displacement curve corresponds to the last interval, in which the geometrical non-linearity has strong impact.

At this stage, we have insufficient experimental data to simultaneously identify elastic moduli and tension, but if we assume the elastic moduli are known, then tension can be identified. To this end we compute the values of central point displacements for the values of tension that corresponds to radial displacements of the boundary in the range $[0,1]\mu$ and pressure varies in the range $[0,300]$ kPa. The resulting dependence is shown on Figure 11 as yellow surface. Then we take the experimental dependence for central point displacements from two series and compute the standard deviation of the experimental data from theoretical estimation for varying value of tension. All that remains is to find such tension values, and for which the standard deviation becomes minimal. These values are

• For the first series of experimental data the stretching in radial direction is $0.698\mu$. The corresponding radial component of stresses is $328.1$ MPa.
• For the second series of experimental data the stretching in radial direction is $0.748\mu$. The corresponding radial component of stresses is $351.9$ MPa.

These values correspond to certain sections of the surface that represents the theoretical distribution for displacements, at the intersection with which the experimental distributions in the form of blue polygons are shown on Figure 11

For a more detailed comparison of theoretical and experimental values, we depict these sections on separate graphs in Figure 12 One can also see on these graphs confidence intervals computed upon experimental data. The slight difference between the theoretical distributions and experimental data can be explained by the fact that the theoretical model takes into account geometric non-linearity, while the constitutive law is still assumed to be linear. It is likely that the use of more suitable hyperelastic potentials will allow us to construct a more accurate theoretical model. We hope that such generalizations will be obtained in subsequent works.

9. Comparison of the proposed mathematical model with standard formulas

Finally, we compare the proposed mathematical model with standard formulas usually used to identify experiments on the deformation of thin films. First we note that three classes of such formulas can be distinguished. The first class includes formulas obtained from solutions of boundary value problems on the bending of elastic plates within the framework of a theory of the Kirchhoff – Love type. Formulas of the second class are obtained from numerical solutions of the Föppl – von Kármán equations, obtained by the Galerkin method. In this case, as a rule, a one-term approximation is used, which corresponds to the assumption that the curved shape of the plate is spherical. The third class contains formulas obtained empirically. For a uniformly loaded circular plate of isotropic material, these three types are represented by the formulae, given in Table 4.

Note that the last, empirical, formula is a formal combination of two different analytical relationships, the first of which determines the bending of the membrane from the solution of the Poisson equation, and the second represents the first term of the Galerkin solution to the bending equation of flexible plates.

To compare the proposed solution with the mentioned approximations, a series of calculations were made for a sequence of thicknesses and pressures. Their values and the results are shown in the Table 5. In this table one can see the vertical displacements of the plate bottom, computed with proposed model, v.s. deflections calculated in the framework of corresponding plate theory with formulae, given in first two rows of Table 4. Differences are shown by shading. From the analysis of the results it is clear that the theory of Kirchhoff-Love plates gives a good approximation for small deflections and thicknesses, while the spherical approximation of the solution to the Föppl – von Kármán equations gives a small error for significant deflections and also small thicknesses. In general, these approximations have a rather narrow range of application. At the same time, they are often used to identify the results of experiments, and therefore we decided to calculate from them the values of theoretically expected displacements of the midpoint of the plate in order to assess the possible error introduced by such approximate calculations. The results are shown in Table 6. The maximal displacements of the plate were compared at a pressure of 150 kPa. The first two columns show the displacements and initial tensions found above from the experimental results. The third column shows the theoretical values of maximal displacements calculated using the proposed model. The last two columns show displacements found using simplified formulae. Errors calculated relative to experimental data are given in parentheses. It can be seen that the errors of the approximate formulas exceed the accuracy required to identify the anisotropy factors. At the same time, we note that the empirical formula (into which, however, we substituted the initial tension found using the proposed method) for large pressure gives rather accurate result. This is apparently due to the fact that at large deflections the membrane stress-strain state prevails and is well determined from the classical Poisson equation (recall that the first part of the empirical formula is precisely the solution to the corresponding Poisson equation).

10. Conclusions

The following conclusions can be drawn from the work.

1. The account of the anisotropic moduli of the crystal in the analysis of single-crystal plate stress-strained state leads to theoretical corrections about 10%. In this regard, to analyze anisotropic properties, the experimental error should be less than 10%. At this level of errors, in the identification of the experimental results the approximation formulas either cannot be used at all, or they can be used in narrow ranges of mechanical parameters. On the contrary, the proposed numerical-analytical solution can be used for the identification in a wide range of parameters.

2. The initial tension has a great influence on the effective flexural rigidity of the plate. For single-crystal circular plates formed on the basis of a silicon-on-insulator structure the initial stresses change the rigidity in several times.

3. The proposed experimental technique makes it possible to quite accurately quantify the initial stresses and, taking them into account, the effective flexural rigidity.

The mathematical modelling and analysis of experimental data allows us to clarify the questions that experimenters need to answer in subsequent works. First, a series of mechanical tests should be carried out in the range of almost linear deformations, from which the anisotropic elastic moduli can be directly determined, independent of the effects caused by geometric non-linearity. To this end, if one use plates of the same thickness as in the experiments described above, then the loading steps should be taken on the order of 100 Pa. It may be technically more convenient to use thicker plates, which will allow the pressure increment to be increased. Secondly, a series of experiments should be carried out with a small step of pressure increment, but adding these increments upon a significant pre-load. This will allow one to feel the effects of superimposing small deformations on finite ones and distinguish the impact of geometric non-linearity and possibly physical non-linearity. The last factor is of particular interest for the finite deformations modelling of anisotropic crystals. Thirdly, a series of experiments, similar to those described in present work, should be carried out, but with plates obtained at different deposition temperatures. Changing the manufacturing conditions will affect the initial tensions and, as the analysis in this paper shows, will significantly change the effective flexural stiffness of the plates.