# Updating Stiffness and Hysteretic Damping Matrices Using Measured Modal Data

Authors: Jiashang Jiang; Yongxin Yuan
Journal: Shock and Vibration
 (2018)

---

## Abstract

A new direct method for the finite element (FE) matrix updating problem in a hysteretic (or material) damping model based on measured incomplete vibration modal data is presented. With this method, the optimally approximated stiffness and hysteretic damping matrices can be easily constructed. The physical connectivity of the original model is preserved and the measured modal data are embedded in the updated model. The numerical results show that the proposed method works well.


---

## Body

## 1. Introduction

The dynamic behavior of a mechanical system withn degrees of freedom is modeled by the following set of second-order ordinary differential equations:(1)Maq¨t+Ka+iHaqt=0, where Ma,Ka∈Rn×n and Ha∈Rn×n are analytical mass, stiffness, and hysteretic damping matrices, respectively. q(t)∈Rn×1 is a displacement vector depending on time t. In general, Ma is symmetric and positive definite, and Ka and Ha are symmetric. Equation (1) is usually modeled by FE techniques and thus a FE model is a numerical model. Assuming a fundamental solution q(t)=xeαt, one obtains the structural eigenproblem(2)λMax+Ka+iHax=0, where λ=α2.The FE method is used in many applications in engineering practice such as structural response prediction, structural control, structural health monitoring, damage detection, and reliability and risk assessment [1–6]. However, the accuracy of the FE model may be adversely affected by the inaccuracies in the model often related to material properties, modeling of joints, boundary conditions, and damping and simplifications made. As a result, a significant discrepancy may exist between the modal properties calculated by the constructed FE model and those identified from the vibration measurements of the actual structure. Many investigations show that the differences between the numerical and experimental frequencies may exceed 10% [7, 8]. The problem of how to modify a numerical model from the dynamic measurements is known as model updating in structural dynamics. Generally speaking, FE model updating involves two major applications. The first aspect is the FE model tuning, in which the goal is to update a numerical model to characterize the behavior of a real structure; the second one is damage detection, where the discrepancies in the structural dynamic properties before and after damage are identified to help locate and quantify structural damage.FE model updating has been an active area of research for the last four decades and a lot of approaches have been established for updating structural dynamic models [5, 6]. The direct matrix updating methods were first introduced by Baruch and Bar-Itzhack [9], Baruch [10], and Berman and Nagy [11]. An optimal solution was obtained using Lagrange multipliers for minimizing the changes in the matrices subject to the orthogonality properties of the modes, the eigenvalue equation, and the symmetry of the updated matrices. The matrix mixing methods were developed by Caesar [12] and Link et al. [13]. This approach utilized experimental modal data and analytical ones to construct the inverses of the mass and stiffness matrices. An alternative approach suggested by Wei [14, 15] was to update the mass and stiffness matrices simultaneously using a measured eigenvector matrix as the reference. The control-based eigenstructure assignment techniques for FE model updating were proposed by Zimmerman and Widengren [16] and Inman and Minas [17]. These early methods are simple and computationally efficient but do not generally respect the structural connectivity in the initial FE model. To deal with this difficulty, Kabe [18] presented a stiffness matrix adjustment technique using the constrained minimization theory. The adjusted stiffness matrix can predict the measured modal data accurately and the connectivity of the original stiffness matrix is preserved. Kammer [19] proposed the projector matrix (PM) method which uses the projector matrix theory and the Moore-Penrose inverse, resulting in a more computationally efficient solution. Halevi and Bucher [20] also established a direct matrix updating method by including the connectivity constraints.It is observed that model updating problems for linear viscously damped elastic systems (or gyroscopic systems) have been considered by many authors [21–29]. However, problems for updating hysteretic damping models have not got enough attention in these years. It is known that the principal difference between viscous and hysteretic damping models is that, for the viscous system, the energy dissipation per cycle depends linearly upon the frequency of oscillation, while for the hysteretic case it is independent of frequency. Very few authors pay attention to hysteretic damping model updating problems because the free vibration response for a system with hysteretic damping is necessarily complex, whereas the modal data of viscously damped elastic systems are closed under complex conjugation. In [30], the authors provided an extended application of the constrained eigenstructure assignment method (CEAM), which was first introduced in [31], to finite element model updating. The existing formulation was modified to accommodate larger systems by developing a quadratic optimization procedure which is unconditionally stable.The present paper proposes an efficient direct updating method for the FE model with hysteretic (or material) damping which preserves the connectivity of the original coefficient matrices. Assume thatKa and Ha are real-valued symmetric (2r+1)-diagonal matrices; that is, Ka and Ha are band matrices and the nonzero elements are n(2r+1)-r(r+1). The problem of updating stiffness and hysteretic damping matrices simultaneously can be stated following the inverse eigenvalue problem and an associated optimal approximation problem, leading to Problems 1 and 2.Problem 1.
LetΓ=diag⁡{γ1,…,γp}∈Cp×p and Z=[z1,…,zp]∈Cn×p represent the diagonal matrix of p measured eigenvalues and the matrix of corresponding p measured eigenvectors, where p≪n. Find real-valued symmetric (2r+1)-diagonal matrices K and H such that(3)MaZΓ+K+iHZ=0.Problem 2.
Find(K^,H^)∈SKH such that(4)K^-Ka2+H^-Ha2=minK,H∈SKH⁡K-Ka2+H-Ha2, where SKH is the solution set of Problem 1.This paper is divided into five sections. In Section2, we give an explicit formula for the solution set SKH of Problem 1 using the Kronecker product and vec operator. In Section 3, we show that the solution of Problem 2 is unique and present the expression of the unique solution (K^,H^) of this problem. Section 4 describes and discusses the numerical results which indicate the accuracy and efficiency of the proposed algorithm. Some concluding remarks will be drawn in Section 5.In this paper, we use the following notations.Cm×n and  Rm×n denote the set of all m×n complex and real matrices,  and SRn×n denotes the set of all symmetric matrices in Rn×n. A⊤ and A+ denote the transpose and the Moore-Penrose inverse of a real matrix A, respectively. In represents the identity matrix of size n. In space Rm×n, we define an inner product as (A,B)=trace(B⊤A), for all A,B∈Rm×n; then, Rm×n is a Hilbert space and the matrix norm ‖·‖ induced by the inner product is the Frobenius norm. Given A=[aij]∈Rm×n,B=[bij]∈Rp×q, the Kronecker product of A and B is defined by A⊗B=[aijB]∈Rmp×nq. Also, vec(·) represents the vec operator; that is, vec⁡(A)=[a1⊤,a2⊤,…,an⊤]⊤, for the matrix A=[a1,a2,…,an]∈Rm×n, where aj,j=1,…,n, is the jth column vector of A. For convenience, the symbols EA and FA stand for the two orthogonal projectors: EA=Im-AA+ and  FA=In-A+A, where A∈Rm×n.


## 2. The Solution of Problem1

To solve Problem1, we need the following lemmas.Lemma 3 (see [32]).
Assume thatL∈Rm×q and b∈Rm. Then, the equation of Ly=b has a solution y∈Rq if and only if LL+b=b. In this case, the set of all solutions of the equation is y=L+b+FLz, where z∈Rq is an arbitrary vector.Lemma 4 (see [33]).
Suppose thatD∈Rm×n,H∈Rn×l, and J∈Rl×s. Then, (5)vec⁡DHJ=J⊤⊗Dvec⁡H.Lemma 5 (see [34]).
Let[A,B] be an arbitrary partitioned matrix, where A∈Rm×k,B∈Rm×l, and let C=EAB,D=A+BFC. Then, the generalized inverse of the matrix [A,B] is (6)A,B+=Ik+DD⊤-1A+-A+BC+C++D⊤Ik+DD⊤-1A+-A+BC+.If we denote the set of alln×n real symmetric (2r+1)-diagonal matrices by S0, then we can easily see that S0 is a linear subspace of SRn×n, and the dimension of S0 is N=1/2(2n-r)(r+1). Define(7)Yij=22eiej⊤+ejei⊤,i=1,…,n-1;j=i+1,…,ti,eiei⊤,i=j=1,…,n, where ti=min⁡{i+r,n} and ej,j=1,…,n, is the jth column vector of In. Clearly, {Yij} forms an orthonormal basis of the subspace S0; that is,(8)Yij,Ykl=0,i≠korj≠l,1,i=k,j=l. Now, if K,H∈SRn×n are (2r+1)-diagonal matrices, then K,H can be expressed as(9)K=∑i,jαijYij,H=∑i,jβijYij, where the real numbers αij,βij,i=1,…,n;j=i,…,ti,ti=min⁡{i+r,n}, need to be determined.Separating (3) into its real and imaginary parts yields(10)KZR-HZI=MaZIΓI-ZRΓR,KZI+HZR=-MaZRΓI+ZIΓR, where Z=ZR+iZI(ZR,ZI∈Rn×p),Γ=ΓR+iΓI(ΓR,ΓI∈Rp×p). Substituting (9) into (10), we have(11)∑i,jαijYijZR-∑i,jβijYijZI=MaZIΓI-ZRΓR,∑i,jαijYijZI+∑i,jβijYijZR=-MaZRΓI+ZIΓR. Let (12)α=α11,…,α1,r+1,…,αn-r,n-r,…,αn-r,n,…,αn-1,n-1,αn-1,n,αn,n⊤,β=β11,…,β1,r+1,…,βn-r,n-r,…,βn-r,n,…,βn-1,n-1,βn-1,n,βn,n⊤,G=vec⁡Y11,…,vec⁡Y1,r+1,…,vec⁡Yn-r,n-r,…,vec⁡Yn-r,n,…,vec⁡Yn-1,n-1,vec⁡Yn-1,n,vec⁡Yn,n∈Rn2×N,(13)A=ZR⊤⊗InG,B=ZI⊤⊗InG,(14)d=vec⁡MaZIΓI-ZRΓR,f=vec⁡-MaZRΓI+ZIΓR.By Lemma 4, (11) are equivalent to(15)Jξ=d,Lξ=f, where J=[A,-B],L=[B,A], and ξ=[α⊤,β⊤]⊤. Applying Lemma 3, the first equation of (15) with unknown vector ξ has a solution if and only if(16)JJ+d=d. In this case, the general solution of the equation is(17)ξ=J+d+FJu, where u∈R2N is an arbitrary vector. Substituting (17) into the second equation of (15) yields(18)LFJu=f-LJ+d. Using Lemma 3 again, the equation of (18) with unknown vector u has a solution if and only if(19)ETf-LJ+d=0, where T=LFJ. When condition (19) is satisfied, the general solution of the equation of (18) is given by(20)u=T+f-LJ+d+FTv, where v∈R2N is an arbitrary vector. Substituting (20) into (17) results in(21)ξ=J+d+FJT+f-LJ+d+FJFTv. By Lemma 5, we have(22)J+=PQ, where(23)C=-EAB,D=-A+BFC,P=IN+DD⊤-1A++A+BC+,Q=C++D⊤IN+DD⊤-1A++A+BC+. Using (22), we have(24)T=LFJ=A1,B1, where(25)A1=B-BPA-AQA,B1=BPB+A+AQB. Using Lemma 5 again, we can get(26)T+=UV, where(27)C1=EA1B1,D1=A1+B1FC1,U=IN+D1D1⊤-1A1+-A1+B1C1+,V=C1++D1⊤IN+D1D1⊤-1A1+-A1+B1C1+. By (22), (24), and (26), (21) can be equivalently written as(28)α=Pd+U-PAU+PBVf-BPd-AQd+IN-PAIN-UA1s-PBVA1s-IN-PAUB1t+PBIN-VB1t,(29)β=Qd+V-QAU+QBVf-BPd-AQd-QAIN-UA1s-IN+QBVA1s+QAUB1t+IN+QBIN-VB1t, where v=[s⊤,t⊤]⊤ and s,t∈RN are arbitrary vectors.To sum up the above discussion, we can get the following result.Theorem 6.
Assume thatΓ=diag⁡{γ1,…,γp}∈Cp×p and Z=[z1,…,zp]∈Cn×p are the measured eigenvalue and eigenvector matrices, where p≪n. Let Z=ZR+iZI(ZR,ZI∈Rn×p),Γ=ΓR+iΓI(ΓR,ΓI∈Rp×p), and G,A,B,d,f be given by (12), (13), and (14), respectively. Also, let N=1/2(2n-r)(r+1),J=[A,-B],L=[B,A], and T=LFJ=[A1,B1]. J+ and T+ are given by (22) and (26), respectively. If conditions (16) and (19) are satisfied, then the solution set SKH of Problem 1 is(30)SKH=K,H∈SRn×n×SRn×n∣K=Sα⊗In,H=Sβ⊗In, where(31)S=Y11,…,Y1,r+1,…,Yn-r,n-r,…,Yn-r,n,…,Yn-1,n-1,Yn-1,n,Yn,n∈Rn×nN,where C,D,P,Q,A1,B1,C1,D1,U,V,α, and β are, respectively, given by (23), (25), (27), (28), and (29), with s,t∈RN being arbitrary vectors.


## 3. The Solution of Problem2

When the setSKH is nonempty, it is easy to check from (30) that SKH is a closed convex subset of SRn×n×SRn×n. It follows from the best approximation theorem (see, e.g., [35]) that there exists a unique solution (K^,H^) in SKH such that (4) holds.Now, we set out to find the unique solution(K^,H^) in SKH. If Ka and Ha are real symmetric (2r+1)-diagonal matrices, then Ka,Ha can be expressed as the linear combinations of the orthonormal basis {Yij}; namely,(32)Ka=∑i,jγijYij,Ha=∑i,jδijYij, where γij,δij,i=1,…,n;j=i,…,ti,ti=min⁡{i+r,n}, are uniquely determined by the elements of Ka and Ha. Let(33)γ=γ11,…,γ1,r+1,…,γn-r,n-r,…,γn-r,n,…,γn-1,n-1,γn-1,n,γn,n⊤,(34)δ=δ11,…,δ1,r+1,…,δn-r,n-r,…,δn-r,n,…,δn-1,n-1,δn-1,n,δn,n⊤. Then, for any (K,H)∈SKH in (30), by (8), (9), (32), (33), and (34), we get (35)g=K-Ka2+H-Ha2=∑i,jαij-γijYij2+∑i,jβij-δijYij2=∑i,jαij-γijYij,∑i,jαij-γijYij+∑i,jβij-δijYij,∑i,jβij-δijYij=∑i,jαij-γij2+∑i,jβij-δij2=α-γ2+β-δ2=αβ-γδ2. Substituting (21) into the relation of g, we get (36)g=J+d+FJT+f-LJ+d+FJFTv-ψ2, where ψ=[γ⊤,δ⊤]⊤. Notice that (37)FTFJ=I2N-LFJ+LFJFJ=FJ-LFJ+LFJ=FJ-T+T, which implies that(38)FTFJ=FJFT. Therefore, (39)g=v⊤FTFJFTv-2v⊤FTFJψ-2ψ⊤J+d-2ψ⊤FJT+f-LJ+d+d⊤JJ⊤+d+f-LJ+d⊤T+⊤FJT+f-LJ+d+ψ⊤ψ,∂g∂v=2FTFJv-2FTFJψ. Clearly, the function g=‖K-Ka‖2+‖H-Ha‖2=min attains the smallest value at ∂g/∂v=0, which yields(40)FTFJv=FTFJψ. By substituting (40) into (28) and (29), we obtain(41)α^=Pd+U-PAU+PBVf-BPd-AQd+IN-PAIN-UA1γ-PBVA1γ-IN-PAUB1δ+PBIN-VB1δ,(42)β^=Qd+V-QAU+QBVf-BPd-AQd-QAIN-UA1γ-IN+QBVA1γ+QAUB1δ+IN+QBIN-VB1δ. By now, we proved the following theorem.Theorem 7.
Suppose that the real symmetric(2r+1)-diagonal matrices Ka and Ha are given. If conditions (16) and (19) hold, then Problem 2 has a unique solution and the unique solution of Problem 2 can be expressed as(43)K^=Sα^⊗In,H^=Sβ^⊗In, where α^,β^ are given by (41) and (42), respectively.


## 4. A Numerical Example

Based on Theorems6 and 7, we can state an algorithm for solving Problems 1 and 2 as follows.Algorithm 8.
(1)
InputMa,Ka,Ha,Γ,Z. (2)
Form the orthonormal basis{Yij} by (7). (3)
SeparatingΓ and Z into real and imaginary parts yields ΓR,ΓI,ZR, and ZI. (4)
ComputeG,A,B,d,f according to (12), (13), and (14), respectively. (5)
ComputeJ=[A,-B],L=[B,A]. (6)
ComputeC,D,P,Q by (23) and compute J+ by (22). (7)
If condition (16) is not satisfied, then stop. (8)
ComputeA1,B1 by (25) and form the matrix T by (24). (9)
ComputeC1,D1,U,V by (27) and compute T+ by (26). (10)
If condition (19) is not satisfied, then stop. (11)
Form vectorsγ,δ by (33) and (34). (12)
ComputeS,α^,β^ by (31), (41), and (42), respectively. (13)
Compute the unique solution(K^,H^) of Problem 2 by (43).Example 1.
Consider a five-DOF system with analytical mass, damping, and stiffness matrices as follows:(44)Ma=diag⁡1,2,5,4,3,Ha=11.0-8.0000-8.014.0-3.5000-3.513.0-7.8000-7.813.5-9.0000-9.015.4,Ka=100-20000-20120-35000-3580-12000-1295-40000-40124.That is, Ha,Ka are symmetric 3-diagonal matrices. The model used to simulate the consistent experimental data is given by M=Ma,H=Ha, and (45)K=100-8000-8120-30000-3080-12000-12120-50000-50124.Note that the differences between Ka and K are (1,2),(2,3),(4,4), and (4,5) elements. The eigensolution of the model is used to create the experimental modal data. Assume that the measured eigenvalue matrix Γ=diag⁡{γ1,γ2} and corresponding eigenvector matrix Z=[z1,z2] are given by (46)Γ=diag⁡-13.553-1.3302i,-20.603-2.4171i,Z=0.032335+0.023541i-0.0090358-0.017342i0.30721-0.013746i-0.12644-0.055369i0.94805-0.05195i-0.3335-0.13811i0.3359+0.17504i0.85743-0.14257i0.19824+0.11415i0.69922-0.082153i. All the tests are performed on an Intel Core 3.39 GHz PC with a main memory of 2.99 GB running MATLAB 6.5. It is easy to check that conditions (16) and (19) hold (‖JJ+d-d‖=1.3277×10-13,  ‖ET(f-LJ+d)‖=0.10614). According to Algorithm 8 with 0.1710 CPU time (in seconds), we solve the optimal approximation solution of Problem 2 as follows: (47)K^=100-8000-8120.04-29.995000-29.99579.986-12.002000-12.002120.04-50.049000-50.049123.99,H^=10.999-8000-813.999-3.5013000-3.501313.009-7.7795000-7.779513.491-9.0188000-9.018815.388. We can see that Algorithm 8 works well and the updated matrices K^ and H^ are very close to the accurate stiffness matrix K and hysteretic damping matrix Ha with absolute errors (48)H^-Ha=0.043152,K^-K=0.09293. Define (49)res⁡γi,zi=γiMazi+K^+iH^zi. By Algorithm 8, we can get the results in Table 1.
Therefore, the measured eigenvalues and eigenvectors are reproduced in the new modelMaZΓ+(K^+iH^)Z=0, and the updated matrices K^,H^ are also symmetric 3-diagonal matrices, which means that the structural connectivity of the original model is preserved.Table 1
Results forres(γi,zi).
( γ i , z i ) ( γ 1 , z 1 ) ( γ 2 , z 2 ) r e s ( γ i , z i ) 0.039544 0.066018


## 5. Concluding Remarks

Updating a structural FE model to match measured modal data has been an important task for engineers. In this paper, a new direct method for updating symmetric and banded mechanical systems has been introduced, which uses measured eigenvalues and eigenvectors and the additional constraint related to structural connectivity to adjust stiffness and hysteretic damping matrices simultaneously. The method was developed based on the constrained optimization theory. The error function was minimized such that the discrepancies between the estimated and actual stiffness and hysteretic damping matrices are minimal. The method can maintain the physical connectivity of the numerical model, and the adjusted model is able to reproduce the measured modal data accurately when conditions (16) and (19) are satisfied. A numerical example confirms the effectiveness of the proposed approach.



---

*Source: 3727021-2018-02-19.xml*