Boundary condition identification of a clamped honeycomb sandwich panel based on thinlayer element
Yu Xu^{1} , Zhifu Cao^{2} , Dong Jiang^{3}
^{1, 2, 3}School of Mechanical and Electronic Engineering, Nanjing Forestry University, Nanjing, 210037, China
^{3}Corresponding author
Journal of Vibroengineering, Vol. 21, Issue 8, 2019, p. 22862295.
https://doi.org/10.21595/jve.2019.20725
Received 15 April 2019; received in revised form 6 September 2019; accepted 26 October 2019; published 31 December 2019
JVE Conferences
The stiffness of boundary conditions in mechanical structures is difficult to represent. An approach on high fidelity modeling of the clamped boundary condition is proposed in this paper. Firstly, the normal and tangential stiffness of the contact surface of the clamped boundary condition is parameterized by using thinlayer element with isotropic material. Secondly, the material parameter is identified by minimizing the discrepancy between the calculated results and the experimental data, and parameter identification can be treated as an optimization problem. Experimental investigation is undertaken to verify the proposed method by employing an aluminum honeycomb panel, the numerical model of which is constructed by using the equivalent theory. Thinlayer elements with different properties are used to simulate the mechanical properties in different area of the boundary conditions, and the experimental modal data is adopted to identify the material parameters. Results show that the weighted matrix is a crucial option in the parameter identification procedure, and the width to thickness ratio of the thin layer element has a great influence on the identification results. After parameter identification, the error of the first three order of the modal frequencies is less than 2.7 %, the thinlayer element can accurately reflect the mechanical performance of clamped boundary condition.
Keywords: boundary conditions, honeycomb panel, thinlayer element, stiffness identification, modal data.
1. Introduction
Structural dynamic performance is strongly affected by boundary conditions, which should be considered appropriately in both theoretical and numerical analysis. The effects of the boundary conditions on the performance of the mechanical structures are always underestimated and difficult to be determined [1]. The exist inappropriate methods for modeling of boundary stiffness, such as idealizing constraints into completely rigid and modelling ostensibly similar boundaries refer to the previous experiences, can be applied to characterize the boundary conditions in an oversimplified way, which will inevitably lead to errors in structural dynamics modeling.
Due to the crucial importance of mechanical behavior of boundary condition in engineering structures [2, 3], investigations on modeling of which have attracted interests of researchers [25]. The intuitive approach to determine the stiffness of clamped boundary conditions by using contact model between rough surfaces [68], Ahmadian and Zamani [9] studied the local nonlinear mechanism by experimentally observing the behavior of a beam structure with microslip/slap at the clamped end, these models can reflect the contact mechanism through analysis the asperities behavior of the real contact area, and the qualitative nonlinear contact stiffness can be derived. However, the contact analysis for boundary conditions will be very expensive in time consuming when apply in modeling of complex structures. Simplified equivalent modeling is the alternative method to represent the contact mechanics relationship and the boundary stiffness.
Identification of the boundary conditions from experimental data is an effective solution to improve the modeling accuracy. Pabst and Hagedorn [10] proposed an inverse method to identify the boundary stiffness and damping parameters of a viscoelastically clamped beam by using modal data, in which the boundary conditions was parametrized with a torsional springdamper, and a translational springdamper. Yang et al. [11] constructed a numerical model to solve multivariables inverse problems for a beam with elastic boundary support, in which the static responses were employed and the elastic supports and constitutive parameters can be identified by using LevenbergMarquardt algorithm. Wang [12] extended the crossmodel crossmode method for model updating by incorporating the connection flexibility and boundary conditions, the flexible joints and boundary conditions can be simultaneously updated. Yang et al. [13] proposed an approach using frequency response functions to identify the joint parameters including the rotational stiffness, which the accuracy of the identified results is strongly sensitive to the parameters. Baruh and Boka [14] modeled boundary conditions using axial and torsional springs, the orthogonality condition and the system eigenvectors were adopted to identify the spring constants. Boundary stiffness are simplified as spring models in the above researches, modeling difficulties will arise when apply spring models to complex structures, and the deficiency is obvious for modeling large contact surface.
Parameterization plays an important role in the stiffness identification of boundary conditions. Mottershead et al. [1] indicated that parameterization is critical in updating of boundary conditions, the boundary stiffness in a cantilever was modeled as the effective length of elements closest to the joint, but in some cases, the corrected matrices will lose connection to the physical significance while matching the experimental modal data. Noll et al. [15] inversely identified an elastic metal beam endsupported by elastomeric joints by revising nondiagonal terms in the stiffness matrix. The above stiffness characterization methods can be approximately parametrized the boundary stiffness, however, modeling complexities are still exit. In order to solve this problem, thinlayer element has been introduced to model the dynamics of mechanical joints [1618], The concept of which was firstly proposed by Desai [19] who applied to the contact analysis in geotechnical engineering, and has been widely applied in equivalent modeling of joints, The normal and tangential stiffness of the contact surface can be accurately represented by the thinlayer elements.
Based on the basic theory of thinlayer element, this paper parameterized the contact surface of the clamped boundary. The elastic parameters of the thinlayer element are identified by using the experimental modal data. Besides, parameters such as the weight matrix used in the parameter identification, the width to thickness ration of the thinlayer element should be considered in the identification procedure, because of these parameters control the identification efficiency and the accuracy.
2. Basic theory
2.1. Thinlayer element
Thinlayer element is used to simulate the contact interface by defining a layer of elements. The size of thinlayer element is ${l}_{1}\times {l}_{2}\times d$, according to the principle of virtual displacement, the virtual work of the element can be expressed as:
in which ${l}_{1}$ and ${l}_{2}$ are the inplane lengths of the thinlayer elements in the local coordinate system, and $d$ is the outofplane thickness. Fig. 1 shows isoparametric transformation of thinlayer element. Stiffness matrix [$K$] of thinlayer element can be expressed as:
where $\xi $, $\eta $, $\zeta $ are natural coordinate symbols, [$J$] is the Jacobian matrix, and when the coordinates of the natural and local coordinates are consistent, which can be written as:
By using of twonode Gaussian integration method, the following numerical expression for stiffness matrix [$K$] of thinlayer element can be obtained:
where ${w}_{\xi}$, ${w}_{\eta}$, and ${w}_{\zeta}$ are the Gaussian integral weight function.
Fig. 1. Isoparametric change of thinlayer element
It is assumed that the thickness $d$ of the thinlayer element is much smaller than the feature sizes ${l}_{1}$ and ${l}_{2}$. In this case, the inplane strain components $\left({\epsilon}_{x},{\epsilon}_{y},{\gamma}_{xy}\right)$ and the stress components $\left({\sigma}_{x},{\sigma}_{y},{\tau}_{xy}\right)$ of the element are ignored. $\partial {N}_{i}/\partial z$ is larger than $\partial Ni/\partial x$ and $\partial Ni/\partial y$ and can be analyzed by shape functions of element. Under this circumstances, $\partial Ni/\partial x=\partial Ni/\partial y\approx $0, which leads to the strain component ${\epsilon}_{x}={\epsilon}_{y}={\gamma}_{xy}\approx $0. Therefore, only three strain components at the Gaussian point of the thinlayer element are not zero. The strain component can be simplified to $\epsilon ={\left[{\epsilon}_{z}{\gamma}_{yz}{\gamma}_{zx}\right]}^{T}$. The normal direction ($n$) and two tangential direction ($t$) of the contact surface are defined as the local coordinate system $\left(z,x,y\right)$ of thinlayer element. According to the previous analysis, assuming that the normal and tangential contact properties of the interface are independent and the two tangential contact properties are consistent. The constitutive equation of the thinlayer element that characterizes the contact properties of the interface is given by the following equation:
${E}_{n}$ and ${G}_{t}$ are the normal elastic modulus and tangential shear modulus, respectively, which are determined by the contact surface properties of the connecting structure. If the contact performance is coupled to the normal and tangent, then the coupling term can be added to the constitutive relationship Eq. (5). In the finite element analysis, the thin layer element integrated isotropic material can be used to simulate the contact interface, the constitutive relationship is expressed as:
in which $\lambda $ is lame constant, $G=E/2\left(1+\nu \right)$ is shear modulus. the isotropic material has only 2 independent material parameters. According to the basic theory of thinlayer element, the constitutive equation of the material will be degenerated to the following form:
In this condition, the normal elastic constant and the tangent elastic constant are correlated.
The selection of different thickness of thin layer element will have a great influence on the calculation results. When the thickness is too large, the element has six strain components, and it is difficult to accurately reflect the mechanical characteristics of the contact interface. When the thickness is too small, the value of the determinant of Jacobian matrix approaches to zero, the matrix is ill conditioned and difficult to solve, and the displacement strain relation can’t be calculated by finite element method. A ratio coefficient is defined for the selection of thickness of a thinlayer element as follows:
Desai et al. [20] recommended that $R$ should be selected in the range of 10100.
The constitutive relation of the material can be integrated into the thinlayer element and the finite element dynamic analysis of the connection structure can be carried out. Vibration test data are used to identify the parameters of the thinlayer element. And the influence of weighted matrix and the ratio coefficient of thinlayer element $R$ on the identification results is researched.
2.2. Parameter identification
In this paper, material parameter identification of thinlayer elements is transformed into an optimization problem. The minimum value of the objective function is as follows:
in which ${\epsilon}_{f}={z}^{m}\uff0d{z}^{a}\left(p\right)$ is the residual term represents the discrepancies between the experimental modal data and the calculated results, the modal data including the modal frequencies and mode shapes. $W$ is the weighted matrix, which reflects the relative weight of modal parameters. In general, $W$ is the identity matrix $\mathbf{I}$ or [diag(${z}^{m}$)]^{2}.
The real value of the elastic parameters is set to $p$ and the initial value is ${P}_{0}$. The first order Taylor expansion of the eigenvalues $\mathrm{\Lambda}$ and the mode matrix $\varphi $ at the ${p}_{0}$ point is as follows:
where ${p}_{j}$ is the updating elastic parameters. When objective function is used for iterative solution, the problem of the $j$th iterations is described as followed:
where ${S}_{j}$ represents the weighted Jacobian matrix of modal parameters to the updating parameters which is generally treated as the sensitivity matrix. ${P}_{j+1}$ can be obtained by iterative solution of numerical analysis method until the identification parameter $p$ converges. At the same time, the calculation of modal parameters satisfies the requirement of precision. Finally, the accurate contact surface material parameters are obtained. The parameter identification procedure is shown as shown in Fig. 2.
Fig. 2. The flow chat of parameter identification
3. Case studies: a clamped honeycomb panel
The one end clamped aluminum honeycomb sandwich panel which can be seen in Fig. 3, The contact surface is modeled by thinlayer elements of three different properties. The honeycomb sandwich panel is shown in Fig. 4, the dimensions is 400 mm×400 mm×16.3 mm. elastic Modulus ${E}_{s}$ is 68 GPa and density is 2700 kg/m^{3}, Adhesive material is epoxy resin, the elastic Modulus ${E}_{s}$ is 5 GPa and density is 1500 kg/m^{3}.
3.1. Equivalent modeling
The honeycomb core is generally considered to resist transverse shear deformation, neglecting lateral shearing stress, the honeycomb core is equivalent to an orthotropic material of uniform thickness and the same height as the honeycomb core. As a result, the honeycomb sandwich composite is a threelayer structure, an upper and a lower panels, and the middle is an equivalent core sandwich structure. Face sheet and the core layer are bonded through the adhesive layer, the mechanical property of the adhesive layer is worse than the panel and the core material. The interlaminar shear effect caused by the weak layer has an important influence on the macroscopic mechanical properties of the sandwich material. Therefore, the influence of the adhesive layer on the dynamic characteristics of the honeycomb sandwich composite can’t be neglected. Fig. 5 shows the finite element model of honeycomb plate.
Fig. 3. The experimental arrangement of the clamped honeycomb panel
Fig. 4. The geometry of honeycomb sandwich panel
The honeycomb sandwich panel is modeled using the equivalent method. Shear deformation, axial tension and compression deformation analysis equivalent elastic parameters are ignored. The cell wall of the core is simplified as a planar EulerBernoulli beam. The equivalent elastic parameters in the $x$$y$ plane can be obtained. In particular, the unit cell of honeycomb core is a regular hexagon that the $l=h$, $\theta =$30°. the inplane elastic parameters of the core are expressed as:
where ${E}_{cx}$ and ${E}_{cy}$ are the equivalent elastic modulus, ${G}_{cxy}$ is equivalent shear modulus, ${E}_{s}$ is elastic modulus of honeycomb core material. Elastic modulus ${E}_{cz}$ in the z direction can be got by tension and compression stiffness equivalent:
Equivalent material parameters which core equivalent shear modulus ${G}_{cxz}$ and ${G}_{cyz}$ have the most impact on the dynamic performance of the equivalent model. The stress distribution is complicated When the honeycomb core considers the pure shear of three planes. The minimum potential energy principle and the principle of the minimum residual energy can be used to estimate the range of the two outofplane shear modulus:
where ${G}_{s}$ in the formula is the shear modulus of the honeycomb core material; $\gamma $ is the updating factor, which theoretical value is equal to 1, generally take 0.40.6 in engineering applications. As the following formula, the equivalent density of honeycomb core can be obtained by mass equivalent:
In summary, aluminum honeycomb sandwich panel and glue layer are constructed as an laminate. Orthotropic honeycomb core layers are modeled using solid elements. Clamped boundary contact layer is modeled by isotropic thinlayer elements. The finite element model’s boundary conditions are achieved by constraining six degrees of freedom.
Fig. 5. FEM of honeycomb plate
3.2. Parameter identification
3.2.1. Weight matrices
The influence of different weighting matrices on the results of parameter identification is studied by establishing a thinlayer element with $R=$100. The frequency convergent curves when weighted matrices are identity matrix and [diag(${z}^{m}$)]^{2} are shown in Fig. 6, Table 1 is the results comparison of natural parameter identification.
Table 1. ${W}_{\epsilon}=\mathbf{I}$ Comparison of natural frequency
Mode

${f}^{m}$ /Hz

Initial

${W}_{\epsilon}=\mathbf{I}$

${W}_{\epsilon}={\left[\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({z}^{m}\right)\right]}^{2}$


${f}^{a}$ /Hz

Error (%)

${f}^{a}$ /Hz

Error (%)

${f}^{a}$ /Hz

Error (%)


1

49.59

71.94

45.07

47.13

–4.96

49.27

–0.64

2

199.8

218.09

9.15

203.84

2.03

205.05

2.63

3

411.93

473.83

15.3

412.04

0.03

416.16

1.03

Fig. 6. Frequency iteration convergence curve: a) ${W}_{\epsilon}=\mathbf{I}$; (b) ${W}_{\epsilon}={\left[\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({z}^{m}\right)\right]}^{2}$
a)
b)
After the identification, the accuracy of the inherent frequency is greatly improved. The error of the first three order calculation natural frequency is less than 2.63 %, when [diag(${z}^{m}$)]^{2} is selected as the weighted matrix ${W}_{\epsilon}$.
3.2.2. Width to thickness ratio
Based on the above results, the weighted matrix is determined as ${W}_{\epsilon}={\left[\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({z}^{m}\right)\right]}^{2}$. Then researches of parameter identification are carried out when the $R$ of thin layer elements is 10, 50, 100, 500, respectively.
Frequency iteration convergence curve under different $R$ is shown in Fig. 7. The convergence error is reduced with the increase of the width and thickness ratio. However, more iterations are need for convergence.
Table 2. Comparisons of calculated frequency errors under different $R$
Mode

${f}^{m}$/Hz

$R=$10

$R=$50

$R=$100

$R=$500


${f}^{a}$ /Hz

Error (%)

${f}^{a}$ /Hz

Error (%)

${f}^{a}$ /Hz

Error (%)

${f}^{a}$ /Hz

Error (%)


1

49.59

42.53

–14.23

49.30

–0.56

49.27

–0.64

49.32

–0.54

2

199.8

200.53

0.36

205.10

2.65

205.05

2.63

205.14

2.67

3

411.93

414.37

0.59

414.70

0.68

416.16

1.03

414.47

0.62

Table 2 is the comparison of calculation natural frequency and parameter identification results of thinlayer element under different $R$. In this case, the $R$ is increased by changing the direction of the thickness. The stiffness properties of the contact surface are determined by the elastic parameters and the thickness of the thin layer element. The ratio coefficient $R$ of the thin layer element has a great influence on the identification results. The error of convergence result is improved with the increase of $R$ but more iterations are needed.
Fig. 7. Frequency iteration convergence curve under different $R$: a) $R=$10; b) $R=$50; c) $R=$ 100; d) $R=$500
4. Conclusions
An approach on identifying the stiffness of the contact surface of the clamped boundary condition is proposed in this paper. The contact surface is modeled based on the theory of thinlayer element, and modal parameters are used to identify the material parameters of thinlayer element which represents the stiffness of the contact surface. The equivalent modeling method are employed to establish the simulation model of aluminum honeycomb panel. Three kinds of thin layer elements with different properties are used to simulate the mechanical change of the boundary. The effect of the weighted matrix and the proportion of $R$ on the convergence results are discussed. It is concluded that the convergence error of the parameter can be reduced effectively when weighted matrix is [diag(${z}^{m}$)]^{2}.
References
 Mottershead J. E., Friswell M. I., Ng G. H. T., et al. Geometric parameters for finite element model updating of joints and constraints. Mechanical Systems and Signal Processing, Vol. 10, Issue 2, 1996, p. 171182. [Publisher]
 Demir C., Izmirli S. B. The effects of support size on the vibration of the point supported plate. International Journal of Physical Sciences, Vol. 6, Issue 8, 2011, p. 19201928. [CrossRef]
 Friswell M. I., Wang D. The minimum support stiffness required to raise the fundamental natural frequency of plate structures. Journal of Sound and Vibration, Vol. 301, Issues 35, 2007, p. 665677. [Publisher]
 Ritto T. G., Sampaio R., Aguiar R. R. Uncertain boundary condition Bayesian identification from experimental data: A case study on a cantilever beam. Mechanical Systems and Signal Processing, Vol. 68, Issue 69, 2016, p. 176188. [Publisher]
 Pietrzyk M., Madej L., Rauch L., et al. Identification of material models and boundary conditions. Computational Materials Engineering, Achieving High Accuracy and Efficiency in Metals Processing Simulations, 2015, p. 153208, https://doi.org/10.1016/B9780124167070.000041. [CrossRef]
 Goedecke A., Jackson R. L., Mock R. A fractal expansion of a three dimensional elasticplastic multiscale rough surface contact model. Tribology International, Vol. 59, 2013, p. 230239. [Publisher]
 Belghith S., Mezlini S., Belhadjsalah H., et al. Modeling of contact between rough surfaces using homogenisation technique. Comptes Rendus Mecanique, Vol. 338, Issue 1, 2010, p. 4861. [Publisher]
 Pei L., Hyun S., Molinari J., et al. Finite element modeling of elastoplastic contact between rough surfaces. Journal of the Mechanics and Physics of Solids, Vol. 53, Issue 11, 2005, p. 23852409. [Publisher]
 Ahmadian H., Zamani A. Identification of nonlinear boundary effects using nonlinear normal modes. Mechanical Systems and Signal Processing, Vol. 23, Issue 6, 2009, p. 20082018. [Publisher]
 Pabst U., Hagedorn P. Identification of boundary conditions as a part of model correction. Journal of Sound and Vibration, Vol. 182, Issue 4, 1995, p. 565575. [Publisher]
 Yang H. T., Yang B., Li H. T. Numerical solution of multivariables inverse problem for a beam with elastic boundary supports. Dalian Ligong Daxue Xuebao/Journal of Dalian University of Technology, Vol. 51, Issue 4, 2011, p. 469472, (in Chinese). [CrossRef]
 Wang S. Model updating and parameters estimation incorporating flexible joints and boundary conditions. Inverse Problems in Science and Engineering, Vol. 22, Issue 5, 2013, p. 727745. [Publisher]
 Yang T. C., Fan S. H., Lin C. S. Joint stiffness identification using FRF measurements. Computers & Structures, Vol. 81, Issues 2829, 2003, p. 25492556. [Publisher]
 Baruh H., Boka J. B. Modeling and identification of boundary conditions in flexible structures. International Journal of Analytical and Experimental Modal Analysis, Vol. 8, 1993, p. 107117. [CrossRef]
 Noll S., Dreyer J. T., Singh R. Identification of dynamic stiffness matrices of elastomeric joints using direct and inverse methods. Mechanical Systems and Signal Processing, Vol. 39, Issues 12, 2013, p. 227244. [Publisher]
 Schmidt A., Bograd S., Gaul L. Measurement of join patch properties and their integration into finiteelement calculations of assembled structures. Shock and Vibration, Vol. 19, Issue 5, 2012, p. 11251133. [Publisher]
 Jalali H., Hedayati A. Ahmadian H. Modelling mechanical interfaces experiencing microslip/slap. Inverse Problems in Science and Engineering, Vol. 19, Issue 6, 2011, p. 751764. [Publisher]
 Bograd S., Reuss P., Schmidt A., et al. Modeling the dynamics of mechanical joints. Mechanical Systems and Signal Processing, Vol. 25, Issue 8, 2011, p. 28012826. [Publisher]
 Desai C. S., Zaman M. M., Lightner J. G., et al. Thinlayer element for interfaces and joints. International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 8, Issue 1, 1984, p. 1943. [Publisher]
 Sharma K. G., Desai C. S. Analysis and implementation of thinlayer element for interfaces and joints. Journal of Engineering Mechanics, Vol. 118, Issue 12, 1992, p. 24422462. [Publisher]