Nonlinear bending models for beams and plates

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ϕ/ds=M needs to be integrated. Here, ϕ(s) is the deflection angle, dϕ/ds 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 2n−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₁,…,P_n. The loads P_j are applied at some points A_j (j=1,…,n), A₀ 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−1A_j (x is the horizontal coordinate), D_j=E_jI_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−1A_j as

Here, L_j is the length of the arc A₀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ϕ/ds, we have another system of equations On differentiating these equations with respect to s, using the relations (2.1) and $\mathrm{d}x/\mathrm{d}s=\cos \varphi$, multiplying the equations obtained by dϕ/ds and integrating them afterwards, we have Here, ${\hat{\varphi }}_j\in [0,\pi /2)$ is a free constant due to integration, $\sin {\hat{\varphi }}_j>\sin \varphi$ for s∈A_j−1A_j, and Let l_j be the length of the arc A_j−1A_j and denote ${\alpha }_j=l_j\sqrt{Q_j/D_j}.$ As the bending does not alter the beam length, l_j and α_j are constants, and from (2.3) we have where ϕ_j∈[0,π/2) is the slope angle at the point A_j. To determine the unknowns ${\hat{\varphi }}_j$ and ϕ_j, we introduce new parameters k_j and θ_j and a function θ, As the curvature dϕ/ds at the point A_n vanishes, we have ${\hat{\varphi }}_n={\varphi }_n$ and θ_n=π/2. In the new notations, the integrals α_j in (2.5) can be written as This brings us the first n equations for the 2n−1 unknowns θ_j ( j=1,…,n−1) and k_j ( j=1,…,n) or, equivalently, by inverting the elliptic integrals, Here, sn(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,

In view of (2.3), they are equivalent to the equations These equations can conveniently be rewritten as We assert that the system of 2n−1 equations we derived form recurrence relations for the parameters k_j and θ_j. Given θ_n=π/2 define and as functions of k_n. Then continue this by determining sequently for j=n−1,n−2,…,2, The last step, j=1, yields Note that F(θ_n,k_n)=K(k_n) is the complete elliptic integral of the first kind. Now, as θ₀ is defined by (2.6), and all the parameters, k₁,…,k_n−1 and θ₁,…,θ_n−1, are expressed through the single parameter k_n, we deduce a single transcendental equation with respect to k_n. It reads On having determined 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₁=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₁=k has the form where $\alpha =L\sqrt{P/D}$. Here, we used the property of the elliptic sine sn(−α+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

on the interval [ϕ_j−1,ϕ_j], we determine This can be expressed through the elliptic integrals of the second kind By summing up these quantities, we find the deflection δ of the free end A_n with respect to the origin The horizontal displacements ${\tilde{\mathrm{\Delta }}}_j$ of the point A_j relatively to the point A_j−1 can be determined by integrating the equation on the interval [ϕ_j−1,ϕ_j]. This ultimately gives and the horizontal displacement of the free end is Finally, we write down the bending moment at the origin, M₁=M, As $\sin {\hat{\varphi }}_1=2k_1^2-1$, we obtain

Figures 1 and 2 show the results of calculations of the deflection and the horizontal displacement of the free end versus the parameter ${\alpha }_0^2$. In figure 1, the beam is uniform, D_j=D₀, loading is modelled by the set of concentrated loads P_j=P₀ (j=1,2,…,20) and ${\alpha }_0=L\sqrt{P_0/D_0}$. The linear EB model is governed by the boundary value problem

where δ(⋅) is the Dirac function. Its solution yields the deflection of the free end x=L and bending moment at the clamped end The horizontal displacement identically equals zero in the linear theory.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

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).

• Download figure
• Open in new tab
• Download PowerPoint

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, $\tilde{L}\ge L$, is prescribed. We thus aim to find the function w(x) that solves the following boundary-value problem:

and satisfies the condition where $\tilde{x}=L-\mathrm{\Delta }$ is the abscissa of the free end of the bent beam, and Δ is a parameter to be determined. In the particular case, when bending does not alter the length of the beam, $\tilde{L}=L$. The solution to the problem (3.1) can be easily found, where G(x,ξ) is the Green function defined by The classical beam theory assumes that $\tilde{x}=L$ and therefore deals with a beam whose bending results in the extension of the beam. As $\tilde{x}=L$, the altered length $\tilde{L}$ attains its maximum, L_*, which can be computed by the formula If the parameter $\tilde{L}$ is prescribed and $L\le \tilde{L}<L_{\ast }$, then the parameter $\tilde{x}$ is unknown and has to be determined from condition (3.2) which can be specified to read where In the particular case q(x)≡0, condition (3.6) is equivalent to the inversion of the elliptic integral On solving numerically the transcendental equation (3.6) (equation (3.8) if q≡0) for $\hat{x}$, we can compute the actual deflection of the free end by employing formula (3.3) for $x=\hat{x}L$. The bending moment at the clamped end due to the shortening of the moment arm is obviously smaller than in the classical EB theory. In the case when q≡0, $M=P_n\hat{x}L$.

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

After we have found the abscissas ${\tilde{x}}_i$, we employ them to recover the deflections $w({\tilde{x}}_i)$ by formula (3.3).

To compare the nonlinear EB model with the elastica model previously analysed, we consider the case when the beam is inextensible, $\tilde{L}=L$, the load P_n is applied to the free end only, P₁=⋯=P_n−1=0 and q(x)≡0. In this case, all the parameters k_i of the elastica model are the same, k₁=⋯=k_n=k. To recover the coordinates $(x_i^{\ast },{\delta }_i^{\ast })$ of the deflected points 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 ${\tilde{\mathrm{\Delta }}}_i$ by (2.21) and (2.24). Then

In the case of the nonlinear EB model, the coordinates $({\tilde{x}}_i,w({\tilde{x}}_i))$ of the deflected points (x_i,0) are determined by inversion of the elliptic integral and by the formula

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 $\alpha =L\sqrt{P_n/D}$. Figure 5 shows the ratio of the bending moments M_e and M_EB at the clamped end associated with the elastica and our nonlinear model, respectively. It is seen that for 0<α²<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.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

(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

subject to the boundary conditions the initial conditions and the condition of inextensibility Here, m is the mass per unit length, and $\tilde{x}(t)$ is to be determined. To avoid unnecessary technicalities and preserve the main feature of the model, the condition of the inextensibility, we assume that q(x,t)≡0 and $P(t)=P_n\sin \omega t$ (ω is the frequency of vibration). We next split $w(x,t)=v(x)\sin \omega t$ and analyse the following nonlinear harmonic vibration model Here, β=(mω²/D)^1/4. The function v(x) can be easily found Upon substituting the derivative of the function v(t) into condition (3.17) and solving the transcendental equation with respect to $\tilde{x}(t)$, we determine the actual moment arm $\tilde{L}(t)=L-\tilde{x}(t)$. The quantity of interest, the bending moment at the clamped end $M(t)=D{v}^{″}(0)\sin \omega t$, therefore has the form Note that the vibration frequencies, the zeros of the function $\cos \beta L\cosh \beta L+1$, coincide with those defined by the classical liner EB model. What is different is the moment 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₁ and x, M₂,

the shear forces N₁ and N₂ acting on the sides x=L and y=±b, and the torsion moments applied on the sides x=L and y=±b, Here, D=Eh³[12(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₂=0 and N₂+∂H₂/∂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₂=0, N₂=0 and H₂=0, y=±b. For the same reason, on the side x=L, we have M₁=0 and N₁+∂H₁/∂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

subject to the boundary conditions

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

The function $\tilde{x}(y)\in (0,L)$ is to be determined by inversion of the integral in (4.6).

(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 It is possible to establish a priori that To do this, we consider two auxiliary bending problems for a quarter plane, problems 1 and 2. In the former problem, one side is clamped, and the second one is free of forces and moments. In problem 2, both sides are free of forces and moments. These two problems admit closed form solutions by the Mellin transform, which imply the relations (4.9). Using this result, we can show that the integral transform (4.7) maps problems (4.4) and (4.5) into the following one-dimensional boundary value problem: where Note that on applying the cosine transform, we integrated by parts and used both boundary conditions on the side y=0 in (4.5) and only the second boundary condition on the side y=b. The first condition, M₂=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 Here, With the Green function at hand, we can now write down the solution to the boundary value problem (4.10) Integration by parts transforms this formula into a new form more efficient for our purposes where The inverse cosine transform (4.7) of the function (4.15) defines the deflection 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 Then the integral equation can be written as where We wish to split the function ${\hat{G}}_n(x,\xi )={\beta }_n^2{\tilde{G}}_n(x,\xi )-\nu {\mathrm{\partial }}^2{\tilde{G}}_n/\mathrm{\partial }{x}^2(x,\xi )$ into two parts It will be shown that the function Φ_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 and To understand the structure of the kernel of the integral equation (4.18), we summarize the series. It is helpful to have the following formulae: For the first part of the kernel, we discover It is seen that up to a regular function, it coincides with the Cauchy kernel, Analyse next the second part of the kernel generated by the functions ${\hat{F}}_n(x,\xi )$ (n=0,1,…). For all x and ξ not simultaneously close to L, where $\hat{F}(x,\xi )$ is an exponentially convergent series for 0≤x≤L−ϵ₁, 0≤ξ≤L−ϵ₂, ϵ₁ and ϵ₂ are positive numbers. Hence, in a neighbourhood of the point x=ξ=0, A similar analysis can be implemented to recover the fixed singularity of the kernel at the point x=ξ=L. We have the representation valid for all x and ξ not simultaneously close to 0. Here, ξ₁=L−ξ, x₁=L−x and $\tilde{F}(x,\xi )$ is an exponentially convergent series for ε₁≤x≤L, ε₂≤ξ≤L and ε_j>0, j=1,2, respectively. In a neighbourhood of the point x=ξ=L, we derive The three relations (4.25), (4.27) and (4.29) when combined with (4.18) and (4.20) amount to the singular integral equation with fixed singularities at both ends of the interval [0,L] where 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 §4b 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

as x→0⁺ and as x→L⁻. It is an easy matter to deduce from equation (4.30) that α and β are not integer. For α≠0,±1,…, On employing these relations, we obtain Here and in (4.34), ϕ_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₀x^α as x→0⁺, and α is the root of the transcendental equation whose real part is positive and the smallest among the real parts of all roots of the equation lying in the half-plane Re α>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), α=α₁+iα₂, and $0<{\alpha }_1<\frac{1}{2}$ (table 1).

A similar analysis implemented in a neighbourhood of the point x=L shows that χ′(x)∼A₁(L−x)^β, x→L⁻ and β cannot be integer. For the integral I₁(x), we have

where ϕ₃(x) is a function bounded in a neighbourhood of the point x=L, and ϕ₃(x)=o((L−x)^β), x→L⁻. The parameter β solves the equation Re β>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 $\beta \in (\frac{1}{2},1)$ (table 1).

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

where $P_m^{\alpha ,\beta }(x)$ are the Jacobi polynomials, and 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 where The integral 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 and the integral $H_m^{\alpha ,\beta ,\gamma }(x)$ through the generalized hypergeometric functions To compute the coefficients J_m, we need $H_m^{\alpha ,\beta ,\gamma }$ for γ=0,1,2 only. For integer γ, it is possible to express the generalized hypergeometric functions ₃F₂ in (4.43) through the function F. In view of formulae 7.512(12) and 9.111 in [16], we deduce To compute the hypergeometric function in (4.44) for x close to 1, we can employ formula 9.131(1) in [16] that yields We next discretize the functional equation (4.40) using a uniform mesh where 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) The final step in the procedure is to fulfil condition (4.6) which guarantees that the plate does not extend in the y-direction and determines the function $\tilde{x}(y)$. The derivative 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 2n−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

where $\tilde{L}$ is a prescribed parameter, $L\le \tilde{L}<L_{\ast }$, and L_* is the deflected beam length computed by the linear EB model. If bending does not alter the beam length, then $\tilde{L}=L$. The solution of this nonlinear problem has been found explicitly for any load, and the final transcendental equation is equivalent to inversion of an integral. In the case when the load is applied to the free end only, one needs to invert an elliptic integral. The numerical results for the deflection, bending moment and the shape of the deflected beam corresponding to the nonlinear EB model are in good agreement with the results of the elastica model, and both models are more accurate than the elementary linear EB model. Based on our static nonlinear model, we have proposed its dynamic version and derived its exact solution. In addition, we have analysed bending of a rectangular Kirchhoff plate. One of its sides is clamped and the others are free of forces and moments. Although the governing equation, the Germain equation and the boundary conditions are linear, the model is nonlinear due to the nonlinear inextensibility condition. The problem has been reduced to a singular integral equation with two fixed singularities. We have examined the singularities of the solution of the equation, χ′(x)=(∂²/∂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).

Footnotes