Full-field dynamic strain reconstruction of an aero-engine blade from limited displacement responses

Abstract

Blade strain distribution and its change with time are crucial for reliability analysis and residual life evaluation in blade vibration tests. Traditional strain measurements are achieved by strain gauges (SGs) in a contact manner at discrete positions on the blades. This study proposes a method of full-field and real-time strain reconstruction of an aero-engine blade based on limited displacement responses. Limited optical measured displacement responses are utilized to reconstruct the full-field strain. The full-field strain distribution is in-time visualized. A displacement-to-strain transformation matrix is derived on the basis of the blade mode shapes in the modal coordinate. The proposed method is validated on an aero-engine blade in numerical and experimental cases. Three discrete vibrational displacement responses measured by laser triangulation sensors are used to reconstruct the full-field strain over the whole operating time. The reconstructed strain responses are compared with the results measured by SGs and numerical simulation. The high consistency between the reconstructed and measured results demonstrates the accurate strain reconstructed by the method. This paper provides a low-cost, real-time, and visualized measurement of blade full-field dynamic strain using displacement response, where the traditional SGs would fail.

Graphical abstract

Keywords

aero-engine blade / displacement response / dynamic strain reconstruction / mode shape / strain gauge

Introduction

Topology optimization (TO), which has been extensively studied over the last decades, is a process of determining optimal layout of materials inside a given design domain. TO has been applied to various structural optimization problems, such as minimum compliance [1,2] vibration [3,4], and thermal issues [5,6], after Bendsøe and Kikuchi [7] proposed the homogenization method. Homogenization is a material distribution method in which a design domain is discretized into small rectangular elements, and each element contains an artificial composite material with microscopic voids. The proposal of the homogenization method was followed by a parallel exploration of the solid isotropic material with penalization (SIMP) method [8,9], which uses an artificial isotropic material whose physical properties are expressed as a function of continuous penalized material density (design variables). The phase-field method [10,11], which is also a material distribution method, is based on the theory of phase transitions. A different type of method called evolutionary structural optimization (ESO) [12,13] has also been proposed. This method eliminates elements with the lowest criterion value on the basis of certain heuristic criteria. ESO is computationally expensive because it requires a much larger number of iterations with an enormous number of intuitively generated solutions compared with material-based methods.

However, these conventional TO methods, which are based on element-wise design variables, suffer from numerical instability problems, such as checkerboards and mesh dependency. Accordingly, several studies have proposed prevention methods. The use of high-order elements has been proven to be an effective means to prevent checkerboards [14,15], but this method entails a considerable increase in computation time. Various filter techniques have been utilized to mitigate checkerboards and mesh dependency because these techniques require only a small amount of extra computation time and are simpler to implement than other methods [16,17]. Notably, filter schemes are purely heuristic. Other prevention schemes, such as perimeter control and gradient constraint, which often make the optimization procedure unstable, are yet to be improved.

A new type of TO approach is the level set-based method (LSM), which was developed by Sethian and Wiegmann [18] to numerically track the motion of structural boundaries. In LSM, boundaries are represented as zero level set contour of an implicit high-dimensional function (level set function or LSF), in which boundary motion, merging, and introduction of new holes are performed. The evolution of a structural interface with respect to time is tracked by solving a Hamilton-Jacobi (H-J) partial differential equation (PDE), where transporting LSF along its outward normal direction is equivalent to moving the boundaries along the descent gradient direction. The conventional level set-based TO approach uses shape derivatives coupled with the original LSM for boundary tracking [19,20]. In this approach, regularization that reinitializes LSF to be signed distance to a zero level set is employed to control the slope of LSF. This conventional approach is updated by solving the H-J equation via an explicit up-wind scheme [21,22]. Variations of the conventional approach include parameterizing LSF using various basis functions, such as finite element method (FEM) basis functions [23], radial basis functions (RBFs) [24,25], and spectral parameterization [26], and corresponding methods for solving the H-J equation. By defining the interfaces between two material phases via the iso-contour of LSF, LSM can handle shape and topology changes during the optimization procedure and provides optimal structures with clear boundaries that are free of checkerboard patterns. Notably, most LSMs rely on finite elements wherein boundaries are still represented by discretized mesh in the analysis field unless alternative techniques are utilized to map the geometry to the analysis model.

Most of these TOs are performed in a fixed domain of finite elements where FEM is used to solve optimization problems. Currently used FEMs are often based on Lagrange polynomials for analysis while the geometrical representation of structures relies on non-uniform rational B-splines (NURBS), which are the criteria in computer aided design (CAD) systems. Thus, conversion of NURBS-based representation into one that is compatible with Lagrange polynomials, that is, mesh generation, is required in structural analysis. The disadvantages of FEM are as follows. First, the geometry approximation inherent in the FEM mesh may generate an approximate error. Second, frequent data interaction between geometry description and the computational mesh, which can be found in several calculations (e.g., fluid, large deformation, and shape optimization problems), is cumbersome and error-prone. An integrating method, namely, isogeometric analysis (IGA) [27], for unifying analysis and CAD processes has been proposed to overcome these disadvantages. This method employs the same basis functions as a technique for describing and analyzing the geometric model, which features the IGA method and CAD-based parameterization of field variables in an isoparametric manner. The first work on isogeometric approximation dates back to 1982 [28]; however, this method is considerably different from the IGA method. Several methods have been devised to help alleviate the difficulties faced by IGA. Special parameterization techniques, such as variational harmonic-based methods [29,30] and analysis-aware parameterization methods for single [31] and multi-domain geometries [32], have been proposed for the computational domain. Alternatives to NURBS, such as T-splines [33,34], polynomial splines over hierarchical T-meshes (PHT-splines) [35–37], and Powell-Sabin splines [38], have been studied for local refinement in IGA due to the limitation of the tensor product form of NURBS in computation refinement. Methods of parameterization of the interior domain while retaining the geometry exactness from the CAD model have been devised [39,40], and isogeometric collocation method is one of the most important among these methods [39]. With regard to interior discretization obstacles, the isogeometric boundary element method is a suitable candidate [41] because only boundary data are required for analysis, and it enables stress analysis [42], fracture analysis [43,44], acoustic analysis [45], and shape optimization [46,47]. Considering that the integral efficiency of IGA is limited by the tensor product structure of NURBS, an efficient quadrature rule, which is more suitable for NURBS-based IGA compared with the Gaussian quadrature rule, has been proposed in Ref. [48]. IGA has been applied to a wide range of problems, such as structural vibrations [49], fluid-structure interaction [50,51], heat conduction analyses [52], shape optimization [53,54], shell analyses [55], TO [56], and electromagnetics [57].

IGA-based shape optimization has been extensively investigated because remeshing is eliminated during the optimization process. IGA has also been recently applied to TO. The most commonly used TO approach is material based, and the most commonly solved TO problem is the minimum compliance case. In Ref. [58], TO was proposed based on isogeometric shape optimization. B-spline curves were introduced to represent the material boundary, and the coordinates of their control points were considered design variables. In Ref. [59], the trimmed spline surface technique was used for spline-based TO. In Refs. [60,61], optimality criteria (OC) and the method of moving asymptotes were implemented in the isogeometric-based SIMP method. In Ref. [62], the TO problem was solved by using a phase-field model, and IGA was utilized for the exact representation of the design domain. In Ref. [56], IGA was introduced to a level set TO method, and NURBS basis functions were used for geometry description and LSF parameterization.

Most studies on TO were restricted to compliance optimization, and the number of studies on TO of dynamic problems is limited. Dynamic problems, such as vibration and noise, are critical in many engineering fields, such as aeronautical and automotive industries. As a typical dynamic problem, structure vibration is controlled by the structure’s dynamic characteristics, which are usually represented by the eigenfrequencies of the structure. Thus, eigenfrequency optimization plays an important role in improving dynamic characteristics.

As an important research topic, eigenfrequency optimization has been studied by many scholars. Díaaz and Kikuchi [63] extended TO to eigenfrequency optimization within the framework of the homogenization approach. Tenek and Hagiwara [64] introduced the SIMP method to structural shape optimization and TO in eigenfrequency problems. Xie and Steven [65] applied the ESO scheme to TO, and the specified eigenfrequency of the structure was used as a constraint. Additional research has also been conducted on frequency optimization [66–70]. However, standard density-based TO methods are unsuitable for eigenfrequency optimization due to localized modes in low-density areas [71]. Low-density areas are much more flexible than areas with full densities; hence, they control the lowest eigenmodes of the entire structure. By changing the penalization of stiffness in the SIMP method, a modified algorithm has been proposed and applied to circumvent this numerical instability problem [3]. LSM that employs a crisp description of structure boundaries has advantages over the density-based approach. The LSM approach can avoid artificial modes localized in the weak phase, which makes LSM a choice for eigenfrequency optimization. Many studies have been performed on level set TO for eigenfrequency problems [72–75].

IGA-based TO has been extensively studied. However, research on the combination of the level set approach and IGA is limited, and that on eigenfrequency problems is absent. In this study, we develop a new optimization method to formulate the TO problem for cases with maximum fundamental eigenfrequency by using LSM under the framework of IGA instead of conventional finite element analysis (FEA). Given that the OC algorithm is particularly efficient for problems with many design variables and few constraints, we consider the OC method for the solution of the optimization problem and conduct a sensitivity analysis. The rest of this paper is organized as follows. In Section 2, IGA and NURBS-based TO are summarized. The eigenvalue optimization problem is described in Section 3, and the TO model is proposed in Section 4. Section 5 presents numerical examples to demonstrate the validity of the proposed approach. The conclusions and discussions are shown in Section 6.

IGA for level set-based TO

In IGA, geometric modeling and analysis are integrated by using NURBS, where the basis function is a bridge of the parameter domain, physical field, and numerical solution. In the proposed method, we consider the NURBS basis function as a bound between IGA and parameterized LSM. We provide a brief review of IGA and NURBS-parameterized LSM.

NURBS-based IGA

We assume that geometrical mapping $\mathrm{\Psi }$ maps the parameter domain $\widehat{\mathrm{\Omega }}$ into the physical domain $\mathrm{\Omega }$. Given a knot vector $\mathbf{\Xi }=({\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\dots ,{\xi }_s)$ with a non-decreasing sequence of values lying in parameter space, the mapping between two domains can be expressed as

(1)$\mathrm{\Psi }:\widehat{\mathrm{\Omega }}\to \mathrm{\Omega },\mathbf{\Xi }\to \mathrm{\Psi }(\mathbf{\Xi }).$

The NURBS curve is constructed by linear combinations of its basis functions, in which the coefficients are a given set of control points. A NURBS curve of p-degree is defined as

(2)$\mathrm{\Psi }(\xi )={\displaystyle \sum _{i=\mathrm{1}}^n{R}_{i,p}}(\xi )P_i,$

where $n=s-p-\mathrm{1}$ is the number of control points, $P_i\in {ℝ}^d$ is the ith control points in the physical domain, and $R_{i,p}(\xi )$ is the ith univariate NURBS basis function defined in the parameter space $\widehat{\mathrm{\Omega }}$ as follows:

(3)$R_{i,p}(\xi )={\displaystyle \frac{N_{i,p}(\xi ){\omega }_i}{{\displaystyle \sum _{i=\mathrm{1}}^n{N}_{i,p}}(\xi ){\omega }_i}},$

where ${\omega }_i\in (\mathrm{0},\mathrm{1})$ are non-decreasing weights associated with control points and $N_{i,p}(\xi )$ represents the ith B-spline basis functions of p-degree; it is defined by the following Cox-de Boor recursion formula [27].

(4)$\{\begin{array}{l}{N}_{i,\mathrm{0}}(\xi )=\{\begin{array}{lll}\mathrm{1}& & \mathrm{if} {\xi }_i\le \xi <{\xi }_{i+\mathrm{1}}\\ \mathrm{0}& & \mathrm{otherwise}\end{array}\\ N_{i,p}(\xi )={\displaystyle \frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}}{N}_{i,p-\mathrm{1}}(\xi )+{\displaystyle \frac{{\xi }_{i+p+\mathrm{1}}-\xi }{{\xi }_{i+p+\mathrm{1}}-{\xi }_{i+\mathrm{1}}}}{N}_{i+\mathrm{1},p-\mathrm{1}}(\xi )\end{array},$

where basis function $N_{i,p}(\xi )$ has its own support domain $[{\xi }_i,\dots ,{\xi }_{i+p+\mathrm{1}}]$ in which $N_{i,p}(\xi )$ is non-zero. A knot vector is deemed open when the knots are repeated $p+\mathrm{1}$ times at the ends of the vector. In IGA, the open-knot vector is generally used to satisfy the Kronecker delta property at boundary points.

A NURBS surface is a tensor product of bivariate NURBS curves in $\mathbf{\Xi }$ and $\mathrm{H}$ directions with p- and p-degrees, respectively, where a knot vector $\mathrm{H}=({\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\dots ,{\eta }_t)$ is given in $\mathrm{H}$ direction.

(5)$\mathrm{\Psi }(\xi ,\eta )={\displaystyle \sum _{i=\mathrm{1}}^n{\displaystyle \sum _{j=\mathrm{1}}^m{R}_{i,j}^{p,q}}}(\xi ,\eta )P_{i,j},$

where $m=t-q-\mathrm{1}$, $P_{i,j}$ are control points and $R_{i,j}^{p,q}$ are bivariate basis functions of the form:

(6)$R_{i,j}^{p,q}={\displaystyle \frac{N_{i,p}(\xi )N_{j,q}(\eta ){\omega }_{i,j}}{{\displaystyle \sum _{i=\mathrm{1}}^n{\displaystyle \sum _{j=\mathrm{1}}^m{N}_{i,p}}}(\xi )N_{j,q}(\eta ){\omega }_{i,j}}},$

where $N_{i,p}(\xi )$ and $N_{j,q}(\eta )$ are B-spline basis functions defined on the knot vector $\mathbf{\Xi }=({\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\dots ,{\xi }_s)$ and$\mathrm{H}=({\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\dots ,{\eta }_t)$, respectively. The interval $\mathbf{\Xi }\times \mathrm{H}$ forms a patch containing all elements, namely, $[{\xi }_k,{\xi }_{k+\mathrm{1}}]\times [{\eta }_l,{\eta }_{l+\mathrm{1}}]$, $\mathrm{1}\le k\le n+p$, and $\mathrm{1}\le l\le m+q$, which are defined by the two knot vectors. The parameter domain and corresponding physical domain for a surface model are depicted in Fig. 1.

Fig.1 Geometrical mapping $\mathrm{\Psi }$maps the common parameter space $(\xi ,\eta )$ onto the physical space

Full size|PPT slide

On the basis of the isoparametric concept, the IGA approach utilizes the same parameters for geometry and analysis models, and the basis functions used for geometry representation are also employed to approximate the numerical solution of PDEs. With Eq. (2), the numerical solution $u$ can be expressed as

(7)$u={\displaystyle \sum _{i=\mathrm{1}}^n{R}_i}(\xi )u_i,$

where $R_i$ is the ith basis function. $u_i$, which is referred to as a control variable at the ith control point, is the coefficient used to approximate the field variable $u$, which plays the same role as the nodal value in FEA. For each element, the shape function and strain-displacement matrix can be expressed as

(8)$\mathit{R}=[R_{\mathrm{1}}{R}_{\mathrm{2}}\mathrm{...}{R}_n],\mathit{B}=\left[\begin{array}{lllll}{\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial x}}& \mathrm{0}& \mathrm{...}& {\displaystyle \frac{\partial R_n}{\partial x}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial y}}& \mathrm{...}& \mathrm{0}& {\displaystyle \frac{\partial R_n}{\partial y}}\\ {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial y}}& {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial x}}& \mathrm{...}& {\displaystyle \frac{\partial R_n}{\partial y}}& {\displaystyle \frac{\partial R_n}{\partial x}}\end{array}\right].$

The strain matrix is given by

(9)$\mathit{\epsilon }=\mathit{B}\mathit{u}.$

NURBS-based parameterized LSM

LSM is a TO method that implicitly defines the interfaces between material and void phases by iso-contours of a high-dimensional LSF. Thus, this method allows a crisp description of the material boundary and helps avoid mesh-dependent problems.

The shape of the interpolating functions of LSF directly influences the smoothness of LSF and the material domain. In its most general form, LSF is described by a summation of interpolating functions scaled by their degree of freedom (DOF).

(10)$\mathrm{\Phi }(\mathit{x})={\varphi }^{\mathrm{T}}(\mathit{x})\mathit{s}={\displaystyle \sum _i^N{\varphi }_i(\mathit{x})s_i},$

where $\mathit{x}$ denotes the spatial coordinate, ${\varphi }_i$ comprises interpolating functions associated with $N$spatial points, and $s_i$ are time-dependent optimization variables.

The most commonly used interpolation functions in present LSM are FEM shape functions and RBFs, and their corresponding optimization variables are nodal values and expansion coefficients, respectively. We introduce NURBS basis functions, which can be used to approximate a given set of points with smooth polynomial functions, for parameterizing LSF. Thus, the NURBS-based parameterized LSF is constructed as

(11)$\mathrm{\Phi }(x)={\mathit{R}}^{\mathrm{T}}\mathit{s}={\displaystyle \sum _i^N{R}_i(\mathit{x})s_i},$

where $R_i$ is the ith NURBS basis function and $s_i$ is the ith time-dependent expansion coefficient related to the ith control point.

The evolution of LSF is governed by the following H-J equation.

(12)${\displaystyle \frac{\partial \mathrm{\Phi }(\mathit{x},t)}{\partial t}}+\nabla \mathrm{\Phi }{\displaystyle \frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}}=\mathrm{0},$

where $t$ denotes pseudo-time, which represents the evolution of the design in the optimization process. The speed of movement of a point on the level set surface can be expressed by $V=\mathrm{d}\mathit{x}/\mathrm{d}t$. ${\mathit{V}}_{\mathrm{n}}=V\cdot \mathit{n}$ defines the speed of propagation of all level sets along the external normal velocity, where $\mathit{n}=-\nabla \mathrm{\Phi }/|\nabla \mathrm{\Phi }|$. Therefore, Eq. (12) can be rewritten as

(13)${\displaystyle \frac{\partial \mathrm{\Phi }(\mathit{x},t)}{\partial t}}={\mathit{V}}_{\mathrm{n}}|\nabla \mathrm{\Phi }|.$

By substituting Eq. (11) into Eq. (13), the H-J equation can be written as

(14)${\mathit{R}}^{\mathrm{T}}{\displaystyle \frac{\partial \mathit{s}}{\partial t}}+{\mathit{V}}_{\mathrm{n}}\cdot \left|{(\nabla \mathit{R})}^{\mathrm{T}}\mathit{s}\right|=\mathrm{0.}$

The moving speed of the material free boundary during evolution is related to the time derivative of the expansion coefficient as follows:

(15)${\mathit{V}}_{\mathrm{n}}=-{\displaystyle \frac{{\mathit{R}}^{\mathrm{T}}}{\left|{(\nabla \mathit{R})}^{\mathrm{T}}\mathit{s}\right|}}{\displaystyle \frac{\partial \mathit{s}}{\partial t}}.$

Optimization problems of maximizing eigenvalue

Definition of the eigenvalue problem

We describe the eigenvalue problems in the linear elastic continuum to facilitate the computation of vibration frequencies and modes. A linear elastic continuum structure with a constant mass density is defined in domain $\mathrm{\Omega }\subset {ℝ}^d(d=\mathrm{2}\mathrm{o}\mathrm{r}\mathrm{3})$ with the boundary $\mathrm{\Gamma }=\partial \mathrm{\Omega }$. The weak formulation of the undamped free vibration problem can be expressed as

(16)$a(\mathit{u},\mathit{v})-\lambda b(\mathit{u},\mathit{v})=\mathrm{0},\mathit{v}\in U,$

where eigenfrequency $\lambda$ and corresponding eigenvector $\mathit{u}$, that is, the displacement subdomain in $\mathrm{\Omega }$, are solutions of this eigenvalue problem, $\mathit{v}$ is adjoint displacement, which satisfies the kinematic boundary condition, and $U$ is a space of kinematically admissible displacement fields. In LSM, $a(\cdot ,\cdot )$ and $b(\cdot ,\cdot )$ are respectively defined as

(17)$a(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}}{\epsilon }_{kl}(\mathit{u}){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega },$

(18)$b(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}\rho }\mathit{u}\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega },$

where $D_{ijkl}$ stands for the elasticity tensor component, ${\epsilon }_{ij}$ is the strain tensor component, $\rho$ is the density of the material, and $H(\mathit{\Phi })$ is the Heaviside function, which takes $\mathrm{0}$ for $\mathit{\Phi }<\mathrm{0}$ and $\mathrm{1}$ otherwise.

The eigenvalue problem has a family of solutions ${\lambda }_k$ and $u_k$, $k\ge \mathrm{1}$. The first eigenfrequency and its eigenvector are related to each other as

(19)${\lambda }_{\mathrm{1}}=\min {\displaystyle \frac{{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}}{\epsilon }_{kl}(\mathit{u}){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}{{\displaystyle {\int }_{\mathrm{\Omega }}\rho }\mathit{u}\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}}.$

Optimization model

We consider the TO problem by maximizing the first eigenfrequency under a volume constraint. Under the NURBS-based level set framework, eigenfrequency TO can be expressed as

(20)$\begin{array}{l}\mathrm{Maximize}\mathrm{\hspace{1em}\hspace{1em}}J(\mathit{u},\mathrm{\Phi })={\lambda }_{\mathrm{1}}\\ \mathrm{subject}\hspace{0.17em}\mathrm{to}:\hspace{1em}\hspace{0.17em}a(\mathit{u},\mathit{v},\mathrm{\Phi })={\lambda }_{\mathrm{1}}b(\mathit{u},\mathit{v},\mathrm{\Phi }),\\ {\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }\le V_{\max },\\ s_{\min }\le s_i\le s_{\max },\mathrm{\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}} \end{array}$

where $J(\mathit{u},\mathrm{\Phi })$ is the objective function, $V_{\max }$ represents the maximum admissible volume of the design domain, and $s_{\min }$ and $s_{\max }$ stand for the lower and upper bounds of the design variables, respectively.

However, in the eigenfrequency optimization process, the value of higher-order eigenfrequency may decrease whereas that of lower-order target eigenfrequency may increase, which may possibly lead to the repetition and exchange of mode order number. Given that the objective and constraint functions are typically defined based on a fixed modal order, the sensitivities of these functions are discontinuous in the repeated eigenfrequency case. Approaches are often used to maintain the simplicity of eigenfrequency during the entire optimization process and overcome this ill-posed problem. The modal assurance criterion (MAC) method, an efficient and accurate strategy, is introduced to monitor a single target mode, which is the first mode in this study. The definition of MAC is

(21)$\mathrm{MAC}({\mathit{u}}_a,{\mathit{u}}_b)={\displaystyle \frac{|{\mathit{u}}_a^{\mathrm{T}}{\mathit{u}}_b{|}^{\mathrm{2}}}{({\mathit{u}}_a^{\mathrm{T}}{\mathit{u}}_a)({\mathit{u}}_b^{\mathrm{T}}{\mathit{u}}_b)}},$

where ${\mathit{u}}_a$ and ${\mathit{u}}_b$ represent two eigenvectors: One is the reference eigenvector of the current optimization cycle and the other is the objective eigenvector of the previous cycle in the eigenvalue optimization process. The value of MAC varies between 0 and 1. Theoretically, a MAC value of 1 means that the two eigenvectors representing modal shapes are exactly the same. However, this condition is impossible because the structural configuration changes in each iteration, and the modal shapes of adjacent iterations are not orthogonal to each other. By comparing a few reference eigenvectors with an objective eigenvector ${\mathit{u}}_{\mathrm{obj}}^{}$, a new objective eigenvector is obtained in each iteration step of the optimization process, and it can be expressed as

(22)${\mathit{u}}_{\mathrm{obj}}^n={\mathit{u}}_k^n\hspace{1em}\mathrm{that} \underset{k}{\max }[\mathrm{MAC}({\mathit{u}}_{\mathrm{obj}}^{n-\mathrm{1}},{\mathit{u}}_k^n)], k=\mathrm{1},\mathrm{2},\dots ,N_m,$

where $N_m$ is the number of modal eigenvectors that need to be checked, superscript of displacement indicates the number of iteration step, that is, in the nth iteration, objective eigenvector of previous iteration step ${\mathit{u}}_{\mathrm{obj}}^{n-\mathrm{1}}$ is used to calculate the MAC value.

Sensitivity analysis

Establishing the relationship between the optimization function and design variables by using a sensitivity analysis approach is necessary to solve the optimization problem. According to the material derivative and the adjoint method, the Lagrangian function can be defined as

(23)$L(\mathit{u},\mathrm{\Phi })=J(\mathit{u},\mathrm{\Phi })+a(\mathit{u},\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b(\mathit{u},\mathit{v},\mathrm{\Phi })+\mathrm{\Lambda }\left[{\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }-V_{\max }\right],$

where $\mathrm{\Lambda }$ is the Lagrangian multiplier. Assuming that $V(\mathit{u},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }-V_{\max }$ is the volume constraint, the shape derivative of Lagrangian function $L(\mathit{u},\mathrm{\Phi })$ is

(24)${\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}={{\lambda }_{\mathrm{1}}}^{\prime }+a^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}{b}^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })-{{\lambda }_{\mathrm{1}}}^{\prime }b(\mathit{u},\mathit{v},\mathrm{\Phi })+\mathrm{\Lambda }{V}^{\prime }(\mathit{u},\mathrm{\Phi }) ,$

where the material derivatives of $a(\mathit{u},\mathit{v},\mathrm{\Phi })$ and $b(\mathit{u},\mathit{v},\mathrm{\Phi })$ are respectively given by

(25)$\begin{array}{l}{a}^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}\\ +{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}({\mathit{v}}^{\prime })H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}\\ +{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}(\mathit{v})\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega ，}}\end{array}$

(26)$b^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}\rho {\mathit{u}}^{\prime }\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}+{\displaystyle {\int }_{\mathrm{\Omega }}\rho \mathit{u}{\mathit{v}}^{\prime }H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}+{\displaystyle {\int }_{\mathrm{\Omega }}\rho \mathit{u}\mathit{v}\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega }},$

where ${\mathit{u}}^{\prime }$ and ${\mathit{v}}^{\prime }$ are partial derivatives of $\mathit{u}$ and $\mathit{v}$, respectively, with respect to pseudo-time. $\delta (x)$ is the Dirac function.

The adjoint state equation can be obtained by the Kuhn-Tucker condition.

(27)$a({\mathit{u}}^{\prime },\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b({\mathit{u}}^{\prime },\mathit{v},\mathrm{\Phi })=\mathrm{0},$

(28)$a(\mathit{u},{\mathit{v}}^{\prime },\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b(\mathit{u},{\mathit{v}}^{\prime },\mathrm{\Phi })=\mathrm{0},$

(29)$\mathrm{1}-b(\mathit{u},\mathit{v},\mathrm{\Phi })=\mathrm{0.}$

Given that the real mode $\mathit{u}$ is equal to the adjoint mode $\mathit{v}$, Eq. (29) is a normalization condition for the eigenvector. In this case, Eq. (24) can be simplified as

(30)${\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}={\displaystyle {\int }_{\mathrm{\Omega }}(F(\mathit{u})+\mathrm{\Lambda })\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega }},$

where $F(\mathit{u})=D_{ijkl}\mathit{\epsilon }(\mathit{u})\mathit{\epsilon }(\mathit{u})-{\lambda }_{\mathrm{1}}\rho \mathit{u}\mathit{u}$. By substituting Eqs. (11) and (15) into the preceding shape derivative equation, we obtain

(31)${\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}=-{\displaystyle {\int }_{\mathrm{\Omega }}(F(\mathit{u})+\Lambda )\delta (\mathrm{\Phi })\mathit{R}{\displaystyle \frac{\partial s}{\partial t}}\mathrm{d}\mathrm{\Omega }}.$

Given that

(32)${\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial s}}={\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s}}+\mathrm{\Lambda }{\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s}}={\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial s}}{\displaystyle \frac{\partial s}{\partial t}},$

the sensitivity of the objective function and volume constraint with respect to the design variables is respectively obtained as follows:

(33)${\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s}}=-{\displaystyle {\int }_{\mathrm{\Omega }}F(\mathit{u})\delta (\mathrm{\Phi })\mathit{R}\mathrm{d}\mathrm{\Omega }},$

(34)${\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s}}=-{\displaystyle {\int }_{\mathrm{\Omega }}\delta (\mathrm{\Phi })\mathit{R}\mathrm{d}\mathrm{\Omega }}.$

Numerical implementation

Many design variables, which correspond to large-scale nonlinear equations in the eigenvalue problem, exist in continuum structural TO. Thus, OC is introduced to solve this eigenvalue TO problem. By properly iterating and updating the design variables, this optimization problem is guaranteed to converge to a final solution. Starting from an initial value, the iterative formula of the design variables is expressed as

(35)$s_i^{(k+\mathrm{1})}=c_i^{(k)}{s}_i^{(k)}.$

Theoretically, the iteration coefficients $c_i^{(k)}$ are obtained by setting Eq. (32) equal to 0, which can be written as

(36)$c_i^{(k)}=-{\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s_i^{(k)}}}/\max \left\{\mu ,{\mathrm{\Lambda }}^{(k)}{\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s_i^{(k)}}}\right\},$

where $\mu$ is a very small number that can avoid singularity when the sensitivity of the volume constraint with respect to the design variables is equal to 0. The Lagrangian multiplier $\mathrm{\Lambda }$ is calculated by the bisection method [76].

A flowchart of the structural TO for maximization of the first eigenfrequency problem is depicted in Fig. 2. Given the condition of constraints, when the relative difference value of the objective function in the current and previous iterations is less than 10⁻³, this optimization process is considered convergent, and the current optimization process is terminated.

Fig.2 Flowchart of the optimization procedure

Full size|PPT slide

Numerical examples

In this section, the proposed IGA-based level set TO framework is applied to two 2D optimization problems. For all examples, the properties of the isotropic material are set as follows: Young’s modulus $E=\mathrm{210}\mathrm{GPa}$and mass density $\rho =\mathrm{7.8}\times {\mathrm{10}}^{\mathrm{3}}{\mathrm{kg/m}}^3$. The properties of the artificial weak material are $E_{\mathrm{0}}=\mathrm{210}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{GPa}$ and mass density $\rho =\mathrm{7.8}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{kg/m}}^3$. Poisson’s ratio $\nu =\mathrm{0.3}$ and plane stress state are assumed for all the materials. The two examples are TO of maximizing the fundamental eigenfrequency of the plane structure with a unit thickness of $\mathrm{0.001}\mathrm{m}$ and a prescribed material volume fraction of $\alpha =\mathrm{50}\%$. In the initial design, the available material is uniformly distributed over the entire admissible design domain. In the following examples, the boundary conditions are imposed by the collocation method, which enforces these conditions to be satisfied pointwise. Given that NURBS basis functions associated with the interior control points vanish at the structural boundary when open-knot vectors are employed, the displacement boundary condition applied on the left and right sides of the beam is imposed by setting the displacement values at left and right boundary control points to zero. All results are produced with programs developed in the MATLAB R2018a environment on a computer with an Intel Core i3-3240 CPU, 3.4 GHz clock speed, and 6 GB RAM. Additional details on the implementation of IGA on the MATLAB platform were presented in Ref. [77].

Cantilever beam

The first numerical example of a short cantilever beam for maximizing the first eigenfrequency optimization problem is shown in Fig. 3. The entire design domain is a rectangle with a size of 0.2 m×0.1 m, a Dirichlet boundary, and fixed displacement at the left edge of the design domain. A concentrated nonstructural mass $M=\mathrm{15.6}$ kg that is one-tenth of the total structural mass of the plate is placed at the center of the right side. Notably, the structure disappears without the nonstructural mass because no structure leads to the highest eigenfrequencies.

Fig.3 Design domain of a cantilever beam structure

Full size|PPT slide

In this example, IGA- and FEA-based TO approaches are applied to a similar problem for comparison. For this 2D example, equally spaced open-knot vectors are used for $x$ and $y$ directions, and the degrees of NURBS basis functions on the two directions are the same, i.e., $p=q=\mathrm{2}$. Nine-node quadratic rectangle elements are used for FEA, and the element number of both methods is the same, which facilitates a comparison under identical conditions.

The computation consumption times of the two methods are shown in Table 1. Notably, for simplicity in the proceeding table and figures, IGA-L and FEM-L represent IGA- and FEM-based LSM, respectively. The solution time of the system equation includes the time spent on assembling the stiffness matrix and solving the equation to obtain an objective function. In this case, when the numbers of Lagrange and NURBS elements are $a\times b$ ($a$ and $b$ represent the number of elements in $x$ and $y$ directions, respectively), the number of DOFs of the IGA-based method is $(a+\mathrm{2})(b+\mathrm{2})=ab+\mathrm{2}a+\mathrm{2}b+\mathrm{4}$ and that of FEM is $(\mathrm{2}a+\mathrm{1})(\mathrm{2}b+\mathrm{1})=\mathrm{4}ab+\mathrm{2}a+\mathrm{2}b+\mathrm{1}$. The DOF of IGA is much less than that of FEM (nearly a quarter of FEM) when the element number is sufficiently large. Table 1 shows that the computation efficiency of IGA-based LSM is higher than that of FEM-based LSM in TO due to the fewer DOFs or smaller size of equations in the IGA-based approach. However, because of the extra calculations of the basis functions and their derivative in the IGA-based approach, the ratio between the iteration time and DOFs of the IGA-based optimization method is higher than that of the FEM-based method.

Tab.1 Comparison of IGA- and FEA-based LSM TO

Method	Number of elements	Number of DOFs	Time of each iteration/s	Time of solution of system equation/s
IGA-L	256×128	67080	873.26	713.68
FEA-L	256×128	263682	1580.42	1237.16
IGA-L	128×64	17160	30.06	25.31
FEA-L	128×64	66306	54.17	43.97
IGA-L	64×32	4488	6.41	5.18
FEA-L	64×32	16770	8.55	6.85

The optimal layouts obtained by using IGA and FEA with 128×64 elements are shown in Fig. 4. In this case, the results are similar, and the crisp boundary is obtained due to the level set method. Adopting B-spline basis functions to parameterize the level set function together with IGA for calculation instead of FEA leads to the rapid convergence of the IGA-based optimization process.

Fig.4 Optimized results obtained by using IGA and FEA methods with 128 × 64 elements

Full size|PPT slide

The convergence history of optimization using IGA and FEA with 128 × 64 elements is shown in Figs. 5 and 6. Figure 5 shows the convergence history of the first eigenfrequency and volume ratio of the structure by using IGA- and FEA-based optimization frameworks, respectively. The initial designs and resultant structures of both methods are basically the same. In optimization with the IGA method, the ${\omega }_i$ of the initial design and the resultant optimum are 173.5 and 188.7, respectively. The volume ratio of the initial design and the resultant optimum are 0.79 and 0.5, respectively. Thus, the fundamental eigenfrequency increased by 8.8% and the volume decreased by 36.7%. Figure 6 shows the iteration history of the first three eigenfrequencies by using both methods. The first eigenfrequency always remains simple, whereas the second and third eigenfrequencies tend to oscillate and overlap. By using the MAC method, problems regarding the orders of eigenfrequency exchange during the optimization process are avoided. The first eigenfrequency generally increases, and the second and third eigenfrequencies decrease as the volume ratio decreases.

Fig.5 Comparison of IGA and FEA in terms of convergence history

Full size|PPT slide

Fig.6 Comparison of IGA and FEA in terms of eigenfrequency history

Full size|PPT slide

Beam with clamped ends

In this section, we present an example of maximizing the fundamental frequency of a clamped beam structure shown in Fig. 7. The working domain has a size of 0.4 m×0.1 m. A fixed displacement boundary condition is imposed on both sides, and a concentrated nonstructural mass M=31.2 kg is placed at the center of structure. In this example, the capability of the proposed method to capture the optimum topology and the effect of the number of elements are studied. For this purpose, a clamped beam is solved with three mesh sizes, namely, 64×16, 128×32, and 256×64. The degree of NURBS basis function is 2 in both directions. The resulting layouts shown in Fig. 8 indicate that the mesh dependency problem is avoided in the proposed method. This figure also shows that an inaccurate optimum topology with a rough boundary is obtained as a consequence of coarse meshes. By refining the mesh, the smoothness of the boundaries of the optimal layout is improved. However, when the number of elements is larger than 128×32, the resulting layout slightly changes. Accounting for the computation cost and precision of results, 128×32 meshes are used in this example.

Fig.7 Design domain of a clamped beam structure

Full size|PPT slide

Fig.8 Optimal layouts obtained by using (a) 64×16 meshes, (b) 128×32 meshes, and (c) 256×64 meshes

Full size|PPT slide

The convergence history of the objective function and volume ratio is given in Fig. 9. The first eigenfrequency of the initial design and the resultant topology are 244.3 and 257.4, respectively. The volume ratio of the initial design and the resultant topology are 0.82 and 0.5, respectively. Thus, the fundamental frequency increases by 5.4%, and the volume decreases by 39%. Figure 10 shows the convergence history of the first three eigenfrequencies. The fundamental frequency remains simple throughout the entire optimization process. The value of first eigenfrequency gradually increases after the fifth iteration, whereas that of the second and third eigenfrequencies decline with the increase in iterations.

Fig.9 Convergence history

Full size|PPT slide

Fig.10 Iteration history of the first three eigenfrequencies

Full size|PPT slide

Conclusions

We solved maximum fundamental eigenfrequency TO problems with a level set model based on the IGA technique. IGA combines the fundamental idea of FEM with the spline technique from a computer-aided geometry design for the integration of CAD and CAE. The IGA method was also introduced to TO due to its superiority over currently used FEM in terms of accuracy and efficiency. The feature of the proposed method is the combination of IGA and LSM in eigenfrequency optimization where the same basis functions (NURBS) are used for geometry representation, dynamic analysis, and parameterization of the implicit LSF. High accuracy and smoothness of LSF were achieved by using smooth NURBS basis functions to approximate LSF. In the case of maximizing the fundamental eigenfrequency, a regularization method for the repetition or exchange of eigenfrequencies was employed to guarantee the simple behavior of structural eigenfrequency.

Two benchmark numerical examples of TO for dynamic problems were applied to verify the validity of the proposed approach. The results obtained from the comparison of FEA- and IGA-based level set TO methods demonstrated that the proposed IGA-based optimization method has better computational efficiency and converges faster than the traditional FEA-based optimization method. The results also showed that solving dynamic TO problems by using IGA together with the level set method is possible. Although we only presented examples of 2D structures, no theoretical difficulties will be encountered if the proposed is extended to the optimization of 3D structures.

References

[1] Chen X F , Wang S B , Qiao B J , Chen Q . Basic research on machinery fault diagnostics: past, present, and future trends. Frontiers of Mechanical Engineering, 2018, 13(2): 264–291

[2] Zhu Y P , Wang Y , Luo Z , Han Q K , Wang D Y . Similitude design for the vibration problems of plates and shells: a review. Frontiers of Mechanical Engineering, 2017, 12(2): 253–264

[3] Ma H , Xie F T , Nai H Q , Wen B C . Vibration characteristics analysis of rotating shrouded blades with impacts. Journal of Sound and Vibration, 2016, 378: 92–108

[4] Gao Y , Liu X Y , Xiang J W . FEM simulation-based generative adversarial networks to detect bearing faults. IEEE Transactions on Industrial Informatics, 2020, 16(7): 4961–4971

[5] Yang F M , Yang Y M , Hu H F , Guan F J , Shen G J , Chen S Y , Bian Z F . Extraction of features due to breathing crack from vibration response of rotated blades considering tenon connection and shroud contact. Shock and Vibration, 2019, 2019: 8729620

[6] Guo H T, Duan F J, Zhang J L. Blade resonance parameter identification based on tip-timing method without the once-per revolution sensor. Mechanical Systems and Signal Processing, 2016, 66–67: 625–639

[7] Du Z H , Chen X F , Zhang H , Zi Y Y , Yan R Q . Multiple fault separation and detection by joint subspace learning for the health assessment of wind turbine gearboxes. Frontiers of Mechanical Engineering, 2017, 12(3): 333–347

[8] Bouchain A , Picheral J , Lahalle E , Chardon G , Vercoutter A , Talon A . Blade vibration study by spectral analysis of tip-timing signals with OMP algorithm. Mechanical Systems and Signal Processing, 2019, 130: 108–121

[9] Xu J H , Qiao B J , Liu M R , Yang Z B , Chen X F . Crack propagation monitoring of rotor blades using synchroextracting transform. Journal of Sound and Vibration, 2021, 509: 116253

[10] Battiato G , Firrone C M , Berruti T M . Forced response of rotating bladed disks: blade tip-timing measurements. Mechanical Systems and Signal Processing, 2017, 85: 912–926

[11] Poozesh P, Baqersad J, Niezrecki C, Avitabile P, Harvey E, Yarala R. Large-area photogrammetry based testing of wind turbine blades. Mechanical Systems and Signal Processing, 2017, 86(Part B): 98–115

[12] Qiao B J , Ao C Y , Mao Z , Chen X F . Non-convex sparse regularization for impact force identification. Journal of Sound and Vibration, 2020, 477: 115311

[13] Huh Y H , Kim J , Hong S G . Detection of local bonding failure damage by digital image correlation technique. Reliability Engineering & System Safety, 2019, 184: 21–26

[14] Khadka A , Fick B , Afshar A , Tavakoli M , Baqersad J . Non-contact vibration monitoring of rotating wind turbines using a semi-autonomous UAV. Mechanical Systems and Signal Processing, 2020, 138: 106446

[15] Retze U, Schüssler M. Dynamic Stress and Strain Measurement. Polytec Technical Papers, 2010

[16] Sever I A , Maguire M . Correlation of full-field dynamic strain measurements with reverse engineered finite element model predictions. Experimental Techniques, 2021, 45(3): 377–387

[17] Liu C H , Zang C P , Zhou B . Extension of continuous scanning laser Doppler vibrometry measurement for complex structures with curved surfaces. Chinese Journal of Aeronautics, 2020, 33(12): 3220–3227

[18] Chen Y C , Escalera Mendoza A S , Griffith D T . Experimental and numerical study of high-order complex curvature mode shape and mode coupling on a three-bladed wind turbine assembly. Mechanical Systems and Signal Processing, 2021, 160: 107873

[19] Vuye C , Vanlanduit S , Presezniak F , Steenackers G , Guillaume P . Optical measurement of the dynamic strain field of a fan blade using a 3D scanning vibrometer. Optics and Lasers in Engineering, 2011, 49(7): 988–997

[20] Maguire M , Sever I . Full-field strain measurements on turbomachinery components using 3D SLDV technology. AIP Conference Proceedings, 2016, 1740(1): 080001

[21] Luo W , Li J H , Ma X F , Wei W . A novel static deformation measurement and visualization method for wind turbine blades using home-made LiDAR and processing program. Optics and Lasers in Engineering, 2020, 134: 106206

[22] Chen Y C , Joffre D , Avitabile P . Underwater dynamic response at limited points expanded to full-field strain response. Journal of Vibration and Acoustics, 2018, 140(5): 051016

[23] Sarrafi A , Poozesh P , Niezrecki C , Mao Z . Mode extraction on wind turbine blades via phase-based video motion estimation. Proceeding SPIE 10171, Smart Materials and Nondestructive Evaluation for Energy Systems, 2017, 10171: 101710E

[24] Baqersad J , Bharadwaj K . Strain expansion-reduction approach. Mechanical Systems and Signal Processing, 2018, 101: 156–167

[25] Bharadwaj K , Sheidaei A , Afshar A , Baqersad J . Full-field strain prediction using mode shapes measured with digital image correlation. Measurement, 2019, 139: 326–333

[26] Chen Z S , Sheng H , Xia Y M , Wang W M , He J . A comprehensive review on blade tip timing-based health monitoring: status and future. Mechanical Systems and Signal Processing, 2021, 149: 107330

[27] Zhang X J , Wang Y R , Jiang X H , Gao S M . Blade vibration stress determination method based on blade tip timing simulator and finite element method. Journal of Engineering for Gas Turbines and Power, 2020, 142(3): 031001

[28] Wang W M , Hu D F , Li Q H , Zhang X L . An improved non-contact dynamic stress measurement method for turbomachinery rotating blades based on fundamental mistuning model. Mechanical Systems and Signal Processing, 2020, 144: 106851

[29] Mohamed M , Bonello P , Russhard P . A novel method for the determination of the change in blade tip timing probe sensing position due to steady movements. Mechanical Systems and Signal Processing, 2019, 126: 686–710

[30] Wu S M , Zhao Z B , Yang Z B , Tian S H , Yang L H , Chen X F . Physical constraints fused equiangular tight frame method for blade tip timing sensor arrangement. Measurement, 2019, 145: 841–851

[31] Bornassi S , Berruti T M , Firrone C M , Battiato G . Vibration parameters identification of turbomachinery rotor blades under transient condition using blade tip-timing measurements. Measurement, 2021, 183: 109861

[32] Chakraborty A , Shishkin S , Birnkrant M J . Optimal control of build height utilizing optical profilometry in cold spray deposits. Proceeding SPIE 10168, Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems, 2017, 10168: 101683H

[33] Kochkin D Y , Zaitsev D V , Kabov O A . Thermocapillary rupture and contact line dynamics in the heated liquid layers. Interfacial Phenomena and Heat Transfer, 2020, 8(1): 1–9

[34] Kranjc T , Slavič J , Boltežar M . A comparison of strain and classic experimental modal analysis. Journal of Vibration and Control, 2016, 22(2): 371–381

[35] Rao S S , Atluri S N . The finite element method in engineering. Journal of Applied Mechanics, 1983, 50(4a): 914

[36] HumarJ. Dynamics of Structures. 3rd ed. Boca Baton: CRC Press, 2012

[37] Zhou Y D , Tao J Y . Theoretical and numerical investigation of stress mode shapes in multi-axial random fatigue. Mechanical Systems and Signal Processing, 2019, 127: 499–512

[38] Chen H , Song H W . Analysis of strain modals and the relationship between the strain modals and displacement modals. Noise and Vibration Control, 2016, 36(4): 7–13

[39] Ao C Y , Qiao B J , Chen L , Xu J H , Liu M R , Chen X F . Blade dynamic strain non-intrusive measurement using L_1/2-norm regularization and transmissibility. Measurement, 2022, 190: 110677

[40] Zhang X D. Matrix Analysis and Applications. Cambridge: Cambridge University Press, 2017

[41] Kharyton V , Bladh R . Using tiptiming and strain gauge data for the estimation of consumed life in a compressor blisk subjected to stall-induced loading. In: Proceedings of the ASME Turbo Expo 2014: Turbine Technical Conference and Exposition. Volume 7B: Structures and Dynamics. Düsseldorf: ASME, 2014, V07BT33A028

[42] Wu S M , Russhard P , Yan R Q , Tian S H , Wang S B , Zhao Z B , Chen X F . An Adaptive online blade health monitoring method: from raw data to parameters identification. IEEE Transactions on Instrumentation and Measurement, 2020, 69(5): 2581–2592

[43] Duffy K P, Provenza A J, Bakhle M, Min J B, Abdul-Aziz A. Laser displacement measurements of fan blades in resonance and flutter during the boundary layer ingesting inlet and distortion-tolerant fan test. In: Proceedings of 2018 AIAA Aerospace Sciences Meeting. Kissimmee: AIAA, 2018

[44] Aadnøy B, Looyeh R. Petroleum Rock Mechanics. Gulf Professional Publishing, 2011

[45] Dowling N E. Mechanical Behavior of Materials: Engineering Methods for Deformation, Fracture, and Fatigue. 3rd ed. London: Pearson Prentice Hall Ltd., 2007

[46] Wu Y , Liu G , Liu Z Q , Wang B . Formability and microstructure of Ti22Al24.5Nb0.5Mo rolled sheet within hot gas bulging tests at constant equivalent strain rate. Materials & Design, 2016, 108: 298–307

[47] ANSYSI. Help System: Mechanical APDL Element Reference, 19.0. 2018

[48] Kamaraj A V. Stress prediction in turbine blades under forced excitation. Thesis for the Master’s Degree. San Diego: San Diego State University, 2016

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 52075414) and the National Science and Technology Major Project, China (Grant No. 2017-V-0009).

RIGHTS & PERMISSIONS

2023 Higher Education Press