Closed Solution of the Föppl – Von Kármán Equations for Square Plates

1. Introduction

The Föppl – von Kármán equations, which describe stress-strain state of flexible plates, has attracted attention of physicists and mathematicians since August Föppl introduced their prototype in the treatise "Vorlesungen uber technische Mechanik" [1]. A. Föppl assumed that in the case of finite deflections, both in-plane linear strains and the rotations are small with respect to unity, while the rotations and squares of in-plane strains are infinitesimals of the same order. Furthermore, he supposed that the bending strains can be thought of infinitesimally small with respect to membrane ones. With these assumptions he derived the following system of equations:

$$,.$$

(1)

Here is a plate thickness, denotes Young’s modulus of its material and is the function of transversal load applied to the plate surface. The function $\left(,\right)$ represents the deflections of the plate middle-plane, while $\left(,\right)$ is the Airy stress function related with the components of the stress tensor within it:

$$1,2,.$$

(2)

Symbols and stand for the biharmonic operator and the Monge – Ampère form [2], respectively. In Cartesian coordinates they can be represent as follows:

$$\left[f\right],\left[f\right]g.$$

(3)

Theodor von Kármán improved first equation in (1) by adding an item that accounts for influence of in-plane stresses on bending and derived the system of equations in a well-known almost symmetric form [3]:

$$,.$$

(4)

Here denotes the bending stiffness of the plate, which expressed through $,$ and the Poisson’s ratio : .

The Föppl – von Kármán equations were initially developed as a mathematical tool for buckling analyses of ship hall parts [4]. However, wide application of the structural elements that were able to bend sufficiently without undergoing plastic deformations calls for further development and improvement of the non-linear theories of plates, subjected to various loads, including transverse one. The Föppl – von Kármán equations proved to be a good approach, taking into account non-linear effects but not having the intimidating bulky form of the non-linear elasticity equations. At present in many modern technologies, particularly MEMS, the problem of mathematical modelling of stress-strain state in flexible parts is extremely relevant [5,6].

First solutions of non-linear bending problem were obtained for circular plates subjected to axial symmetric loading, which allowed one to significantly simplify the Föppl – von Kármán equations. Thus, the first attempt to solve them for a clamped circular plate under uniform pressure was undertaken by Nádai [7] in 1925. Another variants of approximate solution were obtained by Timoshenko who used the Rayleigh – Ritz method [8] and by Way through the power series expansions [9]. In the first part of their study [10], Friedrichs and Stoker presented numerical solutions for a simply supported circular plate under compressive forces at the boundary, obtained by three distinct algorithms. Furthermore, in [11] they introduced a change of variables, which reduced the order of the governing system of equations from 8th to 4th. In addition, in the second part of [10], they presented a proof of the existence theorem and demonstrated that, apart from a trivial solution, there is at most another one (accurate to a sign). In present paper, only square plates will be discussed, so the part of review focused on circular plates is exhausted by a few references. A detailed survey can be found in [12].

One of the first approximate solutions of the Föppl – von Kármán equations for a simply supported square plate under uniformly pressure was derived by Kaiser [13] using finite-differences method. An alternative approach to approximate solution was proposed by Way [14]. He obtained solution for clamped rectangular plates with various aspect ratios by Ritz method. The Galerkin method was also successfully applied for this problem by Panov [15]. Levy first derived exact solutions in terms of double Fourier series for a simply supported rectangular plate [16] and a clamped square one [17]. Later, he and Greenman solved the problem for a clamped rectangular plate having an aspect ratio of 3/2 [18]. This method also has been applied in [19] and [20] for long clamped and simply supported plates under combined load. It is important to note that in works [16,17,18,19] it is assumed that clamped edges are movable they are rigidly clamped against normal displacements and rotations while allowing for movement in plane. A numerical solution for a clamped rectangular plate with immovable edges was obtained by Wang [21]. He also considered the important case of a riveted panel [22]. Lastly, it is worthy noting the article by Green and Southwell on the application of the relaxation method to solve the Föppl – von Kármán system [23]. A detailed description of the majority of the aforementioned results, as well as many others, can be found in the remarkable books [4,24].

Briefly discuss the recent work on this topic. In the works [25,26] the solution for a simply supported plate was obtained by the "accelerated" iterative method. The perturbation method was successfully applied in the works [27,28,29,30]. In addition, the work [29] is notable because of application of the spline function method for solving a set of linear equations. A number of authors successfully applied the finite-difference method [31,32,33]. A modification of this approach was proposed by Petit [34], who employed the finite-difference method in combination with the relaxation one. Furthermore this method was also used in the work [35] for a plate under combined load. The finite elements method was widely used [36,37,38,39]. Also, a few attempts to apply the boundary elements method were made [40,41,42,43,44]. The Ritz method was widely used as well. Boresi and Terner [45] applied it to obtain the solution for a simply supported rectangular plate with stress-free edges. In the work [46] the rectangular plates with various boundary conditions were considered using the automated Ritz method procedure. This procedure was also applied to calculate large deflections of orthotropic plates [47] and functionally graded material ones [48]. Enem [49] introduced the Ritz method solutions for rectangular plates with twelve different boundary conditions. Recently, the homotopy analysis method has begun to be successfully applied to solve the Föppl – von Kármán system [50,51,52,53]. A sufficiently detailed review of works managed before the 1970s are presented in the article [54]. Also, a comprehensive list of works managed before the 2000s is given in the book [55], while detailed review containing more recent works can be found, for example, in the PhD thesis by Khoa [56].

An original approach to solve the Föppl – von Kármán system was proposed by H. Berger in his PhD thesis [57]. This approach was based on the hypothesis that the strain tensor second invariant makes a small contribution to the full potential energy in comparison with the contribution made by the first invariant. Although this hypothesis lacked a physical basis, it allowed one to decouple the Föppl – von Kármán system. Berger obtained solutions for circular and rectangular plates with various boundary conditions using this method. The further development of Berger’s technique was continued by a number of other scientists [58,59,60,61]. A similar method for decoupling the Föppl – von Kármán equations was employed in the recent works [62,63,64]. An alternative approach to decoupling of this system system was proposed by Banerjee and Datta [65].

A large amount of results were obtained by the Galerkin method. Yamaki [66] used a modification of the one-term Galerkin method to solve the dynamic analogue of the Föppl – von Kármán equations. He also gave a classification of flexible plate boundary conditions. In their work [67], Iyengar and Naqvi applied biorthogonal expansions of the displacements and the Airy function to obtain a solutions for a simply supported plate and a clamped one. For each boundary condition type two options were considered: all edges are stress-free and all edges are immovable. In the papers [68,69] the authors obtained the Galerkin method based system of cubic equations for Kármán plates in an explicit form. In the work [70] this result was supplemented with the explicit expression of the inverted Galerkin matrix. This procedure allowed to obtain the solution more effectively due to avoiding a laborious calculation of the inverted Galerkin matrix. Recently, the wavelet-Galerkin method has been widely used [71,72]. This method leads to simplification of calculations due to the choice of wavelet functions as the basis.

It is advisable to consider some works devoted to analyse of the Föppl – von Kármán equations in term of pure maths. Initially, they did not have a rigorous mathematical justification and were derived through an engineering approach significantly based on a few auxiliary hypotheses [2]. The justification of this system has been provided by Ciarlet, who has demonstrated [2,73] that this system can be obtained using the asymptotic expansion of the equations of non-linear elasticity theory by a small parameter, assuming certain hypotheses about the asymptotic orders of displacement vector components. This approach was also applied in the recent work [74] to construct the family of Kármán-like plate models.

The existence theorem for the Föppl – von Kármán boundary-value problem with sufficiently smooth boundary and homogeneous boundary conditions was proofed by Morozov [75] using the Leray –- Schauder method [76]. The conditions for the existence of a non-trivial equilibrium state of a Kármán plate under compression forces were formulated by M. Berger and Fife [77]. Knightly [78] demonstrated the existence of a solution for a plate subjected to combined loading, and moreover, he proofed that if the magnitude of the loading is sufficiently small, then this problem has a unique solution. Further, Knightly and Sather [79] have modified this uniqueness result. Hlaváček and Naumann have obtained similar results for the various types of inhomogeneous boundary conditions [80,81]. It is important to highlight that results of Knightly, Hlaváček and Naumann are applicable to plates with the boundary which consists of curves from $C^3$ and may have corners at junctions of the arcs they are valid for the majority of plates used in engineering.

From this brief analysis, it is observed that, to date there is no long-standing methods to solve Föppl – von Kármán equations for square and rectangular plates fixed at the boundary with respect to all displacement components. Moreover, direct 3D FEM modelling of thin plates results in a ill-conditioned algebraic system due to the small parameter, which is the plate thickness. However, it is the full fixation that is usually implemented in experiment and practice. In this context, the development of Föppl – von Kármán equations solution methods for square and rectangular plates, fully fixed at the boundary, is an urgent problem in computational mechanics. This is the subject of present paper. In it a new solution is proposed in the form of decomposition over basis in functional space of square-integrable functions. The coefficients of expansion are determined from infinite-dimensional cubic system via reduction method. This approach could provide a non-linear counterpart to Galerkin method. The principle difference from similar studies [67,70] consist in a strict formulation of the auxiliary boundary problem, allowing for account natural ways of plate fixing.

2. Problem Statement

Consider a clamped square plate of side 2 and thickness , made from a linear-elastic material with moduli , and subjected to an arbitrary transversal load $\left(,\right)$ (see Figure 1). For a closed formulation of the initial boundary problem, we should complement the equations (4) with suitable boundary conditions, which we will now deal with. The boundary value problem will be formulated in the region: ,

with the boundary: ,

where the parts of the boundary $,$ correspond to pairs of opposite edges,

$$,.$$

The boundary conditions for the flexible plates usually are gathered into two sequences associated with plate bending and with plane stress-strain state respectively. Consider the first one. In the case of clamped plates it can be written as: .

Other variants of the boundary conditions related to bending are well known [24], and we will not comment on them as (5) are the ones that will be stated hereinafter. The boundary conditions associated with in-plain stress-strain state are not so widely covered in literature and we will discuss their formulation in more detail. Following [17,21] we consider two cases:

– movable edges: ;

– immovable edges: .

Note that despite the designation, moving edges are not free from fixing at all. The points on them can only be moved along the edges, while their displacement in a perpendicular direction is constrained. In equations (6a), (6b) $\left(,\right)$ and $\left(,\right)$ denote the in-plane plate displacements (the in-plane displacements of the points on the plate middle plane). These functions are not explicitly included in the governing system (4), so in order to resolve the boundary-value problem, it is necessary to express them through the desired functions $,$. The sought relations can be obtained using the Cesàro formula, which express the displacement vector and the rotation vector via the infinitesimal strain tensor [82]: .

Note that the displacement vector is defined with account of kinematic hypothesis of the Kirchhoff-like plate theory: .

Therefore, to relate the displacements to the Airy stress function , it is necessary to express the infinitesimal strain tensor through the stress tensor . The essence of Kármán theory based on a number of simplifications which allow only the use Saint Venant – Kirchhoff constitutive equation [73,83]: ,

where is the reduced Green – Saint Venant strain tensor (also known as the Kármán strains): .

Finally, substituting the relation (9) into the Cesàro formula (7) one can obtain:

$$\begin{array}{c}\hfill .\end{array}$$

This formula can be significantly simplified by choosing an appropriate path of integration. Taking into account that displacements and rotations of whole plate are defined up to an additive constants, any point of it can be fixed. Apart from that, if the transversal load function is odd or even by both coordinates then the centre of a plate will be displaced vertically without changing its in-plane position. For the sake of simplicity, assume that the load function and, as a consequence, sought functions $,$ are even by both coordinates (three other cases of symmetry are derived from reasoning completely similar to the following): .

Here we also restrict the class of considered functions $,,$ only by square-integrable ones. This allows to decompose them over basis and reduce the boundary-value problem to a system of cubic equations.

With account of the transversal loading function symmetry (10) the displacement on can be calculated as follows:

$$1,1.$$

Substitution of equations (8) and the relations (2) into above formulae results in:

$$\begin{array}{c}1=0,\\ 2=0.\end{array}$$

(5)

However, taking into account that a function is identically zero iff it takes value zero at least at one point and its the first derivative is zero also, the latter integral condition can be reduced to a differential relation:

$$2=0.$$

The last terms in the right-hand sides of this relation can be omitted due to the fact that the deflection function is identically zero on the edges (5):

$$2=0.$$

Thus, above relation together with the first integral equality (5) fully describes the boundary conditions for the in-plane displacements of the plate on . The formulae for the boundary conditions on can be written in the similar way up to replacement of the variable by and function by respectively.

Now we can move on to formulation of boundary value problem. To this end introduce the following dimensionless variables: .

In the new variables the Föppl – von Kármán equations (4) and corresponding boundary conditions can be written as:

$$\begin{array}{c},,\end{array}$$

(6)

$$\begin{array}{c},\end{array}$$

(7)

$$\begin{array}{c}1,2;\end{array}$$

(8)

– movable edges: ;

– immovable edges:

$$1,2,$$

(16b)

where , ${\tilde{\nabla }}^2:=a^2{\nabla }^2,\tilde{L}(·,·):=a^{4}L(·,·)$. In the what follows, we will omit the tildes above these operators for brevity. The equations (6) – (13) or (9) constitute two different boundary-value problems for clamped square plate.

To obtain an approximate solution of the problem (6) – (13) or (9) a recursive approach is often used. Briefly describe its essence for comparison with the proposed solution below. The basis of the recursion can be taken as a solution of the Sophie Germain – Lagrange equation:

$$0:={\mathcal{M}}^{-1},$$

while the step of the recursion is defined as follows (here subscript indicates recursion step number):

$$={\mathcal{M}}^{-1}\left(-\left(-1,{\mathcal{M}}^{-1}\left(-1,-1\right)\right)\right).$$

Stress function can be obtained on each step from the relation

$$:=-{\mathcal{M}}^{-1}\left(,\right).$$

Here $\mathcal{M},\mathcal{M}$ stand for differential operators, defined by the same biharmonic differential expression (3).1, but by different boundary conditions, corresponding to deflection () and stress function (), (13), (9). Computational efficiency of recursive process depends on way how to represent inverse linear operators $\mathcal{M}$ and $\mathcal{M}$. It can be done differently, but for square plates the usage of complex potentials method (Goursat – Muskhelishvili – Kalandiya) [84] seems to be most appropriate. Within this method one can define the solution of linear problem $\mathcal{M}=$ as the sum of some particular solution $*$ that satisfy non-homogeneous biharmonic differential equation, but not met boundary conditions, and boundary solution of homogeneous equation $0$, such that $=*+0$ satisfies the conditions on the boundary.

Particular solution $*$ can be obtained as

$$*\left(z\right)={\displaystyle 116}{\int }_0^z{\int }_0^{\overline{z}}(z-\beta )(\overline{z}-\alpha )\left({\displaystyle \alpha +\beta 2},{\displaystyle \beta -\alpha 2i}\right)d\alpha d\beta ,$$

or

$$\begin{array}{cc}\hfill *\left(x,y\right)& ={\displaystyle 116}[{\int }_0^x{\int }_0^x(x-iy-\alpha )(x+iy-\beta )\left({\displaystyle \alpha +\beta 2},{\displaystyle \beta -\alpha 2i}\right)d\alpha d\beta \hfill \\ \hfill & +{\int }_0^x{\int }_0^y(y-\alpha )(x+iy-\beta )\left({\displaystyle x-i\alpha +\beta 2},{\displaystyle -x+i\alpha +\beta 2i}\right)d\alpha d\beta \hfill \\ \hfill & -{\int }_0^x{\int }_0^y(y-\alpha )(x-iy-\beta )\left({\displaystyle x+i\alpha +\beta 2},{\displaystyle x+i\alpha -\beta 2i}\right)d\alpha d\beta \hfill \\ \hfill & +{\int }_0^y{\int }_0^y(y-\alpha )(y-\beta )\left({\displaystyle 2x-i\alpha +i\beta 2},{\displaystyle \alpha +\beta 2}\right)d\alpha d\beta ].\hfill \end{array}$$

while the boundary solution $0$ can be expressed in terms of complex-valued series with respect to $z=x+iy$ as follows:

$$0=2\mathrm{Re}(\gamma \overline{\psi }+\chi ).$$

Here:

$$\chi =\int \psi {\displaystyle \partial \gamma \partial z}dz,\varphi =\sum ak{z}^k,\psi =\sum bk{z}^k,$$

where $ak,bk$ can be obtained from the solution of algebraic counterpart for boundary problem (see [84,85] for details). In above formulae $\overline{z}$ indicates the complex conjugate z while $\gamma$ stands for conformal map between unit disk and rectangular domain. This mapping can be exactly expressed via hypergeometric function, but for computations it is more continent to use polynomial approximation of order n:

$$\gamma =e^{-{\displaystyle i\pi 4}}{\displaystyle \sqrt{2}\Gamma \left(34\right)\sqrt{\pi }\Gamma \left(54\right)}\sum _{k=0}^n{\displaystyle z^{4k+1}(4k+1)k!}\prod _{s=0}^{k-1}\left(s+{\displaystyle 12}\right).$$

The inverse operator $M_{\varphi }^{-1}$ can be obtained in a similar way. Note, that the advantage of complex potential method is that the integral boundary conditions can be expressed directly via potentials with Love relations:

$$={\displaystyle 12\mu }\mathrm{Re}\left({\displaystyle 3-\nu 1+\nu }\phi -\gamma {\displaystyle \partial \overline{\varphi }\partial \gamma }-\overline{\psi }\right),={\displaystyle 12\mu }\mathrm{Im}\left({\displaystyle 3-\nu 1+\nu }\phi -\gamma {\displaystyle \partial \overline{\varphi }\partial \gamma }-\overline{\psi }\right).$$

Before embarking on general procedure of problem (6) – (13), (9) solution let us try to approach it with recursive method mentioned above. In doing so, keep in mind that in order to identify the characteristic features of non-linear deformation of the plate, it is sufficient to perform at least one step of recursion. The calculation results for the model problem $\left(=1,=1,=1,=1/3,=1/150\right)$ are shown in Figure 2. In it the subfigure (a) shows the distribution of deflections $0$ obtained without of account non-linear effects (it is the base or recursion, solution of the Sophie Germain – Lagrange problem). Subfigure (b) shows the distribution of plain deformation incompatibility, induced by bending, namely $\left(0,0\right)$. In subfigure (c) one can see stress function , corresponding to such incompatibility. In the second line of figures the distributions for normal and tangential in-plane stresses are shown, while the last subfigure (f) demonstrates dummy transverse loading $\tilde{q}:=\left(0,\mathcal{M}-1\left(0,0\right)\right)$ that has to be subtracted from actual loading on the first recursive step.

Already from the analysis of the results obtained in the first step of recursion, it is possible to give a qualitative estimate of the effect of the dummy load and, accordingly, of the non-linear effects depending on the expected value of the maximum deflections. In particular, with evenly distributed load the maximal value of dummy load (according to the distribution, shown in Figure 9.(f) for some load $\alpha \tilde{p}$ with dimensionless intensity $\alpha$ can be estimated as $4\times {10}^{-4}{\alpha }^3$. Therefore, the given and dummy loads will be of the same order when $\alpha \approx 20$, that, with account of distribution for deflections (Figure 9.a), is approximately equivalent to maximum deflections, obtained in linear approach, of 20% from the edge length of the plate. This value is a very rough estimate from above, because in fact the displacements will be much lower if we take into account non-linear effects, calculating the following recursion steps. Also, it can be noticed that:

– the influence of the dummy load increases proportionally to the cube of maximum deflection;

– the geometric shape of the bending surface at the final deformation significantly differs from the corresponding shape for small deflections, because the dummy load changes the distribution of RHS of bending equation on successive recurrent steps. We will see this again later when we consider the general procedure of solution.

Mentioned above strategy for solution of non-linear bending problem has computational stability for moderate bending only. When the amount of deflection increases, the dummy load is no longer a small adjustment and may even exceed given transverse load. In this case the sequence of recursively constructed approximations tends to diverge. Although special regularization procedures [5] can be used, diverging cannot be avoided. For this reason, below a different strategy is proposed, implementing the full algebraization of the problem and transforming it to the infinite-dimensional system of cubic equations. Reducing this system to a finite-dimensional one and using any variation of the successive approximation method to find the solution of the reduced system, it is possible to avoid computational instability even for deflections that are comparable with the size of the plate in the plane. The rest of the article is devoted to the construction and analysis of such algebraization.

3. Solution of the Linear Problem

Since the solution of the non-linear problem will be obtained in the form of decomposition over some basis in , generated by a linear self-adjoined operator, we will start with the linear problem. It is advisable to consider the linear bending of a square plate as such prototype for this problem, because, first, it generates a basis in terms of well-known Krylov-Duncan functions, and second, the problem itself can be used for testing a non-linear solution at zero approach.

The boundary-value problem for clamped Sophie Germain’s plate can be stated as follows:

$$\begin{array}{c},\\ .\end{array}$$

(10)

In order to expand the deflection function into a Fourier series, it is necessary to choose a basis in . According to the Hilbert – Schmidt theorem [86], a system of eigenfunctions of any compact self-adjoint operator constitutes a basis of . However, the biharmonic operator has sufficiently complicated system of eigenfunctions even with homogeneous boundary conditions, so it is more preferable to use an another self-adjoint operator. It is appropriate to use the operator constructed by removing a mixed derivative from the biharmonic one: .

The sequence of eigenfunctions for this operator can be obtained as a combination of solutions of an auxiliary Sturm – Liouville problem.

3.1. Auxiliary Sturm – Liouville Problem

Let us consider a Sturm – Liouville problem:

$$,.$$

A fundamental system of solutions for this equation is well-known: .

By substituting a linear combination of these functions into the boundary conditions, we get the spectral equation with respect to eigenvalues : .

Thus, a spectrum of the operator consists from two sequences of eigenvalues 1:

$$=1\infty ={\left\{1\right\}}_{=1}^{\infty }\cup {\left\{2\right\}}_{=1}^{\infty },$$

which are defined as roots of transcendental equations:

$$1,2.$$

These roots can be found numerically with arbitrary precision, but for sufficiently large values of very precise approximations can be obtained from asymptotic formulae:

$$1,2.$$

(11)

The question arises: what does it mean to be large enough? The following computation (Table 1 and Table 2) shows that for standard machine precision it enough to refine only five first values, while all the rest can be calculated directly by the formulae (11). Of course, for calculations with enhanced accuracy the number of refined roots should be increased.

Each sequence of the eigenvalues is associated with a sequence of the auxiliary eigenfunctions:

$$\begin{array}{cc}\hfill & 1,\hfill \\ \hfill & 2.\hfill \end{array}$$

(12)

Such sequences of function are convenient for derivations because their compact view, however for calculations it is more convenient to use normalized functions:

$$1,2.$$

Curves of the first several normalized eigenfunctions are shown in Figure 3.

It is evident that each pairwise combination of functions (12) presents a eigenfunction of two-dimensional operator . Thus, we arrive at four families of eigenfunctions for : .

Although each term of these families is marked by two indexes, it is preferable to use a one-indexed numeration system. Such system can be provided by any bijection $+{\times }^{+}{\to }^{+}$. The optimal numeration system is well-known from classical courses of calculus: .

This bijection is named the Cantor pairing function for the positive integers and was proposed by him in the article [87]. The correspondence between two systems of numeration can be illustrated by a diagram as shown in Figure 3 (c). This way for numeration is the most preferable for square regions because it allows one to construct the reduced sums symmetric with respect to both indices. In order to achieve this, it is sufficient to include complete diagonals in the sum.

Thus the two-indexed system of basis functions (29) can be replaced by the following one: .

One can go back to the functions (29) with reverse function to (30). In order to construct this function it is necessary to relate a number of the current term with the number of diagonal on which this term is located: ,

where the symbols $\lfloor \dots \rfloor$ denotes the floor. From here, it is simple to obtain the inverse function: .

Hereinafter and are brief notations for components of the inverse function.

3.2. Solution of the Linear Boundary-Value Problem

Now it is time to recover the boundary-value problem (), (10). We will seek the solution of this problem in the form of expansion in terms of a series of normalized basis functions (35): .

This solution is guaranteed to satisfy the boundary conditions () due to the appropriate choice of the basis. Consequently, in order to solve the boundary-value problem, it is sufficient to find the Fourier coefficients , with which (37) satisfy inhomogeneous Sophie Germain – Lagrange equation (10).1: .

Here with we denote the eigenvalues of the operator :

$$\begin{array}{c}.\end{array}$$

One can find by acting the projection operators: ,

on both sides of equation (38):

It is essential to note that the basis functions lying in different sequences have a different parity. Consequently, due to the fact that the second derivative of any function has the same parity as the function itself, the double sum in the previous relation can be reduced to the sum:

Thus, the differential equation is reduced to an infinite-dimensional system of algebraic equations. By applying the reduction method [88,89] we arrive at the sequence of finite-dimensional systems of linear algebraic equations which can be written in a matrix form: ,

where are matrices × known as Galerkin matrices. Their components, can be obtained via formula 2: .

Vector contains Fourier coefficients of transverse load function (for brevity, we will refer on it as load vector): .

Recall now that we consider functions only from (10), so it is sufficient to calculate only the first Galerkin matrix $1$ and the first load vector $1$, which correspond to even by both variables functions. In this regard hereinafter we will omit the first index and adopt the following notations:

$$\begin{array}{c}.\end{array}$$

(13)

Thus, the problem has been reduced to one matrix equation. In the case of an arbitrary non-symmetric load function, further steps would be similar, however they will require finding of all matrices and vectors .

Proceed now to calculating of Galerkin matrix and load vector components. It is convenient to represent an expression for them using an auxiliary integral $G$3: where denote norms of the basis functions: .

The components of the load vector are calculated by formula (42) is assume $=1$.

Now the Fourier coefficients can be calculated from the first matrix equation system (40) using, for example, Cramer’s rule: .

Here is the identity matrix and are so-called matrix units which have only one unit component with indexes , while other components are zero.

Finally, the solution can be obtained by summing the basis functions with the calculated Fourier coefficients : .

The above limit exists because the original infinite system is fully regular (see [89]). It should be noted that the formulae (47), based on Cramer’s rule, are merely of academic interest, while in numerical realization it is convenient to use more efficient algorithms of matrix inverting.

To illustrate obtained solution model calculations were performed for two type of load function: uniform load, $\equiv 1$ and bump-function, which is defined within the region as polynomial function $={\left(1{-}^2\right)}^2{\left(1{-}^2\right)}^2$ and is zero outside of it. Trial load function and their approximations by reduced Fourier sums are shown in Figure 4. The purpose of these calculations is to show how the convergence rate of the sequence of partial sums depends on the smoothness of the load function and its belonging to the differential operator’s domain (obviously the first function does not belong to the domain, whereas the second one belongs to it).

Results of these calculations are presented in Figure 5 and Figure 6. Subfigures (a) contain deflected shapes of plate. In subfigures (b) the inaccuracy functions for corresponding load functions are shown, which is defined as distributions of unbalanced transversal forces absolute values:

$$\delta =\left|-\sum =1\right|.$$

Subfigures (c) give integral estimations of the inaccuracy functions versus number of accounted terms in reduced sums:

$$\epsilon =\left|1-\sum =1{\displaystyle \underset{-1}{\overset{1}{\int }}\underset{-}{\overset{1}{\int }}dd\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}dd}\right|.$$

In these subfigures, the orange points correspond to sums that are symmetric with respect to the coordinates with any term $=\tilde{1}\tilde{1}$ they contain a symmetric term $\tilde{1}\tilde{1}$.

4. Solution of the Föppl-von Kármán System

Proceed to solving the Föppl-von Kármán boundary-value problems (6) – (13) and (6) – (9). As above we will represent the unknown functions and as decompositions over basis functions (35). However, in this case such expansions do not allow the boundary conditions to be satisfied automatically. In order to overcome this difficulty, we augment the decompositions by the solution of the specific boundary problem.

4.1. Auxiliary Boundary Problem

Consider a homogenous biharmonic equation with the following boundary conditions:

$$\begin{array}{c},\end{array}$$

(15)

where $1,2,1,2$ are arbitrary square-integrable functions. The solution of this problem can be found in the form of a sum: ,

where meets inhomogeneous boundary conditions () but not satisfy homogeneous biharmonic equation (). Another summand , instead, meets homogeneous boundary condition and inhomogeneous equation:

$$=-.$$

Thus, the sum (63) will represent the solution of the boundary problem (), ().

Taking into account symmetry relation (10) it is possible to simplify boundary conditions ():

$$\begin{array}{c},\end{array}$$

(16)

where is an even function, while is an odd one. In this case, the function may be presented as follows:

$$\begin{array}{c}.\end{array}$$

(17)

Substitution of this function in () results in expansions of functions and into the Fourier series:

$$\begin{array}{c}.\end{array}$$

Thus, this solution allows one to satisfy any boundary conditions that may occur in the case of the loading function from (10). In (17) the sums are represented in the most natural form, however for further derivation it is convenient to separate the first term of the second sum and renumber the coefficients:

$$\begin{array}{c}.\end{array}$$

It is now necessary to define the function in such a way that the full solution (63) satisfies the homogeneous biharmonic equation. In order to do this, this function will be expressed as a series over basis functions : .

Substitution of this expression and the boundary solution (17) into the biharmonic equation leads to the following relation:

$$\begin{array}{c}\hfill .\end{array}$$

As before, with projector operator (39), (13) one can transform above relation to an infinite-dimensional system of algebraic equations with respect and : and reduce it to a finite-dimensional one [88,89]. It is convenient to write the reduced system in a matrix form: ,

where $,,$ are vectors of corresponding Fourier coefficients, $0$, $0$ are auxiliary vectors and $,,,$ are auxiliary matrices. Recall that the components of the Galerkin matrix, , are presented by the relation (14). The components of other matrices can be tersely expressed through the corresponding auxiliary integrals $I0(\dots ),\dots ,I5(\dots )$ (see appendix):

$$\begin{array}{c},\\ .\end{array}$$

Here the expressions components of only two matrices are shown. Matrices and can be written in the similar way as and respectively using the replacement of eigenvalues $1$ by $1$ and vice versa. The components of vectors $0$ and $0$ can be obtained from corresponding components of and by assuming equal to zero:

$$0=0,0=0.$$

Thus the solution (63) is fully defined. Further, for brevity, we will use the following notations:

With this notation the solution (17) takes the form: .

As an example, consider boundary-value problem (), () with following boundary values (see Figure 7):

$$\left(t\right)={\left(1-t^2\right)}^2,\left(t\right)=t{\left(1-t^2\right)}^2.$$

Corresponding Airy function and stress distributions are shown in Figure 8.

4.2. Plate with Movable Edges

Consider now first main boundary-value problem (6) – (14a) for a square plate. We will seek the solution in following form:

$$,$$

(18)

Due to the specificity of boundary conditions (14a) the solution of the auxiliary boundary problem can be sufficiently simplified: .

Expressions for the coefficients of the other part of the auxiliary solution (66) can be written in the following form: .

Although the equations (6) are non-linear, bilinearity of Monge – Ampère form (3) allows us to use the expansions (18) and (76) for reducing them to the system of algebraic equation:

$$\begin{array}{c},\\ .\end{array}$$

(19)

Boundary condition () and (15) are satisfied automatically due to an appropriate choice of basis functions in (18), (45). Apart from that, symmetry of problem allows us to consider boundary conditions () satisfied if at least one of them is met. Then, substituting (18) and (22) into, for example, the first condition from (18) we have 4: .

As before one can act on these equations by projector operator and reduce them to an infinite-dimensional system of algebraic equations. For equations (19) it is possible to use operator (39):

$$\begin{array}{c}\hfill ,\\ \hfill & .\hfill \end{array}$$

(20)

while boundary relation (33) should be subjected by action of one-dimensional projector operator. Since an auxiliary part of the Airy function is presented through cosines it is convenient to use the projector operators with cosine kernel: .

The action of these operators on the boundary relation results in:

$$\begin{array}{c}\hfill .\end{array}$$

(21)

For solving system (20), (21) it is necessary to reduce it to finite-dimensional one. The application of reduction method to infinite system of linear equation is validated by well-developed theory rooted in early works of Koyalovich [88,89]. However, there is no rigorous foundation for its appliance to non-linear systems in present. In spite of this fact, reduction method is widely used [16,67,70] and provide results in good agreement with practice. "Why should I refuse a good dinner simply because I do not understand the digestive processes involved?" (Oliver Heaviside).

So, applying the reduction method one can transform (20), (21) to the following system:

$$\begin{array}{c},\\ ,\\ ,,\end{array}$$

(22)

where $0,0,0,$ are matrices × and are cubic matrices $\times \times$, the expression for the components of which are presented in Table 3.

Turning now to the solution of constructed finite-dimensional algebraic system. With account of (72) exclude the auxiliary vector from last equation of (22): .

Moreover, it is possible to express in turn vectors $,$ and constant $0$ from system (22) – (63):

$$\begin{array}{c},\end{array}$$

and substitute them into the first equation of (22), then the whole system can be reduced to the following system of cubic equations: .

Here is the cubic matrix with known components:

$$\begin{array}{c},\end{array}$$

where is an auxiliary cubic matrix deriving through known matrices: .

Although calculation of even finite numbers of components of cubic matrix requires some rather cumbersome calculations, it does not present any insurmountable difficulties.

To solve the cubic equation (35), it is possible to use Newton’s method linearizing the equation at each iteration. Denote a non-linear algebraic operator in LHS of (12) as :

$$:=-\left[+{\displaystyle 12\left(1-\right)}\left(\left(0\right)\right)\left(0{+}^0\right)\right].$$

Let $0$ be the initial value of the solution obtained from linear problem (48) or as expansion of some approximate solution. Then the cubic equation can be reduced to a linear one in neighbour of $0$ as follows:

$$=\left(0\right)+{\displaystyle \partial \partial }|=0\left(-0\right)+o\left(-0\right).$$

(23)

where:

$$\begin{array}{cc}\hfill {\displaystyle \partial \partial }|=0& =\sum =1-\sum =1\left\{\left[0\right]0+\left[\left(0\right)\right]0+\left[\left(0\right)0\right]\right\}\hfill \\ \hfill & -{\displaystyle 2\left(1-\right)}\sum =1\{\left[\left(0\right)0+\left(00\right)\right]\left(0{+}^0\right)0\hfill \\ \hfill & +\left(\left(00\right)0\right)\left(0{+}^0\right)\}.\hfill \end{array}$$

4.3. Plate with Immovable Edges

Consider second main boundary-value problem, defined by boundary conditions (14b). In this case the additional term appear in the right-hand site of (35): .

Hereafter we denote terms containing in corresponding equations in previous case by symbol $\left[\dots \right]$. The first and the third equations of system (19) are modified as follows:

$$\begin{array}{c},\\ ,\end{array}$$

while the second equation remains unchangeable. Apart from that, this system must be added by one more boundary condition equation:

$$\begin{array}{c}.\end{array}$$

As in the previous case, projecting into basis functions with subsequent reduction allows to obtain a finite system of algebraic equations:

$$\begin{array}{cc}& \\ \hfill ,& \\ \hfill ,& \\ \hfill ,& \\ \hfill ,& \\ \hfill .\end{array}$$

(24)

Here the additional vectors, matrices and cubic matrices are used. The expressions for components of some of it are presented in Table 4. Components of vectors $,1,2,$ and matrices $0,0$ can be found from corresponding matrices $,1,2,$ and cubic matrices if assume their last indexes as zero.

Turn to solving this system. First, scalars $0$ can be expressed through other variables as follows:

$$\begin{array}{cc}\hfill 0=& -\{\left[0{+}^0-1\left(+\right)\right]+\left[0{+}^0-1\left(+\right){-}^0{+}^0\right]\hfill \\ \hfill & +\left[\left(0{-}^0-1\right)\right]\}{\left[{\displaystyle 2\left(5-\right)\pi }{+}^0-1\left(0{+}^0\right)\right]}^{-1}\hfill \\ \hfill & =:\mathcal{S}0+\mathcal{S}0+\left(\mathcal{S}0\right).\hfill \end{array}$$

Further, one can express vector :

$$\begin{array}{cc}\hfill =-& {\left[0+\left(+\right)+\mathcal{S}\right]}^{-1}\hfill \\ \hfill & \{\left[\left({+}^{-1}\left(0{+}^0\right)\right)\otimes \mathcal{S}0{+}^{-1}\left(+\right)\right]\hfill \\ \hfill -& \left({+}^{-1}\left(0{+}^0\right)\right)\left(\left(\mathcal{S}0\right)\right)+\left(-1\right)\}\hfill \\ & =:\mathcal{S}+\left(\mathcal{S}\right)+\mathcal{T}\left(\mathcal{S}0\right),\hfill \end{array}$$

where ⊗ is the Kronecker product for vectors, which assigns a matrix to pair of vectors by rule:

$$a\otimes b:=M,Mij=aibj.$$

Then $0,$:

$$\begin{array}{cc}\hfill =& {\{}^{-1}\left[\left(+\right)+\left(+\right)\mathcal{S}+\left(0{+}^0\right)\otimes \left(\mathcal{S}0\mathcal{S}+\mathcal{S}0\right)\right]+\hfill \\ \hfill +& \left({+}^1{-}^2\right)\mathcal{S}+\left({+}^1{-}^2\right)\otimes \left(\mathcal{S}0\mathcal{S}+\mathcal{S}0\right){\}}^{-1}\hfill \\ \hfill & \left\{\right[{(}^{-1}{-}^{-1}\left(+\right)\mathcal{S}-\left({+}^1{-}^2\right)\mathcal{S}+\hfill \\ \hfill +& \left(-1\left(0{+}^0\right)+\left({+}^1{-}^2\right)\right)\otimes \mathcal{S}0\mathcal{S}\left)\right]\hfill \\ \hfill +& [-1\left(+\right)+\left({+}^1{-}^2\right)\mathcal{T}\hfill \\ \hfill +& \left(-1\left(0{+}^0\right)+\left({+}^1{-}^2\right)\right)\left(1+\mathcal{S}0\mathcal{T}\right)\left]\left(\mathcal{S}0\right)\right\}\hfill \\ \hfill =:& \left(\mathcal{S}\right)+\mathcal{T}\left(\mathcal{S}0\right),\hfill \\ \\ \hfill 0=& {\displaystyle 11-}\left\{\right[({\displaystyle 2\pi }\left(1-\right)\left(\mathcal{S}0\mathcal{S}+\mathcal{S}0\mathcal{S}\mathcal{S}+\mathcal{S}0\mathcal{S}+\mathcal{S}0\right)\hfill \\ \hfill -& \left(0{-}^0\right)\left(\mathcal{S}\mathcal{S}+\mathcal{S}\right){+}^0\left)\right]\hfill \\ \hfill +& [{\displaystyle 2\pi }\left(1-\right)\left(\mathcal{S}0\mathcal{T}+\mathcal{S}0\mathcal{S}\mathcal{T}+\mathcal{S}0\mathcal{T}\right)\hfill \\ \hfill -& \left(0{-}^0\right)\left(\mathcal{S}\mathcal{T}+\mathcal{T}\right)\left)\right]\left(\mathcal{S}\right)\}=:\left(\mathcal{S}0\right)+\mathcal{T}0\left(\mathcal{S}0\right).\hfill \end{array}$$

After introducing a set of auxiliary structure $\mathcal{S}$ and $\mathcal{T}$ one can simply reduce the full system to cubic system: ,

where $˜$ is a cubic matrix:

$$\begin{array}{cc}\hfill ˜=& {-}^{-1}+{\displaystyle 12}\left(0{+}^0\right)\mathcal{S}0+-1\left(+\right)++\mathcal{S}\hfill \\ & +\left[-1\left(+\right)++\right]\left(\mathcal{S}\mathcal{S}+\mathcal{S}\right)\hfill \\ & +-1\left(0{+}^0\right){+}^0{+}^0\left(\mathcal{S}0\mathcal{S}+\mathcal{S}0\mathcal{S}\mathcal{S}+\mathcal{S}0\mathcal{S}\right)\hfill \end{array}$$

and $\tilde{T}$ is a vector:

$$\begin{array}{cc}\hfill \tilde{T}=& {\displaystyle 12}\mathcal{T}0\left(0{+}^0\right)+\left[-1\left(+\right)++\right]\mathcal{T}\hfill \\ & +\left[-1\left(+\right)++\right]\left(\mathcal{S}\mathcal{T}+\mathcal{T}\right)\hfill \\ & +\left[-1\left(0{+}^0\right){+}^0{+}^0\right]\left(\mathcal{S}0\mathcal{T}+\mathcal{S}0\mathcal{S}\mathcal{T}+\mathcal{S}0\mathcal{T}+1\right)\hfill \end{array}$$

As in previous case, obtained cubic equation can be linearized by formula similar with (23) and solved by Newton’s method.

5. Numerical Results and Discussion

In order to verify the proposed solution, mechanical tests with thin square completely clamped plates were carried out. The essence of tests is as follows. A plate initially undergo finite strain under given pressure $0$ and then reload by small increments of additional pressure $\delta$. Corresponding increments of normal deflections $\delta$ are determined with holographic interferometry [90]. To this end the plates, fixed in square thick-wall steel chamber were employed. Initial deformed state considered as a reference and recorded on a reference hologram. Superimposing on it another hologram, recorded after the increase in pressure, results in fringe pattern that identify contour lines of deflection increment field $\delta$. Repeating these experiments for various initial pressures and a number of incremental steps, one can get a comprehensive description of the non-linear deformation at different levels of finite strains.

The physical and geometric parameters of the plate are given in Table 5.

Experimentally obtained fields of deflection increments were compared the results of calculations according to the proposed solution. Intermediate theoretical results for initial pressure $0=3$ kPa are shown in Figure 9. Comparison of theoretical and experimental fields of deflection increments for three different initial pressures $0=0$, $0=3$ kPa, $0=6$ kPa and 9 additional loading with equal increments in 10 Pa is shown on chart in Figure 10. On this chart one can see experimentally and theoretically obtained maximal deflections (in the centre of the plate), represented in terms of contour lines quantities. This representation of the displacements is convenient for comparison with the holographic experimental method, because in the latter we directly obtain contour lines of the displacement with a step equal to half the measuring laser wavelength.

A more detailed juxtaposition of the theoretical and experimental contours for 5th and 8th steps of incremental loading are presented in Figure 11. The above results demonstrate a good agreement between theoretical results, obtained with the Föppl-von Kármán model and proposed solution, and the experimental data.

6. Conclusion

To sum up, we would like to highlight the key advantages of proposed solution:

– Computational implementation of the solution reduces to iterative solving of cubic equation system, with all coefficients of this system being determined analytically via finite formulas.

– The theoretical model accurately reproduces the tendency of the plate effective stiffness increasing with the increasing of in-plane stress state.

– For thin plates with an in-plane characteristic size to thickness ratio of less than 1/100 there is a substantially non-linear load-deflection dependence and this should be taken into account when calculating and designing the elements of micro-mechanical systems.