JPH0921720A - Method for analyzing vibration of structure - Google Patents

Abstract

PROBLEM TO BE SOLVED: To acquire a highly accurate response value in a short time by subjecting a plurality of vectors simultaneously to mass matrix or solving a stiffness matrix to obtain an eigenvector and then adding a vector for compensating the external force component. SOLUTION: Node data of various structures are obtained and the coordinates of nodes or the number of equation are outputted in order to obtain a mass matrix or stiffness matrix for a plurality of components (e.g. springs, trusses and pipes) of a structure. These matrixes are then superposed each other while decreasing the bandwidth thereof and the mass matrixes and stiffness matrixes in the entire structural system are assembled to make up a load spectrum of the entire system. When the load spectrum acts statically, deformation and stress are calculated, for example. The coefficients of load spectrum vary with time and the combined dynamic response is calculated thus obtaining the dynamic response value (vibration) of structure or the like.

Description

Detailed Description of the Invention

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to an analysis method for calculating the dynamic response of various devices and structures using an electronic computer,
Particularly, the present invention relates to a structural vibration analysis method for calculating vibration stress with high accuracy and in a short calculation time.

[0002]

2. Description of the Related Art Various machines and structures are modeled as a linear elastic body, and a response when an external force which changes with time is applied thereto, a dynamic response at the time of an earthquake, etc. are calculated by using a computer. Often. In such cases, first
It is common to calculate the natural vibration mode by eigenvalue analysis and use the result to calculate the dynamic response by a method called modal analysis. In that case, for many large-scale models that require a long calculation time, the calculation time of the eigenvalue analysis accounts for a large proportion of the total calculation time.

A method called a subspace method is used for eigenvalue analysis in structural vibration analysis (for example, 1979 Sep.
Published by Science and Technology Publishing Co., K. J. Bathe, E .;
L. Wilson, translated by Fumio Kikuchi, numerical calculation of finite element method, disclosed in p556 to 582) is considered to be effective, put to practical use, and a program by this method is commercially available.
On the other hand, a method called the Lanczos method is faster, and a paper by Fujikake, who claims to have solved the problems of the conventional Lanczos method and completed it as a practical algorithm (March 1989, published by the Japan Society of Mechanical Engineers, Japan Opportunity Society Proceedings, Volume 55, No. 511, p. 564-572, Fast solution of large-scale eigenvalue problem by Lanczos method). However, the explanation by that is difficult to understand, and it seems that Fujikake thinks that it can be made a tough solution that general users can use with peace of mind, and the effect is very attractive according to Fujikake's paper. Nevertheless, even after several years, ADI
It seems that commercial programs such as NA, ANSYS, and NASTRAN have adopted at least a standard solution method, and have not significantly reduced the calculation time for eigenvalue analysis.

On the other hand, with respect to modal analysis, problems of the conventional method and methods for solving them (Hatatake, Shigeta, correction method for higher-order mode deletion in modal analysis, January 1984, published by The Japan Society of Mechanical Engineers, Japan Opportunity) Academic Papers, Volume 50, Volume 449,
p11-16) have been proposed. However, although this problem is considered to be an important problem in practical use, it seems that it has not been adopted as a standard method or a feasible option in the above-mentioned commercially available programs even after 10 years have passed. . And both have been proposed individually, and no method or program for using both at the same time has been reported. That is, in the eigenvalue analysis by the Lanczos method, the vector to be added in order to execute the modal analysis with high accuracy is completely ignored. In addition, Hatake's paper (August 1985, published by The Japan Society of Opportunities, The Japan Society of Mechanical Engineers, Volume 51, Volume 468, p1897-
1905, the generalized coordinates considering the external force term in the dynamic response analysis (1st report, the basic concept and the calculation method of the proposed generalized coordinates) use the subspace method in the eigenvalue analysis.

As described above, in the conventional analysis method, even if the vector [Ψr] for calculating the vibration stress with high accuracy is used when only the low-order mode is used, the calculation efficiency of the eigenvector is high. However, even if a fast eigenvalue analysis method is used, the calculation accuracy is poor because the vector [Ψr] for calculating the vibration stress with high accuracy is not used. It was

[0006]

SUMMARY OF THE INVENTION An object of the present invention is to model a machine or structure as a linear elastic body and calculate the response when an external force that changes with time is applied, that is, its dynamic response. At this time, it is possible to obtain an accurate response value efficiently, that is, in a short calculation time.

[0007]

[Means for Solving the Problems] A method called modal analysis is often used in linear dynamic response calculation. At that time, eigenvalue analysis and eigenvector [Φ] are calculated, and the result is used. Calculate the dynamic response. In the present invention, a method called the Lanczos method, which is evaluated to be a high-speed method for eigenvalue analysis, or a method similar thereto is used. Subsequently, the vector [Ψr] or [Ψ] required to execute the response calculation with high accuracy.
h] is calculated. The dynamic response calculation is performed using both [Φ] and [Ψr] or [Φ] and [Ψh] to always obtain a highly accurate solution. Obtaining a highly accurate response value among the above objectives is achieved by using both [Φ] and [Ψr] or both [Φ] and [Ψh] in the mode analysis, and a solution can be obtained in a short calculation time. This is achieved by using the Lanczos method or a similar method for the eigenvalue analysis.

[0008]

In the present invention, in the eigenvalue analysis, the process of multiplying one or a small number of vectors by the mass matrix and solving with the stiffness matrix is repeated. Therefore, a method for simultaneously processing a large number of vectors, such as the subspace method, is used. By comparison, eigenvectors can be obtained in a short calculation time in eigenvalue analysis. At the same time, a vector is added to compensate for the component of the external force that would fall out only with the eigenvector, so that a highly accurate response value can be obtained in the dynamic response calculation.

[0009]

FIG. 1 shows an embodiment of the present invention. Step 1
Then, create the nodal data, and output the nodal coordinates, nodal equation numbers, etc. to a file. In step 2, the mass matrix and stiffness matrix of the elements that make up the structure are output to a file. Details of step 2 are shown in FIG.

In step 2, steps 21 to 29 are repeated N times. Here, N is the number of element blocks. In step 21, a character or number indicating which type of element the element block is, and a numerical value indicating the overall property of the element block may be read. If it is a spring element according to the data read in step 21, a rigidity matrix is calculated for the spring element and output to a file in step 22.

If it is a truss element according to the data read in step 21, a mass matrix, a rigidity matrix, etc. are calculated for the truss element in step 23 and output to a file. Similar to step 23, steps 24 to 28 are executed when the type of each element is selected in step 21.

Step 29 is executed when the superelement is selected in step 21. In this step, unlike steps 22 to 28, the superelement mass matrix and rigidity are previously stored in the superelement file in step 16 of FIG. The matrix must be created. And
In this step 29, the posture of the super-element in the coordinate system when the super-element is created and the posture in the coordinate system when the super-element is incorporated into the whole system are generally different. Coordinate conversion.

Returning to FIG. 1, in step 3, the equation numbers of the nodes are changed in order to reduce the band width (band width) of the rigidity matrix and the mass matrix. In step 4, the mass matrix and rigidity matrix of the elements created in step 2 are overlaid to assemble the mass matrix and rigidity matrix of the entire structural system. At that time, the stiffness matrix and the mass matrix are the matrix (band matrix) in which all the components distant from the diagonal to some extent are zero, so only the inside of the band width is used, and since they are symmetrical matrices, only the upper right or lower left part is filed. Suffice it to store in. Further, when the element is the element of step 13-19 in FIG. 2, in a large-scale system, most of the components are zero even within the bandwidth, and the non-zero components are often less than a few percent, so the value of the non-zero component is It suffices to store in the file only the information needed to recover the row and column numbers of the and its components.

In step 5, a load vector for the entire system is created. In the case of static analysis, the deformation and stress when the load vector created here acts statically are calculated in steps 6 and 7. In the case of modal analysis, the coefficient applied to the load vector created here changes with time, and the dynamic response when an external force, which is the sum of them, acts is calculated in Step 10.
Calculate with ~ 12. Further, in the case of creating a super element, the load vector created here becomes the element load vector of the super element, and the response when an external force acts on the inside of the super element can be calculated.

In the case of mode analysis, steps 8 and 9 are executed, and steps 10, 11 and 12 are executed only when it is instructed by the input data.

Details of step 8 are shown in FIG. First, N1 and N3 are determined. N1 is how many times the process of step 35 is performed. This number may be zero. If it is zero, step 35 is not executed. N3 is how many vectors are generated by the processing of steps 37-40. [U] in step 40 becomes a matrix of N3 columns, and [K _w ] in step 41,
[M _w ] becomes a matrix having a size of N3 × N3.
Initially, N2 is set to 1, and the processes of steps 31 to 43 are repeated until all necessary eigenvectors are obtained.

At least initially (j = 1), the stiffness matrix [K] is triangulated. When moving the origin, in step 31, [K] -μ _s [M] is substituted into [K ^* ].
When such [K ^* ] is created and the processes of steps 33 to 43 are performed, an eigenvector whose eigenvalue is close to μ _s is obtained. When the origin is not moved, in step 32, [K] is directly substituted into [K ^* ]. At that time, the eigenvectors with lower natural frequencies are sequentially obtained. Next step 33
Triangularly decompose [K ^* ] with.

Triangular decomposition is the calculation of [L] and [D] such that [K ^* ] is the product of [L], [D] and [L] ^T. [L] is called a lower triangular matrix, and is a matrix in which all the diagonal elements are 1, the upper right half is all zero, and the lower left half has nonzero elements. [D] is a diagonal matrix in which all the components other than the diagonal components are zero. [L] ^T is the transposed matrix of [L].

In the first (j = 1) loop, step 34
The initial vector {v ₁ } is created with. This vector must be defined so as to include all the components of the eigenvector to be obtained. When N1 is 1 or more, the process of step 35 is performed N1 times. In step 36, {v
Substitute _{N1 + 1} ) into {u}. When N1 is zero, {v ₁ } created in step 34 is directly substituted into {u}.

Next, in steps 37 to 40, i is set to N2.
To N3, increment by one and execute. In step 37

[Equation 1] To calculate {v}. The matrix on the right-hand side and the vector multiplication can be calculated easily. After the right side is represented by one vector, {v} can be calculated with a small number of operations if [K ^* ] is triangulated in step 33.

If step 38 is omitted and the processing of step 37 is repeated, {v} gradually converges to the eigenvector closest to μ _s , so that they become nearly parallel vectors. As a result, the eigenvalue analysis in step 42 becomes numerically difficult. Therefore, in step 38, {v} is changed to {u
₁ } to {u _i-1 ) and [M] are orthogonalized. This process is as i increases from N2 to N3,
It need not be performed every time, and may be performed every several times. Here, {u _i } is the vector in the i-th column of [U].

Step 39 is called normalization, where {v} is multiplied by a constant.

[Equation 2] Adjust the size so that If this process is not performed, the absolute value of the numerical value of each component of the vector will be gradually increased or decreased by the process of step 37. This process also need not be performed every time as in step 38, and may be performed every several times. However, since most of the calculation time is often occupied by the processing of step 37, even if step 38 or step 39, especially step 39 is executed every few times, the total calculation time is not saved so much. It may be.

At step 40, {u} obtained by the processing at steps 37 to 39 is substituted into the i-th column of [U]. In step 41, [U] is used to [K ^* ] and [M].
To a subspace to create a reduced stiffness matrix [Kw] and mass matrix [Mw] (N3 × N3 matrix).

At step 42, eigenvalue analysis is performed on [Kw] and [Mw] to obtain all N3 eigenvectors to obtain an eigenmatrix [Q]. In FIG. 3, the expression is such that the general eigenvalue problem is solved, but if the processing of steps 38 and 39 is performed for all i, [Mw] becomes a unit matrix, so in step 42 the standard eigenvalue problem

(Equation 3) You can solve it. At that time, in step 41, [Mw]
The calculation of may be omitted. Although step 42 is a small-scale eigenvalue problem, it is necessary to calculate all eigenvalues and eigenvectors, and therefore the Jacobian method may be used.

At step 43, [U] is multiplied by [Q] to calculate [W]. Error norm as eigenvector for each column vector of [W]

(Equation 4) Is calculated, and if it is sufficiently small, it is regarded as an eigenvector and is substituted into [U] in order from the left. At this time, the error norm calculation may be omitted for the vector that has already been determined to be the eigenvector in the previous j. A value obtained by adding 1 to the obtained number of eigenvectors is set in N2. And
Some of the vectors not assigned to [U] as eigenvectors are added up and assigned to {v ₁ }.

Then, at the end of the loop, if all the required eigenvectors have been obtained, the loop is exited (EXI
T). The above operation is repeated until all eigenvectors are obtained. This is the end of the description of step 8 in FIG. 3 or FIG.

In FIG. 3, in steps 35, 37, etc., one vector is multiplied by [M] and is solved by [K ^* ]. However, in order to be advantageous to the conventional subspace method, the operation is performed on a small number of vectors, and step 41 is performed using a larger number of vectors.
It is necessary to project onto the subspace of. Therefore, the most efficient solution is to repeat steps 35 and 37 for one vector as shown in FIG. 3, but even if it is repeated for two or three vectors, the conventional subspace method is used. Expect to be faster.

Next, the processing in step 9 of FIG. 1 will be described with reference to FIG. In step 51, the component [Fh] that excites only the higher-order mode is extracted from the external force vector [F] by extracting the component that excites the lower-order mode. In step 52, as a result of the operation in step 51, the ones that are close to the zero vector are removed and packed. In step 53

(Equation 5) To solve for [U]. In the Hansteen method, [U] obtained here is used as [Ψh] in the response calculation. According to Hatake's method, in step 54, [U] is used to project [K] and [M] onto the subspace, and [Kw] and [M] are projected.
w] is calculated. In step 55, the eigenvalue analysis in the subspace is performed to obtain the eigenvalue [Ω ² ] and the eigenvector [Q]. In step 56, [U] is multiplied by [Q] to calculate the compensation vector [Ψr]. Then, this is used for response calculation. To calculate [Ψr], in addition to the method shown here, Guram in the paper by Hatake (August 1985)
It can also be calculated by orthogonalization of -Schmidt.

Further, in a structure floating in the air such as an airplane, or in a system having an incomplete constraint and an eigenvalue of zero, the stiffness matrix cannot be calculated as it is as [K] in step 53 and step 54. , Move the origin, find [K] -μ _s [M], and set it to [K ^* ]
And Then, in steps 53 and 54, [K ^* ] is used instead of [K]. At this time, it is preferable to use a value of μ _{s that} is as small as possible so that the calculation is not unstable.

Returning to FIG. 1, in steps 8 and 9, the eigenvector [Φ] and the compensation vector [Ψh] or [Ψr] are used instead of the eigenvector in the case of the conventional modal analysis using the eigenvectors, instead of the eigenfrequency. If the equivalent frequency and the generalized load are obtained, the time history response analysis of step 10 and the frequency response analysis of step 11 can be performed in exactly the same manner as the conventional method. In the response spectrum analysis in step 12, the response to the compensation vector is taken as one vector by taking the algebraic sum of the response values of each vector, and it is combined with the response of the true eigenvector in the same manner as the conventional method as the response of one eigenvector. Good.

Next, steps 13 to 16 of FIG. 1 will be described. Steps 13 to 16 are executed when the analysis mode is the super element creation mode. Regarding the calculation using superelements, a detailed description is given in JP-B-5-46891. Since the present invention is basically the same as the superelement in its contents, the explanation using the formula is omitted. However, in the present invention, the calculation of the standard mode (natural vibration mode when the boundary degrees of freedom are constrained) in step 13 is executed by the method shown in FIG. By using this method,
It can be calculated in a shorter calculation time than in the case of the embodiment of FIG. The calculation of the compensation vector in step 14 is executed by the method shown in FIG. Also for this part, by using this method,
It can be executed in a shorter calculation time than in the case of the illustrated embodiment. In the calculation of the constraint mode in step 15, only the non-zero component of the stiffness matrix is read from the file and the calculation is performed. Also in the superelement definition processing in step 16, only nonzero components in the stiffness matrix and mass matrix are read from the file, and calculations are performed using only nonzero components.

An example of calculation using the program according to the present invention will be shown. The dimensions of the calculation model are shown in FIG. This model is a model created such that the calculation result can be predicted and the effect of the eigenvalue analysis method (FIG. 3) included in the present invention can be sufficiently verified. When the load F as shown in FIG. 5 is applied, it may be calculated only for one rod.
First, the model in FIG. 5 is calculated by static analysis to calculate a reference solution. In this case, the static analysis, that is, steps 1, 2, 4, 5, 6, and 7 in FIG. 1 is executed. As a result, the distribution of the shearing force of the member A is obtained as shown in FIG.

Next, the response when the external force F temporally changes as shown in FIG. 7 is calculated without using the superelement. For that purpose, in FIG. 1, time history response analysis by mode analysis, that is, steps 1, 2, 4, 5, 8, 9, and 10 are executed. As a result, the distribution of the shearing force of the member A at time t = 10 seconds becomes the distribution shown in FIG. In the case of time-dependent change of force as shown in FIG. 7, it has to agree with FIG. 6, and according to the present invention, it can be confirmed that the solution is correctly obtained. On the other hand, when the conventional method, that is, step 9 in FIG. 1 is omitted and only the eigenvector is calculated by the subspace method and the calculation is performed using the eigenvector, the result shown in FIG. 9 is obtained. This is different from the result shown in FIG. 6 and it can be seen that the solution cannot be obtained correctly. Further, in step 8 of FIG. 1, when eigenvalue analysis is performed using FIG. 3 of the present invention, the eigenvector can be obtained in 346 seconds, but when the subspace method of Japanese Patent Publication No. 5-46891 is used, 1500 seconds of calculation time is required.
That is, if the method of Japanese Patent Laid-Open No. 63-231233 is used, the same calculation result as in FIG. 8 can be obtained, but the calculation time of the eigenvalue analysis is long. Therefore, according to the present invention, a highly accurate solution can be obtained in a short calculation time.

Then, using the superelement, the same problem as above is solved. The model of FIG. 5 can be created by assembling three structures shown in FIG. Therefore, the structure of FIG. 10 is defined as one super element, and the model of FIG. 5 is solved using it. First, in FIG. 1, superelement creation, that is, steps 1, 2, 4, 5, 13, 14, 15, and 16 are executed for the structure of FIG. Then, using the result, the mode analysis, that is, steps 1, 2, 4, 5, 8, 9, and 10 of FIG. 1 is executed for the model of FIG. This time,
In step 2, step 29 of FIG. 2 is executed.

If the step 14 of FIG. 1 is not omitted in the present invention, that is, the superelement formation, and the step 9 of FIG. 1 is not omitted in the mode analysis of the entire system, the same shear force distribution as that of FIG. 8 is obtained. However, when both steps 9 and 14 in FIG. 1 are omitted, the shear force distribution is the same as that in FIG. Further, when only step 14 of FIG. 1 is omitted in the super element creation, even if step 9 of FIG. 1 is executed in the calculation of the whole system, the shear force distribution as shown in FIG. I can only get it.

[0036]

Since the present invention is configured as described above, it has the following effects. In the convergence calculation for calculating the eigenvectors, the operation of multiplying one or a small number of vectors by the mass matrix and solving with the stiffness matrix is repeated sequentially, so that a vector including many low-order eigenvector components can be generated with a small number of calculations. Since the eigenvalue analysis is performed in the subspace spanning these vectors, low-order eigenvectors can be calculated efficiently. Furthermore, since a dynamic response calculation is performed by adding a compensation vector for compensating a portion of the external force vector that cannot be represented by the eigenvector, a highly accurate dynamic response calculation result can be obtained.

In addition, as an element constituting the entire structure, 2
Since it is possible to include an element called a super element that combines one or more elements into one element, it is possible to efficiently perform vibration analysis of a large-scale or complex structure.
Furthermore, it is possible to accurately calculate the response when an external force is applied to the inside of the superelement.

[Brief description of drawings]

FIG. 1 is an overall flow of an embodiment of a program according to the present invention.

FIG. 2 is a flow of an element data creation part of an embodiment of a program according to the present invention.

FIG. 3 is a flow of an eigenvalue analysis part of an embodiment of a program according to the present invention.

FIG. 4 is a flow of a compensation vector calculation part of an external force vector of an embodiment of a program according to the present invention.

FIG. 5 is dimensional data of a model of a calculation example.

FIG. 6 is a calculation result of static analysis of a calculation model.

FIG. 7 is a temporal change of an external force in a time history response analysis of a calculation model.

FIG. 8 is a calculation result of time history response analysis according to the present invention.

FIG. 9 is a calculation result when an external force vector compensation vector is not used.

FIG. 10 is a dimensional specification of a super element.

FIG. 11 is a calculation result when the present invention is not used in the super element creation mode.

[Explanation of symbols]

 8,13 Steps for calculating eigenvectors

Claims (2)

[Claims] 1. A method of calculating a low-order eigenvector or the like, and using the result to analyze a dynamic response such as displacement or stress of a structure or a machine when an external force varying with time is applied, , multiplied by the mass matrix [M] to a vector of N present at the same time in the step of calculating the eigenvectors, performed a plurality of times a process of solving a rigid matrix [K] or origin movement to the resulting matrix [K ^*], their When the J eigenvectors [Φ] are obtained using the J vectors obtained from the processing, L is larger than N and J is larger than or equal to L, and the eigenvectors of lower order among the given external force vectors Find a compensation vector to compensate for the part that cannot be expressed by and calculate the dynamic response to the given external force vector using the vector and the low-order eigenvector. A structural vibration analysis method characterized by the following. 2. An element called a super element, which is formed by assembling two or more elements into one element, as an element constituting the entire structure in claim 1, can be included. When calculating the stiffness matrix and mass matrix, as the mode corresponding to the standard mode,
A structural vibration analysis method characterized by using [Φ] and [Ψr] or [Φ] and [Ψh] of claim 1 calculated from a mass matrix, a stiffness matrix and an element load vector.