Inverse Problems for an Euler-Bernoulli Beam: Identification of Bending Rigidity and External Loads

Hiroaki Katori

Department of Vehicle System Engineering, Faculty of Science and Technology, Meijo University, Nagoya, Japan.

DOI: 10.4236/wjm.2018.85014                

Abstract

We present a method for identifying the flexural rigidity and external loads acting on a beam using the finite-element method. We used mixed beam elements possessing transverse deflection and the bending moment as the primary degrees of freedom. The first step is to determine the bending moment from the transverse deflection and boundary conditions. The second step is to substitute the bending moment into the final equations with respect to the unknown parameters (flexural rigidity or external load). The final step solves the resulting system of equations. We apply this method to some inverse beam problems and provide an accurate estimation. Several numerical examples are performed and show that present method gives excellent results for identifying bending stiffness and distributed load of beam.

Keywords

Inverse Problem, Beam, Finite-Element Method, Structural Identification

1. Introduction

In recent years, analyses of inverse problems have been actively promoted in a variety of science and engineering fields. The usefulness of such analysis is that physical quantities or phenomena that are difficult to measure or observe directly can be determined from inverse problems via their outputs.

Inverse problems are diverse but typically include 1) obtaining the shape of a domain being analyzed, 2) obtaining boundary conditions for an entire boundary or the initial conditions, 3) obtaining loads applied to a domain, and 4) assuming the material characteristics of a field [1] . For analysis of inverse problems, the relationship between inputs and outputs is used in a manner similar to that used for direct problems. In the analysis of inverse problems, discretization methods, such as the finite-element method (FEM) and boundary-element method (BEM), have been used by engineers and researchers.

While identifying elastic structures, it is important to identify the location and shape of defects or cracks in the structure, extensive research has been carried out in this regard [2] [3] [4] . Yoshimura et al. [5] proposed a method based on a neural network and computational mechanics for identifying the solid-state properties of structures from their vibration characteristics using a plate. Arai et al. [6] demonstrated a BEM for identifying out-of-plane dynamic pressure distributions applied to a thin elastic plate. Both Ronasi et al. [7] and Rao & Ritti [8] demonstrated methods for estimating the load using a FEM. Jadamba et al. [9] employed output least-squares functional to identify the Lamé parameters in the linear elasticity.

In the present paper, a methodology for inverse problems in Euler-Bernoulli beams using static deflections measurements is proposed. In the first, we detail a methodology for identifying the distribution of flexural rigidity in a beam by using static deflection measurements. In the second, we develop a method for the identification of external loads acting on a beam. We give examples that demonstrate this scheme will yield accurate solutions.

2. Formulation of the Inverse Problem

2.1. Formulation of Mixed Beam Element

Consider the Euler-Bernoulli equations in the following form:

$EI\left(x\right)\frac{{\text{d}}^{2}w\left(x\right)}{\text{d}{x}^2}=-M\left(x\right)$ (1)

and

$\frac{{\text{d}}^{2}M\left(x\right)}{\text{d}{x}^2}=-q\left(x\right)$ (2)

where w denotes the transverse deflection and M(x) represents the bending moment at point x. The coefficient EI(x), called the flexural rigidity, is the product of the modulus of elasticity E and the moment of inertia I of the beam’s cross-section. The function q(x) represents the transversely distributed load (see Figure 1). We will demonstrate the discretization process based on the mixed form of the beam bending equations, Equations (1) and (2). Here, we begin by assuming that w and M are approximated in the usual manner by appropriate shape functions and unknown parameters:

$w=N{w}_e,M=N{M}_e$ (3)

where N represents element shape functions and w_e and M_e are the nodal parameters to be determined. Higher order functions can be used for shape functions, but the same linear shape functions were used for w_e and M_e to simplify the calculation. For the Galerkin solution, the original shape functions are simply used as weighting functions. If we consider the Galerkin form for an element of length $\mathcal{l}$ , we obtain:

Figure 1. The beam subjected to a normal load.

$-\frac{1}{EI}{\displaystyle {\int }_0^{\mathcal{l}}{N}^{\text{T}}N\text{d}x{M}_e}+{\displaystyle {\int }_0^{\mathcal{l}}\frac{\text{d}{N}^{\text{T}}}{\text{d}x}\frac{\text{d}N}{\text{d}x}dx{w}_e}={\left[N^{\text{T}}\theta \right]}_0^{\mathcal{l}}$ (4)

and

${\displaystyle {\int }_0^{\mathcal{l}}\frac{\text{d}{N}^{\text{T}}}{\text{d}x}\frac{\text{d}N}{\text{d}x}\text{d}x{M}_e}={\displaystyle {\int }_0^{\mathcal{l}}{N}^{\text{T}}q\text{d}x}+{\left[N^{\text{T}}V\right]}_0^{\mathcal{l}}$ (5)

where θ and V denote the slope and shear force, respectively. Here N is the linear shape function, which is:

$N_1=1-\xi ,N_2=\xi ,\xi =x/\mathcal{l}$ (6)

With the above approximations, Equations (4) and (5) yields the following

$\left[\begin{array}{cc}{k}_{11}& k_{12}\\ k_{21}& 0\end{array}\right]\left\{\begin{array}{c}{M}_e\\ w_e\end{array}\right\}=\left\{\begin{array}{c}\theta \\ F\end{array}\right\}$ (7)

where:

$k_{11}=-\frac{1}{EI}{\displaystyle {\int }_0^{\mathcal{l}}{N}^{\text{T}}N\text{d}x},k_{12}=k_{21}^{\text{T}}={\displaystyle {\int }_0^{\mathcal{l}}\frac{\text{d}{N}^{\text{T}}}{\text{d}x}\frac{\text{d}N}{\text{d}x}\text{d}x}$ (8)

$\left\{\begin{array}{c}{M}_e\\ w_e\end{array}\right\}=\left\{\begin{array}{c}\begin{array}{c}{M}_1\\ M_2\end{array}\\ \begin{array}{c}{w}_1\\ w_2\end{array}\end{array}\right\},\left\{\begin{array}{c}\theta \\ F\end{array}\right\}=\left\{\begin{array}{c}\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\\ \begin{array}{c}{V}_1+Q_1\\ V_2+Q_2\end{array}\end{array}\right\}$

2.2. Identification of the Flexural Rigidity of a Beam

We consider the inverse analysis of the flexural rigidity of a beam as applied to a finite-element solution of Equation (7). Suppose that the force acting on the beam domain, the boundary conditions, and the deflection at each point on the beam are given. Under these conditions, the bending moment at each point can be easily obtained from the second equation of Equation (7). The first equation of Equation (7) is written as follows:

$-\frac{\mathcal{l}}{6EI}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]\left\{\begin{array}{c}{M}_1\\ M_2\end{array}\right\}+\frac{1}{\mathcal{l}}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\left\{\begin{array}{c}{w}_1\\ w_2\end{array}\right\}=\left\{\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right\}$ (9)

By transforming Equation (9) into an equation for unknown quantity 1/EI_e of flexural rigidity, the following equation is obtained.

${\left[A\right]}_e{\left\{X\right\}}_e={\left\{B\right\}}_e$ (10)

with:

${\left[A\right]}_e=\left[\begin{array}{c}\frac{\mathcal{l}}{3}{M}_1+\frac{\mathcal{l}}{6}{M}_2\\ \frac{\mathcal{l}}{6}{M}_1+\frac{\mathcal{l}}{3}{M}_2\end{array}\right],{\left\{X\right\}}_e=\left\{\frac{1}{E{I}_e}\right\},{\left\{B\right\}}_e=-\left\{\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right\}+\left[\begin{array}{cc}\frac{1}{\mathcal{l}}& -\frac{1}{\mathcal{l}}\\ -\frac{1}{\mathcal{l}}& \frac{1}{\mathcal{l}}\end{array}\right]\left\{\begin{array}{c}{w}_1\\ w_2\end{array}\right\}$ (11)

When a beam is divided into m elements, the number of flexural rigidity values to be considered as unknown quantities will become m, i.e., the same as the number of elements. The next step is the assembly of the final equations of the type given by Equation (12). This is accomplished, according to the rule of Equation (10), by simple addition of all the numbers in the appropriate space of the global matrix:

$\left[A\right]\left\{X\right\}=\left\{B\right\}$ (12)

in which [A] is an $n\times m$ matrix, for which m is the total number of unknowns and n is the number of equations. For this problem, singular value decomposition is applied to [A], which yields:

$\left[A\right]=\left[U\right]\left[\Lambda \right]{\left[D\right]}^{\text{T}}$ (13)

where $\left[U\right],\left[\Lambda \right]$ and $\left[D\right]$ each denote the presence of $n\times n,n\times m$ and $m\times m$ matrices, respectively. Matrix $\left[\Lambda \right]$ is given as follows, while the matrix rank $\left[A\right]$ is expressed by p:

$\left[\Lambda \right]=\left[\begin{array}{ccc}\begin{array}{cc}{\lambda }_1& \\ & {\lambda }_2\end{array}& \begin{array}{cc}& \\ & \end{array}& \begin{array}{cc}& 0\\ & \end{array}\\ \begin{array}{cc}& \\ & \end{array}& \begin{array}{cc}\ddots & \\ & {\lambda }_p\end{array}& \begin{array}{cc}& \\ & \end{array}\\ \begin{array}{cc}& \\ 0& \end{array}& \begin{array}{cc}& \\ & \end{array}& \begin{array}{cc}0& \\ & \ddots \end{array}\end{array}\right],{\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_p>0,p\le \min \left(m,n\right)$ (14)

where ${\lambda }_i$ is the singular value of $\left[\Lambda \right]$ . This immediately yields:

$\left\{X\right\}=\left[D\right]{\left[\Lambda \right]}^{-1}{\left[U\right]}^{\text{T}}\left\{B\right\}$ (15)

2.3. Identification of Loads Acting on a Beam

If we consider the system of equations in Equation (7) to be a mixed system, in which M and w, the solution can proceed by eliminating M from the first equation and substituting it into the second to obtain:

$\left\{F\right\}=-\left[K_{21}\right]{\left[K_{11}\right]}^{-1}\left[K_{12}\right]\left\{w\right\}+\left[K_{21}\right]{\left[K_{11}\right]}^{-1}\left\{\theta \right\}$ (16)

As the initial slope $\theta \to 0$ , the equation above changes to:

$\left\{F\right\}=-\left[K_{21}\right]{\left[K_{11}\right]}^{-1}\left[K_{12}\right]\left\{w\right\}$ (17)

The vector {F} is termed the equivalent nodal force. The distributed forces q using interpolations of the form:

$q=N{q}_e$ (18)

where $q_e$ represents the distributed force parameters to be determined. For linearly distributed forces q over the element, the force vector by

$\frac{\mathcal{l}}{6}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]\left\{\begin{array}{c}{q}_1\\ q_2\end{array}\right\}=\left\{\begin{array}{c}{Q}_1\\ Q_2\end{array}\right\}\text{or}{\left[H\right]}_e{\left\{q\right\}}_e={\left\{Q\right\}}_e$ (19)

The next step is the assembly of the final equation of the type given by:

$\left[H\right]\left\{q\right\}=\left\{F\right\}$ (20)

Finally, the unknown distributed forces vector q can be obtained by solving the resulting equation system.

3. Numerical Examples

To verify the proposed method, we attempted to identify the flexural rigidity of a beam and the forces acting on it. In practice, the beam deflection is found through direct measurements. In this calculation, however, the beam deflection is calculated using beam theory prior to the inverse analysis. Since it is a model calculation, it is sufficient to use a unified unit, so the unit is assumed to be a dimensionless quantity.

3.1. Identification of the Flexural Rigidity of a Beam

We consider a simply supported beam subject to a concentrated load p applied at the center of the beam. Since any unit system may be used, provided it is unified throughout, dimensionless quantities are used here. Thus the length of the beam is l = 10.

$EI\left(x\right)=\{\begin{array}{l}10,\left(0\le x\le 3,7\le x\le 10\right)\\ 20,\left(3\le x\le 7\right)\end{array}$

In Figure 2, the solid line shows the distribution of flexural rigidity assumed in solving the problem, and the circles show the results of identification obtained using the proposed method. For the monitoring point, at which deflection is determined, the beam is divided into ten equal parts in Figure 2(a) and forty equal parts in Figure 2(b). In Figure 2, this method reproduces the behavior of a stepped beam almost exactly.

The next example is a beam with small defects in close proximity, as shown in Figure 3, where the flexural rigidity is expressed as:

$EI\left(x\right)=\{\begin{array}{l}10,\left(0\le x\le 3,3.2\le x\le 4.2,4.8\le x\le 10\right)\\ 9,\left(3\le x\le 3.2\right)\\ 8,\left(4.2\le x\le 4.8\right)\end{array}$

In this example, it is assumed that the reduction in flexural rigidity might be dependent on the level of damage. The beam is divided into fifty elements for identification. The results obtained are shown by plotting of the circles. It can be observed in Figure 3 that reductions in rigidity are identified with high accuracy. Figure 3 indicates that the flexural rigidity was evaluated precisely, except around the simple support, where the beam deflection was very small.

(a)(b)

Figure 2. Identification for the stepped beam. (a) 10 Elements (11 Nodes); (b) 40 Elements (41 Nodes).

Figure 3. Identification of damage in the beam.

3.2. Identification of Loads Acting on a Beam

We show the numerical results for the first example the simply supported beam subject to two concentrated loads, p = 1. In this example, we assume that the deflection is measured at 21 points. The circles at x = 0 and x = 10 in Figure 4 denote the vertical reactions at the supporting point. In Figure 4, the method reproduces the loads acting on a beam almost exactly.

In the next example, we show the numerical results for identification of a simply supported beam subjected to a uniform load acting over a portion of the beam, as indicated in Figure 5. Although a significant overshoot occurred around x = 3.5 and 8.5, on the whole, good identification results were obtained. In this figure, circles at x = 0 and x = 10 denote the vertical reactions at the supporting point. Although there was a part disturbed by overshoot, its amount was small and good one was obtained as an estimated value.

Figure 4. Identification of two concentrated loads.

Figure 5. Identification of a uniformly distributed load.

4. Conclusion

In this paper, we evaluated the efficacy of analysis techniques for identifying the flexural rigidity and distribution of loads acting on a beam. We applied an inverse analysis technique based on the finite-element method and singular value decomposition to the identification of the flexural rigidity of a simply supported beam subject to a concentrated load applied at the center of the beam. We further investigated the identification of forces acting on the beam in the form of a concentrated load or distributed load. The results of the identification of flexural rigidity showed that our method was effective for identification of distributed beam stiffness values. Although 10% overshoot occurred in the identification of distributed forces, which changed rapidly in a stepwise fashion, good results for the identification of loads were obtained on the whole.

Acknowledgements

Author would like to thank to anonymous referee for rigorous review, constructive comments and valuable suggestions.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Kubo S. ,et al. (1988)Inverse Problems Related to the Mechanics and Fracture of Solids and Structures. JSME International Journal Series 1 31, 157-166.
[2]	Burczynski, T. and Belich, W. (2001) The Identification of Cracks Using Boundary Elements and Evolutionary Algorithms. Engineering Analysis with Boundary Elements, 25, 313-322. https://doi.org/10.1016/S0955-7997(01)00027-3
[3]	Engelhardt, M., Stavroulakis, G.E. and Antes, H. (2006) Crack and Flaw Identification in Elastodynamics Using Kalman Filter Techniques. Computational Mechanics, 37, 249-265.
[4]	Teidj, S., Khamlichi, A. and Driouach, A. (2016) Identification of Beam Cracks by Solution of an Inverse Problem. Procedia Technology, 22, 86-93. https://doi.org/10.1016/j.protcy.2016.01.014
[5]	Yoshimura, S. and Yagawa, G. (1993) Inverse Analysis by Means of Neural Network and Computational Mechanics: Its Application to Structural Identification of Vibrating Plate. In: Kubo, S., Ed., Inverse Problems, Atlanta Technology Publications, Georgia, 184-193.
[6]	Arai, M., Nishida, T. and Adachi, T. (2000) Identification of Dynamic Pressure Distribution Applied to the Elastic Thin Plate. Inverse Problems in Engineering Mechanics II. In: Tanaka, M. and Dulikravich, G.S., Eds., International Symposium on Inverse Problems in Engineering Mechanics 2000, Elsevier, Nagano, 129-138.
[7]	Ronasi, H., Johansson, H. and Larsson, F. (2012) Load Identification for a Rolling Disc: Finite Element Discretization and Virtual Calibration. Computational Mechanics, 49, 137-147.
[8]	Rao, R.R. and Ritti, L.K. (2015) Estimation of Load Acting on a Crane Hook by Inverse Finite Element Method. International Journal of Engineering Research & Technology, 4, 126-129.
[9]	Jadamba, B., Khan, A.A. and Raciti, F. (2008) On the Inverse Problem of Identifying Lamé Coefficients in Linear Elasticity. International Journal Computers and Mathematics with Applications, 56, 431-443. https://doi.org/10.1016/j.camwa.2007.12.016