This work has been supported by the grant 23-06220S of the Czech Science Foundation within
institutional support RVO:61388998

Cimrman, R.

Kolman, R.

Kotek, V.

Kylar, J.

Musil, L.

DSpace at University of West Bohemia

Dynamics of a cantilever beam with piezoelectric sensor:Parameter identificationR. Cimrmana,b, R. Kolmana,c, L. Musila, V. Koteka, J. KylaraaInstitute of Thermomechanics, Czech Academy of Sciences, Dolejškova 1402/5, 18200, Praha 8, Czech RepublicbNew Technologies Research Centre, University of West Bohemia, Univerzitnı́ 8, 301 14 Plzeň, Czech RepubliccCollege of Polytechnics Jihlava, Tolstého 16, 586 01 Jihlava, Czech RepublicThe piezoelectric materials are electroactive materials often applied for real-time sensing orstructural health monitoring. Mathematical models of such structures contain several materialparameters that need to be identified from experiments. The (proportional) damping coefficientsare an example of such parameters, difficult to obtain in another way.Our aim was to develop a computer model of a simple experiment performed in our labora-tory [6] involving dynamics of a cantilever beam with an attached piezoelectric sensor excitedby a suddenly removed weight, see Figs. 1–2. The sensor was connected to an oscilloscopethat measured the voltage on the top side of the sensor, the bottom side was grounded. Wehad been interested in correctly simulating the experiment for various materials of the beamand/or in identifying parameters that were unknown or uncertain in advance. Preliminary simu-lations revealed that the oscilloscope had a finite resistance, motivating us to augment the modelpresented in [4]. The new model is briefly summarized below.Let us denote by ΩE the elastic part and by ΩP the piezoelectric part of a body with Ω ⊂ R3,see Fig. 1. The weak formulation of the model is as follows. Let V u0 (Ω) = {u ∈ [H1(Ω)]3,u =0 on Γu}, V ϕ0 = {ϕ ∈ H1(Ω), ϕ = 0 on Γp0}. We seek u(t), ϕ(t), ϕ̄(t) such that∫Ωρv · ü+∫Ωε(v)TCε(u) + α∫Ωρv · u̇+ β∫Ωε(v)TCε(u̇)−∫ΩPε(v)TeT∇ϕ−∫Ωv · b = 0 ∀v ∈ V u0 (Ω) , (1)∫ΩP(∇ψ)Teε(u) +∫ΩP(∇ψ)Tκ∇ϕ−∫ΓpQ(κ∇ϕ) · nψ+∫ΓpQ(κ∇ψ) · n(ϕ− ϕ̄) = 0 ∀ψ ∈ V ϕ0 (Ω) , (2)∫ΓpQ(κ∇ϕ̇) · n+ ϕ̄/R = 0 , (3)u = 0 on Γu × [0, T ] , (4)ϕ = 0 on Γp0 × [0, T ] , (5)u(0) = u0, u̇(0) = 0, ϕ(0) = ϕ0 in Ω , (6)where u is the mechanical displacement vector inducing the Cauchy strain ε, ϕ is the electricpotential, ϕ̄ is the unknown potential on the top side of the sensor ΓpQ, ρ is the density. Cis the matrix of elastic properties (under constant electric field intensity in ΩP ), α, β are the21W↓ LΩE ΩPΓuFig. 1. The computational mesh.A weight is attached to the bottom sideat the pointW at time t = 0 s. Dynamicquantities are recorded in the point L0.250.000.25A3.655 3.660 3.665 3.670 3.6750.000.02BA: velocity, B: displacement2.040 2.045 2.050 2.055 2.0604020020 AA: voltageFig. 2. The experimental data: (left) laser vibrometer, (right)oscilloscopeproportional damping coefficients, e the piezoelectric modulus, κ the permittivity under con-stant deformation and n is the unit outward normal. The structure is loaded by the self-weightvolume forces b. The initial conditions u0, ϕ0 correspond to a static loading of the body by alocalized surface force of a weight (280 g) attached, using a nylon string, at point W in Fig. 1.The external circuit is modeled by (3) with R being the oscilloscope resistance. The last twoterms in (2) correspond to the weak enforcement of the Dirichlet boundary condition ϕ = ϕ̄(t)on ΓpQ using the non-symmetric Nitsches method [8] without the penalty term [2].The finite element discretization of (1)–(3) usesu(ξ) = Nu(ξ)u , ϕ(ξ) = Nϕ(ξ)p , ϕ̄ = 1p̄ , (7)where Nu(ξ), Nϕ(ξ) are the displacement and the potential shape functions, respectively. Theimplicit second order Newmark method [7] is used to discretize the dynamical part (1), whilethe central difference scheme is used for (3). Because we want to perform a sensitivity analysisof the time-dependent system on several problem parameters, we present the discrete equationsin the extended form with the unknowns vectors u, u̇, ü, p, ṗ, p̄ of the time step n, the n − 1time step quantities indicated by the superscript 0Ü : Mü + Cu̇ + Ku−BTp− b = 0 , (8)P : Bu + (D− F + FT )p− FT1p̄ = 0 , (9)P̄ : 1TFṗ + 12R(p̄+ p̄0) = 0 , (10)U̇ : u̇− u̇0 − (1− γN)∆tü0 − γN∆tü = 0 , (11)U : u− u0 −∆tu̇0 − (12− βN)∆t2ü0 − βN∆t2ü = 0 , (12)Ṗ : ṗ + 1∆tp0 − 1∆tp = 0 , (13)where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, B is the piezoe-lastic coupling matrix, D is the electrostatic potential matrix and b is the volume forces vector.The surface flux matrix F =∫ΓpQNTϕnκN′ϕ is used to impose weakly the Dirichlet boundarycondition ϕ = ϕ̄ on ΓpQ, 1, the matrix of ones, is used to sum the rows of F, performing thusthe integration of (3). The Newmark scheme is expressed by (11), (12), βN , γN are its param-eters, the central difference scheme by (13) and the average in the second term of (10). Whencomputing the direct problem, u̇, u, ṗ expressed from (11)–(13) are substituted into (8)–(10),resulting in three primary unknowns ü, p, p̄ to compute in each time step by solving a linearsystem.22The parameter identification is done by fitting a time series of measured data – in our casethe velocity component vz(t) measured by a laser vibrometer in the point L – to the corre-sponding simulated quantity u̇z(L, t) by a nonlinear least-squares solver (least squares()function from SciPy [10]). This solver can compute the Jacobian matrix of the objective func-tion F numerically or use a user-supplied function returning partial sensitivities. For our linearproblem, we do not expect to achieve a faster elapsed time with a semi-analytical Jacobianmatrix (see below) than with the default 2-point finite difference scheme, because the numberof linear systems to solve is the same and computing the Jacobian analytically involves someadditional calculations. Nevertheless it allows us to test our automatically differentiated termsbased on JAX [1] and provides a proof-of-concept for future nonlinear extensions of the model.The state problem (8)–(13) in an abstract form (see also [9]) can be written asΦ(α,y,y0) = 0 , y ≡ [u, u̇, ü,p, ṗ, p̄]T , (14)where α are the parameters to be identified. Our aim is to calculate the partial sensitivities∂u̇z(L,t)∂α. By differentiating (14) w.r.t. α a recurrent relation∂Φ∂y∂y∂α= −(∂Φ∂α+∂Φ∂y0∂y0∂α) (15)is obtained allowing us to compute ∂y∂αin all time steps, initialized by using the initial conditions.Then ∂u̇z(L,t)∂αare simply a component in ∂y∂α.To briefly demonstrate the parameter identification, below we identify the parameters β (thestiffness proportional damping parameter), E (the Young’s modulus of the beam) and u0 (ascalar multiplier of the initial condition u0), using preliminary data for t = [0, 0.02] s with thefollowing (initial) material parameters:• Steel elastic beam: ρ = 7800 kg/m3, E = 210 GPa, ν = 0.3.• Piezoelectric disc: ρ = 7800 kg/m3, the vacuum permittivity ε0 = 8.8541878128 ·10−12 F/m andin Voigt notation: CP =127.2050 80.2122 84.6702 0.0000 0.0000 0.000080.2122 127.2050 84.6702 0.0000 0.0000 0.000084.6702 84.6702 117.4360 0.0000 0.0000 0.00000.0000 0.0000 0.0000 22.9885 0.0000 0.00000.0000 0.0000 0.0000 0.0000 22.9885 0.00000.0000 0.0000 0.0000 0.0000 0.0000 23.4742GPa,e =0.00000 0.00000 0.0000 0.0000 17.0345 0.00.00000 0.00000 0.0000 17.0345 0.0000 0.0−6.62281 −6.62281 23.2403 0.0000 0.0000 0.0 C/m2, κ = ε01704.4 0.0 0.00.0 1704.4 0.00.0 0.0 1433.6 F/m .All results were obtained using the Open Source finite element software SfePy [5, 3]. Theconvergence of the identification procedure is shown in Fig. 3 (left), where the evolution of theparameters and the l2 norm of the objective function are depicted. The initial and identified timehistories of u̇z(L, t) are shown in Fig. 3 (right), together with the experimental data. The initialparameters correspond to reference values (a steel beam) and the attached weight loading. Theresults indicate some discrepancies in the experiment: the initial deflection of the beam shouldhave been higher, see u0, to explain the amplitude of the data, and the elastic beam Young’smodulus E lower to match the principal frequency of the data. Note that the decreased Young’smodulus E does not explain the initial amplitude by itself. Based on that, new experimentsare being undertaken, where e.g. the influence on the applied initial load of cutting the nylonstring of the weight is decreased, and the material parameters of the beam are determined exper-imentally instead of using reference material data from a data sheet. Finally, the proportionaldamping coefficient was determined to be about 1.92× 10−5 s.231.0 · 1001.2 · 1001.4 · 1001.6 · 1001.9 · 100u0 [1]1.5 · 10111.7 · 10111.8 · 10112.0 · 10112.2 · 1011E [Pa]0 10 20iteration1.0 · 10−51.7 · 10−52.8 · 10−54.7 · 10−57.8 · 10−5β [m]0 10 20iteration1.8 · 10−14.0 · 10−19.2 · 10−12.1 · 1004.8 · 100|F |0.000 0.005 0.010 0.015 0.020t [s]−0.4−0.20.00.20.4v z[ms−1]identifiedinitialexperimentFig. 3. Parameter identification: (left) the evolution of parameters and objective function during theobjective function calls; (right) experiment, original and identified u̇z(L, t)AcknowledgementThis work has been supported by the grant 23-06220S of the Czech Science Foundation withininstitutional support RVO:61388998.References[1] Bradbury, J., Frostig, R., et al., JAX: Composable transformations of Python+NumPy programs,https://github.com/google/jax, 2023. Ver. 0.4.14.[2] Burman, E., A penalty-free nonsymmetric nitsche-type method for the weak imposition of bound-ary conditions, SIAM Journal on Numerical Analysis 50 (4) (2012) 1959-1981.[3] Cimrman, R., Fast evaluation of finite element weak forms using python tensor contraction pack-ages, Advances in Engineering Software 159 (2021) No. 103033.[4] Cimrman, R., Kolman, R., Musil, L., Kotek, V., Kylar, J., Dynamics of a cantilever beam withpiezoelectric sensor: Finite element modeling, Proceedings of the 29th International ConferenceEngineering Mechanics 2023, Milovy, May 9-11, 2023, (editors) V. Radolf, I. Zolotarev, Instituteof Thermomechanics of the Czech Academy of Sciences, 2023, pp. 51-54.[5] Cimrman, R., Lukeš, V., Rohan, E., Multiscale finite element calculations in Python using SfePy,Advances in Computational Mathematics 45 (2019) 1897-1921.[6] Kolman, R., Kylar, J., Kotek, V., Cimrman, R., Musil, L., Dynamics of a cantilever beam withpiezoelectric sensor: Experimental study, Proceedings of the 29th International Conference Engi-neering Mechanics 2023, Milovy, May 9-11, 2023, (editors) V. Radolf, I. Zolotarev, Institute ofThermomechanics of the Czech Academy of Sciences, 2023, pp. 119–122.[7] Newmark, N. M., A method of computation for structural dynamics, Journal of the EngineeringMechanics Division 85 (3) (1959) 67-94.[8] Nitsche, J., About a variation principle for solving the Dirichlet problem when using subspacesthat are not subject to boundary conditions, Abhandlungen aus dem Mathematischen Seminar derUniversität Hamburg 36 (1) (1971) 9-15. (in German)[9] Rohan, E., Cimrman, R., Sensitivity analysis and material identification for activated smooth mus-cle, Computer Assisted Mechanics and Engineering Sciences Vol. 9, No. 4 (2002) 519–541.[10] Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski,E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, I., Feng,Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero,E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., SciPy1.0: fundamental algorithms for scientific computing in Python, Nature Methods 17 (33) (2020)261-272.24

Dynamics of a cantilever beam with piezoelectric sensor: Parameter identification

https://dspace5.zcu.cz/bitstream/11025/54834/2/CM2023_Conference_Proceedings-30-33.pdf

Dynamics of a cantilever beam with piezoelectric sensor: Parameter identification

Abstract

Similar works

Full text

Available Versions

DSpace at University of West Bohemia