CN105260563A - Finite element prestressed mode analysis method for three-dimensional entity blade and two-dimensional axisymmetric wheel disc variable-dimensionality model of impeller structure - Google Patents

Abstract

The invention belongs to the technical field of structural dynamics, and particularly relates to a finite element prestressed mode analysis method for a three-dimensional entity blade and two-dimensional axisymmetric wheel disc variable-dimensionality model of an impeller structure. In the method, the deformation of blades with complex structures is represented by a three-dimensional displacement mode, and a blade model is subjected to freedom-degree condensation polymerization by adopting a dynamic substructure method; the deformation of a revolution body wheel disc is simplified into a two-dimensional axisymmetric mode comprising a circumferential displacement component, the degree of freedom of a wheel disc model is reduced, and the displacement coordination between the blades and the wheel disc is kept. Through the two operations, the degree of freedom of a blade and wheel disc coupled dynamic model is reduced by multiple orders of magnitude, so that the analysis efficiency of intrinsic modes (comprising a prestressed mode) of the impeller structure is improved, and the computation results are accurate. In the method, the shape of the impeller structure is not limited by period of revolution; the method is also suitable for a plurality of groups of impellers which are provided with different quantities of blades in the circumferential direction and mounted on a wheel disc, and can solve the difficulty that a complex impeller structure cannot be subjected to mode analysis by using a sector model.

Description

A kind of 3D solid blade of blade wheel structure and two-dimensional axial symmetric wheel disc become dimensional model finite element prestressed modal analysis method

Technical field

The invention belongs to the technical field such as Structural Dynamics and FEM (finite element) calculation, the 3D solid blade and the two-dimensional axial symmetric wheel disc that are specially a kind of blade wheel structure become dimensional model finite element prestressed modal analysis method.

Background technology

Blade wheel structure is extensively present in rotating machinery, as electric organ, steamer device and aeromotor etc.The important step that model analysis is Impeller Design is carried out to blade wheel structure, significant to stablizing of impeller system under guarantee High Rotation Speed.For the list group blade impeller structure with characteristic gyration period, current more employing sector models carries out model analysis, but sector models is often larger, and computing time is long, and cannot calculate the blade wheel structure not meeting characteristic gyration period.When being optimized design to blade wheel structure, need the model analysis that model carries out repeatedly be calculated, the bottleneck becoming optimal design analysis time that sector models is longer.

The inventive method proposes the prestressed modal analysis method that the 3D solid blade of blade wheel structure and two-dimensional axial symmetric wheel disc become dimensional model, its meaning is the prestressed modal analysis efficiency greatly improving impeller class model gyration period, does not reduce its model analysis precision simultaneously.The inventive method three-D displacement pattern characterizes the distortion of complex configuration blade meticulously, carries out degree of freedom polycondensation on this basis with dynamic sub-structure methods to 3D solid leaf model; The distortion of solid of revolution wheel disc is reduced to the two-dimensional axial symmetric pattern comprising circumferential displacement component, reduces the degree of freedom of wheel disc model, and keeps the displacement coordination between blade and wheel disc.By above two operations, the degree of freedom of blade and wheel disc Coupling Dynamic Model can be reduced some magnitudes, improve the efficiency that blade wheel structure natural mode of vibration (comprising pre-stressed mode) is analyzed, and modal analysis result is close with complete three-dimensional blades finite element model result of calculation.In this method, the shape of blade wheel structure is not by the restriction of gyration period, for being provided with many groups, the impeller of the different blade of circumferential number on wheel disc, also can use this method analysis, thus complicated blade wheel structure cannot carry out a model analysis difficult problem with sector models can be solved.

Summary of the invention

The 3D solid blade and the two-dimensional axial symmetric wheel disc that the object of this invention is to provide a kind of blade wheel structure become dimensional model finite element prestressed modal analysis method.The method has universality for impeller class rotationally periodic structure, it becomes dimension finite element model by setting up with the 3D solid blade of Complete three-dimensional impeller pattern equivalence and two-dimensional axial symmetric wheel disc, reduce the number of degrees of freedom of model, thus improve the model analysis efficiency of this class model.

The 3D solid blade of a kind of blade wheel structure provided by the invention and two-dimensional axial symmetric wheel disc become dimensional model finite element prestressed modal analysis method, and concrete steps comprise:

(1) set up the blade wheel structure be made up of blade and wheel disc and become dimension finite element model; First blade 3D solid unit represents by it, sets up 3D solid leaf model; Wheel disc two-dimensional axial symmetric unit is represented, sets up two-dimensional axial symmetric wheel disc model; Again by arranging one group of connected node respectively in 3D solid leaf model and two-dimensional axial symmetric wheel disc model, realize the combination that impeller becomes dimensional model, and between the partial interior node and connected node of 3D solid leaf model and two-dimensional axial symmetric wheel disc model, set up the linear restriction relation of one group of displacement component respectively, make each connected node all as static determinacy restraint joint;

(2) utilize dynamic sub-structure methods to obtain the stiffness matrix after the polycondensation of 3D solid leaf model and mass matrix, the coupling process proposed by this method, be coupled with the stiffness matrix of two-dimensional axial symmetric wheel disc model and mass matrix; The kinetic model after coupling is utilized to carry out the model analysis of change dimensional model.

In the present invention, utilize the Coupling Deformation computing method of 3D solid leaf model and two-dimensional axial symmetric wheel disc model, obtain the quiet distortion of impeller under centrifugal force preloads effect, should be the tangent stiffness matrix of GEOMETRICALLY NONLINEAR simultaneously with stiffness matrix, analyze the pre-stressed mode under impeller change dimensional model centrifugal force field.

Technical scheme of the present invention specifically describes as follows.

One, the foundation of two-dimensional axial symmetric wheel disc model in dimensional model is become

Two-dimensional axial symmetric wheel disc model defines in cylindrical coordinate Oxr θ, and three axles of this coordinate system are called x-axis (axis), and r axle (radial direction) and θ axle (circumference), the transform mode of its coordinate and standard straight angle coordinate system is

$\left\{\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right\}=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\mathrm{\θ}& -\sin \mathrm{\θ}\\ 0& sin\mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left\{\begin{array}{c}{u}_x\\ u_r\\ u_{\mathrm{\θ}}\end{array}\right\}---\left(1\right)$

Two-dimensional axial symmetric wheel disc finite element model defines in the unilateral Oxr of θ=0, is made up of two-dimensional axial symmetric unit, and rotating a circle around x-axis to obtain the on all four solid of revolution wheel disc model with three-dimensional model.

In cylindrical coordinates, be the distortion of n solid of revolution wheel disc for circumferential wave number, be called n pitch diameter displacement model.The upper axial displacement u in its optional position (x, r, θ) _x, radial displacement u _rwith circumferential displacement components u _θthe cosine displacement function of θ=0 meridian ellipse (hereinafter referred to as 0-meridian ellipse) can be defined in one group with sinusoidal displacement function be expressed as

$\left\{\begin{array}{c}{u}_x(x,r,\mathrm{\θ})\\ u_r(x,r,\mathrm{\θ})\\ u_{\mathrm{\θ}}(x,r,\mathrm{\θ})\end{array}\right\}\stackrel{\mathrm{\Δ}}{=}\left\{\begin{array}{c}{u}_x^{\left(c\right)}(x,r)\\ u_x^{\left(c\right)}(x,r)\\ u_{\mathrm{\θ}}^{\left(c\right)}(x,r)\end{array}\right\}\cos n\mathrm{\θ}\left\{\begin{array}{c}{u}_x^{\left(s\right)}(x,r)\\ u_x^{\left(s\right)}(x,r)\\ u_{\mathrm{\θ}}^{\left(s\right)}(x,r)\end{array}\right\}\sin n\mathrm{\θ}---\left(2\right)$

For zero pitch diameter displacement model, because sinn θ perseverance is zero, then only comprise cosine displacement function.

Under cylindrical coordinate, the components of strain of solid of revolution wheel disc n pitch diameter displacement are

$\mathrm{\ϵ}=\left\{\begin{array}{c}{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_r\\ {\mathrm{\ϵ}}_{\mathrm{\θ}}\\ {\mathrm{\γ}}_{xr}\\ {\mathrm{\γ}}_{r\mathrm{\θ}}\\ {\mathrm{\γ}}_{\mathrm{\θ}x}\end{array}\right\}=\left[\begin{array}{ccc}\frac{\∂}{\∂x}& 0& 0\\ 0& \frac{\∂}{\∂r}& 0\\ 0& \frac{1}{r}& \frac{1}{r}\frac{\∂}{\∂\mathrm{\θ}}\\ \frac{\∂}{\∂r}& \frac{\∂}{\∂x}& 0\\ 0& \frac{1}{r}\frac{\∂}{\∂\mathrm{\θ}}& \frac{\∂}{\∂r}-\frac{1}{r}\\ \frac{1}{r}\frac{\∂}{\∂\mathrm{\θ}}& 0& \frac{\∂}{\∂x}\end{array}\right]\left\{\begin{array}{c}{u}_x(x,r,\mathrm{\θ})\\ u_r(x,r,\mathrm{\θ})\\ u_{\mathrm{\θ}}(x,r,\mathrm{\θ})\end{array}\right\}---\left(3\right)$

Potential energy and the kinetic energy of the n pitch diameter displacement of solid of revolution wheel disc are

$\begin{array}{c}U=\frac{1}{2}\∫\∫\∫{\mathrm{\ϵ}}^T{D}_0\mathrm{\ϵ}dxdrd\mathrm{\θ}\\ E=\frac{1}{2}\mathrm{\ρ}\∫\∫\∫{\stackrel{\·}{u}}^T\stackrel{\·}{u}dxdrd\mathrm{\θ}\end{array}---\left(4\right)$

Wherein D ₀for elastic matrix, ρ is density.

Wheel disc hatch region on 0-meridian ellipse is divided into two-dimensional axial symmetric finite element grid, total n _dindividual node.In this finite element model, arrange one group of node be connected with blade, connecting joint is counted as n _j.So-called wheel disc tie point, refers to that the circular arc of this group node after x-axis revolution (being called for short revolution arc) is all connected with the corresponding node of in blade finite element model.

The front n of two-dimensional axial symmetric wheel disc model might as well be supposed _jindividual node is connected node, its node cosine motion vector with node sinusoidal displacement vector respectively there are six components, comprise three translation displacements components and three rotational displacement components.By defining linear restriction relation between wheel disc connected node partial interior adjacent thereto nodal displacement component, each connected node is made all to can be used as static determinacy restraint joint, for the coupling of leaf dish provides condition.The cosine motion vector of other nodes with sinusoidal displacement vector without particular restriction, with reference to conventional finite element modeling method.

Utilize finite element interpolation, the cosine displacement function of wheel disc and sinusoidal displacement function can be separated into

$\left\{\begin{array}{c}{u}_x^{\left(c\right)}(x,r)\\ u_r^{\left(c\right)}(x,r)\\ u_{\mathrm{\θ}}^{\left(c\right)}(x,r)\end{array}\right\}=N^{\left(D\right)}(x,r)q^{(D,c)}---\left(5\right)$

$\left\{\begin{array}{c}{u}_x^{\left(s\right)}(x,r)\\ u_r^{\left(s\right)}(x,r)\\ u_{\mathrm{\θ}}^{\left(s\right)}(x,r)\end{array}\right\}=N^{\left(D\right)}(x,r)q^{(D,s)}---\left(6\right)$

Wherein, N ^(D)(x, r) is shape function

$\begin{array}{c}{q}^{(D,c)}=\left\{\begin{array}{c}{q}_I^{(D,c)}\\ q_J^{(D,c)}\end{array}\right\},q^{(D,s)}=\left\{\begin{array}{c}{q}_I^{(D,s)}\\ q_J^{(D,s)}\end{array}\right\}\\ q_I^{(D,c)}={\mathrm{col}}_{i=n_J+1}^{n_D}\left\{q_i^{(D,c)}\right\},q_J^{(D,c)}={\mathrm{col}}_{j=1}^{n_J}\left\{q_j^{(D,c)}\right\}\\ q_I^{(D,s)}={\mathrm{col}}_{i=n_J+1}^{n_D}\left\{q_i^{(D,s)}\right\},q_J^{(D,s)}={\mathrm{col}}_{j=1}^{n_J}\left\{q_j^{(D,s)}\right\}\end{array}---\left(7\right)$

Like this, the displacement function of whole wheel disc all can by q ^(D)={ q ^{(D, c) T}, q ^{(D, s) T}} ^tbe expressed as

$\left\{\begin{array}{c}{u}_x(x,r,\mathrm{\θ})\\ u_r(x,r,\mathrm{\θ})\\ u_{\mathrm{\θ}}(x,r,\mathrm{\θ})\end{array}\right\}\stackrel{\mathrm{\Δ}}{=}\left[\begin{array}{cc}{N}^{\left(D\right)}(x,r)\cos n\mathrm{\θ}& N^{\left(D\right)}(x,r)\sin n\mathrm{\θ}\end{array}\right]\left\{\begin{array}{c}{q}^{(D,c)}\\ q^{(D,s)}\end{array}\right\}---\left(8\right)$

For zero pitch diameter displacement model, then motion vector q ^(D)only include q ^{(D, c)}.For dissimilar two-dimensional axial symmetric unit, formula (8) is substituted into potential energy and kinetic energy expression (4), element stiffness matrix and the element mass matrix of two-dimensional axial symmetric unit n pitch diameter displacement model can be obtained.

Element stiffness matrix is

$K_e=\left\{\begin{array}{cc}2\mathrm{\π}\∫\∫B^{\left(c\right)T}{D}_0{B}^{\left(c\right)}rdrdx& (n=0)\\ \mathrm{\π}\∫\∫B^{T}DBrdrdx& (n>0)\end{array}\right.---\left(9\right)$

$B=\left[\begin{array}{cc}{B}^{\left(c\right)}& -B^{\left(s\right)}\\ B^{\left(s\right)T}& B^{\left(c\right)}\end{array}\right],D=\left[\begin{array}{cc}{D}_0& 0\\ 0& D_0\end{array}\right]---\left(10\right)$

$\begin{array}{c}{B}^{\left(c\right)}=\[\begin{array}{cccc}{B}_1^{\left(c\right)}& B_2^{\left(c\right)}& ...& B_{n_P}^{\left(c\right)}\end{array}\]\\ B^{\left(s\right)}=\[\begin{array}{cccc}{B}_1^{\left(s\right)}& B_2^{\left(s\right)}& ...& B_{n_P}^{\left(s\right)}\end{array}\]\end{array}---\left(11\right)$

$B_i^{\left(c\right)}=\left[\begin{array}{ccc}\frac{\∂N_i}{\∂x}& 0& 0\\ 0& \frac{\∂N_i}{\∂x}& 0\\ 0& \frac{N_i}{r}& 0\\ \frac{\∂N_i}{\∂r}& \frac{\∂N_i}{\∂x}& 0\\ 0& 0& \frac{\∂N_i}{\∂r}-\frac{N_i}{r}\\ 0& 0& \frac{\∂N_i}{\∂x}\end{array}\right],B_i^{\left(s\right)}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& -\frac{n}{r}{N}_i\\ 0& 0& 0\\ 0& -\frac{n}{r}{N}_i& 0\\ -\frac{n}{r}{N}_i& 0& 0\end{array}\right]---\left(12\right)$

Wherein D ₀for elastic matrix, n _pfor cell node number.

Element mass matrix is

$M_e=\left\{\begin{array}{cc}2\mathrm{\π}\∫\∫N^{T}Nrdrdx& (n=0)\\ \mathrm{\π}\∫\∫\left[\begin{array}{cc}{N}^{T}N& 0\\ 0& N^{T}N\end{array}\right]rdrdx& (n>0)\end{array}\right.---\left(13\right)$

THE TANGENTIAL STIFFNESS MATRICES is derived as follows:

Cylindrical coordinates intrinsic displacement gradient formula is such as formula shown in (14)

${\mathrm{\θ}}_x=\frac{\∂}{\∂x}\left[\begin{array}{c}{u}_x\\ u_r\\ u_{\mathrm{\θ}}\end{array}\right],{\mathrm{\θ}}_r=\frac{\∂}{\∂r}\left[\begin{array}{c}{u}_x\\ u_r\\ u_{\mathrm{\θ}}\end{array}\right],{\mathrm{\θ}}_{\mathrm{\θ}}=\frac{1}{r}\left[\begin{array}{c}\frac{\∂u_x}{\∂\mathrm{\θ}}\\ \frac{\∂u_r}{\∂\mathrm{\θ}}-u_{\mathrm{\θ}}\\ \frac{\∂u_{\mathrm{\θ}}}{\∂\mathrm{\θ}}+u_r\end{array}\right]---\left(14\right)$

Formula (8) is substituted into formula (14) can obtain

${\mathrm{\θ}}_x=\[\begin{array}{cc}\frac{\∂N}{\∂x}\cos n\mathrm{\θ}& \frac{\∂N}{\∂x}\sin n\mathrm{\θ}\end{array}\]\left\{\begin{array}{c}{q}^{(D,c)}\\ q^{(D,s)}\end{array}\right\}---\left(15\right)$

${\mathrm{\θ}}_r=\[\begin{array}{cc}\frac{\∂N}{\∂r}\cos n\mathrm{\θ}& \frac{\∂N}{\∂r}\sin n\mathrm{\θ}\end{array}\]\left\{\begin{array}{c}{q}^{(D,c)}\\ q^{(D,s)}\end{array}\right\}---\left(16\right)$

${\mathrm{\θ}}_{\mathrm{\θ}}=\frac{1}{r}\[\begin{array}{cc}-nN\sin n\mathrm{\θ}+NS\cos n\mathrm{\θ}& nN\cos n\mathrm{\θ}+NS\sin n\mathrm{\θ}\end{array}\]\left\{\begin{array}{c}{q}^{(D,c)}\\ q^{(D,s)}\end{array}\right\}---\left(17\right)$

Wherein

$S_{(3*nP)*(3*nP)}=diag\left(\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -1\\ 0& 1& 0\end{array}\right]\right)---\left(18\right)$

Formula (15), formula (16) and formula (17) is utilized to obtain

$\left[\begin{array}{c}{\mathrm{\θ}}_x\\ {\mathrm{\θ}}_r\\ {\mathrm{\θ}}_{\mathrm{\θ}}\end{array}\right]\stackrel{\mathrm{\Δ}}{=}G\left\{\begin{array}{c}{q}^{(D,c)}\\ q^{(D,s)}\end{array}\right\}=(G^{\left(c\right)}\cos n\mathrm{\θ}+G^{\left(s\right)}\sin n\mathrm{\θ})\left\{\begin{array}{c}{q}^{(D,c)}\\ q^{(D,s)}\end{array}\right\}---\left(19\right)$

Wherein G matrix can be write as cosine component G ^(c)with sinusoidal component G ^(s)form.

Definition matrix M _σfor

$M_{\mathrm{\σ}}=\left[\begin{array}{ccc}{\mathrm{\σ}}_x{I}_3& {\mathrm{\τ}}_{xr}{I}_3& {\mathrm{\τ}}_{x\mathrm{\θ}}{I}_3\\ {\mathrm{\τ}}_{xr}{I}_3& {\mathrm{\σ}}_r{I}_3& {\mathrm{\τ}}_{r\mathrm{\θ}}{I}_3\\ {\mathrm{\τ}}_{x\mathrm{\θ}}{I}_3& {\mathrm{\τ}}_{r\mathrm{\θ}}{I}_3& {\mathrm{\σ}}_{\mathrm{\θ}}{I}_3\end{array}\right]---\left(20\right)$

Wherein

$\left[\begin{array}{c}{\mathrm{\σ}}_x\\ {\mathrm{\σ}}_r\\ {\mathrm{\σ}}_{r\mathrm{\θ}}\\ {\mathrm{\τ}}_{xr}\\ {\mathrm{\τ}}_{r\mathrm{\θ}}\\ {\mathrm{\τ}}_{\mathrm{\θ}x}\end{array}\right]=D\left[\begin{array}{c}{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_r\\ {\mathrm{\ϵ}}_{\mathrm{\θ}}\\ {\mathrm{\γ}}_{xr}\\ {\mathrm{\γ}}_{r\mathrm{\θ}}\\ {\mathrm{\γ}}_{\mathrm{\θ}x}\end{array}\right]---\left(21\right)$

$\begin{array}{c}{\mathrm{\ϵ}}_x=\frac{\∂u_x}{\∂x},{\mathrm{\ϵ}}_r=\frac{\∂u_r}{\∂r},{\mathrm{\ϵ}}_{\mathrm{\θ}}=\frac{u_r}{r}+\frac{1}{r}\frac{\∂u_{\mathrm{\θ}}}{\∂\mathrm{\θ}}\\ {\mathrm{\γ}}_{xr}=\frac{\∂u_x}{\∂r}+\frac{\∂u_r}{\∂x},{\mathrm{\γ}}_{r\mathrm{\θ}}=\frac{1}{r}\frac{\∂u_r}{\∂\mathrm{\θ}}+\frac{\∂u_{\mathrm{\θ}}}{\∂r}-\frac{u_{\mathrm{\θ}}}{r},{\mathrm{\γ}}_{\mathrm{\θ}x}=\frac{\∂u_{\mathrm{\θ}}}{\∂x}+\frac{1}{r}\frac{\∂u_x}{\∂\mathrm{\θ}}\end{array}---\left(22\right)$

Formula (8) is substituted into formula (22) each components of strain can be obtained, then substitute into formula (21) each components of stress can be obtained.Then the prestress tangent stiffness matrix of unit is

K _Te＝π∫∫(G ^(c)TM _σG ^(c)+G ^(s)TM _σG ^(s))rdrdx(23)

Stiffness matrix and the mass matrix of two-dimensional axial symmetric model can be obtained by assembling.Owing to describing the nodal displacement vector q of solid of revolution wheel disc distortion ^{(D, c)}and q ^{(D, s)}all be defined on meridian ellipse Oxr, vector dimension is much smaller than complete three-dimensional wheel disc model.

Two, the foundation of three dimendional blade model in dimensional model is become

Suppose to only have one group of blade in complete impeller pattern, blade circumference number is N, and by it around the shaft counterclockwise successively from 0 to N-1 numbering, wherein No. 0 blade of starting point is specified arbitrarily.

By No. 0 3D solid blade modeling in rectangular coordinate system Oxyz, its true origin and x-axis with require that in 2, two-dimensional axial symmetric wheel disc model is identical, and y-axis is identical with requiring the r axle in 2 in the plane of θ=0.The position of No. 0 blade should rotate with two-dimensional axial symmetric wheel disc model the solid of revolution wheel disc obtained and coordinate mutually.

No. 0 blade defines cosine displacement function with sinusoidal displacement function so for n pitch diameter displacement model, kth blade (k=0 ..., N-1) on displacement function be

$\left\{\begin{array}{c}{u}_x(x,y,z)\\ u_y(x,y,z)\\ u_z(x,y,z)\end{array}\right\}=\left\{\begin{array}{c}{u}_x^{(B,c)}(x,y,z)\\ u_y^{(B,c)}(x,y,z)\\ u_z^{(B,c)}(x,y,z)\end{array}\right\}\cos nk\mathrm{\α}+\left\{\begin{array}{c}{u}_x^{(B,s)}(x,y,z)\\ u_y^{(B,s)}(x,y,z)\\ u_z^{(B,s)}(x,y,z)\end{array}\right\}\sin nk\mathrm{\α}---\left(24\right)$

Wherein α=2 π/N is the differential seat angle of adjacent two blades.For zero pitch diameter displacement model, because sinnk α perseverance is zero, then only comprise cosine displacement function.

No. 0 blade is divided into three-dimensional finite element mesh, total n _bindividual node.One group of node be connected with two-dimensional axial symmetric wheel disc is set bottom it, connected node number n in connected node number and two-dimensional axial symmetric wheel disc should be ensured _jidentical, and the x-axis coordinate of both corresponding point is identical.So-called blade tie point, be exactly the corresponding tie point of a wheel disc that this group node all drops on revolution arc on.

The front n of No. 0 blade three-dimensional finite element model might as well be supposed _jindividual node is connected node, its node cosine motion vector with node sinusoidal displacement vector respectively there are six components, comprise three translation displacements components and three rotational displacement components.By defining linear restriction relation between blade connected node partial interior adjacent thereto nodal displacement component, each connected node is made all to can be used as static determinacy restraint joint, for the coupling of leaf dish provides condition.The cosine motion vector of other nodes with sinusoidal displacement vector without particular restriction, with reference to conventional finite element modeling method.

Utilize finite element interpolation, No. 0 blade cosine displacement function and sinusoidal displacement function can be separated into

$\left\{\begin{array}{c}{u}_x^{(B,c)}(x,y,z)\\ u_y^{(B,c)}(x,y,z)\\ u_z^{(B,c)}(x,y,z)\end{array}\right\}=N^{\left(B\right)}(x,y,z)q^{(B,c)}---\left(25\right)$

$\left\{\begin{array}{c}{u}_x^{(B,s)}(x,y,z)\\ u_y^{(B,s)}(x,y,z)\\ u_z^{(B,s)}(x,y,z)\end{array}\right\}=N^{\left(B\right)}(x,y,z)q^{(B,s)}---\left(26\right)$

Wherein, N ^(B)(x, y, z) is shape function

$\begin{array}{c}{q}^{(B,c)}=\left\{\begin{array}{c}{q}_I^{(B,c)}\\ q_J^{(B,c)}\end{array}\right\},q^{(B,s)}=\left\{\begin{array}{c}{q}_I^{(B,s)}\\ q_J^{(B,s)}\end{array}\right\}\\ q_I^{(B,c)}={\mathrm{col}}_{i=n_J+1}^{n_B}\left\{q_i^{(B,c)}\right\},q_J^{(B,c)}={\mathrm{col}}_{j=1}^{n_J}\left\{q_j^{(B,c)}\right\}\\ q_I^{(B,s)}={\mathrm{col}}_{i=n_J+1}^{n_B}\left\{q_i^{(B,s)}\right\},q_J^{(B,s)}={\mathrm{col}}_{j=1}^{n_J}\left\{q_j^{(B,s)}\right\}\end{array}---\left(27\right)$

For the displacement of n pitch diameter, the connected node on No. 0 blade and the corresponding connection node on wheel disc meet following displacement coordination relation

$\begin{array}{c}\left\{\begin{array}{c}{q}_i^{(B,c)}\\ q_i^{(B,s)}\end{array}\right\}=\left[\begin{array}{cc}{T}_i{\mathrm{cosn\β}}_i& T_i{\mathrm{sinn\β}}_i\\ -T_i{\mathrm{sinn\β}}_i& T_i{\mathrm{cosn\β}}_i\end{array}\right]\left\{\begin{array}{c}{q}_i^{(D,c)}\\ q_i^{(D,s)}\end{array}\right\},i=1,...,n_J\\ T_i=\left[\begin{array}{cc}{R}_i& 0\\ 0& R_i\end{array}\right],R_i=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\mathrm{cos\β}}_i& -{\mathrm{sin\β}}_i\\ 0& {\mathrm{sin\β}}_i& {\mathrm{cos\β}}_i\end{array}\right]\end{array}---\left(28\right)$

Wherein β _ibe the differential seat angle of i-th tie point tie point corresponding to wheel disc in cylindrical coordinate in No. 0 blade.

Thus, kth bar blade (k=0 ..., N-1) displacement function can by q ^(B)={ q ^{(B, c) T}q ^{(B, s) T}} ^tbe expressed as

$\left\{\begin{array}{c}{u}_x^{\left(k\right)}(x,y,z)\\ u_y^{\left(k\right)}(x,y,z)\\ u_y^{\left(k\right)}(x,y,z)\end{array}\right\}=\[\begin{array}{cc}{N}^{\left(B\right)}\cos nk\mathrm{\α}& N^{\left(B\right)}\sin nk\mathrm{\α}\end{array}\]\left\{\begin{array}{c}{q}^{(B,c)}\\ q^{(B,s)}\end{array}\right\}---\left(29\right)$

For zero pitch diameter displacement model, motion vector q ^(B)only include q ^{(B, c)}.

Total potential energy under the n pitch diameter displacement model of N bar blade is

$\begin{array}{c}U=\frac{1}{2}\[q^{(B,c)T},q^{(B,s)T}\]\underset{k=0}{\overset{N-1}{\Σ}}\left[\begin{array}{cc}{K}_{b0}{\cos }^{2}nk\mathrm{\α}& K_{b0}\cos nk\mathrm{\α}\sin nk\mathrm{\α}\\ K_{b0}\cos nk\mathrm{\α}\sin nk\mathrm{\α}& K_{b0}{\sin }^{2}nk\mathrm{\α}\end{array}\right]\left\{\begin{array}{c}{q}^{(B,c)}\\ q^{(B,s)}\end{array}\right\}\\ =\left\{\begin{array}{cc}\frac{N}{2}{q}^{(B,c)T}{K}_{b0}{q}^{(B,c)}& n=0\\ \frac{N}{4}\[q^{(B,c)T},q^{(B,s)T}\]\left[\begin{array}{cc}{K}_{b0}& 0\\ 0& K_{b0}\end{array}\right]\left\{\begin{array}{c}{q}^{(B,c)}\\ q^{(B,s)}\end{array}\right\}& n>0\end{array}\right.\end{array}---\left(30\right)$

K _b0＝∫∫∫B ^TD ₀Bdxdydz(31)

Wherein B is the strain-transposed matrix represented by shape function, D ₀for elastic matrix.Use herein

$\begin{array}{c}\underset{k=0}{\overset{N-1}{\Σ}}{\cos }^2\left(\frac{2\mathrm{\π}nk}{N}\right)=\underset{k=0}{\overset{N-1}{\Σ}}{\sin }^2\left(\frac{2\mathrm{\π}nk}{N}\right)=\frac{N}{2}\\ \underset{k=0}{\overset{N-1}{\Σ}}cos\left(\frac{2\mathrm{\π}nk}{N}\right)\sin \left(\frac{2\mathrm{\π}nk}{N}\right)=0\end{array}---\left(32\right)$

Total kinetic energy under the n pitch diameter displacement model of N bar blade is

$\begin{array}{c}E=\frac{1}{2}\[{\stackrel{\·}{q}}^{(B,c)T},{\stackrel{\·}{q}}^{(B,s)T}\]\underset{k=0}{\overset{N-1}{\Σ}}\left[\begin{array}{cc}{M}_{b0}{\cos }^{2}nk\mathrm{\α}& M_{b0}\cos nk\mathrm{\α}\sin nk\mathrm{\α}\\ M_{b0}\cos nk\mathrm{\α}\sin nk\mathrm{\α}& M_{b0}{\sin }^{2}nk\mathrm{\α}\end{array}\right]\left\{\begin{array}{c}{\stackrel{\·}{q}}^{(B,c)}\\ {\stackrel{\·}{q}}^{(B,s)}\end{array}\right\}\\ =\left\{\begin{array}{cc}\frac{N}{2}{\stackrel{\·}{q}}^{(B,c)T}{M}_{b0}{\stackrel{\·}{q}}^{(B,c)}& n=0\\ \frac{N}{4}\[{\stackrel{\·}{q}}^{(B,c)T},{\stackrel{\·}{q}}^{(B,s)T}\]\left[\begin{array}{cc}{M}_{b0}& 0\\ 0& M_{b0}\end{array}\right]\left\{\begin{array}{c}{\stackrel{\·}{q}}^{(B,c)}\\ {\stackrel{\·}{q}}^{(B,s)}\end{array}\right\}& n>0\end{array}\right.\end{array}---\left(33\right)$

M _b0＝ρ∫∫∫N ^(B)TN ^(B)dxdydz(34)

Wherein ρ is density.

N pitch diameter displacement model global stiffness matrix and the mass matrix that can calculate N bar blade are further

$K_B=\left\{\begin{array}{cc}{\mathrm{NK}}_{b0}& n=0\\ \frac{N}{2}\left[\begin{array}{cc}{K}_{b0}& 0\\ 0& K_{b0}\end{array}\right]& n>0\end{array}\right.---\left(35\right)$

$M_B=\left\{\begin{array}{cc}{\mathrm{NM}}_{b0}& n=0\\ \frac{N}{2}\left[\begin{array}{cc}{M}_{b0}& 0\\ 0& M_{b0}\end{array}\right]& n>0\end{array}\right.---\left(36\right)$

Be not difficult to find, K _b0and M _b0the stiffness matrix calculated for single 3 D entity blade conventional finite element method and mass matrix, be called single blade stiffness matrix and single leaf quality matrix.If calculated prestressing force mode, K _b0should be the tangent stiffness matrix of GEOMETRICALLY NONLINEAR.

For the blade wheel structure of many group blades, then respectively organize No. 0 independent modeling of blade of blade, stiffness matrix and the mass matrix of each group of blade can be obtained.

Three, become dimensional model Coupling Deformation in centrifugal force field to calculate

For the Coupling Deformation problem preloaded, this method refers in particular to impeller by preloading problem under centrifugal action, and correspond to the zero pitch diameter displacement model becoming dimensional model, therefore No. 0 blade displacement vector only comprises q ^{(B, c)}.

Consider the blade wheel structure of single group blade, using the interface point of the tie point of No. 0 leaf three-dimensional model as static(al) polycondensation, utilize static condensation method to carry out static(al) polycondensation [1].Suppose that blade is subject to centrifugal action, the balance equation of leaf model is

$\left[\begin{array}{cc}{K}_{bii}& K_{bij}\\ K_{bji}& K_{bjj}\end{array}\right]\left\{\begin{array}{c}{q}_I^{(B,c)}\\ q_J^{(B,c)}\end{array}\right\}=\left\{\begin{array}{c}{p}_I^{(B,c)}\\ p_J^{(B,c)}\end{array}\right\}---\left(37\right)$

Wherein for blade interior nodal displacement; for the displacement of blade interface node; with the vector that centrifugal force produces suffered by node is known quantity.Then internal node displacement can be gone out by blade interface node offset table, namely

$q_I^{(B,c)}=K_{bii}^{-1}{p}_I^{(B,c)}-K_{bii}^{-1}{K}_{bij}{q}_J^{(B,c)}---\left(38\right)$

(38) are substituted into (37) second equations, the nodal force vector of static(al) polycondensation rear interface point under centrifugal action can be obtained and Stiffness Matrix

$\begin{array}{c}{\tilde{P}}_b=p_J^{(B,c)}-K_{bji}{K}_{bii}^{-1}{p}_I^{(B,c)}\\ {\tilde{K}}_b=K_{bjj}-K_{bij}{K}_{bii}^{-1}{K}_{bij}\end{array}---\left(39\right)$

By interface node motion vector from the rectangular coordinate system of definition 3D solid leaf model be transformed into two-dimensional axial symmetric wheel disc model cylindrical coordinate, namely

$\begin{array}{c}{q}_J^{(B,c)}=\tilde{T}{\hat{q}}_J^{(B,c)}\\ \tilde{T}={\mathrm{diag}}_{i=1}^{n_J}\left(T_i\right)\end{array}---\left(40\right)$

Wherein T _isee formula (28), obtain convert after polycondensation stiffness matrix and force vector be

$\begin{array}{c}{K}_b={\tilde{T}}^T{\tilde{K}}_b\tilde{T}\\ P_b={\tilde{T}}^T{\tilde{P}}_b\end{array}---\left(41\right)$

For two-dimensional axial symmetric wheel disc, equally using the interface node of tie point as polycondensation, the balance equation that can coil under centrifugal action is

$\left[\begin{array}{cc}{K}_{dii}& K_{dij}\\ K_{dji}& K_{djj}\end{array}\right]\left\{\begin{array}{c}{q}_I^{(D,c)}\\ q_I^{(D,c)}\end{array}\right\}=\left\{\begin{array}{c}{p}_I^{(D,c)}\\ p_I^{(D,c)}\end{array}\right\}---\left(42\right)$

Wherein for the displacement of dish internal node; for the displacement of dish interface node; with the vector that centrifugal force produces suffered by node is known quantity. by show out for

$q_I^{(D,c)}=K_{dii}^{-1}{p}_J^{(D,c)}-K_{dii}^{-1}{K}_{dij}{q}_J^{(D,c)}---\left(43\right)$

In like manner obtain the force vector P of polycondensation rear interface point _dwith stiffness matrix K _d

$\begin{array}{c}{P}_d=p_J^{(D,c)}-K_{dij}{K}_{dii}^{-1}{p}_I^{(D,c)}\\ K_d=K_{djj}-K_{dji}-K_{dii}^{-1}{K}_{dij}\end{array}---\left(44\right)$

Then the rigidity of unitized construction and force vector are

K＝K _d+NK _b

(45)

P＝P _d+NP _b

Wherein N is the circumferential number of blade.

Utilize the stiffness matrix after assembling and force vector, the displacement of placing interface node can be separated

$q_J^{(D,c)}=K^{-1}P---\left(46\right)$

With the displacement of blade interface node

$q_J^{(B,c)}=\tilde{T}{q}_J^{(D,c)}---\left(47\right)$

Further formula (47) and formula (46) are substituted into respectively the internal node motion vector that formula (38) and formula (43) can obtain blade and wheel disc.

For 3D solid leaf model, strain and stress can be obtained by nodal displacement vector by general finite Meta algorithm and distribute; For two-dimensional axial symmetric wheel disc model, strain and stress can be obtained by nodal displacement vector substitution formula (22) and formula (21) and distribute.

For the situation of multiple blade, then the rigidity of unitized construction and force vector need to superpose more, and computation process is similar.

Four, the Dynamics Coupling of dimensional model Leaf and wheel disc is become

For n pitch diameter displacement model (n > 0), the Free Vibration Equations of single 3 D entity blade cosine motion vector is

$M_{b0}{\stackrel{\·\·}{q}}^{(B,c)}+K_{b0}{q}^{(B,c)}=0---\left(48\right)$

The form being write as partitioned matrix is

$\left[\begin{array}{cc}{M}_{bii}& M_{bij}\\ M_{bji}& M_{bjj}\end{array}\right]\left\{\begin{array}{c}{\stackrel{\·\·}{q}}_I^{(B,c)}\\ {\stackrel{\·\·}{q}}_J^{(B,c)}\end{array}\right\}+\left[\begin{array}{cc}{K}_{bii}& K_{bij}\\ K_{bji}& K_{bjj}\end{array}\right]\left\{\begin{array}{c}{q}_I^{(B,c)}\\ q_J^{(B,c)}\end{array}\right\}=0---\left(49\right)$

Wherein for blade interior nodal displacement; for the displacement of blade interface node.

Utilize the Craig-Bampton method [2] of Dynamic Substructure, by single 3 D entity leaf model cosine motion vector q ^{(B, c)}with part cosine immobile interface cylindrical coordinates ξ ^(c)and whole tie point cosine motion vector be expressed as

$q^{(B,c)}=\[\begin{array}{cc}{\mathrm{\Φ}}_b& {\mathrm{\Ψ}}_b\end{array}\]\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(c\right)}\\ q_J^{(B,c)}\end{array}\right\}=T_b\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(c\right)}\\ q_J^{(B,c)}\end{array}\right\}---\left(50\right)$

In formula (50), Φ _bbe called immobile interface master mode, meet

$\begin{array}{c}{\mathrm{\Φ}}_b=\left[\begin{array}{c}{\widetilde{\mathrm{\Φ}}}_b\\ 0\end{array}\right]\\ {\widetilde{\mathrm{\Φ}}}_b^T{M}_{bii}{\widetilde{\mathrm{\Φ}}}_b=I_{n_K}\\ {\widetilde{\mathrm{\Φ}}}_b^T{K}_{bii}{\widetilde{\mathrm{\Φ}}}_b=diag\left({\mathrm{\ω}}_{b1}^2,{\mathrm{\ω}}_{b2}^2,...,{\mathrm{\ω}}_{b{n}_K}^2\right)\end{array}---\left(51\right)$

Namely for vibration equation (52) corresponding about M _biinormalized front n _krank modal vector, ω _bifor corresponding circle frequency.

$M_{bii}{\stackrel{\·\·}{q}}_I^{(B,c)}+K_{bii}{q}_I^{(B,c)}=0---\left(52\right)$

In formula (50), Ψ _bbe called Constrained mode, meet

${\mathrm{\Ψ}}_b=\left[\begin{array}{c}C\\ I\end{array}\right]=\left[\begin{array}{c}-K_{bii}^{-1}{K}_{bij}\\ I\end{array}\right]---\left(53\right)$

Formula (50) is substituted into formula (49), and simultaneously premultiplication can obtain the stiffness matrix after polycondensation is

$\begin{array}{c}{\stackrel{\‾}{K}}_b=\left[\begin{array}{cc}{\stackrel{\‾}{K}}_{\mathrm{\ξ}\mathrm{\ξ}}& {\stackrel{\‾}{K}}_{{\mathrm{\ξq}}_J}\\ {\stackrel{\‾}{K}}_{{\mathrm{\ξq}}_J}& {\stackrel{\‾}{K}}_{q_J{q}_J}\end{array}\right]\\ {\stackrel{\‾}{K}}_{\mathrm{\ξ}\mathrm{\ξ}}=diag\left({\mathrm{\ω}}_{b1}^2,{\mathrm{\ω}}_{b2}^2,...,{\mathrm{\ω}}_{b{n}_K}^2\right)\\ {\stackrel{\‾}{K}}_{{\mathrm{\ξq}}_J}=0\\ {\stackrel{\‾}{K}}_{q_J{q}_J}=K_{bjj}-K_{bij}^T{K}_{bii}^{-1}{K}_{bij}\end{array}---\left(54\right)$

Mass matrix after polycondensation is

$\begin{array}{c}{\stackrel{\‾}{M}}_b=\left[\begin{array}{cc}{\stackrel{\‾}{M}}_{\mathrm{\ξ}\mathrm{\ξ}}& {\stackrel{\‾}{M}}_{{\mathrm{\ξq}}_J}\\ {\stackrel{\‾}{M}}_{{\mathrm{\ξq}}_J}^T& {\stackrel{\‾}{M}}_{q_J{q}_J}\end{array}\right]\\ {\stackrel{\‾}{M}}_{\mathrm{\ξ}\mathrm{\ξ}}=I_{n_K}\\ {\stackrel{\‾}{M}}_{{\mathrm{\ξq}}_J}={{\widetilde{\mathrm{\Φ}}}_b}^T{M}_{bij}+{{\widetilde{\mathrm{\Φ}}}_b}^T{M}_{bii}C\\ {\stackrel{\‾}{M}}_{q_J{q}_J}=M_{bjj}+C^T{M}_{bii}C+C^T{M}_{bij}+M_{bij}^{T}C\end{array}---\left(55\right)$

The wherein corresponding modal coordinate of ξ, q _jeach connected node component corresponding.

Example is transformed to, by connected node motion vector with stiffness matrix transform in Oxr θ cylindrical coordinate

$q_J^{(B,c)}=\tilde{T}{\hat{q}}_J^{(B,c)}---\left(56\right)$

Then the matrixing of blade reduced model list blade stiffness is

$\begin{array}{c}{\hat{K}}_{b0}={\stackrel{\‾}{T}}^T{\stackrel{\‾}{K}}_{b0}\stackrel{\‾}{T}\left[\begin{array}{cc}{\hat{K}}_{\mathrm{\ξ}\mathrm{\ξ}}& {\hat{K}}_{{\mathrm{\ξq}}_J}\\ {\hat{K}}_{{\mathrm{\ξq}}_J}& {\hat{K}}_{q_J{q}_J}\end{array}\right]\\ \stackrel{\‾}{T}=\left[\begin{array}{cc}{I}_{n_K}& 0\\ 0& \tilde{T}\end{array}\right]\end{array}---\left(57\right)$

In like manner offset of sinusoidal motion vector carries out polycondensation conversion, can obtain displacement relation

$\left\{\begin{array}{c}{q}_I^{(B,s)}\\ q_J^{(B,s)}\end{array}\right\}=\[\begin{array}{cc}{\mathrm{\Φ}}_b& {\mathrm{\Ψ}}_b\end{array}\]\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(s\right)}\\ q_J^{(B,s)}\end{array}\right\}=\left[\begin{array}{cc}{\widetilde{\mathrm{\Φ}}}_b& C\\ 0& I\end{array}\right]\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(s\right)}\\ q_J^{(B,s)}\end{array}\right\}=T_b\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(s\right)}\\ q_J^{(B,s)}\end{array}\right\}---\left(58\right)$

Its immobile interface master mode is identical with cosine motion vector polycondensation result of calculation with Constrained mode, and coordinate conversion mode is identical, thus gained polycondensation matrix with identical, then n pitch diameter displacement model (n > 0) the global stiffness matrix polycondensation of N bar blade is transformed to

${\hat{K}}_b=\frac{N}{2}\left[\begin{array}{cc}{\hat{K}}_{b0}& 0\\ 0& {\hat{K}}_{b0}\end{array}\right]---\left(59\right)$

By formula (28), in cylindrical coordinate, the interfacial displacement vector of No. 0 blade reduced model can be expressed as by wheel disc tie point motion vector

$\left\{\begin{array}{c}{\hat{q}}_J^{(B,c)}=diag\left(I_6{\mathrm{cosn\β}}_i\right)q_J^{(D,c)}+diag\left(I_6{\mathrm{sinn\β}}_i\right)q_J^{(D,s)}\\ {\hat{q}}_J^{(B,s)}=diag(-I_6{\mathrm{sinn\β}}_i)q_J^{(D,c)}+diag\left(I_6{\mathrm{cosn\β}}_i\right)q_J^{(D,s)}\end{array}\right.---\left(60\right)$

Note

$T=\left[\begin{array}{cccc}{I}_{n_K}& & 0& \\ & diag\left(I_6{\mathrm{cosn\β}}_i\right)& & diag\left(I_6{\mathrm{sinn\β}}_i\right)\\ 0& & I_{n_K}& \\ & diag(-I_6{\mathrm{sinn\β}}_i)& & diag\left(I_6{\mathrm{cosn\β}}_i\right)\end{array}\right]---\left(61\right)$

Then have

$\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(c\right)}\\ {\hat{q}}_J^{(B,c)}\\ {\mathrm{\ξ}}^{\left(c\right)}\\ {\hat{q}}_J^{(B,s)}\end{array}\right\}=T\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(c\right)}\\ q_J^{(D,c)}\\ {\mathrm{\ξ}}^{\left(s\right)}\\ q_J^{(D,s)}\end{array}\right\}---\left(62\right)$

Global stiffness matrix then after all blade polycondensations is transformed to further

$K=T^T{\hat{K}}_{b}T=\frac{N}{2}\left[\begin{array}{cc}{K}_c& K_s\\ K_s^T& K_c\end{array}\right]---\left(63\right)$

Wherein

For mass matrix, the oeverall quality matrix after all blade polycondensations can be obtained by identical method, namely

$M_b=\frac{N}{2}\left[\begin{array}{cc}{M}_c& M_s\\ M_s^T& M_c\end{array}\right]---\left(66\right)$

Wherein

${\hat{M}}_{b0}={\stackrel{\‾}{T}}^T{\stackrel{\‾}{M}}_{b0}\stackrel{\‾}{T}---\left(67\right)$

Finally by conversion after stiffness matrix and mass matrix general assembly in the stiffness matrix of two-dimensional axial symmetric wheel disc model and the correspondence position of mass matrix, obtain global stiffness matrix and the mass matrix of model, just can carry out the n pitch diameter model analysis of two-dimensional axial symmetric wheel disc and 3D solid blade change dimensional model.Suppose that the coupling Free Vibration Equations after general assembly is

$M\left\{\begin{array}{c}{\stackrel{\·\·}{\mathrm{\ξ}}}^{\left(c\right)}\\ {\stackrel{\·\·}{q}}_J^{(D,c)}\\ {\stackrel{\·\·}{q}}_I^{(D,c)}\\ {\stackrel{\·\·}{\mathrm{\ξ}}}^{\left(s\right)}\\ {\stackrel{\·\·}{q}}_J^{(D,s)}\\ {\stackrel{\·\·}{q}}_I^{(D,s)}\end{array}\right\}+K\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(c\right)}\\ q_J^{(D,c)}\\ q_I^{(D,c)}\\ {\mathrm{\ξ}}^{\left(s\right)}\\ q_J^{(D,s)}\\ q_I^{(D,s)}\end{array}\right\}=0---\left(70\right)$

The generalized eigenvalue problem that calculating formula (70) is corresponding, can obtain becoming the n pitch diameter natural frequency of dimensional model and corresponding mode thereof.Suppose that certain the rank Mode Shape vector obtained is

Six components are corresponding blade modal coordinates respectively, dish on connected node displacement with coil the cosine component of internal node displacement and sinusoidal component to should the component of rank mode.

For two-dimensional axial symmetric wheel disc, can utilize substitution formula (8), obtains the modal displacement vector coiling each node upper, thus recovers the mode of whole dish.

For kth bar blade, utilize recover the motion vector that it is corresponding in rectangular coordinate system namely

Wherein for formula (57) homography, T is formula (61) homography.Formula (72) is substituted into formula (50) and formula (58), and acquired results substitutes into formula (29), then the node modal displacement vector of kth bar blade can be expressed as

For zero pitch diameter displacement model, the global stiffness after becoming dimensional model polycondensation and mass matrix can be obtained by similar method, carry out zero pitch diameter model analysis.

Beneficial effect of the present invention is:

1, its method can realize the prestressed modal analysis that the 3D solid blade of blade wheel structure and two-dimensional axial symmetric wheel disc become dimensional model, under the prerequisite ensureing precision, improves the model analysis counting yield of impeller class rotationally periodic structure.

2, its method can specify the pitch diameter number of model analysis, targeted for particular problem.Such as in engine speed model analysis, only have impeller a pitch diameter mode can with rotating shaft generation coupled vibrations, then utilize this method directly to analyze a pitch diameter mode, avoid calculating other unnecessary mode.

3, in this method, the shape of blade wheel structure is not subject to the restriction of gyration period.At present, model analysis is carried out for impeller class formation many employings sector models, makes analyzable blade wheel structure be limited to characteristic gyration period.For wheel disc being provided with many groups, the impeller of the different blade of circumferential number, often identical sector strucre cannot be marked off.The blade of this method and wheel sections can Independent modelings, and the shape of blade wheel structure by the restriction of gyration period, thus solves complicated impeller pattern cannot carry out a model analysis difficult problem with sector models.

Accompanying drawing explanation

Fig. 1 is the Complete three-dimensional finite element model of blade wheel structure.

Fig. 2 is the finite element model of two-dimensional axial symmetric wheel disc.

Fig. 3 is that the multi-point constraint at two-dimensional axial symmetric wheel disc connected node place adds situation schematic diagram.

Fig. 4 is the finite element model of 3D solid blade.

Fig. 5 is that the multi-point constraint at 3D solid blade bottom connected node place adds situation schematic diagram.

Fig. 6 is that 3D solid blade and two-dimensional axial symmetric wheel disc become dimensional model schematic diagram.

Fig. 7 is the sector finite element that blade wheel structure is corresponding.

Embodiment

With example, the present invention is elaborated further by reference to the accompanying drawings.

Technology of the present invention realizes the model analysis that the 3D solid blade of blade wheel structure and two-dimensional axial symmetric wheel disc become dimensional model, under the prerequisite ensureing precision, improves the model analysis counting yield of impeller class rotationally periodic structure.Specific implementation process comprises the steps:

(1) set up two-dimensional axial symmetric wheel disc model and 3D solid leaf model, mate corresponding tie point.

As shown in Figure 1, it has 17 blades to example three-dimensional blades model used.Impeller material parameter is: Young modulus 1.18e+11Pa; Poisson ratio 0.3; Density 4400kg/m ³.

The modeling method proposed according to the present invention and requirement, set up two-dimensional axial symmetric wheel disc model and 3D solid leaf model respectively, and add certain multi-point constraint constraint in each comfortable model.Wherein two-dimensional axial symmetric wheel disc model as shown in Figure 2, and its multi-point constraint adds situation as shown in Figure 3; As shown in Figure 4, its multi-point constraint adds situation as shown in Figure 5 to 3D solid leaf model.

3D solid leaf model and two-dimensional axial symmetric wheel disc model are placed in same view, the schematic diagram that equivalence becomes dimensional model can be obtained, as shown in Figure 6.As can be seen from Figure 6, the position of 3D solid blade is coordinated mutually with two-dimensional axial symmetric wheel disc model.

The corresponding connection node information of 3D solid leaf model and two-dimensional axial symmetric wheel disc model is as shown in table 1.Can see from table 1, two model connected node numbers are identical, and the x-axis coordinate of corresponding connection node is substantially identical, and error is all within 2/10000ths millimeters.Therefore the 3D solid blade of this example and two-dimensional axial symmetric wheel disc model meet the modeling demand of the inventive method.

Table 1 corresponding connection node information

(2) the leaf dish coupling scheme utilizing the present invention to propose, calculate and become dimensional model mode.

The pitch diameter mode that this example becomes dimensional model for impeller calculates.The reference object of result of calculation is the FEM modal analysis and modal of business software Nastran to same Impeller fan section model.

Nastran sector models used as shown in Figure 7, is 1/17th of complete impeller pattern.

The pitch diameter FEM modal analysis and modal that impeller central is fixed is as shown in table 2.In general, the error of sector models result that calculates of the change dimensional model FEM modal analysis and modal that proposes of the present invention and Nastran is within 1%.Become computing time of dimensional model much smaller than Nastran three-dimensional sectors model, within shortening to 1/10th of Nastran computing time.

This example shows, the 3D solid blade of blade wheel structure that the inventive method proposes and two-dimensional axial symmetric wheel disc become dimensional model and modal analysis method effective, it is close that it calculates three-dimensional sectors model result in result and existing business software Nastran, computing time much smaller than the Nastran time used, the analysis efficiency of raising impeller class formation that can be larger and do not lose analysis precision.

Table 2 one pitch diameter FEM modal analysis and modal

List of references

[1]Wilson,EdwardL.Thestaticcondensationalgorithm[J].InternationalJournalforNumericalMethodsinEngineering,1974,8(1):198-203.

[2]M.C.C.Bampton,R.R.CraigJr.Couplingofsubstructuresfordynamicanalysis.AIAAJ[J].AIAAJournal,1968,6:1313-1319。

Claims (6)

1. the 3D solid blade of blade wheel structure and two-dimensional axial symmetric wheel disc become a dimensional model finite element prestressed modal analysis method, and it is characterized in that, concrete steps are as follows:

(1) set up the blade wheel structure be made up of blade and wheel disc and become dimension finite element model; First blade 3D solid unit represents by it, sets up 3D solid leaf model; Wheel disc two-dimensional axial symmetric unit is represented, sets up two-dimensional axial symmetric wheel disc model; Again by arranging one group of connected node respectively in 3D solid leaf model and two-dimensional axial symmetric wheel disc model, and between the partial interior node and connected node of 3D solid leaf model and two-dimensional axial symmetric wheel disc model, set up the linear restriction relation of one group of displacement component respectively, make each connected node all as static determinacy restraint joint, realize the combination that impeller becomes dimensional model;

(2) dynamic sub-structure methods is utilized to obtain the stiffness matrix after the polycondensation of 3D solid leaf model and mass matrix, be coupled with the stiffness matrix of two-dimensional axial symmetric wheel disc model and mass matrix, utilize the kinetic model after coupling to carry out the model analysis of change dimensional model.

2. change dimensional model finite element prestressed modal analysis method according to claim 1, is characterized in that:

Utilize the Coupling Deformation computing method of 3D solid leaf model and two-dimensional axial symmetric wheel disc model, obtain the quiet distortion of impeller under centrifugal force preloads effect, should be the tangent stiffness matrix of GEOMETRICALLY NONLINEAR simultaneously with stiffness matrix, analyze the pre-stressed mode under impeller change dimensional model centrifugal force field.

3. change dimensional model finite element prestressed modal analysis method according to claim 1, is characterized in that:

In step (1), the modeling method of two-dimensional axial symmetric wheel disc model is as follows:

Two-dimensional axial symmetric wheel disc model defines in cylindrical coordinate Oxr θ, and three axles of this coordinate system are respectively by the x-axis of axis, and the θ axle of radial r axle and circumference forms; Two-dimensional axial symmetric wheel disc finite element model defines in the section Oxr of θ=0, is made up of two-dimensional axial symmetric unit, and rotating a circle around x-axis obtains the on all four solid of revolution wheel disc model with three-dimensional model;

For circumferential wave number be n solid of revolution wheel disc distortion, be called n pitch diameter displacement model; In cylindrical coordinates, the axial displacement u on its optional position (x, r, θ) _x, radial displacement u _rwith circumferential displacement components u _θθ=0 meridian ellipse is defined in, i.e. the cosine displacement function of 0-meridian ellipse with one group with sinusoidal displacement function be expressed as:

$\left\{\begin{array}{c}{u}_x(x,r,\mathrm{\θ})\\ u_r(x,r,\mathrm{\θ})\\ u_{\mathrm{\θ}}(x,r,\mathrm{\θ})\end{array}\right\}\stackrel{\mathrm{\Δ}}{=}\left\{\begin{array}{c}{u}_x^{\left(c\right)}(x,r)\\ u_r^{\left(c\right)}(x,r)\\ u_{\mathrm{\θ}}^{\left(c\right)}(x,r)\end{array}\right\}\cos n\mathrm{\θ}+\left\{\begin{array}{c}{u}_x^{\left(s\right)}(x,r)\\ u_r^{\left(s\right)}(x,r)\\ u_{\mathrm{\θ}}^{\left(s\right)}(x,r)\end{array}\right\}\sin n\mathrm{\θ}$

For zero pitch diameter displacement model, because sinn θ perseverance is zero, then only comprise cosine displacement function;

The components of strain under cylindrical coordinate are expressed as:

$\mathrm{\ϵ}=\left\{\begin{array}{c}{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_r\\ {\mathrm{\ϵ}}_{\mathrm{\θ}}\\ {\mathrm{\γ}}_{xr}\\ {\mathrm{\γ}}_{r\mathrm{\θ}}\\ {\mathrm{\γ}}_{\mathrm{\θ}x}\end{array}\right\}=\left[\begin{array}{ccc}\frac{\∂}{\∂x}& 0& 0\\ 0& \frac{\∂}{\∂r}& 0\\ 0& \frac{1}{r}& \frac{1}{r}\frac{\∂}{\∂\mathrm{\θ}}\\ \frac{\∂}{\∂r}& \frac{\∂}{\∂x}& 0\\ 0& \frac{1}{r}\frac{\∂}{\∂\mathrm{\θ}}& \frac{\∂}{\∂r}-\frac{1}{r}\\ \frac{1}{r}\frac{\∂}{\∂\mathrm{\θ}}& 0& \frac{\∂}{\∂x}\end{array}\right]\left\{\begin{array}{c}{u}_x(x,r,\mathrm{\θ})\\ u_r(x,r,\mathrm{\θ})\\ u_{\mathrm{\θ}}(x,r,\mathrm{\θ})\end{array}\right\}$

Potential energy and the kinetic energy of the n pitch diameter displacement of solid of revolution wheel disc are expressed as:

$U=\frac{1}{2}\∫\∫\∫{\mathrm{\ϵ}}^{T}D\mathrm{\ϵ}dxdrd\mathrm{\θ}$

$E=\frac{1}{2}\mathrm{\ρ}\∫\∫\∫{\stackrel{\·}{u}}^T\stackrel{\·}{u}dxdrd\mathrm{\θ}$

Wherein D is elastic matrix, and ρ is density;

Wheel disc hatch region on 0-meridian ellipse is divided into two-dimensional axial symmetric finite element grid, total n _dindividual node; In two-dimensional axial symmetric wheel disc finite element model, arrange one group of node be connected with blade, connecting joint is counted as n _j, so-called wheel disc tie point, refers to the circular arc of this group node after x-axis revolution, namely turns round arc and is all connected with the corresponding node of in blade finite element model;

Suppose the front n of two-dimensional axial symmetric wheel disc model _jindividual node is connected node, its node cosine motion vector with node sinusoidal displacement vector respectively there are six components, comprise three translation displacements components and three rotational displacement components; By defining linear restriction relation between wheel disc connected node partial interior adjacent thereto nodal displacement component, each connected node is made all to can be used as static determinacy restraint joint, for the coupling of leaf dish provides condition; The cosine motion vector of other nodes with sinusoidal displacement vector without particular restriction, with reference to conventional finite element modeling method.

Utilize finite element interpolation, the cosine displacement function of wheel disc and sinusoidal displacement function be separated into:

$\left\{\begin{array}{c}{u}_x^{\left(c\right)}(x,r)\\ u_r^{\left(c\right)}(x,r)\\ u_{\mathrm{\θ}}^{\left(c\right)}(x,r)\end{array}\right\}=N^{\left(D\right)}(x,r)q^{(D,c)}$

$\left\{\begin{array}{c}{u}_x^{\left(s\right)}(x,r)\\ u_r^{\left(s\right)}(x,r)\\ u_{\mathrm{\θ}}^{\left(s\right)}(x,r)\end{array}\right\}=N^{\left(D\right)}(x,r)q^{(D,s)}$

Wherein N ^(D)(x, r) is shape function;

$q^{(D,c)}=\left\{\begin{array}{c}{q}_I^{(D,c)}\\ q_J^{(D,c)}\end{array}\right\},q^{(D,s)}=\left\{\begin{array}{c}{q}_I^{(D,s)}\\ q_J^{(D,s)}\end{array}\right\}$

$q_I^{(D,c)}={\mathrm{col}}_{i=n_J+1}^{n_D}\left\{q_i^{(D,c)}\right\},q_J^{(D,c)}={\mathrm{col}}_{j=1}^{n_J}\left\{q_j^{(D,c)}\right\}$

$q_I^{(D,s)}={\mathrm{col}}_{i=n_J+1}^{n_D}\left\{q_i^{(D,s)}\right\},q_J^{(D,s)}={\mathrm{col}}_{j=1}^{n_J}\left\{q_j^{(D,s)}\right\}$

Like this, the displacement function of whole wheel disc is by q ^(D)={ q ^{(D, c) T}, q ^{(D, s) T}} ^trepresent, for zero pitch diameter displacement model, then motion vector q ^(D)only include q ^{(D, c)}; For dissimilar two-dimensional axial symmetric unit, different shape functions is substituted into potential energy and kinetic energy expression, obtains element stiffness matrix and the element mass matrix of two-dimensional axial symmetric unit n pitch diameter displacement model; Stiffness matrix and the mass matrix of two-dimensional axial symmetric wheel disc model is obtained by assembling.

4. change dimensional model finite element prestressed modal analysis method according to claim 1, it is characterized in that, in step (1), the modeling method of 3D solid leaf model is as follows:

(1) for the impeller pattern only having one group of blade, blade circumference number is N, and by it around the shaft counterclockwise successively from 0 to N-1 numbering, wherein No. 0 blade of starting point is specified arbitrarily;

By No. 0 blade modeling in rectangular coordinate system Oxyz, its true origin is identical with two-dimensional axial symmetric wheel disc model with x-axis, and y-axis is identical with the r axle in θ=0 plane in two-dimensional axial symmetric wheel disc model; Position and the two-dimensional axial symmetric wheel disc model of No. 0 blade rotate the solid of revolution wheel disc obtained to be coordinated mutually;

No. 0 blade defines cosine displacement function with sinusoidal displacement function so for n pitch diameter displacement model, kth blade (k=0 ..., N-1) on displacement function be:

$\left\{\begin{array}{c}{u}_x(x,y,z)\\ u_y(x,y,z)\\ u_z(x,y,z)\end{array}\right\}=\left\{\begin{array}{c}{u}_x^{(B,c)}(x,y,z)\\ u_y^{(B,c)}(x,y,z)\\ u_z^{(B,c)}(x,y,z)\end{array}\right\}\cos nk\mathrm{\α}+\left\{\begin{array}{c}{u}_x^{(B,s)}(x,y,z)\\ u_y^{(B,s)}(x,y,z)\\ u_z^{(B,s)}(x,y,z)\end{array}\right\}\sin nk\mathrm{\α}$

Wherein α=2 π/N is the differential seat angle of adjacent two blades; For zero pitch diameter displacement model, because sinnk α perseverance is zero, then only comprise cosine displacement function.

No. 0 blade is divided into three-dimensional finite element mesh, total n _bindividual node; One group of node be connected with two-dimensional axial symmetric wheel disc model is set bottom it, ensures connected node number n in connected node number and two-dimensional axial symmetric wheel disc _jidentical, and the x-axis coordinate of both corresponding point is identical; So-called blade tie point, be exactly the corresponding tie point of a wheel disc that this group node all drops on revolution arc on;

Suppose the front n of No. 0 blade three-dimensional finite element model _jindividual node is connected node, its node cosine motion vector with node sinusoidal displacement vector respectively there are six components, comprise three translation displacements components and three rotational displacement components; By defining linear restriction relation between blade connected node partial interior adjacent thereto nodal displacement component, each connected node is made all to can be used as static determinacy restraint joint, for the coupling of leaf dish provides condition; The cosine motion vector of other nodes with sinusoidal displacement vector without particular restriction, with reference to conventional finite element modeling method;

Utilize finite element interpolation, No. 0 blade cosine displacement function and sinusoidal displacement function can be separated into:

$\left\{\begin{array}{c}{u}_x^{(B,c)}(x,y,z)\\ u_y^{(B,c)}(x,y,z)\\ u_z^{(B,c)}(x,y,z)\end{array}\right\}=N^{\left(B\right)}(x,y,z)q^{(B,c)}$

$\left\{\begin{array}{c}{u}_x^{(B,s)}(x,y,z)\\ u_y^{(B,s)}(x,y,z)\\ u_z^{(B,s)}(x,y,z)\end{array}\right\}=N^{\left(B\right)}(x,y,z)q^{(B,s)}$

Wherein: N ^(B)(x, y, z) is shape function;

$q^{(B,c)}=\left\{\begin{array}{c}{q}_I^{(B,c)}\\ q_J^{(B,c)}\end{array}\right\},$ $q^{(B,s)}=\left\{\begin{array}{c}{q}_I^{(B,s)}\\ q_J^{(B,s)}\end{array}\right\}$

$q_I^{(B,c)}={\mathrm{col}}_{i=n_J+1}^{n_B}\left\{q_i^{(B,c)}\right\},$ $q_J^{(B,c)}={\mathrm{col}}_{j=1}^{n_J}\left\{q_j^{(B,c)}\right\}$

$q_I^{(B,s)}={\mathrm{col}}_{i=n_J+1}^{n_B}\left\{q_i^{(B,s)}\right\},$ $q_J^{(B,s)}={\mathrm{col}}_{j=1}^{n_J}\left\{q_j^{(B,s)}\right\}$

For the displacement of n pitch diameter, the connected node on No. 0 blade and the corresponding connection node on wheel disc meet following displacement coordination relation:

$\left\{\begin{array}{c}{q}_i^{(B,c)}\\ q_i^{(B,s)}\end{array}\right\}=\left[\begin{array}{cc}{T}_i{\mathrm{cosn\β}}_i& T_i{\mathrm{sinn\β}}_i\\ -T_i{\mathrm{sinn\β}}_i& T_i{\mathrm{cosn\β}}_i\end{array}\right]\left\{\begin{array}{c}{q}_i^{(D,c)}\\ q_i^{(D,s)}\end{array}\right\},i=1,\mathrm{...},n_J$

$T_i=\left[\begin{array}{cc}{R}_i& 0\\ 0& R_i\end{array}\right],$ $R_i=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\mathrm{cos\β}}_i& -{\mathrm{sin\β}}_i\\ 0& {\mathrm{sin\β}}_i& {\mathrm{cos\β}}_i\end{array}\right]$

Wherein β _ibe i-th tie point and the differential seat angle of corresponding tie point in cylindrical coordinate in two-dimensional axial symmetric wheel disc in No. 0 blade;

Thus, the displacement function of all N bar blades is by q ^(B)={ q ^{(B, c) T}, q ^{(B, s) T}} ^trepresent, for zero pitch diameter displacement model, motion vector q ^(B)only include q ^{(B, c)}; Potential energy under the n pitch diameter displacement model of then N bar blade and kinetic energy are by q ^(B)quadratic form table go out, the n pitch diameter displacement model global stiffness matrix and the mass matrix that calculate N bar blade are further:

$K_B=\left\{\begin{array}{cc}{\mathrm{NK}}_{b0}& n=0\\ \frac{N}{2}\left[\begin{array}{cc}{K}_{b0}& 0\\ 0& K_{b0}\end{array}\right]& n>0\end{array}\right.$

$M_B=\left\{\begin{array}{cc}{\mathrm{NM}}_{b0}& n=0\\ \frac{N}{2}\left[\begin{array}{cc}{M}_{b0}& 0\\ 0& M_{b0}\end{array}\right]& n>0\end{array}\right.$

Wherein K _b0and M _b0the stiffness matrix calculated for single 3 D entity blade conventional finite element method and mass matrix, be called single blade stiffness matrix and single leaf quality matrix.If calculated prestressing force mode, K _b0should be the tangent stiffness matrix of GEOMETRICALLY NONLINEAR.

(2) for the blade wheel structure of many group blades, then respectively organize No. 0 independent modeling of blade of blade, obtain stiffness matrix and the mass matrix of each group of blade.

5. change dimensional model finite element prestressed modal analysis method according to claim 1, it is characterized in that, in step (2), the method that the stiffness matrix in model analysis and mass matrix carry out being coupled is as follows:

For n pitch diameter displacement model (n>0), utilize the Craig-Bampton method of Dynamic Substructure, by single 3 D entity leaf model cosine motion vector q ^{(B, c)}with part cosine immobile interface cylindrical coordinates ξ ^(c)and whole tie point cosine motion vector represent, gained K _b0polycondensation stiffness matrix be:

${\stackrel{\‾}{K}}_{b0}=\left[\begin{array}{cc}{\stackrel{\‾}{K}}_{\mathrm{\ξ}\mathrm{\ξ}}& {\stackrel{\‾}{K}}_{{\mathrm{\ξq}}_J}\\ {\stackrel{\‾}{K}}_{{\mathrm{\ξq}}_J}& {\stackrel{\‾}{K}}_{q_J{q}_J}\end{array}\right]$

The wherein corresponding immobile interface cylindrical coordinates of ξ, q _jeach connected node component corresponding;

By connected node motion vector transform in Oxr θ cylindrical coordinate,

$q_J^{(B,c)}=\tilde{T}{\hat{q}}_J^{(B,c)}$

Then the matrixing of blade reduced model list blade stiffness is:

${\hat{K}}_{b0}={\stackrel{\‾}{T}}^T{\stackrel{\‾}{K}}_{b0}\stackrel{\‾}{T}=\left[\begin{array}{cc}{\hat{K}}_{\mathrm{\ξ}\mathrm{\ξ}}& {\hat{K}}_{{\mathrm{\ξq}}_J}\\ {\hat{K}}_{{\mathrm{\ξq}}_J}& {\hat{K}}_{q_J{q}_J}\end{array}\right]$

$\stackrel{\‾}{T}=\left[\begin{array}{cc}{I}_{n_K}& 0\\ 0& \tilde{T}\end{array}\right]$

Wherein n _kfor blade retains mode number, sized by be n _kunit matrix; In like manner offset of sinusoidal motion vector carries out polycondensation conversion, gained polycondensation matrix with identical; Then n pitch diameter displacement model (n>0) the global stiffness matrix polycondensation of N bar blade is transformed to:

${\hat{K}}_b=\frac{N}{2}\left[\begin{array}{cc}{\hat{K}}_{b0}& \\ & {\hat{K}}_{b0}\end{array}\right]$

Interfacial displacement vector in No. 0 blade reduced model can be expressed as by wheel disc tie point motion vector:

$\left\{\begin{array}{c}{\hat{q}}_J^{(B,c)}=diag\left(I_6{\mathrm{cosn\β}}_i\right)q_J^{(D,c)}+diag\left(I_6{\mathrm{sinn\β}}_i\right)q_J^{(D,s)}\\ {\hat{q}}_J^{(B,s)}=diag(-I_6{\mathrm{sinn\β}}_i)q_J^{(D,c)}+diag\left(I_6{\mathrm{cosn\β}}_i\right)q_J^{(D,s)}\end{array}\right.$

Be designated as:

$T=\left[\begin{array}{cccc}{I}_{n_K}& & 0& \\ & diag\left(I_6{\mathrm{cosn\β}}_i\right)& & diag\left(I_6{\mathrm{sinn\β}}_i\right)\\ 0& & I_{n_K}& \\ & diag(-I_6{\mathrm{sinn\β}}_i)& & diag\left(I_6{\mathrm{cosn\β}}_i\right)\end{array}\right]$

Then have

$\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(c\right)}\\ {\hat{q}}_J^{(B,c)}\\ {\mathrm{\ξ}}^{\left(s\right)}\\ {\hat{q}}_J^{(B,s)}\end{array}\right\}=T\left\{\begin{array}{c}{\mathrm{\ξ}}^{\left(c\right)}\\ q_J^{(D,c)}\\ {\mathrm{\ξ}}^{\left(s\right)}\\ q_J^{(D,s)}\end{array}\right\}$

Global stiffness matrix then after all blade polycondensations is transformed to further:

$K=T^T{\hat{K}}_{b}T=\frac{N}{2}\left[\begin{array}{cc}{K}_c& K_s\\ K_s^T& K_c\end{array}\right]$

Wherein:

The mass matrix after 3D solid leaf model polycondensation conversion is obtained by identical method;

Finally by the stiffness matrix after conversion and mass matrix general assembly to the stiffness matrix and mass matrix of two-dimensional axial symmetric dish model, obtain the overall stiffness matrix of model and mass matrix; Utilize this dynamic coupling model corresponding to two matrixes, n pitch diameter displacement model 3D solid blade and two-dimensional axial symmetric wheel disc being become to dimensional model carries out model analysis; Recovered by mode, obtain 3D solid blade and the modal displacement of two-dimensional axial symmetric wheel disc under the n pitch diameter frequency of each rank;

For zero pitch diameter displacement model, use the same method and obtain the global stiffness after becoming dimensional model polycondensation and mass matrix, carry out zero pitch diameter model analysis.

6. change dimensional model finite element prestressed modal analysis method according to claim 2, it is characterized in that, in centrifugal force field, the computing method of 3D solid blade and two-dimensional axial symmetric wheel disc Coupling Deformation are as follows:

For the Coupling Deformation problem preloaded, refer in particular to impeller by preloading problem under centrifugal action, correspond to the zero pitch diameter displacement model becoming dimensional model, therefore No. 0 blade displacement vector only comprises q ^{(B, c)};

Consider the blade wheel structure of single group blade, using the interface point of the tie point of No. 0 leaf three-dimensional model as static(al) polycondensation, static(al) polycondensation is carried out to a blade static condensation method by centrifugal force, obtain the nodal force vector of polycondensation rear interface point and Stiffness Matrix by interface node motion vector be transformed into the cylindrical coordinate Oxr θ of two-dimensional axial symmetric wheel disc model modeling from the rectangular coordinate system of definition 3D solid leaf model, namely

$q_J^{(B,c)}=\tilde{T}{\hat{q}}_J^{(B,c)}$

$\tilde{T}={\mathrm{diag}}_{i=1}^{n_J}\left(T_i\right)$

Obtain convert after polycondensation stiffness matrix and force vector be:

$K_b={\tilde{T}}^T{\tilde{K}}_b\tilde{T}$

$P_b={\tilde{T}}^T{\tilde{P}}_b$

For two-dimensional axial symmetric wheel disc, carry out the force vector P that static(al) polycondensation obtains polycondensation rear interface point equally _dwith stiffness matrix K _d;

Then blade wheel structure becomes the interfacial displacement vector of dimensional model from following equation solution:

(K _d+NK _b)q _J＝P _d+NP _b

Wherein N is the circumferential number of blade, the internal displacement of blade and wheel disc then uses static condensation method according to q _jcalculate;

For the situation of multiple blade, then the rigidity and the force vector that superpose more unitized constructions calculate.