Computers, Materials & Continua
DOI:10.32604/cmc.2022.020245
images

Article

Identification of Composite-Metal Bolted Structures with Nonlinear Contact Effect

Mohammad Ghalandari1, Ibrahim Mahariq2, Majid Pourghasem3, Hasan Mulki2 and Fahd Jarad4,5,*

1Department of Aerospace Engineering, Sharif University of Technology, Tehran, Iran
2College of Engineering and Technology, American University of the Middle East, Kuwait
3Faculty of Mechanical and Mechatronics Engineering, IHU, Tehran, Iran
4Department of Mathematics, Cankaya University, Etimesgut, Ankara, Turkey
5Department of Medical Research, China Medical University Hospital, China Medical University, Taichung, Taiwan
*Corresponding Author: Fahd Jarad. Email: fahd@cankaya.edu.tr
Received: 17 May 2021; Accepted: 02 July 2021

Abstract: The middle layer model has been used in recent years to better describe the connection behavior in composite structures. The influencing parameters including low pre-screw and high preload have the main effects on nonlinear behavior of the connection as well as the amplitude of the excitation force applied to the structure. Therefore, in this study, the effects of connection behavior on the general structure in two sections of increasing damping and reducing the stiffness of the structures that lead to non-linear phenomena have been investigated. Due to the fact that in composite structure we are faced to the limitation of increasing screw preload which tend to structural damage, so the investigation on the hybrid connection (metal-composite) behavior is conducted. In this research, using the two-dimensional middle layer theory, the stiffness properties of the connection are modeled by normal stiffness and the connection damping is modeled using the structural damping in the shear direction. Nonlinear frequency response diagrams have been extracted twice for two different excitation forces and then proposed by a high-order multitasking approximation according to the response range of the nonlinear finite element model for stiffness and damping of the connection. The effect of increasing the amplitude of the excitation force and decreasing the preload of the screw on the nonlinear behavior of the component has been extracted. The results show that the limited presented novel component model has been accurately verified on the model obtained from the vibration experimental test and the reduction of nonlinear model updating based on that is represented. The comparison results show good agreement with a maximum of 1.33% error.

Keywords: Nonlinear equation; hybrid joint behavior; thin laminated element; nonlinear frequency response function

1  Introduction

Nowadays, the hybrid structures is a common material which is widely employed in design [1], and analyzing [2–8] of various industrial engineering component with its specific bottlenecks and uncertainties. Since the advent of composite or metal-composite materials in aerospace structures, their applications have been steadily increasing. In addition to their primary applications in non-structural parts or secondary structural components, the application of these materials in the main parts of the components of structures is inevitably increasing.

The finite element method is a common method in the industry for estimating [6] the dynamic response of structures. But predicting a structure's behavior, using the finite element model, is always different from the actual behavior of the structure. The vibrational response of composite structures has complex mechanical behavior than metal structures. For example, composites are relatively more brittle compared to metals.

So far, many models have been proposed to describe the levels of screw connection and model their dynamic behavior. Spring-damper models [9,10], compensatory elements [9], generic elements [10], intermediate layer elements [11], porch models and more recently intermediate elements [12] are examples of these efforts.

Two types of middle layer elements are zero-layer thickness middle elements and thin middle layer element [13,14]. In the middle layer element with zero thickness, the compatibility relationships of the element usually include fixed values for both shear stiffness (tangent to the connection surface plate) and normal stiffness (perpendicular to the connection surface plate) [15,16]. However, in the thin middle element, the properties of the bonding surface are controlled by a thin band around the surface with different characteristics from the adjacent materials [17].

The studies on thin layer element were conducted in order to update and extract the stiffness coefficients of the connection using linear elastic compatibility relationships and test results [17,18]. The investigations prove the linear behavior of contact surface under the condition of low force value, isotropic material, or high preload value. The studies reveal that the behavior of connection can be modeled in the compatibility equation of thin middle layer element using the normal and shear stiffness with acceptable accuracy of nearly 99.8% percentage of reliability [18].

For nonlinear behavior of a structure that begins with the increase in excitation force, a careful modeling of the connection is needed, so that the nonlinear behavior of the contact surface can be described. Previous experimental studies have shown that the connection behavior depends on the domain of the structural response [19]. Since the structure weight decreases due to increase of stiffness and damping of the structure, by using the thin middle layer model and its compatibility relationships, the approximation of the second order is presented according to the response range [20]. In the non-linear model, the stiffness and damping is represented in the directions of normal stiffness and shear stiffness, respectively [21,22].

The main purpose of this paper is to provide a novel non-linear thin intermediate layer of four ordered nonlinear equation according to the response range for modeling the metal-composite behavior with two screws. Here, it is assumed that the total damping of the structure, including the composite damping and the screws connections damping [23], is supposed in the contact surface damping model and the composite structures and bonding layers are undamaged. So, the linear connection model and its linear frequency response are firstly extracted. Finally, the linear model obtained by the sensitivity method has been supplied with special values.

2  Connection Modeling Formulations

So far, a variety of connection models have been proposed for two identical metal beams. According to Fig. 1 in this study, we have a different connection including a metal beam and an isotropic quasi-composite beam with completely different elastic modulus in the direction of width.

Figure 1: Composite-metal single lap bolted joint

The most important thing in providing a connection model is that the model presented is compatible to the physical view governing the connection. In the middle layer theory proposed by Desai [16], it is assumed that the stiffness in the normal direction and the shear are separate. In this model, the stiffness changes in X direction can be neglected according to Fig. 2, and the middle layer at distance H from the neutral axis.

Figure 2: The 2D thin-layer element and degrees of freedom

2.1 Linear Modeling of Connection

The middle layer element used to model the connection is assumed to be a rectangular membrane element with degrees of rotational freedom [18]. In order to calculate the elemental stiffness matrix, the compatibility (stress-strain) relations of the matrix [D]must be determined.

{σ}=[D]{ε}(1)

{σxσyτxy}=[ε000Ec000Gc]{εxεyγxy}(2)

As mentioned in the middle layer theory, elements outside the original diameter have a very small value due to the separation of the stiffness in the normal direction and a cut equal to zero, as well as D11, and are assumed to reduce numerical errors. It suffices to place the matrix [D] in relation 3, where L is the length of the element and H is the height and t is the thickness of the middle layer (Fig. 2).

[Ke]=t∫0h⁡∫0L⁡[B]T[D][B]dxdy(3)

Matrix B is in the form of [B]=[L][N] in which [N] is the matrix of shape functions and [L] is the derivative matrix in two-dimensional mode [18]. It is only necessary to identify the parameters Ec and Gcusing the test results to model the linear behavior of the connection correctly.

2.2 Nonlinear Connection Modeling

As shown in previous research, the onset of a weak nonlinear area requires a sufficient amount of stimulation or looseness of the preload screw [19,20]. We know that the relationship between stress and strain (Eq. (3).) can be assumed that as follow for nonlinear mode:

σ(t)=f{ε(t),ε˙(t)}(4)

Considering that the stimulation is mono-harmonic, which is explained in the following test conditions, so the answer can be assumed to be mono-harmonic and therefore strain is represented as shown in Eq. (5). The linear stress-strain relationship is obtained by describing the concept of function as follows [21].

ε(t)=Xεsin⁡(ωt+θ)(5)

σ(t)=~(Er(Xε,ω)+jEj(Xε,ω))ε(t)=Eeq(Xε,ω)ε(t)(6)

Eeq(Xε,ω)=1Xεπ∫02π⁡f(Xεsin⁡(β),ωXεcos⁡(β))∗(sin⁡(β)+jcos⁡(β))dβwithβ=(ωt+θ)(7)

As shown in Eqs. (5)–(7), the strain-stress relationship is a function of displacement and frequency. In this study, it is assumed that because nonlinear behavior is around the first frequency and in the range of small changes in frequency, nonlinear behavior is only a polynomial function of the response range.

The following relationships are applied to model the stiffness and damping of the non-linear connection.

Ec(X)=Ec0−F1(X)(8)

Gc(X)=Gc0+j×F2(X)(9)

where Ec0 and Gc0 are the values obtained from the linear solution and X is the size of the response range of the structure to the sinusoidal excitation force with frequency ω. In the case of stiffness modeling, the value of the F1 function decreases from the linear solution due to the softening phenomenon. For damping modeling, structural stiffness, i.e., complex function form F2, with stiffness in the shear direction has been used to model it. According to empirical observations, the cause is the stiffening phenomenon of damping [21,22].

3  The Sample Under Study

The metal-composite hybrid structure used in this experimental study (See Fig. 3) includes an AL7075-T651 aluminum beam with dimensions of 360 × 39 × 11 mm3 and a carbon/epoxy composite beam with dimensions of 353 × 39 × 8.85 mm3, which is connected by two M8 bolts. Properties of aluminum beam is equal to 69.2 e 10 (Pa) elastic modulus and density of 2850 kg/m3, which have been extracted and employed based on the data available in the references.

Figure 3: Schematic of test setup

The composite beam is made of double-sided carbon fibers with 32 layers and symmetrical arrangement [[(±45)4/(±45)4]2/(0/90)2]s with Epoxy adhesives cyanoacrylate bonding. The density of composite with isotropic arrangement is equal to 1200\; Kgm3. Its elastic modulus is obtained by modal testing and model updating using Euler-Bernoulli beam theory. For this purpose, the desired beam has been subjected to a modal test in the free-free mode according to Fig. 4. The beam is triggered by a modal hammer and the sensor as illustrated in Fig. 4.

Figure 4: Composite beam test setups

The frequency response diagram is shown in Fig. 5. The results of model updating of the composite beam model by analyzing the sensitivity of the frequency response and its error percentage are given in Tab. 1. Based on the data in this table, the equivalent elastic modulus for the composite beam is obtained with a quasi-isotropic arrangement equivalent to 38.566 Gpa.

Figure 5: FRF of a composite beam with simple free-free boundary condition

4  Modal Test of Linear and Nonlinear System

The position of the hybrid structure is shown in Fig. 6 during the vibration modal test. As can be seen in Fig. 3, the hybrid structure in the aluminum section (according to industrial applications) is fixed to the boundary condition as cantilevered form. The shaker is connected to the beam at x = 125 mm distance from the retaining edge.

Figure 6: FRF of two bolted lap joint for composite-metal

To get the linear response, the structure is excited by a random signal and with a force of 5 N. The applied force to the structure is measured by a power sensor and its response is measured by one accelerator installed in positions x = 699 mm. The screws preload are set to 9 N.m and 14 N.m by one torque meter. Fig. 6 shows the frequency response chart for the 9 N.m preloads. The results of the resulting frequencies are shown in Tab. 2.

Figure 7: Composite-metal single lap bolted joint modal test setup

To obtain the nonlinear frequency response, we then extract the nonlinear frequency response with the help of vibration in the two amplitude ranges of 5 and 15 N in near range of the first frequency. In this way, we first determine the excitation force and excitation frequency and stimulate the structure with this harmonic force by shaker.

In each frequency, by obtaining the amplitude of the response and the amplitude of the stimulus, its value can be determined in the frequency response diagram in this specific frequency, and at the end, by connecting these points; the nonlinear frequency response graph can be reached. As shown in Fig. 8, the displacement of the resonant point and the reduction of the amplitude are observed with the increase of the excitation force.

Figure 8: nonlinear FRF of joint with 9 & 14 N.m preload

5  Identify the Linear Parameters of the Middle Layer

One of the most applicable approaches in identifies parameters and reduces the computational error in finite element method is the models based on sensitivity eigenvalues method. The matrix [S1] is the sensitivity matrix of eigenvalues to identification parameters, and {ε0} is the residual vector of the difference between the eigenvalues of experimental and what predicted values. To calculate {Δp}, it is sufficient to calculate the inverse of the sensitivity matrix [S1] and multiply by {ε}. Therefore, relationship 10 can only be solved by being in a repetitive cycle as follows. The goal in each step is to minimize the objective function as ∑i=1m⁡(((ωjb)2−(ωja)2)(ωjb)2)2 where ωjb are the test frequencies and ωja is the analytical frequencies.

{ΔP}=[S1]+{ε0}(10)

{P}new={P}old+{ΔP}(11)

According to the frequency response diagram of the metal-composite bonded with two bolt of Fig. 7 and the values of the frequencies obtained from the test in two pre-loads of 9 and 14 N.m, the identification of the parameters of the middle layer of E1 and G1 has been conducted. As can be seen in Tab. 3, the initial and final values of the identified parameters and the percentage of errors of the updated frequencies are given. The error nom L2of updated frequencies and test for corresponding modes is represented in Tab. 3.

6  Nonlinear Parameters Identification of the Middle Layer

In this section, discussion of nonlinear modeling is firstly completed, then a brief overview of the discussion of updating by frequency response diagram and parameter adjustment method during the convergence process is provided, and finally, the results of identifying nonlinear parameters are presented.

6.1 Nonlinear Middle Layer Element Relationships

Based on the mentioned theories in Section 2 and given that the nonlinear behavior of the joint is investigated in its first mode, it is assumed that the joint stiffness behavior is in the direction of the normal stiffness of the middle layer and the damping behavior of the joint occurs in the direction of its shear stiffness. As described in the previous section2, we use the following equations to model the stiffness and damping of our connection.

E1(X)=E10−F1(X)(12)

G1(X)=G10+j×F2(X)(13)

where, Fi(X) are functions to describe the stiffness and nonlinear damping of the joints. If we consider the Fi(X) function to be dependent on the response domain, we can assume that it consists of a polynomial of order 4 in terms of X as follows. The reason for using 4 ordered equations is to be able to provide a more accurate update existing study [15].

Fi(X)=δi1X4+δi2X3+δi3X2+δi4X+δi5(14)

i=1,2

Therefore, we can identify δij coefficients using the update method and nonlinear frequency response diagrams obtained from the test in Fig. 8.

In fact, by placing the functions Fi(X) in the relations (12) to (13) and obtaining the nonlinear Ec and Gc it is enough to form a new matrix D and integrate it around the element using Eq. (3). The obtained stiffness matrix has a real part and an imaginary part to describe the attenuation. Using the update method of the frequency response function to update the analytical model has the advantages. The frequency response sensitivity matrix to the P-update parameter is used to update the analytical model using the frequency response functions and is calculated as follows: where [Z] is the dynamic stiffness matrix of the system [23].

∂[H]∂p=−[H]T∂[Z]∂p[H](15)

[Z(ω)]=[−ω2[M]+jω[C]+[K]](16)

[H(ω)]=[Z(ω)]−1(17)

where [M], [C], [K] and ω are the matrices of mass, damping, stiffness and frequency, respectively. Using Eq. (15), the ratio update formula that the parameter P in the frequency ωK with the excitation and response points j and i can be calculated, respectively.

Hije(ωk)−Hija(ωk)=(−{Hia(ωk)}T∂[Z(ωk)]∂p{Hja(ωk)})Δp(18)

{εk}=[Sk]{Δp}

where {Hia(ωk)}Tand {Hja(ωk)}where the r-line vectors i and the column j-th are the frequency response matrix [Hij(ωk)], respectively. The e and a parameters represent the experimental and analytical FRFs. Arruda and Duarte [19] and later Balmes [20] defined the objective function as follows.

J({p})=‖20log10|{Hija(ωk,{p})}|−20log10|{Hije(ωk)}|‖(19)

By deriving from 20log10|Hij|| With respect to P, the derivative of the dynamic stiffness matrix and the frequency response matrix can be obtained.

‖20log10|{Hija(ωk,{p})}|−20log10|{Hije(ωk)}|‖=8.6859(Re(Hij)∂Re(Hij)∂p+Im(Hij)∂Im(Hij)∂pRe(Hij)2+Im(Hij)2)(20)

Using Eq. (20) at each excitation frequency, the sensitivity matrix [Sk] and the remaining vector {Δp}resulting from the difference of the values obtained from the test and analysis of the upgrade operation can be performed.

6.2 Investigation on the Results of Nonlinear Parameters of the Middle Layer Element

Examining the results of identifying the nonlinear parameters of the middle layer element, the examination of the updated nonlinear frequency response diagrams is conducted using the 4th order approximation according to the response amplitude and based on the representing result in Figs. 9 and 10. As can be seen, the frequency response diagrams for both preload and excitation force have been updated with high accuracy. The Figs. 9 and 10 are obtained by using the relation 18 the coefficients δij\;  for each of the frequencies ωk using the update method by the frequency response diagram and finally from the relations (12) to (13). The stiffness and damping matrices are obtained and using them the predicted frequency response diagram of the solution is obtained. Tab. 4 is calculated using the remaining 15% criterion for each of the graphs, indicating the high accuracy of the update process.

Figure 9: calculated and experimental results of FRF for torque 9 N.m

Figure 10: calculated and experimental results of FRF for torque 14 N.m

The changing diagrams of the Fi(X) functions are given in Fig. 11. Firstly, it is considered that the identification of nonlinear parameters located in the range 5 N. The linear functions F1(X), F2(X) are depicted on both preload forces of 9 N.m, 14 N.m. Parameters F1(X), F2(X) have been updated in the iterative process until convergence is achieved.

When the screws preload is 14 Nm, F1(X) is linear and F2(X) is a double polynomial, increasing and decreasing along the response amplitude, but at 9 Nm where the screw is looser, both functions are F1(X), F2(X) are polynomials of the second order and also the maximum amount of damping and stiffness of these functions is higher than the case where the clamp screws is closed. Using the 4th order polynomial allows the parameters F1(X), F2(X) to find the best way to identify the connection parameters.‏ Finally, the nonlinear parameters of the thin middle layer element for the excitation force of 10 N are identified using its frequency response diagram, both functions of F1(X), F2(X) are polynomial. In the case that the excitation force is increased in the same preloads, the values of the mentioned functions are increased compared to the previous case.‏ The results of Fig. 11 show that the Fi(X) diagrams increase with respect to increasing the excitation force range or decreasing the preload load of the polynomial degree, which is the nonlinear damping and stiffness factor.

Figure 11: The F1(X), F2(X) which updated at preloads 9 N.m, 14 N.m and loads 5 N, 15 N

7  Conclusion

In this research, the use of a thin middle layer model for a metal-composite joint has been studied. This modeling has been evaluated in both linear and non-linear modes by coding in MATLAB software. First, the flexural modulus of a composite beam has been identified due to its quasi-isotropic behavior. Then, with the help of linear test results, the parameters of connection, normal stiffness and shear were identified using frequency sensitivity analysis with undamaged composite structures and bonding layers. The results of this identification indicate that the joint stiffness is reduced by 103 order compared to the side beam stiffness. In the following, by describing the finite element model of the fourth order relative to the response amplitude, the description of the nonlinear connection behavior is investigated. The nonlinear modal connection test was performed and the nonlinear properties of the connection were identified by analyzing the sensitivity of the frequency response and using parameter adjustment to reduce the dissolution time and better convergence for first optimal mode. Since composite parts have high damping and have a limitation in applying pre-screws preload due to damage to the part, so we have no choice but to do nonlinear analysis. The high-order model presented in this paper has high accuracy (nearly about 99 percent) and sufficient efficiency (maximum 1.3% differences) to use it in applying similar analyzes. Also, considering that the obtained results satisfy the analysis correctly according to the previous observations, and the stiffness and damping parameters of the connection are polynomial and have been identified with high accuracy.

Funding Statement: This work was supported by College of Engineering and Technology, American University of the Middle East, Kuwait.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

 1.  Z. Wei, J. Feng, M. Ghalandari, A. Maleki and Z. Abdelmalek, “Numerical modeling of sloshing frequencies in tanks with structure using New presented DQM-BEM technique,” Symmetry, vol. 12, no. 4, pp. 655–667, 2020.

 2.  I. Mahariq, M. Kavyanpoor, M. Ghalandari, M. A. Nazari and D. T. Bui, “Identification of nonlinear model for rotary high aspect ratio flexible blade using free vibration response,” Alexandria Engineering Journal, vol. 3, no. 11, pp. 1073–1086, 2020.

 3.  M. Ghalandari, A. Ziamolki, A. Mosavi, S. Shamshirband, K. Chau et al., “Aeromechanical optimization of first row compressor test stand blades using a hybrid machine learning model of genetic algorithm, artificial neural networks and design of experiments,” Engineering Applications of Computational Fluid Mechanics, vol. 13, no. 1, pp. 892–904, 2019.

 4.  T. Abdeljawad, I. Mahariq, M. Kavyanpoor, M. Ghalandari and N. Nabipour, “Identification of nonlinear normal modes for a highly flexible beam,” Alexandria Engineering Journal, vol. 59, no. 4, pp. 1–11, 2020.

 5.  I. Mahariq and A. Erciyas, “A spectral element method for the solution of magnetostatic fields,” Turkish Journal of Electrical Engineering & Computer Sciences, vol. 25, no. 4, pp. 2922–2932, 2017.

 6.  I. Mahariq, M. Kuzuoğlu and H. I. Tarman, “On the attenuation of perfectly matched layer in electromagnetic scattering problems with spectral element method,” Applied Computational Electromagnetics Society Journal, vol. 29, pp. 701–710, 2014.

 7.  I. Mahariq, “On the application of the spectral element method in electromagnetic problems involving domain decomposition,” Turkish Journal of Electrical Engineering & Computer Sciences, vol. 25, pp. 1059–1069, 2017.

 8.  I. Mahariq, H. Kurt, H. I. Tarman and M. Kuzuoglu, “Photonic nanojet analysis by spectral element method,” IEEE Photonics Journal, vol. 6, pp. 1–14, 2014.

 9.  M. Ghalandari, S. Bornassi, S. Shamshirband, A. Mosavi and K. W. Chau, “Investigation of submerged structures flexibility on sloshing frequency using a boundary element method and finite element analysis,” Engineering Applications of Computational Fluid Mechanics, vol. 13, no. 1, pp. 519–528, 2019.

10. R. Omar, M. N. A. Rani, M. A. Yunus and M. H. Othman, “Finite element modelling and updating for bolted lap joints,” Journal of Mechanical Engineering, vol. 4, no. 3, pp. 202–222, 2017.

11. H. Ahmadian and H. Jalali, “Identification of bolted lap joints parameters in assembled structures,” Mechanical Systems and Signal Processing, vol. 21, no. 2, pp. 1041–1050, 2007.

12. C. Peter, M. N. A. Rani, H. Hashim, M. A. Yunus, W. Mirza et al., “Improving finite element modelling of the dynamic behaviour of a car trunk lid with modal based updating method,” In IOP Conf. Series: Materials Science and Engineering, vol. 788, no. 1, pp. 12002–12002, 2020.

13. E. Shahabi and H. Ahmadian, “Nonlinear lap joint interface modeling and updating strategies for assembled structures,” Journal of Theoretical and Applied Vibration and Acoustics, vol. 5, no. 2, pp. 131–152, 2019.

14. H. Ahmadian, J. E. Mottershead, S. James, M. I. Friswell and C. A. Reece, “Modelling and updating of large surface-to-surface joints in the AWE-MACE structure,” Mechanical Systems and Signal Processing, vol. 20, no. 4, pp. 868–880, 2006.

15. G. Bfer, “An isoparametric joint/interface element for finite element analysis,” International Journal for Numerical Methods in Engineering, vol. 21, no. 4, pp. 585–600, 1985.

16. C. S. Desai, M. M. Zaman, J. G. Lightner and H. J. Siriwardane, “Thin-layer element for interfaces and joints,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 8, no. 1, pp. 19–43, 1984.

17. D. Wang, “Identification for joint interfaces with correlation analysis of instantaneous dynamics,” Archive of Applied Mechanics, vol. 90, no. 1, pp. 187–198, 2020.

18. Y. Q. Long, S. Cen and Z. F. Long, Advanced Finite Element Method in Structural Engineering, Berlin, Heidelberg: Springer, pp. 100–112, 2009.

19. W. E. Velde, Multiple-input Describing Functions and Nonlinear System Design, New York: Citeseer, pp. 80–85, 1968.

20. E. Ferhatoguglu, E. Cigugeroguglu and H. N. Özgüven, “A modal superposition method for the analysis of nonlinear systems,” in Special Topics in Structural Dynamics, vol. 6, Springer, vol. 12, no. 20, pp. 269–278, 2020.

21. H. Jalali, H. Ahmadian and J. E. Mottershead, “Identification of nonlinear bolted lap-joint parameters by force-state mapping,” International Journal of Solids and Structures, vol. 44, no. 25, pp. 8087–8105, 2007.

22. H. G. Natke, “On regularization methods within system identification,” in Inverse Problems in Engineering Mechanics, Berlin, Heidelberg: Springer, pp. 3–20, 1993.

23. M. Friswell and J. E. Mottershead, in Finite Element Model Updating in Structural Dynamics, University of Wales Swansea, UK: Springer Science & Business Media, pp. 35–40, 2013.

images
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.