JOURNAL OF THEORETICAL AND APPLIED MECHANICS 44, 1, pp. 91-105, Warsaw 2006 IDENTIFICATION OF BEAM BOUNDARY CONDITIONS IN ILL-POSED PROBLEM Leszek Majkut Faculty of Mechanical Engineering and Robotics, AGH University of Science and Technology e-mail: majkut@agh.edu.pl The purpose of the study is demonstration of the possibility of identifi- cation of boundary conditions in an ill-posed problem. The problem is understood as determination of four constants necessary for a descrip- tion of functions of forced vibration amplitudes from three equations.To this end, the Singular Value Decomposition (SVD) algorithm is used. After determining the amplitudes of forced vibration, elasticity coeffi- cients of supports can be calculated from the equations describing the boundary conditions. Verification of the obtained mathematical model (elastically supported Bernoulli-Euler’s beam) was done by comparing natural frequencies ob- tained from analytical and numerical models, and analysing the cor- relation of forced vibration amplitude vectors for different excitation frequencies. Key words: Singular Value Decomposition, identification, inversemodel 1. Introduction Analysis of dynamic processes of real objects can be expensive, time- consumingand, in certain cases, impossible,whereas experiments canbe easily carried out on models which can be used to simulate dynamic responses. For this purpose, physical and mathematical models of an object should be built and followed by estimation ofmodel parameters and verification. This process is called the identification of mechanical systems (Giergiel and Uhl, 1990). The bibliography includes many definitions of the identification, such as one given byBellman (1965), which corresponds best to the analysed problem of identification of boundary conditions: 92 L.Majkut ”The identification is a process, in which based on certain given information on the system structure and certain information on inputs, outputsandoperationof theobject themissing information on the structure, inputs and outputs can be obtained.” Here, the analysed system is a beam, described byBernoulli-Euler’smodel with unknownboundary conditions,modelled by elastic supports.Themathe- matical model of the boundary conditions is described by equations coupling respectively the bendingmoment and angle of rotation of the cross-section as well as the lateral force and amplitude of vibrations in cross sections, inwhich the beam is supported. In order to determine support elasticities, an inversemodel of thebeamhas been created. The definition of the inverse model was given by Engel (Engel and Engel, 2005): ”Inverse modelling is done using the current results of several me- asurements of visible parameters in order to infer about actual values of the model parameters.” In thepaper, the inversemodel is understoodas determination of functions of amplitudes of forced vibration based on measurements of amplitudes in several points, and then as calculation of support elasticities from equations describing boundary conditions. The functionofvibrationamplitudes causedbyagiven force isdescribedby an equation containing four sought constants. The simplestmethod of finding the constants is tomeasure the amplitudes of vibrations caused by a force of a knownamplitude and frequency in fourpoints. Inorder tominimizemeasuring errors,measurements shouldbe taken in a larger number of points, afterwhich one of standard statistical methods, e.g. regression analysis, can be used. The biggest problem is encountered when the available number of measurement values is lower than the number of constants to be determined the so-called under-determined problem, see Lanczos (1961). This paper concerns such a problem of identification in an ill-posed problem. In such cases, four constants from three equations can be determined by using decomposition of the main matrix by the Singular Value Decomposition (SVD) algorithm (Golub and Van Loan, 1989). Decomposition by the SVD algorithm is also used in diagnostics of over- determined systems (with excess information about the system) (Cempel and Tobaszewski, 2005; Cholewa and Kiciński, 1997, 2001; Żółtowski and Cem- pel, 2004), or for determination of inverse models, e.g. in order to find power of acoustic wave sources (Engel, 2004; Engel et al., 2002; Moorhouse, 2003; Stryczniewicz, 2004). Identification of beam boundary conditions... 93 2. The inverse model of a beam In order to identify boundary conditions of a beam, an inverse model was used. In the paper, it is understood as determination of functions of forced vibrationamplitudesbasedonmeasurementsofvibrationamplitudes in several points, and then on calculation of support elasticities from relations describing the boundary conditions. Themodel of the beam on elastic supports is shown in Fig.1. Fig. 1. A beamwith general boundary conditions The differential equation of motion has form EI ∂4y(x,t) ∂x4 +ρA ∂2y(x,t) ∂t2 = q(x,t) (2.1) The function of load distribution is q(x,t)=Fδ(x,xf)e iωwt Equation (2.1) can be solved by separating the variables, i.e.: y(x,t) = =X(x)T(t). In the steady-state, the ”time” equation may be expressed as T(t)= eiωwt In this case, the differential equation of ”space” variable takes form X(4)(x)−λ4X(x)= F EI δ(x,xf) (2.2) where λ4 =ω2w ρA EI 94 L.Majkut and its solution is function (2.3) X(x)=P coshλx+Qsinhλx+Rcosλx+S sinλx+ (2.3) + F 2EIλ3 [ sinhλ(x−xf)− sinλ(x−xf) ] H(x,xf) where EI – bending stiffness A – cross-section area ρ – material density δ(x,xf) – Dirac delta function H(x,xf) – Heaviside step function at x=xf. Relation (2.3) describes the function (vector) of amplitudes of steady-state vibrations caused by a force with the amplitude F and frequency ωw applied to the beam at the point with the coordinate x = xf. The constants P , Q, R,S can be determined based on themeasurement of vibration amplitudes in several points of the beam. The procedure of determination of the constants is shortly described in the next section of the paper. After establishing the integration constants P ,Q,R, S, the sought values of support elasticity coefficients can be calculated from equations describing the boundary conditions. For the position x=0 EIX′′′(0)=−kT0X(0) −EIX ′′(0)=−kR0X ′(0) (2.4) hence, the lateral elasticity coefficient is kT0 =EIλ 3S−Q P +R (2.5) and the rotational elasticity coefficient kR0 =EIλ P −R Q+S (2.6) The boundary conditions at x= l are described by EIX′′′(l)= kTlX(l) EIX ′′(l)=−kRlX ′(l) (2.7) hence, the lateral elasticity coefficient is kTl =EIλ 3P sinhλl+Qcoshλl−Rsinλl−S cosλl−f1 P coshλl+Qsinhλl+Rcosλl+S sinλl−f2 (2.8) Identification of beam boundary conditions... 95 where f1 = F 2EIλ3 [coshλ(l−xf)+cosλ(l−xf)] f2 = F 2EIλ3 [sinhλ(l−xf)− sinλ(l−xf)] and the rotational one kRl =−EIλ P coshλl+Qsinhλl−Rcosλl−S sinλl−f3 P sinhλl+Qcoshλl−Rsinλl+S cosλl−f4 (2.9) where f3 = F 2EIλ2 [sinhλ(l−xf)+sinλ(l−xf)] f4 = F 2EIλ2 [coshλ(l−xf)− cosλ(l−xf)] The model has been developed on the assumption that the support ela- sticities are constant, i.e. they do not depend on the amplitude or vibration frequency. 3. The method of determination of the integration constants The simplest method of finding the constants P , Q, R, S is to measure vibration amplitudes caused by a force of a known amplitude and frequency at four points. In that case, the four constants can be determined from fo- ur equations (2.3) describing the vibration amplitudes at the four measuring points. In order to minimize measuring errors, themeasurements should be taken ata largernumberofpoints, andthenoneof statisticalmethods, e.g. regression analysis, can be used (equation (2.3) is a linear equation due to P , Q, R, S constants). In many cases of diagnostics or identification, it is not possible to obtain full information on the analysed object. In the case analysed here, the ”in- complete information” is to be understood that themeasurements of vibration amplitudes were taken only in three measuring points. By assuming that the measuring points are points with the coordinates x = a, x = b and x = c (values of the forced vibration amplitudes in these 96 L.Majkut points aremarked respectively as X(a),X(b),X(c)), weobtain three algebraic equations in form (2.3), which can be written in matrix form MC =B hence    coshλa sinhλa cosλa sinλa coshλb sinhλb cosλb sinλb coshλc sinhλc cosλc sinλc         P Q R S      =    b1 b2 b3    (3.1) where b1 =X(a)− F 2EIλ3 [sinhλ(a−xf)− sinλ(a−xf)]H(a,xf) b2 =X(b)− F 2EIλ3 [sinhλ(b−xf)− sinλ(b−xf)]H(b,xf) b3 =X(c)− F 2EIλ3 [sinhλ(c−xf)− sinλ(c−xf)]H(c,xf) Equation (2.10) is amatrix equationa solution towhichcanbeobtainedby inversion of themainmatrix, i.e. determination of thematrix M−1. In the case of a rectangular matrix, the inverse matrix can be calculated by decomposing the main matrix according to the Singular Value Decomposition algorithm (Golub andVan Loan, 1989) M=UWV> where U – square matrix of the 3×3 rank, having 3 orthogonal columns corresponding to 3 singular values given in thematrix W, such that U>U=1 W – pseudo-diagonal matrix of the rank 3×4, having non-negative singular values of the matrix M on its diagonal W=    w1,1 0 0 0 0 w2,2 0 0 0 0 w3,3 0    V – squarematrix of thee rank 4×4, having 4 orthogonal columns, such that V>V=1 Thus, the matrix M−1 can be calculated M −1 =(UWV>)−1 =(V>)−1W−1U−1 =VW−1U> Identification of beam boundary conditions... 97 and sought constant vector can be determined based on matrix obtained as result of decomposition C =VW−1U>B where W −1 =         1 w1,1 0 0 0 0 1 w2,2 0 0 0 0 1 w3,3 0         After computing the integration constants P , Q, R, S, the sought values of the support elasticities can be calculated from relations (2.5), (2.6), (2.8) and (2.9). 4. Comparative criteria for the models Verification of the model based on data obtained from an experiment is one of the main problems of the identification. The model obtained from the identification of the systemwith ”incomplete information” (under-determined problem) is an approximation of the real system. The verification stage of the identification involves checking of the obtained approximation for sufficiency of the objective for which the model was created. The identification objective adopted in the paper is the determination of amplitude vectors of steady state vibrations causedby a force of any frequency below the second eigenfrequency. The first basic criterion for comparison of the mechanical models is the comparison of their natural frequency (Uhl, 1997). Due to limitations of exci- tation frequency, the first two free vibration frequencieswill be comparedhere. The second criterion used in the analysis of the correlation between the models is the visual comparison of amplitude vectors of vibrations caused by forces of different excitation frequencies. Very commonly, a diagram in the Cartesian coordinate system is created, where amplitudes obtained from the experimental model are placed on the x axis, and the same amplitudes from the analyticalmodel are placed on the y axis. If themodels are coincident, the correspondingpoints shouldbe located on the straight line inclinedbyanangle of 45◦. Themathematical notation of this type of comparison (described and used inmodal analysis for finding the correlation between the eigenvectors) is 98 L.Majkut done using the MAC (Modal Assurance Criterion) (Uhl, 1997) coefficient MAC(x,y)= |x∗>Wgy| 2 (y∗>Wgy)(x∗>Wgx) where x∗,x and y∗,y are twovectors of forcedvibration amplitudes obtained from the analytical and experimentalmodel; Wg is aweightmatrix indicating which coordinates of the vector are themost importantduring the comparison. The analysis is performed on the assumption that the matrix Wg is a unit matrix. i.e. amplitudes in all points of the beam are identically important. 5. Numerical examples The subject of the analysis is an elastically supported beam shown in Fig.1 with the following material data: Young’s modulus E = 2.1 · 1011Pa; material density ρ = 7860kg/m3 and geometric data: cross-section b×h = 0.03×0.03m; beam length l=1.3m. The ”experimental” data required for the identification and verification come fromvibrationanalysisusing thefinite elementmethod.For thispurpose, amplitudes were computed for vibrations excited by a force of the amplitude F =100N applied to the beamat the pointwith a coordinate xf =0.9mand frequency ω=2πf (for different frequencies f). The elasticity constants of supports were determined based on themeasu- rement (obtained from FEM analysis) of the vibration amplitudes at three points, where two of them are located at the beam ends (a=0, c= l). After determining the elasticity coefficients and decomposing the main matrix ac- cording to the Singular Value Decomposition algorithm, the analytical model was verified against the criteria specified above. The free vibration frequencies were compared by finding deviation defined by the following formula δi = |ωie−ωia| ωie ·100% i=1,2 (5.1) where ωie denotes the ith natural frequency obtained from the experimental model (fromFEManalysishere), ωia – ithnatural frequencyobtained fromthe analytical model, where i=1,2 due to limitations of the excitation frequency (below the second eigenfrequency), i.e. only the first two natural frequencies will be compared. Identification of beam boundary conditions... 99 Figure 2 shows the above defineddeviation in thedetermination of the first and second natural frequency as a function of location of the third measure- ment point (other two points at the beam ends) for the excitation frequency ω=2π ·25=157.1rad/s. Fig. 2. Deviation in the determination of the first two natural frequencies as a function of location of the third sensor; ω=157.1rad/s Natural frequencies from FEM analysis are: ω1 = 262.1rad/s and ω2 = 1039.4rad/s. Figure 3 shows the above defined deviation in the determination of the first and second natural frequency as a function of location of the third me- asurement point for the excitation frequency ω=2π ·50=314.2rad/s (higher than the first eigenfrequency). Fig. 3. Deviation in the determination of the first two natural frequencies as a function of location of the third sensor; ω=314.2rad/s 100 L.Majkut According to the analysis of the results shown in Figs 2 and 3, smaller identification errors are made when the system is loaded by a frequency lower than the first frequency of free vibration. In such a case, location of the central sensor (other two are located on the beam ends) has no significant effect on uncertainty of the determination of the first two eigenfrequencies (deviations δ1 and δ2 below 5%). In the case of identificationwith excitation by a force of a frequencyhigher than the first natural frequency, it is essential (due to the identification error) to find the ”appropriate” position of the measuring element (deviation of the frequency determination varies from 0 to 33%). Afterwords, we will determine the correlation coefficients for the vectors of amplitudes of forced vibration obtained from the experiment Xe and the analytical model Xa MAC(Xa,Xe)= |X>aXe| 2 (X>aXa)(X > eXe) Table 1 summarizes the MAC coefficients computed by comparing the forced vibration vectors obtained for 7 different frequencies of the excitation force. Theboundary conditions necessary to calculate the vectors of vibrations from the analytical model are obtained from the identification measurements with the excitation frequency ω = 2π · 25 = 157.1rad/s and ω = 2π · 50 = =314.2rad/s. Table 1. Correlation coefficients MAC for vibration amplitude vectors obtained from the analytical and experimental models Identification Verification for ω=2πf with ω=2πf f =10 f =25 f =50 f =75 f =100 f =125 f =150 f =25Hz 0.9983 0.9982 0.9975 0.9989 0.9930 0.9965 0.9953 f =50Hz 0.9998 0.9997 0.9976 0.9992 0.9992 0.9907 0.9811 The identification and verification measurements were performed for the case in which the central measuring point was in the beam center (b= l/2). The verification results for the analyticalmodel indicate the correct identi- fication of the boundary conditions of the beam, at least in the assumed band of thr excitation frequency (below the second eigenfrequency). Another example of identification in the ill-posed problem can be a system shown in Fig.4 characterised by the following parameters: Young’s modulus E = 2.1 ·1011Pa; material density ρ=7860kg/m3; cross-section dimensions b×h=0.03×0.03m; beam length l=1.3m. Identification of beam boundary conditions... 101 Fig. 4. A cantilever beamwith an elastic support Thepurpose of the identificationmeasurements is the determination of the lateral elasticity coefficient kTl bymeasuring forced vibration amplitudes. As thr identification quantities, two boundary conditions for the beam left end X(0)= 0 X′(0)= 0 and forced vibration amplitudes in one point of the beam have been adopted. Using the above mentioned procedure (SVD), we can determine the inte- gration constants P , Q, R, S, which describe the amplitude of forced vibra- tions (2.3). Having determined the integration constants, the sought value of the elastic coefficient is calculated from the relation kTl =EIλ 3P sinhλl+Qcoshλl−Rsinλl−S cosλl−f1 P coshλl+Qsinhλl+Rcosλl+S sinλl−f2 (5.2) where f1 = F 2EIλ3 [coshλ(l−xf)+cosλ(l−xf)] f2 = F 2EIλ3 [sinhλ(l−xf)− sinλ(l−xf)] Figure 5 shows the above defined deviation, (5.1), in the determination of the first and second natural vibration frequency as a function of location of the measurement point (x = a) for the excitation frequency ω = 2π · 20 = =125.7rad/s (below the first eigenfrequency). The natural frequencies found by FEManalysis are: ω1 =214.0rad/s and ω2 =622.9rad/s. Figure 6 shows the deviation in the determination of the first and second natural frequency as a function of location of the measurement point (x= a) for excitation frequency ω=2π ·60=377.0rad/s (above the first eigenfrequ- ency). 102 L.Majkut Fig. 5. Deviation in the determination of the first two natural frequencies as a function of the sensor location; ω=125.7rad/s Fig. 6. Deviation in the determination of the first two natural frequencies as a function of the sensor location; ω=377.0rad/s Similarly as in the previous example, smaller identification errors aremade when the system is loadedbya frequency lower than thefirst eigenfrequency of the system.However, it is essential to find such a position of the sensor, which will ensure the minimal identification error (deviation varies from 0 to 20%). Table 2 summarizes the MAC coefficients computed by comparison of the forced vibration vectors obtained for 7 different frequencies of the excitation force. The boundary conditions necessary to calculate the vibration vectors from the analytical model are obtained from the identification measurements with the excitation frequency ω = 2π · 20 = 125.7rad/s and ω = 2π · 60 = =377.0rad/s. Identification of beam boundary conditions... 103 Table 2. Correlation coefficients MAC for vibration amplitude vectors obtained from analytical and experimental models Identification Verification for ω=2πf with ω=2πf f =10 f =20 f =40 f =50 f =60 f =80 f =90 f =20Hz 1.0 1.0 1.0 1.0 1.0 1.0 0.9999 f =60Hz 0.956 0.9538 0.9441 0.935 0.9213 0.8731 0.8713 Verification results for the analytical model indicate correct identification of the boundaryconditions of thebeam(it is assumed that good coincidence of vectors is for MAC > 0.8 (Uhl, 1997). In all calculations, it has been adopted that the measuring point is located in the beam center. 6. Conclusions The purpose of the study was to demonstrate the possibility of identifi- cation of beam boundary conditions in an ill-posed problem. In the analysed cases, the problem was to be understood as determination of four constants necessary fordescriptionof functions of forcedvibration amplitudes fromthree equations. In the first case, these equations described amplitudes of vibrations inpointswheremeasuring elements (sensors)were located, or in the second ca- se, two equationswere available for the description of the boundary conditions and vibration amplitude at the measuring point. Verification of the obtained mathematical model (elastically supported Bernoulli-Euler’s beam) was done by comparing natural frequencies obtained from the analytical model and the numerical experiment, and analysing the correlation of forced vibration amplitude vectors for different excitation fre- quencies. In both analysed cases, the deviation in the determination of the first and second free vibration frequencywas computed as a function of location of one measuring element. Identification of boundary conditions was done for two frequencies of the excitation force: below and above the first eigenfrequency. According to the results of analysis presented in Figs 2, 3, 5 and 6, smaller identification errors are made when the system is loaded by a frequency lower than the first eigenfrequency. However, it is essential to find such a position of the measuring element, which will ensure the minimal identification error. Another criterion used in the analysis for studying the correlation between the models (described and used in modal analysis for finding the correlation 104 L.Majkut between eigenvectors) is the MAC (Modal Assurance Criterion) coefficient. MAC coefficients for vectors of forced vibration obtained from the experiment and the analytical model have been summarized in Tables 1 and 2. All MAC coefficients are greater than 0.8, above which, good coincidence of vectors is assumed. The MAC coefficient defines only the similarity of forms between vectors – its value is not affected by vibration amplitudes, i.e. multiple vectors (e.g. amplitudes of vibrations obtained from the experiment can be n-times greater than those obtained from the analytical model in each cross-section of the beam) have the same shape (i.e. MAC =1). Therefore, it is onlynecessary to analyse thevibrationamplitudeofany sin- gle cross-section of the beam. The amplification factor (for a given excitation frequency) depends on the location of resonant frequencies. So selection of the position of measuring sensors which ensures the minimal error of determina- tion of the natural frequency, also ensures theminimal error of determination of the amplitudes of forced vibrations. References 1. BellmanR., 1965,Mathematical aspects of the theory of systems,Symposium on System Theory, NewYork 2. Cempel C., Tabaszewski M., 2005, Skalowanie obserwacji w wielowymia- rowej diagnostycemaszyn,XXXII Ogólnopolskie Sympozjum Diagnostyka Ma- szyn, WęgierskaGórka 3. Cholewa W., Kiciński J., 1997, Diagnostyka techniczna. Odwrotne modele diagnostyczne, Wyd. Pol. Ślaskiej, Gliwice 4. CholewaW.,Kiciński J., 2001,Diagnostyka techniczna. Metody odwracania modeli obiektów, Wyd. Pol. Ślaskiej, Gliwice 5. Engel Z., 2000, Zasada wzajemności, Kraków,Wyd. AGH 6. Engel Z., 2004, Współczesne metody badania procesów wibroakustycznych, Kwartalnik AGH Mechanika, 23, 23-31 7. Engel Z., Engel J., 2005, Metody inwersji i ich zastosowanie w mechanice, XXXII Ogólnopolskie Sympozjum Diagnostyka Maszyn, WęgierskaGórka 8. Engel Z., Pleban D., Stryczniewicz L., 2002, Investigation of emission sound pressure using inversivemethod,Archives of Acoustics, 27, 277-290 9. Giergiel J.,UhlT., 1990, Identyfikacja układówmechanicznych, PWN,War- szawa Identification of beam boundary conditions... 105 10. GolubG.H.,VanLoanC.F., 1989,Matrix computations,TheJonasHopkins University Press, Baltimor and London 11. Lanczos C., 1961, Linear Differential Operators, Van Nostad, London 12. MoorhouseA.T., 2003,Compensation for discarded singular values in vibro- acoustic inversemethods, Journal of Sound and Vibrations, 267, 245-252 13. Stryczniewicz L., 2004, Identification of zones of increased vibroacoustic emission bymeans of the inversemethod, Archives of Acoustics, 29, 563-575 14. UhlT., 1997,Komputerowo wspomagana identyfikacja modeli konstrukcji me- chanicznych, WNT,Warszawa 15. Żółtowski B., Cempel C., 2004, Inżynieria diagnostyki maszyn, Bibl. Pro- blemówEksploatacji, PTDT, ITEks.-PIB,Warszawa, Bydgoszcz, Radom Identyfikacja warunków brzegowych belki w układach o niepełnej informacji Streszczenie Praca dotyczy identyfikacji warunków brzegowych belki w przypadkach, w któ- rych nie ma możliwości uzyskania pełnej informacji o układzie. Niepełna informacja wynika, w rozważanym w pracy przypadku, z problemu wyznaczenia czterech sta- łych całkowania,niezbędnychdo opisania funkcji amplituddrgańwymuszonychbelki, z trzech równań. Do tego celu wykorzystano algorytm rozkładu macierzy względem wartości szczegónych (Singular Value Decomposition). Powyznaczeniu stałychcałkowania i funkcji amplituddrgańwymuszonych,uogól- nionewspółczynniki sprężystości podparciawyznaczono z równańopisującychwarun- ki brzegowe. Weryfikacji tak uzyskanegomodelumatematycznego dokonanopoprzez porówna- nie częstości drgańwłasnych i wyznaczeniewspólczynnikówkorelacji wektorówdrgań wymuszonych. Manuscript received July 28, 2005; accepted for print October 19, 2005