Solution of Solid Mechanic Equilibrium Problems by Power Series

Equilibrium problems are defined by Crandall [1] as

with

where, $L_{2 \mu}(\psi)$ is a partial differential operator of order $2 \mu$, containing only even partial derivatives, and $B_{2 \mu-1}^\zeta(\psi)$ are $\mu$ differential operators of order up to $2 \mu-1$.

A solution of some equilibrium problems is presented in this paper consisting of power series $\psi$ that satisfy identically the homogeneous differential eqn (1) and approximately the boundary conditions in eqn (2). The collocation technique is used at the boundary in such a way that the boundary conditions are satisfied at selected points.

The approximation to the exact solution will depend on the degree $n$ of the elements of the power series used.

Let $\psi$ be a complete power series of degree $n$. Let also $\psi_i, i=0 \ldots \kappa$ be some incomplete power series of degree $\leq n$, such that

If all series $\psi_i$ satisfy a given differential equation, then $\psi$ will also satisfy the said equation. Therefore, if the series $\psi_i$ are mutually independent and their sum forms a complete polynomial expression of degree $n$, then the elements $\psi_i$ can be considered as a base for the solution of the differential equation. The number of elements in the base will be $m$, so that $\kappa=m-1$.

General explicit formulas for $\psi_i$ can be developed for defining the base of solutions for each differential equation. Then the solution of a particular problem is reduced to finding the parameters $\lambda_i$ that satisfy the boundary conditions.

Let us consider the homogeneous differential equation

The elements of the base that satisfies differential equation (4) for power series of degree $n$ are

Note that functions $w_i$ as defined by eqn (5) are mutually independent. This property can be verified by noting that there are no two terms with equal exponents for $x$. Also, the base contains all powers of $x$ and $y$ up to degree $n$. The number of elements of the base is $v$ for all $n$, i.e. $m=v \forall n$.

Let us consider a beam on linear elastic foundation, the equilibrium equation of which is given by Hetényi [2]:

where $w$ : vertical displacement

$k$ : coefficient of soil reaction

$E$ : Young's modulus

$I$ : moment of inertia of beam section.

Introducing $v=4$ and $\rho=k / E I$ into eqn (5), the elements of the solutions' base are

According to eqn (7), the number of parameters $\lambda_i$ that have to be determined is $m=4$. As the differential equation is of order $4(\mu=2$ in eqn (1)), there are two boundary conditions of order up to 3 .

Then the complete solution is given by

Let us consider the beam of length $L$ in Figure 1. Figure 1a shows the displacements $w(x)$ and rotations $\theta(x)=d w(x) / d x$ and Figure 1b shows the shear forces $Q(x)=E I d^3 w(x) / d x^3$ and bending moments $M(x)=E I d^2 w(x) / d x^2$.

Then, from eqn (8) by derivation and evaluation at x = 0 and x = L, relations (9) and (10) are obtained as follows:

or

$W=B \lambda$ (9a)

or

$F=D \lambda$ (10a)

$\lambda=D^{-1} F$ (10b)

From eqns (9a) and (10b), the solution for displacements and rotations is obtained:

Alternatively, the solution for displacements is obtained by introducing eqn (10b) in eqn (8).

Note that this is the exact solution of the differential equation and that it satisfies all boundary conditions. Numerically, the approximation depends on the degree $n$ chosen.

The stiffness matrix K is defined by eqn (12):

From eqn (11) it follows

Given K, the problem is solved using the typical procedure of finite element method (FEM). In this case, the solution is exact.

Let $\phi(x, y)$ be a function that satisfies the homogeneous harmonic equation

The complete base of power series for the harmonic equation is

$\phi_1(x, y)=x$ (15b)

$\phi_2(x, y)=y $ (15c)

$\phi_3(x, y)=x y $ (15d)

$\phi_{3+i}(x, y)=\sum_{j=0}^i\binom{2 i}{2 j}(-1)^j x^{2 i-2 j} y^{2 j} \quad i=1 \ldots\left\lfloor\frac{n}{2}\right\rfloor $ (15e)

$\phi_{3+\left\lfloor\frac{n}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+1}{2 j}(-1)^j x^{2 i-2 j+1} y^{2 j} \quad i=1 \ldots\left\lfloor\frac{n-1}{2}\right\rfloor$ (15f)

$\phi_{3+\left\lfloor\frac{n}{2}\right\rfloor+\left\lfloor\frac{n-1}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+1}{2 j}(-1)^j x^{2 j} y^{2 i-2 j+1} \quad i=1 \ldots\left\lfloor\frac{n-1}{2}\right\rfloor$ (15g)

$\phi_{3+\left\lfloor\frac{n}{2}\right\rfloor+2\left\lfloor\frac{n-1}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i \frac{\binom{2 i+1}{2 j+1}}{(2 i-2 j+1)}(-1)^j x^{2 i-2 j+1} y^{2 j+1} \quad i=1 \ldots\left\lfloor\frac{n-2}{2}\right\rfloor$. (15h)

These expressions can be easily obtained by expansion of the binomial series $(x+i y)^n$ and then by considering that the real and imaginary parts are solutions of the harmonic equation.

The number of elements in the base is

It can be verified that each element of the base satisfies the harmonic equation. It can also be verified by inspection that all powers of the bases are different.

The base contains all powers of x and y up to degree n. In effect, the number of terms in a complete series of maximum exponent n is

From (15), the total number of terms in the base is

Expressions in eqns. (16) and (18) can be verified by taking separately the case $n=$ even and $n=$ odd.

According to eqn (16), the number of parameters $\lambda_i$ that have to be determined is $m=2 n+1$. As the differential equation is of order $2(\mu=1$ in eqn (1)), there is only one boundary condition of order up to 1.

The parameters $\lambda_i$ can be obtained by forcing the solution to satisfy the boundary condition at $m$ points distributed along the boundary.

As can be seen in Connor [3], the classic Saint-Venant's torsion problem can be solved by the Prandtl stress function $\Phi$ in the form

$\tau_{y z}(x, y)=-G \chi \frac{\partial \Phi(x, y)}{\partial x}$. (19b)

The equilibrium equations are then satisfied identically and the compatibility condition gives

At the boundary, the conditions are simply

The twist $\chi$ is related to the twisting moment M by

Let us consider the complete solution as the sum of a particular solution and the homogeneous solution:

The solution for a circular bar (24) can be used as a particular solution:

and the complete series in x and y (25) can be used for the homogeneous solution:

at the boundary

Considering m boundary points so that q = 0…(m − 1), the following m x m system of equations can be formulated:

In matrix notation

and

Finally

with

$a_{1,0}=\lambda_1$ (31b)

$a_{0,1}=\lambda_2$ (31c)

$\alpha_{1,1}=\lambda_3$ (31d)

$a_{2 i-2 j, 2 j}=\lambda_{3+i}\binom{2 i}{2 j}(-1)^j \quad i=1 \ldots\left\lfloor\frac{n}{2}\right\rfloor j=0 \ldots i$ (31e)

$a_{2 i-2 j+1,2 j}=\lambda_{3+\left\lfloor\frac{n}{2}\right\rfloor+i}\binom{2 i+1}{2 j}(-1)^j \quad \quad i=1 \ldots\left\lfloor\frac{n-1}{2}\right\rfloor j=0 \ldots i$ (31f)

$a_{2 j, 2 i-2 j+1}=\lambda_{3+\left\lfloor\frac{n}{2}\right\rfloor+\left\lfloor\frac{n-1}{2}\right\rfloor+i}\binom{2 i+1}{2 j}(-1)^j \quad i=1 \ldots\left\lfloor\frac{n-1}{2}\right\rfloor j=0 \ldots i$ (31g)

$a_{2 i-2 j+1,2 j+1}=\lambda_{3+\left\lfloor\frac{n}{2}\right\rfloor+2\left\lfloor\frac{n-1}{2}\right\rfloor+i}^{\frac{\binom{2 i+1}{2 j+1}}{(2 i-2 j+1)}(-1)^j} \quad i=1 \ldots\left\lfloor\frac{n-2}{2}\right\rfloor j=0 \ldots i$. (31h)

4.2.1.1 Equilateral triangular section

The exact homogeneous solution f of Saint-Venant’s problems for an equilateral triangle (Figure 2) is given in eqn (32); see Sokolnikoff [4]:

Exact homogeneous solution (32) is a polynomial of degree $n=3$; so if the general power series approximation is used, any value of $n \geq 3$ will give the exact solution.

Let us consider the case $n=4$. Then the number of independent parameters to be determined is $m=2 n+1=9$. Therefore this same number of conditions has to be imposed at the boundary. As is shown in Figure 3, every triangle side has been divided into three parts thus obtaining nine boundary points in total. As expected, the exact solution is obtained. The values of $\alpha$ given by (31) associated with $n=4$ are null.

4.2.1.2 Cylindrical bar with a circular slot

Another classical torsion problem, the exact solution of which is not polynomial, is shown in Figure 4; see Sokolnikoff [4].

As is well known, the exact homogeneous solution $\phi$ for this case is

An approximation by power series can be obtained using, for example, $n=20$. Then $m= 2 n+1=41$. Figure 5 shows the shape of functions $\Phi$ obtained by the approximation.

Let us consider a segment of unitary length of the pure torsion problem. Then, the principle of virtual work (PVW) gives the relation (34):

Considering (19) and (23), the following equation is obtained:

Rewriting compactly, we have

$\iint_A\left\langle\begin{array}{ll}\delta \phi_{, x}^H & \left.\delta \phi_{, y}^H\right\rangle\left\{\begin{array}{l}\phi_{, x}^H \\ \phi_{, y}^H\end{array}\right\} d A=\iint_A\left\langle\delta \phi_{, x}^H \quad \delta \phi_{, y}^H\right\rangle\left\{\begin{array}{l}x \\ y\end{array}\right\} d A .\end{array}\right\}$. (35a)

Considering $m-1$ boundary points, so that $q=1 \ldots(m-1)$, the system (36) of ($m-1$ ) $x(m$ - 1) algebraic equations can be formulated (remember that $\Phi^H(x, y)=\phi(x, y)$ according to eqn (25)). Furthermore, because the terms are affected by the first derivatives of $\phi$, the base (15a) does not collaborate in the relation. As the variation corresponds to parameters $\lambda_i$, the following equation can be obtained:

Rewriting compactly, we have

$\iint_A M \phi(x, y) M \phi(x, y)^T d A\left\{\begin{array}{c}\lambda_1 \\ \vdots \\ \lambda_{m-1}\end{array}\right\}=\iint_A M \phi(x, y)\left\{\begin{array}{l}x \\ y\end{array}\right\} d A$ (36a)

or

$D \lambda=F$. (36b)

Let us now assume that we want to develop a finite element of $\eta$ nodes located at its boundary. Evaluating eqn (25) at the $\eta$ boundary point, the values of $\phi$ at those points will be

Then

or, after introducing eqn (38) into (36b), the following relation is obtained:

that is

The relation in eqn (40) has the classical form of the finite element theory, where the unknown parameters are the values of the function at the nodes. Note that complete polynomial expressions of any degree that satisfy the differential equation are used. The maximum degree to be considered will depend, of course, on the number of nodes to be defined in the element.

In the case of an equilateral triangle, the solution by FEM also replicates the exact solution for $\Phi$ shown in Figure 3.

In the case of a cylindrical bar with a circular slot, the solution by FEM provides a better approximation than the solution with boundary conditions.

Figure 6 shows the results by FEM obtained with 24 elements of 6 nodes and 61 nodal points.

The nondimensional warping function $\varphi$ is the harmonic conjugate of the homogeneous stress function, that is,

$\frac{\partial \varphi(x, y)}{\partial y}=-\frac{\partial \phi(x, y)}{\partial x}$. (41b)

Let us consider $\varphi(x, y)$ as a complete power series of degree $n$ :

By integration of eqn (41a) with respect to $x$ and considering eqn (30), the following relation is obtained for $\varphi(x, y)$ :

This can be transformed into

and by integration of eqn (41b) with respect to y and considering eqn (30):

Expressions (44) and (45) should be identical. It can be observed that they have redundant information in the first terms and that the term $f_x(y)$ of expression (43) can be obtained from the redundant terms of expression (45). In the same way, the term $f_y(x)$ can be obtained from the redundant terms of relation (43).

For $\varphi(x, y)$ to be a complete power series of degree $n$, it is necessary to have

Finally, $\varphi(x, y)$ takes the form

with

$\beta_{g, h}=a_{(g-1),(h+1)} \frac{(h+1)}{g}=-a_{(g+1),(h-1)} \frac{(g+1)}{h} \quad g, h=1 \ldots(n-1) \quad(g+h) \leq(n-1)$.

It can be verified that eqn (49) satisfies the harmonic equation.

The warping functions of Figure 7 were obtained by applying eqn (49) to the cases of the triangular and the cylindrical bar with a slot, respectively.

Let $\phi(x, y)$ be a function that satisfies the homogeneous biharmonic equation

The complete base of power series for the biharmonic equation is as follows:

$\phi_1(x, y)=x$ (51b)

$\phi_2(x, y)=y$ (51c)

$\phi_3(x, y)=x^2$ (51d)

$\phi_4(x, y)=x y$ (51e)

$\phi_5(x, y)=y^2$ (51f)

$\phi_6(x, y)=x^3$ (51g)

$\phi_7(x, y)=x^2 y$ (51h)

$\phi_8(x, y)=x y^2$ (51i)

$\phi_9(x, y)=y^3$ (51j)

$\phi_{10}(x, y)=x^3 y$ (51k)

$\phi_{11}(x, y)=x y^3$ (51l)

$\phi_{11+i}(x, y)=\sum_{j=0}^i\binom{2 i+1}{2 j}(-1)^j x^{2 i-2 j+2} y^{2 j} \quad i=1 \ldots\left\lfloor\frac{n-2}{2}\right\rfloor$ (51m)

$\phi_{11+\left\lfloor\frac{n-2}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+1}{2 j}(-1)^j x^{2 j} y^{2 i-2 j+2} \quad i=1 \ldots\left\lfloor\frac{n-2}{2}\right\rfloor$ (51n)

$\phi_{11+2\left\lfloor\frac{n-2}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+3}{2 j} \frac{(i-j+1)}{(i+1)}(-1)^j x^{2 i-2 j+3} y^{2 j} \quad i=1 \ldots\left\lfloor\frac{n-3}{2}\right\rfloor$. (51o)

$\phi_{11+2\left\lfloor\frac{n-2}{2}\right\rfloor+\left\lfloor\frac{n-3}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+1}{2 j} \frac{1}{(2 j+1)}(-1)^j x^{2 j+1} y^{2 i-2 j+2} \quad i=1 \ldots\left\lfloor\frac{n-3}{2}\right\rfloor$ (51p)

$\phi_{11+2\left\lfloor\frac{n-2}{2}\right\rfloor+2\left\lfloor\frac{n-3}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+3}{2 j} \frac{(i-j+1)}{(i+1)}(-1)^j x^{2 j} y^{2 i-2 j+3} \quad i=1 \ldots\left\lfloor\frac{n-3}{2}\right\rfloor$ (51q)

$\phi_{11+2\left\lfloor\frac{n-2}{2}\right\rfloor+3\left\lfloor\frac{n-3}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+1}{2 j} \frac{1}{(2 j+1)}(-1)^j x^{2 i-2 j+2} y^{2 j+1} \quad i=1 \ldots\left\lfloor\frac{n-3}{2}\right\rfloor$ (51r)

$\phi_{11+2\left[\frac{n-2}{2}\right]+4\left\lfloor\frac{n-3}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+1}{2 j} \frac{(2 i+3)}{(2 i-2 j+3)(2 j+1)}(-1)^j x^{2 i-2 j+3} y^{2 j+1} \quad i=1 \ldots\left\lfloor\frac{n-4}{2}\right\rfloor$ (51s)

$\begin{array}{r}\phi_{11+2\left\lfloor\frac{n-2}{2}\right\rfloor+4\left\lfloor\frac{n-3}{2}\right\rfloor+\left\lfloor\frac{n-4}{2}\right\rfloor+i}(x, y)=\sum_{j=0}^i\binom{2 i+1}{2 j} \frac{(2 i+3)}{(2 i-2 j+3)(2 j+1)}(-1)^j x^{2 j+1} y^{2 i-2 j+3} \\ i=1 \ldots\left\lfloor\frac{n-4}{2}\right\rfloor .\end{array}$ (51t)

The number of elements in the base is

It can be verified that each element of the base satisfies the biharmonic equation. Unlike the case of the harmonic equation, in this case there are some bases that repeat powers (compare (51m) with (51n), (51o) with (51p), (51q) with (51r), and (51s) with (51t)).

The base contains all powers of x and y up to degree n. Similarly to the case of the harmonic equation, this is demonstrated by the following identity:

Expressions (52) and (53) can be verified by taking separately the case $n=$ even and $n=$ odd.

According to eqn (52), the number of parameters $\lambda_i$ that have to be determined is $m= 4 n-2$. As the differential equation is of order $4(\mu=2$ in eqn (1)), there are two boundary conditions of order up to 3.

The $\lambda_i$ parameters can be obtained by forcing the solution to satisfy the boundary condition at $N$ points distributed along the boundary.

Using Kirchhoff’s theory, the differential equation governing the bending of thin plates is given by (54):

where, $\quad q$: applied load

$w$: vertical displacement

$D$: plate stiffness ($D=E h^3 / 12\left(1-\sigma^2\right)$).

The exact solution for a simple supported rectangular plate is given by Navier's method; see Timoshenko and Woinowsky-Krieger [5]:

Figure 8 shows the displacements of a simple supported rectangular plate subjected to a uniform load.

Defining as a particular solution of eqn (54)

and considering for the homogeneous solution the bases defined by (51), a solution similar to the one described in section 4.2.1 is possible.

The simple supported conditions are

and the bending moments normal to the boundary = 0 are shown in Figure 9:

$M_y(x, b)=0, \quad \forall 0 \leq x \leq a$ (58b)

$M_x(0, y)=0, \quad \forall 0 \leq y \leq b$ (58c)

$M_x(a, y)=0, \quad \forall 0 \leq y \leq b$, (58d)

where

$M_y(x, y)=-D\left(\frac{\partial^2 w(x, y)}{\partial y^2}+\sigma \frac{\partial^2 w(x, y)}{\partial x^2}\right)$. (59b)

Figure 10 shows the approximate solution considering nine points at the boundary.

A solution of some equilibrium problems in solid mechanics has been presented, consisting of a complete base power series that satisfies identically the homogeneous differential equation that governs the phenomena. The boundary conditions are satisfied only approximately at selected points. This technique is applied to the pure torsion problem giving excellent results for particular cases. The approximation to the exact solution will depend on the number of elements in the base of power series used, which in turn determines the number of points where the boundary conditions have to be satisfied.

The method is extended to the construction of finite elements of high order by selecting boundary nodes in the element and then establishing a relationship between the values of the unknown function at those nodes. In this way the interpolation functions are complete polynomialex pressions that satisfy the differential equation and the approximation will depend on the number of nodes to be used.

Finally, the complete base series for the biharmonic differential equation is presented. The theory is then applied to the plate-bending Kirchhoff’s problem and, in particular, to the case of a rectangular simple supported plate subjected to a uniform load. For other applications, see González [6].