## Abstract

A new nonlinear model for large deflections of a beam is proposed. It comprises the Euler–Bernoulli boundary value problem for the deflection and a nonlinear integral condition. When bending does not alter the beam length, this condition guarantees that the deflected beam has the original length and fixes the horizontal displacement of the free end. The numerical results are in good agreement with the ones provided by the elastica model. Dynamic and two-dimensional generalizations of this nonlinear one-dimensional static model are also discussed. The model problem for an inextensible rectangular Kirchhoff plate, when one side is clamped, the opposite one is subjected to a shear force, and the others are free of moments and forces, is reduced to a singular integral equation with two fixed singularities. The singularities of the unknown function are examined, and a series-form solution is derived by the collocation method in terms of the associated Jacobi polynomials. The procedure requires solving an infinite system of linear algebraic equations for the expansion coefficients subject to the inextensibility condition.

## 1. Introduction

The background motivation of this study is provided by modelling of bioinspired microsensors as hair-like flexible micropillars and pillar arrays. The sensors serve for the wall shear stress measurements in the boundary layer for laminar and turbulent flow. For these measurements, some techniques based on optical imaging and image processing need the values of the displacements of the hair-like sensor tip [1,2]. The procedure requires the solution of the dynamic linear Euler–Bernoulli (EB) beam equation with allowance for damping. Some other measurement techniques take the sensor's output as the bending moment at the base and employ a linear [3] or nonlinear [4] viscoelastic Kelvin–Voigt model of an EB beam coupled to an unsteady non-uniform flow environment. Both models, the former linear one and the latter nonlinear model, have a good accuracy for small deflections only.

It is known that the solution for large deflections of a beam cannot be obtained from the elementary linear EB theory. This theory neglects the square of the curvature derivative and disregards shortening of the moment arm due to the deflection. If the material of the beam remains linear, and the deflections are large, then the exact differential equation *D* d*ϕ*/d*s*=*M* needs to be integrated. Here, *ϕ*(*s*) is the deflection angle, d*ϕ*/d*s* is the curvature, *M* is the moment and *D* is the flexural rigidity. The exact slope of the elastic curve recovered from this equation is called the elastica [5]. When one of the beam ends is clamped, the other is free and a load is applied to the free end, the exact solution was found by Barten [6,7] and Bisshopf & Drucker [8]. The case of a uniformly distributed load was approximately treated by Rohde [9]. Frisch-Fay [10] analysed the case when loading is modelled by a system of *n* concentrated loads and reduced the problem to a system of 2*n*−1 transcendental equations. Wang *et al.* [11] presented a series-form solution for the case when a load is applied to the free end (the Barten problem). To the author's knowledge, exact solutions of the elastica model for an arbitrary continuous load and its generalizations for the dynamic and two-dimensional cases remain out of reach.

The main goal of this paper is to propose and analyse an alternative nonlinear model which is applicable for large deflections and may be used not only in the static but also in the dynamic case and for an arbitrary load. In addition, as some sensors are manufactured as rectangular plates [12], we aim to develop an accurate nonlinear plate bending model when bending does not alter the plate size.

In §2, we consider the Frisch-Fay elastica model problem for *n* concentrated loads, reduce it to a single transcendental equation and solve it numerically. In §3, we propose the static and dynamic nonlinear models for an EB beam. The problems are solved exactly for any continuous load. We show that the solution for an EB beam when bending does not alter the beam length provides the same level of accuracy for large deflections as the elastica model does. However, the solution procedure is simpler, and the model admits the dynamic generalization. In §4, we analyse bending of a Kirchhoff rectangular plate when one side is clamped, the opposite one is subjected to a shear force and the others are free of moments and forces. On employing the method of finite integral transforms, we derive a governing singular integral equation with two fixed singularities. The singularities of the solution are studied, and an approximate series-form solution is obtained in terms of the Jacobi polynomials. The series coefficients solve an infinite linear algebraic system. The abscissas of the deflected points of the plate are fixed by inversion of the integral inextensibility condition.

## 2. Large deflection of a beam: the elastica model

A beam of length *L* is subjected to an arbitrary distributed load modelled by a system of *n* vertical concentrated loads *P*_{1},…,*P*_{n}. The loads *P*_{j} are applied at some points *A*_{j} (*j*=1,…,*n*), *A*_{0} and *A*_{n} are the clamped and free ends of the beam, and the flexural rigidity of the beam is a piecewise constant function, *D*(*x*)=*D*_{j}, *x*∈*A*_{j−1}*A*_{j} (*x* is the horizontal coordinate), *D*_{j}=*E*_{j}*I*_{j}, *E*_{j} is the Young modulus, and *I*_{j} is the moment of inertia. The main assumptions of the model are (i) the curvature is proportional to the bending moment (the Bernoulli–Euler theorem), and (ii) bending does not alter the beam length.

We choose downward deflections to be positive and denote the slope angle by *ϕ*(*s*) and the arc length by *s*. Let *R* and *M* be the total reaction and the moment at the clamped end. Taking into account the shortening of the moment arm as the points *A*_{j} (*j*=1,…,*n*) deflect, we can define the bending moments *M*_{j} on the segment *A*_{j−1}*A*_{j} as
*L*_{j} is the length of the arc *A*_{0}*A*_{j}, and Δ_{j} is the horizontal displacement of the point *A*_{j} (the displacement to the left is positive). On the other hand, as the bending moment is proportional to the beam curvature d*ϕ*/d*s*, we have another system of equations
*s*, using the relations (2.1) and *ϕ*/d*s* and integrating them afterwards, we have
*s*∈*A*_{j−1}*A*_{j}, and
*l*_{j} be the length of the arc *A*_{j−1}*A*_{j} and denote *l*_{j} and *α*_{j} are constants, and from (2.3) we have
*ϕ*_{j}∈[0,*π*/2) is the slope angle at the point *A*_{j}. To determine the unknowns *ϕ*_{j}, we introduce new parameters *k*_{j} and *θ*_{j} and a function *θ*,
*ϕ*/d*s* at the point *A*_{n} vanishes, we have *θ*_{n}=*π*/2. In the new notations, the integrals *α*_{j} in (2.5) can be written as
*n* equations for the 2*n*−1 unknowns *θ*_{j} ( *j*=1,…,*n*−1) and *k*_{j} ( *j*=1,…,*n*)
*x*) and *F*(*x*,*k*) are the elliptic sine and elliptic integral of the first kind, respectively.

The continuity of curvature at the points *A*_{j} (*j*=1,…,*n*−1) yields us *n*−1 extra conditions,
*n*−1 equations we derived form recurrence relations for the parameters *k*_{j} and *θ*_{j}. Given *θ*_{n}=*π*/2 define
*k*_{n}. Then continue this by determining sequently for *j*=*n*−1,*n*−2,…,2,
*j*=1, yields
*F*(*θ*_{n},*k*_{n})=**K**(*k*_{n}) is the complete elliptic integral of the first kind. Now, as *θ*_{0} is defined by (2.6), and all the parameters, *k*_{1},…,*k*_{n−1} and *θ*_{1},…,*θ*_{n−1}, are expressed through the single parameter *k*_{n}, we deduce a single transcendental equation with respect to *k*_{n}. It reads
*k*_{n} from this equation numerically, we can recover the other unknown parameters by equations (2.13)–(2.16). In the particular case *n*=1, when a concentrated vertical load *P*_{1}=*P* is applied to the free end of the beam, and the rigidity is a constant *D*, the transcendental equation with respect to the single parameter *k*_{1}=*k* has the form
*α*+**K**)=sn(*α*+**K**).

Denote next the deflection of the point *A*_{j} with respect to the previous point *A*_{j−1} by *δ*_{j}. Then by integrating the equation
*ϕ*_{j−1},*ϕ*_{j}], we determine
*δ* of the free end *A*_{n} with respect to the origin
*A*_{j} relatively to the point *A*_{j−1} can be determined by integrating the equation
*ϕ*_{j−1},*ϕ*_{j}]. This ultimately gives
*M*_{1}=*M*,

Figures 1 and 2 show the results of calculations of the deflection and the horizontal displacement of the free end versus the parameter *D*_{j}=*D*_{0}, loading is modelled by the set of concentrated loads *P*_{j}=*P*_{0} (*j*=1,2,…,20) and *δ*(⋅) is the Dirac function. Its solution yields the deflection of the free end *x*=*L*

It is observed from figures 1 and 3 that the difference between the deflection, the horizontal displacement of the free end and the moment at *x*=0 computed according to the elastica and the linear theory drastically increase as the parameter *α* grows. Figure 2 shows that if the rigidity of the beam, *D*(*x*), is a decreasing function, then the deflection and horizontal displacement of the free end are greater than the ones for a beam with the constant rigidity *D*(0).

## 3. Nonlinear model for an Euler–Bernoulli beam

### (a) Static model

Assume that the left end of a beam 0<*x*<*L* is clamped, its right end is subjected to a concentrated vertical load *P*_{n}, while the beam itself is subjected to a distributed normal load *p*(*x*). The bending rigidity is constant, and the deflection *w*(*x*) is governed by the EB equation. The main assumption of this model is that bending may alter the beam length, and its new length, *w*(*x*) that solves the following boundary-value problem:
*G*(*x*,*ξ*) is the Green function defined by
*L*_{*}, which can be computed by the formula
*q*(*x*)≡0, condition (3.6) is equivalent to the inversion of the elliptic integral
*q*≡0) for *q*≡0,

The profile of the deflected beam can be reconstructed as follows. We split the beam into *n* equal segments *x*_{i−1}*x*_{i}, *i*=1,…,*n*, *x*_{0}=0 and *x*_{i}=*iL*/*n*. Assume that bending does not alter the length *l*_{i} of each segment *x*_{i−1}*x*_{i}. The abscissa, *x*_{i} is determined by solving the transcendental equation

To compare the nonlinear EB model with the elastica model previously analysed, we consider the case when the beam is inextensible, *P*_{n} is applied to the free end only, *P*_{1}=⋯=*P*_{n−1}=0 and *q*(*x*)≡0. In this case, all the parameters *k*_{i} of the elastica model are the same, *k*_{1}=⋯=*k*_{n}=*k*. To recover the coordinates *A*_{i}, we need to solve the system of recurrence relations (2.13) for *θ*_{i}, the transcendental equation (2.17) for *k* and determine *δ*_{i} and *x*_{i},0) are determined by inversion of the elliptic integral

Our calculations implemented for the elastica and the nonlinear EB models for a uniform beam when the load is applied to the free end only are presented in figures 4–6. The discrepancy between the deflection of the free end (figure 4) becomes notable only for large values of the parameter *M*_{e} and *M*_{EB} at the clamped end associated with the elastica and our nonlinear model, respectively. It is seen that for 0<*α*^{2}<20, *M*_{e}/*M*_{EB}∈(0.985,1.005). The profile of the bending beam is presented in figure 6 for the case *α*=1. The solid line corresponds to the elastica curve. It was recovered by computing the coordinates of the points *A*_{j} by formulae (3.10). The dashed line in Figure 6 is recovered by the nonlinear EB model with the horizontal and vertical coordinates being computed from (3.11) and (3.12), respectively. Again, the nonlinear theory is in good agreement with the elastica model.

### (b) Nonlinear dynamic model

A dynamic analogue of the nonlinear model for an inextensible EB beam 0<*x*<*L* reads as follows.

Find a function *w*(*x*,*t*) which solves the EB equation
*m* is the mass per unit length, and *q*(*x*,*t*)≡0 and *ω* is the frequency of vibration). We next split *β*=(*mω*^{2}/*D*)^{1/4}. The function *v*(*x*) can be easily found
*v*(*t*) into condition (3.17) and solving the transcendental equation with respect to *M*(*t*), the horizontal and vertical displacements, and therefore the profile *w*(*x*,*t*) of the beam.

## 4. Bending of a plate

### (a) Formulation

The upper surface of a thin plate *Π*={0<*x*<*L*,−*b*<*y*<*b*,−*h*/2<*z*<*h*/2} is subjected to a normal load *q*(*x*,*y*), |∂*q*/∂*y*|≪|*q*|, (*x*,*y*)∈*Π*. The side {*x*=0,−*b*<*y*<*b*} is clamped, a load *P*(*y*) is applied to the side {*x*=*L*,−*b*<*y*<*b*}, whereas the other two sides of the plate are free of moments and forces.

We assume that (i) the normal stress *σ*_{z} in the middle surface is small enough to be neglected, (ii) plane cross sections normal to the middle surface remain plane and normal to the middle surface after deflection (the Kirchhoff hypothesis), (iii) the deformations in the middle surface are negligible, (iv) the deflection variation in the *y*-direction is much smaller than that in the *x*-direction, and (v) bending does not alter the plate length *L*.

To formulate the boundary conditions, we need the expressions for the bending moments with respect to the axes *y*, *M*_{1} and *x*, *M*_{2},
*N*_{1} and *N*_{2} acting on the sides *x*=*L* and *y*=±*b*,
*x*=*L* and *y*=±*b*,
*D*=*Eh*^{3}[12(1−*ν*^{2})]^{−1} is the cylindrical rigidity of the plate, *E* is the Young modulus and *ν* is the Poisson ratio. The free edge boundary conditions on the sides *y*=±*b* generally excepted in the Kirchhoff theory of plates require *M*_{2}=0 and *N*_{2}+∂*H*_{2}/∂*x*=0. The second boundary condition takes into account the shear force due to the presence of the stresses *τ*_{yz} and the torsion moment due to the tangential stresses *τ*_{xy}. Also, it allows to avoid dealing with the three boundary conditions *M*_{2}=0, *N*_{2}=0 and *H*_{2}=0, *y*=±*b*. For the same reason, on the side *x*=*L*, we have *M*_{1}=0 and *N*_{1}+∂*H*_{1}/∂*y*=*P*(*y*). For simplicity, we assume that *q*(*x*,*y*)=*q*(*x*,−*y*) for all *x* and *y*. Because of this assumption, it is sufficient to consider a half of the plate, say the one with positive *y*. Then the boundary conditions on the line *y*=0 are the standard symmetry conditions. On the clamped edge, the deflection and the deflection slope vanish.

The deflection *w*(*x*,*y*) is the solution of the Germain equation

To satisfy the fifth assumption of the model, we impose the inextensibility condition

### (b) Singular integral equation with fixed singularities

The approach we aim to apply is based on the method of integral transforms and eventually leads to a single integral equation. It is convenient to apply the finite cosine integral transform with respect to *y*
*β*_{n}=*πn*/*b*, and introduce a new function
*a priori* that
*y*=0 in (4.5) and only the second boundary condition on the side *y*=*b*. The first condition, *M*_{2}=0, will be later used for deriving a governing integral equation for the unknown function *χ*(*x*). Using the fundamental function of the differential operator in (4.10) and the boundary conditions, we derive the Green function
*w*(*x*,*y*). It is directly verified that the solution found solves equation (4.4) and meets all the boundary conditions (4.5) except for *w*_{yy}+*νw*_{xx}=0, *y*=*b*,0<*x*<*L*. Employing formulae (4.15) and (4.7) maps this boundary condition into an integral equation. Denote
*Φ*_{n}(*x*−*ξ*) generates a Cauchy-type singular kernel, whereas the second function, *F*_{n}(*x*,*ξ*), generates a kernel with fixed singularities at the points *x*=*ξ*=0 and *x*=*ξ*=*L*. These functions take the form
*n*=0,1,…). For all *x* and *ξ* not simultaneously close to *L*,
*x*≤*L*−*ϵ*_{1}, 0≤*ξ*≤*L*−*ϵ*_{2}, *ϵ*_{1} and *ϵ*_{2} are positive numbers. Hence, in a neighbourhood of the point *x*=*ξ*=0,
*x*=*ξ*=*L*. We have the representation
*x* and *ξ* not simultaneously close to 0. Here, *ξ*_{1}=*L*−*ξ*, *x*_{1}=*L*−*x* and *ε*_{1}≤*x*≤*L*, *ε*_{2}≤*ξ*≤*L* and *ε*_{j}>0, *j*=1,2, respectively. In a neighbourhood of the point *x*=*ξ*=*L*, we derive
*L*]
*K*(*x*,*ξ*) is a regular kernel having a series representation whose rate of convergence is exponential,

### (c) Solution of the integral equation

The integral equation derived in §4*b* does not admit a closed form solution. An approximation solution can be derived by the method of collocation or the method of orthogonal polynomials. For both methods, it is essential to study the behaviour of the unknown function *χ*′(*x*) at the ends *x*=0 and *x*=*L*. First, we analyse the two integrals
*x*→0^{+} and
*x*→*L*^{−}. It is an easy matter to deduce from equation (4.30) that *α* and *β* are not integer. For *α*≠0,±1,…,
*ϕ*_{j}(*x*) (*j*=0,1,2) are functions bounded in a neighbourhood of the point *x*=0, and *ϕ*_{j}(*x*)=*o*(*x*^{α}), *x*→0^{+}. It follows from the integral equation (4.30) and the asymptotic relation (4.35) that *χ*′(*x*)∼*A*_{0}*x*^{α} as *x*→0^{+}, and *α* is the root of the transcendental equation
*α*>0. By denoting *p*=*α*+1, we deduce the equation that coincides with the equation for a wedge plate with free-clamped mixed boundary conditions analysed by the separation of variables method by Uflyand [13]. It turns out that in the strip 0<Re *α*<1, there is only one root of equation (4.36), *α*=*α*_{1}+*iα*_{2}, and

A similar analysis implemented in a neighbourhood of the point *x*=*L* shows that *χ*′(*x*)∼*A*_{1}(*L*−*x*)^{β}, *x*→*L*^{−} and *β* cannot be integer. For the integral *I*_{1}(*x*), we have
*ϕ*_{3}(*x*) is a function bounded in a neighbourhood of the point *x*=*L*, and *ϕ*_{3}(*x*)=*o*((*L*−*x*)^{β}), *x*→*L*^{−}. The parameter *β* solves the equation
*β*>0 and Re *β* is the smallest among the real parts of all the roots in the half-plane Re *β*>0. This equation is consistent with the one obtained by Uflyand [13] for a wedge plate whose both sides are free of moments and forces. Our numerical results show that in the strip 0<Re *β*<1, there is only one root. It is real and

We next construct an approximate solution to the integral equation (4.30). Because of the asymptotics at the ends, the function *χ*′(*Lx*) admits the expansion
*c*_{m} are complex coefficients to be determined. To find the coefficients *c*_{m}, we employ the collocation method. Upon substituting expansion (4.39) into the integral equation (4.30), we derive
*K*_{m}(*x*) is regular and can be evaluated by standard quadrature formulae. The other two integrals, *S*_{m}(*x*) and *J*_{m}(*x*), are singular. On using the Mellin convolution theorem and the theory of residues [14,15], we express the integral *S*_{m}(*x*) through the hypergeometric function
*J*_{m}, we need *γ*=0,1,2 only. For integer *γ*, it is possible to express the generalized hypergeometric functions _{3}*F*_{2} in (4.43) through the function *F*. In view of formulae 7.512(12) and 9.111 in [16], we deduce
*x* close to 1, we can employ formula 9.131(1) in [16] that yields
*x*_{n}=*n*/*N*. The solution of this system of linear algebraic equations approximately defines the coefficients *c*_{m} and therefore approximates the function *χ*′(*Lx*)
*y*-direction and determines the function *w*_{x}(*x*,*y*) in (4.6) is expressed through the function *χ*_{N}′(*ξ*) as

## 5. Conclusion

We have revised the Frisch-Fay problem for the static nonlinear elastica model when bending does not alter the beam length, *L*, the beam rigidity is a piecewise constant function and loading is a set of normal concentrated loads. We have converted the problem into a system of 2*n*−1 recurrence relations involving elliptic functions and then to a single transcendental equation solved numerically. It is known that the elastica model cannot be generalized to the case of a continuous load, needless to say to the dynamic case or bending of a bounded plate. We have proposed a nonlinear model for bending of a beam that comprises the standard boundary value problem for the linear EB equation and the nonlinear condition
*L*_{*} is the deflected beam length computed by the linear EB model. If bending does not alter the beam length, then *χ*′(*x*)=(∂^{2}/∂*x*∂*y*)*w*(*x*,*b*) (*w* is the deflection) and shown that in a neighbourhood of the clamped edge *x*=0 the solution oscillates and vanishes, while in the vicinity of the free end *x*=*b* the function *χ*′(*x*) is monotonically decaying to zero. An approximate solution has been derived by expanding *χ*′(*x*) in terms of the Jacobi polynomials with a weight and employing the collocation method. As in the one-dimensional models, the final step of the procedure is to verify that bending does not alter the plate length, *L*. This condition is the two-dimensional counterpart of the transcendental equation in the nonlinear EB beam model.

## Funding statement

This work was partly funded by American Society for Engineering Education (Air Force Summer Faculty Fellowship Programme).

- Received January 24, 2014.
- Accepted July 23, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.