JP4154271B2 - Structure optimum design method, program, and storage medium - Google Patents

Description

【００４０】
上式中
ａ_i,j = (λ+2μ)＜∂_xe_i,∂_xe_j＞ + μ＜∂_ye_i,∂_ye_j＞ …（式８）
ｂ_i,j = λ＜∂_xe_i,∂_ye_j＞ + μ＜∂_xe_j,∂_ye_i＞ …（式９）
ｃ_i,j = λ＜∂_xe_j,∂_xe_i＞ + μ＜∂_xe_i,∂_ye_j＞ …（式１０）
ｄ_i,j = (λ+2μ)＜∂_ye_i,∂_ye_j＞ + μ＜∂_xe_i,∂_xe_j＞ …（式１１）
ここで、ｅ_jはノードｊに対する基底ベクトル、＜ , ＞はベクトルの内積を表す。∂_xはｘに関する偏微分、∂_yはｙに関する偏微分である。ξは前述の（非特許文献４）に従い、２とする。また、λ、μはラメ定数と呼ばれる材料定数であり、ヤング率Ｅ及びポアソン比νより次式で計算される。
【００４１】
λ = Eν/(1+ν)(1-2ν) …（式１２）
μ = E/{2(1+ν)} …（式１３）
（式７）より容易に分かるように、ある要素の存在率が０となると要素剛性行列Ａ_jのすべての成分が０となる。全体剛性行列は要素剛性行列の重ね合わせで構成される。従って、ノードjを含む要素の設計変数ｆ_jが０になると、全体剛性行列のノードｊに対応する行及び列の全ての成分が０になる。このように、構造最適化問題において全体剛性行列は一般に特異な行列となる。
【００４２】
設計変数に関する拘束条件は、構造最適設計の場合、総重量の上限値Ｗ₀とする。
【００４３】
【数２】

【００４４】
また、設計変数には値域に対する拘束条件が以下のように定められている。
【００４５】
0 ≦ ｆ_j ≦ 1 j = 1, … , n …（式１５）
以上の表記より、最適設計問題は、拘束条件（式１４）及び（式１５）の下で、（式５）を極小化する問題として定式化される。
【００４６】
また、状態変数ｘについては、拘束条件（式３）の下で（式４）を極小化する解として得られる。
【００４７】
＜本実施形態の構造最適設計装置の動作例＞
上記定式化に基づき、本実施形態の処理手順例を説明する。
【００４８】
図３には、本実施形態の流れ図が示されている。
【００４９】
図中、ステップＳ１０１はシミュレーション対象となる系の諸元を読み込み、変数の初期化を行う処理である。諸元の読み込みは、入力装置２０３或いは通信装置２０６からの入力データで行っても良いし、予め２次記憶装置２０５にファイルとして格納しておいたデータを読み出して利用することもできる。系の諸元にはｘ、ｆの初期値としてｘ⁽⁰⁾、f⁽⁰⁾、境界条件Ｂ、評価関数Ｌ₁、Ｌ₂が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保して、値を設定する。感度係数の初期値Ｃ_j ⁽⁰⁾は次式で計算する。
【００５０】
Ｃ_j ⁽⁰⁾ = - (ｘ_j ⁽⁰⁾)^T(ξ(ｆ_j ⁽⁰⁾)^{ξ -1}Ａ_j)ｘ_j ⁽⁰⁾ …（式１６）
ステップＳ１０２では、反復回数ｋを１に初期化する。反復回数ｋは、上述の設計変数の求解工程に対する反復処理に対するものであり、予め定められた終了条件（評価関数Ｌ₁，Ｌ₂の条件）を満たしたら終了する。
【００５１】
ステップＳ１０３では、ラグランジェ未定定数λの更新を行う。ステップＳ１０４では、設計変数ベクトルｆの更新を行う。ステップＳ１０５では、状態変数ベクトルｘの更新を行う。ステップＳ１０６では、感度係数の計算を行う。
【００５２】
上記ステップＳ１０３及びＳ１０４は、拘束条件付最適化問題の解法として採用する方法によって異なる。例えば、この設計変数ベクトルｆ^(k)を算出する求解工程は、逐次線形計画法、最適性規準法、逐次凸関数近似法のいずれかにより実現される。以下、解法として最適性規準法を用いる場合と、逐次凸関数近似法を用いる場合について詳述する。
【００５３】
（最適性規準法）
まず、最適性規準法を用いる場合について説明する。構造最適化問題に対する最適性規準法は（非特許文献４）及び（非特許文献５）等によって公知である。
【００５４】
最適性規準法では、（式４）及び（式５）をラグランジェ未定定数λを用いて、１つの式にまとめる。
【００５５】
Ｌ(f,λ) = Ｌ₂(f,ｘ) - λ(Ｗ(f) - Ｗ₀) …（式１７）
ステップＳ１０３で実行されるλの更新は、次式を用いて行われる。
【００５６】
λ^(k+1) = min[0,[(1/Ｗ₀)Ｗ(f^(k))] ^αλ^(k)] …（式１８）
ただし、上付きのｋ，ｋ＋１は反復回数を表す。また、αはべき乗係数で、通常0.85程度の値を設定する。
【００５７】
ステップＳ１０４で実行される設計変数の更新は次式による。
【００５８】

ただし、ζは設計変数更新の変動幅の制約値で、0.3程度の値を取る。また、ｓ_j ^(k)は次式で与えられる。
【００５９】
ｓ_j ^(k) = [(1/λ^(k-1))Ｃ_j ^(k-1)]αf_j ^(k-1) …（式２０）
ここで、Ｃ_j ^(k)は第２の評価関数Ｌ₂の、設計変数ｆ_j ^(k)に関する感度係数と呼ばれる量であり、ステップＳ１０６の処理で計算されるものである。
【００６０】
更に（式１９）で計算された設計変数ｆ_j ^(k)は、予め設定された値より小さくなると強制的に０にしてもよい。この予め設定された値とは、例えば１０^-3程度の値である。
【００６１】
ステップＳ１０５では、状態変数ベクトルｘ^(k)を更新する。
【００６２】
状態変数ベクトルｘ^(k)は、（式４）の評価関数の極小化の過程で得られる。この解法は共役残差法と呼ばれており、後述する。
【００６３】
ステップＳ１０６では、感度係数Ｃ_j ^(k)を次式により計算する。
【００６４】
Ｃ_j ^(k) = - (ｘ_j ^(k))^T(ξ(ｆ_j ^(k))^{ξ -1}Ａ_j)ｘ_j ^(k) …（式２１）
（逐次凸関数近似法）
次に、逐次凸関数近似法を用いた場合の、ステップＳ１０３及びＳ１０４の処理について詳述する。逐次凸関数近似法では設計変数ベクトルの各成分に対してスケーリングを行うことによって、数値計算誤差を避けることができるが、以下では説明を簡単にするため、スケーリングについては記述していない。逐次凸関数近似法の詳細は、（参考文献１）及び（参考文献２）等によって公知である。
【００６５】
（参考文献１）藤井：“パソコンで解く構造デザイン”, 丸善(2002年4月)
（参考文献２）C. Fleury : "CONLIN, an efficient dual optimizer based on convex approximation concepts," Structural Optimization, 1, pp.81-89 (1989)
ステップＳ１０３では、ラグランジェ未定定数を次式により更新する。
【００６６】
λ^(k) = λ^(k-1) + Δλ^(k) …（式２２）
ただし、初期値λ⁽⁰⁾は例えば0.01程度の値を設定する。
【００６７】
Δλ^(k)は次の方程式の解として与えられる。
【００６８】
Δλ^(k) = ２×( 1 - Ｗ₀/Ｗ(f^(k-1)) )λ^(k-1) …（式２３）
ステップＳ１０４では、設計変数ベクトルの更新を行う。
【００６９】

更に、（式２４）で計算された設計変数ｆ_j ^(k)は、予め設定された値より小さくなると強制的に０にしてもよい。この予め設定された値とは、１０^-3程度の値である。
【００７０】
ステップＳ１０５及びＳ１０６の処理は前述のとおりである。
【００７１】
（共役勾配法）
ステップＳ１０５において状態変数ベクトルｘ^(k)を算出する求解工程として、共役勾配法を用いた処理を図４を用いて説明する。
【００７２】
ステップＳ３０１では、シミュレーション対象となる系の諸元を読み込む。系の諸元にはｘ^(k)の初期値ｘ⁽⁰⁾、設計変数ベクトルｆ^(k)の初期値ｆ⁽⁰⁾、変位拘束条件Ｂ₁、評価関数Ｌ₁が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保して、値を設定する。また、残差ベクトルｒ、探索方向ベクトルｐを初期化する。それぞれの初期値をｒ⁽⁰⁾、ｐ⁽⁰⁾と表す。
【００７３】
まず、ｒ⁽⁰⁾は、ｘ⁽⁰⁾と各ノード（節点）に加わる力をベクトル形式で記述した節点力ベクトルｂとを用いて次式で計算される。
【００７４】
ｒ⁽⁰⁾ = ｂ - Ａｘ⁽⁰⁾ …（式２６）
探索方向ベクトルｐの初期値ｐ⁽⁰⁾はｒ⁽⁰⁾と等しくする。
【００７５】
ステップＳ３０２では前処理を行う。ます、全体剛性行列Ａの対角成分のうち、０となる行あるいは列番号を調べる。０となる成分がなければ何もせずにステップＳ３０３に進む。０となる成分があれば、ｂ（節点力ベクトル）のうち一致する番号の成分を０にする。
【００７６】
ステップＳ３０３では、反復回数ｔを１に初期化する。反復回数ｔは、上述の状態変数の求解工程に対する反復処理に対するものであり、以下のステップＳ３０４からＳ３０８の処理の反復回数が設定値を超えるか、残差ベクトルｒ^(t)のノルムの２乗（ｒ^(t)，ｒ^(t)）が予め設定された所定値より小さくなったら、反復を終了する。以下、反復回数ｔ回目の値をｘ^(t)のように書く。
【００７７】
ステップＳ３０４では状態ベクトルの更新係数α^(t)を計算する。
【００７８】
α = (r^(t-1),p^(t-1))/(p^(t-1),Ａp^(t-1)) …（式２７）
ステップＳ３０５では、状態変数ベクトルｘを更新する。
【００７９】
ｘ^(t) = ｘ^(t-1) + αｐ^(t-1) …（式２８）
ステップＳ３０６では、残差ベクトルｒ^(t)を次式で計算する。
【００８０】
ｒ^(t) = b - Ａｘ^(t-1) …（式２９）
ステップＳ３０７では、探索方向ベクトルｐの更新係数βを計算する。
【００８１】
β = (ｒ^(t),ｒ^(t))/(ｒ^(t-1),ｒ^(t-1)) …（式３０）
ただし、ｔ＝１のときはβ＝１とする。
【００８２】
ステップＳ３０８では、探索方向ベクトルｐを更新する。
【００８３】
ｐ^(t) = ｒ^(t) + βｐ^(t-1) …（式３１）
共役勾配法及びその様々な派生法は、（参考文献３）等によって公知である。
【００８４】
（参考文献３） 森，杉原，室田：“線形計算，”岩波講座 応用数学，岩波書店(１９９４年)
（本実施形態の構造最適設計の具体例）
本例では、任意の位置に加重を受ける片持ち梁の最適形状自動設計に、本実施形態を適用するものである。説明を簡単にするために平面歪問題に限定する。
【００８５】
図５に示すように、構造部材の存在を可能とする設計領域は長方形４０２であり、有限要素法に従って、該領域を縦ｎ_y、横ｎ_xに等間隔に分割する。分割された部分領域をセルと呼び、左下及び右上のセルがそれぞれ(１,１)及び(ｎ_y,ｎ_x)となるように番号付けを行う。
【００８６】
同様に格子点をノードと呼び、左下及び右上のノードが(１,１)及び(ｎ_y＋１,ｎ_x＋１)となるように番号付けを行う。
【００８７】
図中、４０１は支持部材、４０３は荷重ベクトルである。
【００８８】
セル(j,k)には特性関数値ｆ(j,k)が対応する。ここで、特性関数値とはセル(j,k)における構造部材の存在確率を示す０から１の正の実数値をとる変数であり、本実施形態における設計変数ベクトルｆの要素である。
ｆ = (f(1,1), f(1,2), …, f(n_y,n_ｘ))^T …（式３２）
同様にノード(j,k)には横方向変位ｕ(j,k)と縦方向変位ｖ(j,k)が対応する。これらは任意の値を取る実数であり、本実施形態における状態変数ベクトルｘの要素である。
ｘ = (u(1,1), v(1,1), u(1,2), v(1,2),…,u(n_y+1,n_x+1), v(n_y+1,n_x+1))^T…（式３３）
図６に、本具体例に置ける計算結果を示す。図中、黒い領域が構造部材が存在する領域である。本具体施例において、構造部材の設計領域のアスペクト比は縦/横が２対１であり、その解析解は、水平方向に対して±４５度の梁が組み合わせられたものであることが知られている。図６に示した計算結果は、このような解析解とよく一致していることがわかる。
【００８９】
尚、本実施形態では、図１のように、バスライン２０７に各装置が接続した１つのシステムとして説明しているが、これら各要素をネットワークなどを介して接続して複数の装置のシステムで実現することも、あるいは汎用のコンピュータを複数接続して並列に多重処理あるいは分担処理により処理するようにしても良い。
【００９０】
又、本発明は、上記構造最適設計方法を実現するプログラム、及び該プログラムをコンピュータ読み出し可能に記憶する記憶媒体をも含むものである。
【００９１】
【発明の効果】
以上説明したように、本発明によれば、全体剛性行列が特異になったときでも、特別な処理を施すこと無しに計算を実行することができ、これによってプログラムが簡単になり、更に計算量の軽減が実現できる。
【図面の簡単な説明】
【図１】本実施形態の構造最適設計装置の構成例を示す図である。
【図２】図１の１次記憶装置及び２次記憶装置の記憶構成例を示す図である。
【図３】本実施形態の処理手順例を示すフローチャートである。
【図４】図３のステップＳ１０５を更に詳細に示すの処理手順のフローチャートである。
【図５】本実施形態に適応する具体例の問題設定の説明図である。
【図６】図５に係る本実施形態による計算結果を示す図である。[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a structure optimum design method及beauty programs, storage medium for shape optimization, including the topology of optimum design, particularly structural members of the structure which is the solution of the optimization problem design parameters.
[0002]
[Prior art]
The optimal design of the structural topology is a problem of determining the optimal topology and geometric dimensions of the structural member under given conditions. Hereinafter, the topology and shape dimensions of the structural member are referred to as a design variable function, and the determination problem is referred to as a design variable function optimization problem. The reason for the variable function is that the topology and shape dimensions are functions of a three-dimensional space. In the design variable function optimization problem, the state variable function optimization problem must be solved for each design variable function value. In this sense, the optimal design of the structural topology can be regarded as a double structure optimization problem having an optimization problem of a state variable function on the inner side and an optimization problem of a design variable function on the outer side.
[0003]
The inner state variable function optimization problem employs the idea of dividing the space into a finite number of elements by accumulating the prior art. In particular, in the problem where the strain energy of the structural member is an evaluation functional, the finite element method is generally used as the analysis method. A direct method for linear equations is adopted as a solution method of the finite element method.
[0004]
On the other hand, the design variable function optimization problem is roughly classified into the following three methods (see Non-Patent Document 1 or Non-Patent Document 2).
1. Evolutionary method (E method)
2. Homogenization method (H method)
3. Material distribution method (MD method) or Density method (D method)
In the E method, each of the partial spaces obtained by dividing the space is referred to as a cell, and generation and disappearance are repeated according to appropriate rules. The structural member is given as a collection of cells that are finally present. By allowing only two states, whether or not a cell is present, a clear structural member is obtained. In addition, since the differential information of the evaluation functional is not used, it is not trapped by the local optimum solution, so that it is effective when the evaluation functional is multimodal. In patent document 1, the optimization design apparatus of the frame structure member using the genetic search method which is a kind of E method is provided. This optimized design device requires trial calculation based on accumulated know-how, and therefore solves the problems of the prior art that could not cope with the actual design problem with a large number of design variables by the following method. ing. That is, an approximate optimization calculation device that uses an approximate expression of discrete design variable data such as a cross-sectional dimension of a frame member and a detailed optimization calculation device that uses the design variable data are provided, and these two calculation devices are combined to provide a framework. It constitutes an optimal design device for the structure.
[0005]
In the H method, it is possible to adopt sensitivity analysis by newly introducing a design variable function that assumes a finer structure and introduces a continuous value in a structural element located in each divided partial region. Here, sensitivity analysis is an analysis method that uses differential information of an evaluation functional related to a design variable function. If sensitivity analysis becomes possible, an iterative solution method such as a gradient method can be used. Compared to the brute force method as described above, at least the calculation time for searching for the local optimum solution is significantly shortened (see Non-Patent Document 3).
[0006]
The MD method or the D method is a method of expressing the topology or shape dimension change of the structural member by assigning each element a real number in the range of 0 to 1 indicating the existence rate of the structural member. This is the same as the H method in that the sensitivity analysis is made possible by replacing the discrete information on whether or not there is a structural member with a continuous value of the presence rate, but with a smaller number of parameters than the H method. Modeling is easy and the application range is wide.
[0007]
Non-Patent Document 4 discloses a structure phase optimization method based on the D method. The features of this method are as follows. Since the voxel finite element method (dividing the space into equal intervals) is used, the element stiffness matrix for all elements is the same. Therefore, if the element stiffness matrix is calculated once in advance, it can be used for subsequent calculations. Furthermore, since the elements are regularly arranged, it is not necessary to store node number information of each element. In order to solve large-scale simultaneous linear equations, the combined gradient method with preprocessing and the element-by-element method are used in combination, so the solution can be obtained without assembling the entire stiffness matrix. Less is enough.
[0008]
The homogenization method requires six design variables (one-dimensional case) for one element. In addition, the element stiffness matrix must be recalculated each time the design variable changes. On the other hand, one design variable may be used for one element by adopting a density method that expresses the abundance of the structural member as a density ratio. Also, changes in design variables do not affect the element stiffness matrix.
[0009]
[Patent Document 1]
JP 11-314631 [Non-Patent Document 1]
S. Bulman, J. Sienz, E. Hinton: "Comparisons between algorithms for structural topology optimization using a series of benchmark studies," Computers and Structures, 79, pp.1203-1218 (2001).
[Non-Patent Document 2]
YL. Hsu, MS. Hsu, CT. Chen: "Interpreting results from topology optimization using density contours," Computers and Structures, 79, pp.1049-1058 (2001).
[Non-Patent Document 3]
Hiroshi Yamakawa: “Optimization Design,” Computational Mechanics and CAE Series 9, Baifukan (1996)
[Non-Patent Document 4]
Fujii, Suzuki, Otsubo: “Topological optimization of structures using voxel finite element method,” Transactions of JSCES, Paper No. 20000010 (2000).
[0010]
[Problems to be solved by the invention]
However, the above prior art has the following problems.
[0011]
In general, the structure optimization problem is formulated as a double optimization problem including a state variable vector optimization problem at each iteration step of the design variable vector optimization problem. When the optimization problem of the design variable vector is called outer optimization and the optimization problem of the state variable vector is called inner optimization, the inner optimization takes the design variable vector as a parameter, that is, the design variable vector is fixed, This is a problem of obtaining a state variable vector. This is usually called structural analysis, and can be solved using a linear equation solving method by the finite element method.
[0012]
However, when the structure is changed and there is no structural member in a certain region, for example, when a hole is formed in the member, the design variable of the corresponding element becomes 0, and as a result, the Young's modulus of the element becomes 0. Then, the coefficient matrix of the linear equation is not full rank, and the inverse matrix cannot be calculated, so it cannot be solved by the direct method.
[0013]
Therefore, in the prior art, when the design variable vector becomes 0, the following measures have been taken.
[0014]
The equations are reconfigured so that the coefficient matrix is full rank, that is, the number of ranks and the number of equations are equal.
[0015]
The value of the design variable vector is not exactly zero, but is replaced with a minute value close to zero.
[0016]
However, in the former method, there is a possibility that the equation is updated every time the design variable vector is updated, which requires a calculation time.
[0017]
Further, the latter method, that is, setting the element design variable to a minute value corresponds to the fact that a physically thin film or a weak member exists. In other words, the conventional method does not accurately represent a portion where there is no substance, and thus the accuracy of the obtained calculation result remains a question.
[0018]
In the present invention, in view of the above-described conventional problems, even when the overall stiffness matrix becomes singular, the calculation can be executed without performing special processing, thereby simplifying the program and further reducing the calculation amount. structural Optimization method及beauty program to reduce, to a storage medium.
[0019]
[Means for Solving the Problems]
In order to solve the above-described problems, the structural optimization design method of the present invention includes first and second storage means for storing variables, third storage means for storing a program, and execution means for executing the program. And the execution means realized by executing a program stored in the third storage means, first and second solution finding means, design variable updating means, first and second In the structural optimization design apparatus having the state variable update means, the preprocessing means, the constant update means, the sensitivity coefficient update means, the residual update means, and the search direction update means, the first solving means includes: A first solution step (a) for solving a first evaluation functional optimization problem with respect to the state variable vector and the design variable vector; and the second solution solving means includes a second solution for the state variable vector and the design variable vector. Of the evaluation functional A structural optimization design method for finding a solution of a structural optimization design problem formulated as a double optimization problem having a second solution step (b) for solving the optimization problem, wherein the structure to be optimized is divided vertically and horizontally Each part as an element, the vertex of each element as a node, the state variable vector as a displacement at each node, the design variable vector as a presence rate of a structural member in each element, and the first solving step (a) (A-1) a constant updating step in which the constant updating unit updates an undetermined constant; and (a-2) a design variable vector and a state variable vector stored in the first storage unit by the design variable updating unit. , Updates the design variable vector based on the undetermined constant and a sensitivity coefficient for the design variable vector of the first evaluation functional, and stores the updated design variable vector in the first storage unit Variable update (A-3) The first state variable updating means reads the design variable vector and the state variable vector stored in the second storage means, updates the state variable vector, and the updated state variable vector. A first state variable update step of storing in the second storage means, and (a-4) sensitivity at which the sensitivity coefficient update means updates the sensitivity coefficient based on the design variable vector and the design variable vector. The coefficient updating step is repeated until a predetermined end condition is satisfied, and the first state variable updating step (a-3) includes the second solution solving step (b), and the second solution solving step. (b) is configured by (b-1) the conjugate gradient method, and (b-2) the preprocessing means sets the nodal force vector component corresponding to the diagonal component that is 0 of the entire stiffness matrix to 0. (B-3-1) the second state variable updating means includes a preprocessing step for performing preprocessing to A second state variable updating step of updating the state variable vector based on the direction vector and the residual vector, (b-3-2) the residual updating means, in said search direction vector and the global stiffness matrix A residual update step of updating the residual vector based on (b-3-3) a search in which the search direction update means updates the search direction vector based on the residual vector and the overall stiffness matrix The direction updating step is repeated until the square of the norm of the residual vector becomes smaller than a set value or exceeds a predetermined number of iterations.
[0020]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, a configuration example and an operation example of a structure optimum design apparatus for realizing the structure optimum design method of the present invention will be described.
[0021]
<Configuration example of the structure optimum design apparatus of this embodiment>
In FIG. 1, the structural example of the structure optimal design apparatus of this embodiment is shown.
[0022]
In the figure, 201 is a CPU for calculation / control, 202 is a display device, 203 is an input device, 204 is a primary storage device, 205 is a secondary storage device, 206 is a communication device, and 207 is a bus line. The components in this embodiment are stored in advance in the secondary storage device 205 as a program, loaded into the primary storage device by a command input from the input device 203 or the communication device 206, and executed by the CPU 201. is there.
[0023]
FIG. 2 shows a storage example of data and programs stored in the primary storage device 204 and the secondary storage device 205 in order to realize the structure optimum design method of the present embodiment.
[0024]
The primary storage device 204 includes a RAM and a ROM, and has a data storage area 214 and a program storage area 224. FIG. 2 does not show an OS or BIOS program that is resident without change, an OS or BIOS parameter that is not changed, and the like. In the data storage area 214, various data for realizing the structure optimum design of the present embodiment is temporarily stored. It is essential that data used again / data used repeatedly is temporarily stored. Temporary storage is not essential for data that does not need to be saved or retained.
[0025]
For example, in the data storage area 214, x (state variable vector) 214a, f (design variable vector) 214b, L ₁ (state variable vector evaluation function) 214c, L ₂ (design variable vector evaluation function) 214d, B1 (Boundary conditions for state variables) 214e, W ₀ (Weight conditions for design variables) 214f, X ⁽⁰⁾ and f ⁽⁰⁾ (initial value) 214g, λ (Lagrange undetermined constant) 218h, A ( (Element stiffness matrix) 214i, k (number of design variable change iterations) 214j, C (sensitivity coefficient) 214k, r (residual vector) 124m, p (search direction vector) 215n, r ⁽⁰⁾ and p ⁽⁰⁾ ( Initial value) 214p, t (number of state variable change iterations) 214q, b (nodal force vector) 214r, other data 214s, and the like are included. On the other hand, a program loaded from the secondary storage device 205 or the communication device 206 is stored in the program storage area 224 and executed by the CPU 201.
[0026]
The secondary storage device 205 is, for example, an external large-capacity memory such as a floppy disk or a CD, and it is more preferable that the secondary storage device 205 is removable. The secondary storage device 205 also has a data storage area 215 and a program storage area 225.
[0027]
Although not shown in detail, the data storage area 215 stores various initial values 215a, various conditions 215b, or a database 215c of double optimal value problems. The database 215c may store and use a standardized matrix or the like.
[0028]
In the program storage area 225, the optimality criterion method module 225a, the sequential convex function approximation method module 225b, the sequential linear programming module 225c, and the other method module 225d, which are used as a solution for the design variable vector in this embodiment. Is stored and used as a solution in the state variable vector, conjugate gradient method module 225e, GCR method module 225f, GCR (k) method 225g, Orthomin (k) method 225h, GMREG (k) method 225i, and other methods Module 225j is stored. A module corresponding to a method selected from the programs stored in the program storage area 225 is loaded into the program storage area 224 of the primary storage device 204 and executed by the CPU 201 to realize the structure optimum design method of the present embodiment. .
[0029]
That is, the present structure optimal design apparatus includes a first solution solving means for solving the optimization problem of the first evaluation functional for the state variable vector and the design variable vector, and a second evaluation for the state variable vector and the design variable vector. A structure optimization design apparatus for solving a structure optimization design problem formulated as a double optimization problem, wherein the state variable vector at each node is a second optimization means for solving a functional optimization problem. The displacement, the design variable vector is an abundance ratio of the structural member in each element, and the second solving means is configured by the conjugate gradient method, and before the pre-processing is performed on the nodal force vector based on the entire stiffness matrix Including a processing means, the state variable vector is not initialized when the processing of the second solution finding means is started. The pre-processing means sets the contact force vector component corresponding to the row or column where the diagonal component of the entire stiffness matrix is zero to zero. The first solving means executes one of a sequential linear programming method, an optimality criterion method, and a sequential convex function approximation method.
[0030]
<Formulation of structure optimization problem of this embodiment>
For the following explanation, the structural optimization problem is formulated.
[0031]
Assuming that the variable function is expressed by a finite-dimensional vector by the formulation by the finite element method, the evaluation functional of the variable function is an evaluation function of the variable vector. Hereinafter, it is described as being expressed by a finite dimensional vector.
[0032]
The state variable vector x and the design variable vector f are written as column vectors as follows.
[0033]
x = (x ₁ , x ₂ ,..., x _m ) ^T (Expression 1)
f = (f ₁ , f ₂ ,..., f _n ) ^T (Expression 2)
Here, T represents transposition. x is an m-dimensional vector, and f is an n-dimensional vector. In the structure optimization problem, the state variable vector is composed of the displacement of each node, and the design variable vector is composed of the existence rate of the structure in each element.
[0034]
For example, in the structure optimization problem in a two-dimensional plane, splitting the design area vertical n _y, the horizontal n _x,
The number of elements n is (n _x × n _y ) and the number of nodes is (n _x +1) × (n _y +1).
[0035]
Since the state variable vector is a set of vertical and horizontal displacements of each node, the dimension m of x is 2 × ( _nx + 1) × ( _ny + 1).
[0036]
The boundary condition of x is usually given as a displacement constraint condition. Describe it B ₁ and.
B ₁ (x) = 0 (Formula 3)
When the evaluation functions regarding the state variable vector x and the design variable vector f are L ₁ and L ₂ , respectively, the following equations are defined.
L ₁ : = L ₁ (x, f) = (b−Ax) ^T (b−Ax) (Formula 4)
L ₂ : = L ₂ (f, x) = (1/2) x ^T Ax (Formula 5)
Here, A is the overall stiffness matrix and b is the nodal force vector. Here, b is given in advance, and A is a function of the design variable vector f.
[0037]
The overall stiffness matrix can be configured by superimposing the element stiffness matrix A _{j on} the element j. In the case of a plane strain problem, an element stiffness matrix in consideration of a design variable vector, that is, a structural element existence rate f _j is shown below. However, a four-node element is adopted as a two-dimensional finite element, and the element displacement vector x _j is configured as follows.
x _j = (u ₁ , v ₁ , u ₂ , v ₂ ,..., u ₄ , v ₄ ) (Expression 6)
Here, u _k and v _k are displacements of the node k in the horizontal direction and the vertical direction.
[0038]
The element stiffness matrix corresponding to the element displacement vector of (Expression 6) is given by the following expression.
[0039]
[Expression 1]

[0040]
A _{i, j} = (λ + 2μ) <∂ _x e _i , ∂ _x e _j > + μ <∂ _y e _i , ∂ _y e _j > (Expression 8)
_{b i, j = λ <∂} x e i, ∂ y e j> + μ <∂ x e j, ∂ y e i> ... ( Eq. 9)
c _{i, j} = λ <∂ _x e _j , ∂ _x e _i > + μ <∂ _x e _i , ∂ _y e _j > (Equation 10)
_{d i, j = (λ +} 2μ) <∂ y e i, ∂ y e j> + μ <∂ x e i, ∂ x e j> ... ( Equation 11)
Here, e _j is a basis vector for node j, <,> represents an inner product of vectors. ∂ _x partial differential related to x, ∂ _y are partial derivatives with respect to y. ξ is set to 2 in accordance with the above (Non-Patent Document 4). Further, λ and μ are material constants called lame constants, which are calculated from the Young's modulus E and Poisson's ratio ν according to the following equation.
[0041]
λ = Eν / (1 + ν) (1-2ν) (Formula 12)
μ = E / {2 (1 + ν)} (Formula 13)
As easily understood from (Expression 7), when the existence ratio of a certain element becomes zero, all the components of the element stiffness matrix A _j become zero. The overall stiffness matrix is composed of a superposition of element stiffness matrices. Therefore, when the design variable f _{j of the} element including the node j becomes 0, all the components in the row and column corresponding to the node j of the overall stiffness matrix become 0. Thus, in the structure optimization problem, the overall stiffness matrix is generally a unique matrix.
[0042]
The constraint condition regarding the design variable is the upper limit value W ₀ of the total weight in the case of the structure optimum design.
[0043]
[Expression 2]

[0044]
Further, the constraint conditions for the range of values are determined as follows in the design variables.
[0045]
0 ≦ f _j ≦ 1 j = 1,..., N (Equation 15)
From the above notation, the optimal design problem is formulated as a problem of minimizing (Expression 5) under the constraint conditions (Expression 14) and (Expression 15).
[0046]
The state variable x is obtained as a solution that minimizes (Equation 4) under the constraint condition (Equation 3).
[0047]
<Operation example of the structure optimum design apparatus of this embodiment>
Based on the above formulation, a processing procedure example of the present embodiment will be described.
[0048]
FIG. 3 shows a flowchart of this embodiment.
[0049]
In the figure, step S101 is a process of reading the specifications of the system to be simulated and initializing variables. The specification can be read by input data from the input device 203 or the communication device 206, or data stored as a file in the secondary storage device 205 in advance can be read and used. The system specifications include x ⁽⁰⁾ , f ⁽⁰⁾ , boundary condition B, and evaluation functions L ₁ , L ₂ as initial values of x, f. Based on this information, the program secures a necessary variable area in the primary storage device 204 and sets a value. The initial value C _j ⁽⁰⁾ of the sensitivity coefficient is calculated by the following equation.
[0050]
C _j ⁽⁰⁾ = − (x _j ⁽⁰⁾ ) ^T (ξ (f _j ⁽⁰⁾ ) ^{ξ −1} A _j ) x _j ⁽⁰⁾ (Equation 16)
In step S102, the number of iterations k is initialized to 1. The number of iterations k is for the iteration process for the above-described design variable solution process, and ends when a predetermined termination condition (conditions of the evaluation functions L ₁ and L ₂ ) is satisfied.
[0051]
In step S103, the Lagrange undetermined constant λ is updated. In step S104, the design variable vector f is updated. In step S105, the state variable vector x is updated. In step S106, a sensitivity coefficient is calculated.
[0052]
Steps S103 and S104 differ depending on the method employed as a solution for the constraint-conditional optimization problem. For example, the solving process for calculating the design variable vector f ^(k) is realized by any one of a sequential linear programming method, an optimality criterion method, and a sequential convex function approximation method. Hereinafter, a case where the optimality criterion method is used as a solution method and a case where the successive convex function approximation method is used will be described in detail.
[0053]
(Optimality criteria method)
First, the case where the optimality criterion method is used will be described. The optimality criterion method for the structure optimization problem is known from (Non-Patent Document 4) and (Non-Patent Document 5).
[0054]
In the optimality criterion method, (Expression 4) and (Expression 5) are combined into one expression using a Lagrange undetermined constant λ.
[0055]
L (f, λ) = L ₂ (f, x) −λ (W (f) −W ₀ ) (Expression 17)
The update of λ executed in step S103 is performed using the following equation.
[0056]
λ ^{(k + 1)} = min [0, [(1 / W ₀ ) W (f ^(k) )] ^α λ ^(k) ] (Equation 18)
However, the superscript k, k + 1 represents the number of iterations. Α is a power coefficient, and is usually set to a value of about 0.85.
[0057]
The design variable update executed in step S104 is performed according to the following equation.
[0058]

However, ζ is a constraint value of the fluctuation range of the design variable update, and takes a value of about 0.3. S _j ^(k) is given by the following equation.
[0059]
s _j ^(k) = [(1 / λ ^(k-1) ) C _j ^(k-1) ] αf _j ^(k-1) (Equation 20)
Here, C _j ^(k) is an amount called a sensitivity coefficient related to the design variable f _j ^(k) of the second evaluation function L ₂ and is calculated in the process of step S106.
[0060]
Furthermore, the design variable f _j ^(k) calculated by (Equation 19) may be forcibly set to 0 when it becomes smaller than a preset value. This preset value is, for example, a value of about 10 ⁻³ .
[0061]
In step S105, the state variable vector x ^(k) is updated.
[0062]
The state variable vector x ^(k) is obtained in the process of minimizing the evaluation function of (Equation 4). This solution is called a conjugate residual method and will be described later.
[0063]
In step S106, the sensitivity coefficient C _j ^(k) is calculated by the following equation.
[0064]
C _j ^(k) = − (x _j ^(k) ) ^T (ξ (f _j ^(k) ) ^{ξ −1} A _j ) x _j ^(k) (Equation 21)
(Sequential convex function approximation method)
Next, the processing in steps S103 and S104 when the successive convex function approximation method is used will be described in detail. In the successive convex function approximation method, numerical calculation errors can be avoided by performing scaling on each component of the design variable vector. However, in order to simplify the explanation, scaling is not described below. Details of the successive convex function approximation method are known from (Reference Document 1) and (Reference Document 2).
[0065]
(Reference 1) Fujii: “Structural design solved with a personal computer”, Maruzen (April 2002)
(Reference 2) C. Fleury: "CONLIN, an efficient dual optimizer based on convex approximation concepts," Structural Optimization, 1, pp.81-89 (1989)
In step S103, the Lagrange undetermined constant is updated by the following equation.
[0066]
λ ^(k) = λ ^(k-1) + Δλ ^(k) (Equation 22)
However, the initial value λ ⁽⁰⁾ is set to a value of about 0.01, for example.
[0067]
Δλ ^(k) is given as a solution to the following equation:
[0068]
Δλ ^(k) = 2 × (1 −W ₀ / W (f ^(k−1) )) λ ^(k−1) (Equation 23)
In step S104, the design variable vector is updated.
[0069]

Furthermore, the design variable f _j ^(k) calculated by (Equation 24) may be forcibly set to 0 when it becomes smaller than a preset value. This preset value is about 10 ⁻³ .
[0070]
The processing in steps S105 and S106 is as described above.
[0071]
(Conjugate gradient method)
A process using the conjugate gradient method will be described with reference to FIG. 4 as a solution step for calculating the state variable vector x ^(k) in step S105.
[0072]
In step S301, the specifications of the system to be simulated are read. Initial value x of x ^(k) is the specifications of the system ^(0), the initial value f of the design variable vector f ^{(k) (0),} the displacement constraints B _1, which includes an evaluation function L _1. Based on this information, the program secures a necessary variable area in the primary storage device 204 and sets a value. Also, the residual vector r and the search direction vector p are initialized. Respective initial values are represented as r ⁽⁰⁾ and p ⁽⁰⁾ .
[0073]
First, r ⁽⁰⁾ is calculated by the following equation using x ⁽⁰⁾ and a nodal force vector b describing a force applied to each node (node) in a vector format.
[0074]
r ⁽⁰⁾ = b-Ax ⁽⁰⁾ (Formula 26)
The initial value p ⁽⁰⁾ of the search direction vector p is set equal to r ⁽⁰⁾ .
[0075]
In step S302, preprocessing is performed. First, among the diagonal components of the overall stiffness matrix A, the row or column number that becomes 0 is examined. If there is no component that becomes 0, the process proceeds to step S303 without doing anything. If there is a component that becomes 0, the component with the matching number in b (nodal force vector) is set to 0.
[0076]
In step S303, the number of iterations t is initialized to 1. The number of iterations t is for the iteration processing for the above-described state variable solution process, and the number of iterations of the following steps S304 to S308 exceeds the set value or the square of the norm of the residual vector r ^(t) When (r ^(t) , r ^(t) ) becomes smaller than a predetermined value set in advance, the iteration is terminated. Hereinafter, the value of the tth iteration is written as x ^(t) .
[0077]
In step S304, the state vector update coefficient α ^(t) is calculated.
[0078]
α = (r ^(t-1) , p ^(t-1) ) / (p ^(t-1) , Ap ^(t-1) ) (Equation 27)
In step S305, the state variable vector x is updated.
[0079]
x ^(t) = x ^(t-1) + αp ^(t-1) (Formula 28)
In step S306, the residual vector r ^(t) is calculated by the following equation.
[0080]
r ^(t) = b-Ax ^(t-1) (Formula 29)
In step S307, an update coefficient β of the search direction vector p is calculated.
[0081]
β = (r ^(t) , r ^(t) ) / (r ^(t-1) , r ^(t-1) ) (Equation 30)
However, when t = 1, β = 1.
[0082]
In step S308, the search direction vector p is updated.
[0083]
p ^(t) = r ^(t) + βp ^(t-1) (Formula 31)
The conjugate gradient method and its various derivation methods are known from (Reference 3) and the like.
[0084]
(Reference 3) Mori, Sugihara, Murota: “Linear Computation,” Iwanami Lecture, Applied Mathematics, Iwanami Shoten (1994)
(Specific example of structural optimum design of this embodiment)
In this example, the present embodiment is applied to the optimum shape automatic design of a cantilever beam that receives a load at an arbitrary position. For simplicity of explanation, we will limit it to the plane distortion problem.
[0085]
As shown in FIG. 5, the design area that allows the presence of the structural member is a rectangle 402, and the area is divided into vertical n _y and horizontal _nx at equal intervals according to the finite element method. The divided partial region is called a cell, performing numbered as lower left and upper right of the cell is respectively (1,1) and (n _y, n _x).
[0086]
Similarly called grid points and node performs numbered as lower left and upper right of the node becomes (1,1) and _{(n y + 1, n x} +1).
[0087]
In the figure, 401 is a support member, and 403 is a load vector.
[0088]
The characteristic function value f (j, k) corresponds to the cell (j, k). Here, the characteristic function value is a variable taking a positive real value from 0 to 1 indicating the existence probability of the structural member in the cell (j, k), and is an element of the design variable vector f in the present embodiment.
f = (f (1,1), f (1,2),..., f (n _y , n _x )) ^T (Expression 32)
Similarly, node (j, k) corresponds to lateral displacement u (j, k) and longitudinal displacement v (j, k). These are real numbers taking arbitrary values, and are elements of the state variable vector x in the present embodiment.
x = (u (1,1), v (1,1), u (1,2), v (1,2), ..., u (n _y + 1, n _x +1), v (n _y + 1, n _x +1)) ^T (Formula 33)
FIG. 6 shows the calculation results that can be put in this example. In the figure, black areas are areas where structural members are present. In this specific example, the aspect ratio of the design area of the structural member is 2: 1 in the vertical / horizontal direction, and the analytical solution is known to be a combination of ± 45 degrees beams in the horizontal direction. It has been. It can be seen that the calculation result shown in FIG. 6 is in good agreement with such an analytical solution.
[0089]
In this embodiment, as shown in FIG. 1, the system is described as one system in which each device is connected to the bus line 207. However, these elements are connected via a network or the like to form a system of a plurality of devices. Alternatively, a plurality of general-purpose computers may be connected and processed in parallel by multiple processing or shared processing.
[0090]
The present invention also includes a program that realizes the above-described structure optimization design method and a storage medium that stores the program in a computer-readable manner.
[0091]
【The invention's effect】
As described above, according to the present invention, even when the overall stiffness matrix becomes singular, the calculation can be executed without performing special processing, thereby simplifying the program and further reducing the amount of calculation. Can be reduced.
[Brief description of the drawings]
FIG. 1 is a diagram illustrating a configuration example of a structure optimum design apparatus according to an embodiment.
FIG. 2 is a diagram showing a storage configuration example of the primary storage device and the secondary storage device of FIG. 1;
FIG. 3 is a flowchart illustrating an example of a processing procedure according to the present embodiment.
FIG. 4 is a flowchart of a processing procedure showing step S105 of FIG. 3 in more detail.
FIG. 5 is an explanatory diagram of problem setting of a specific example applied to the present embodiment;
FIG. 6 is a diagram showing a calculation result according to the present embodiment related to FIG. 5;

Claims (5)

1st and 2nd memory | storage means for memorize | storing a variable, The 3rd memory | storage means which memorize | stored the program, The execution means which executes a program, The said execution means memorize | stores in the said 3rd memory | storage means First and second solution finding means, design variable updating means, first and second state variable updating means, preprocessing means, constant updating means, which are realized by executing the program, In the structure optimal design apparatus having a sensitivity coefficient update unit, a residual update unit, and a search direction update unit,
A first solving step (a) in which the first solving means solves an optimization problem of the first evaluation functional with respect to the state variable vector and the design variable vector; and A structure optimum design method for finding a solution of a structure optimum design problem formulated as a double optimization problem having a second solution step (b) for solving an optimization problem of a second evaluation functional for the design variable vector Because
Dividing the structure to be optimized vertically and horizontally, each part as an element, the vertex of each element as a node, the state variable vector as a displacement at each node, the design variable vector as the abundance of structural members in each element,
The first solving step (a) includes:
(a-1) the constant update means, a constant update step of updating an undetermined constant;
(a-2) The design variable update means reads the design variable vector and the state variable vector stored in the first storage means, and the sensitivity coefficient for the design variable vector of the undetermined constant and the first evaluation functional A design variable update step of updating the design variable vector on the basis of and storing the updated design variable vector in the first storage means;
(a-3) The first state variable update unit reads the design variable vector and the state variable vector stored in the second storage unit, updates the state variable vector, and updates the updated state variable vector. A first state variable update step to be stored in the storage means;
(a-4) The sensitivity coefficient update means repeats the sensitivity coefficient update step of updating the sensitivity coefficient based on the design variable vector and the design variable vector until a predetermined end condition is satisfied,
The first state variable updating step (a-3) includes the second solution finding step (b);
The second solution step (b) includes:
(b-1) Constructed by the conjugate gradient method,
(b-2) the pre-processing means includes a pre-processing step of performing pre-processing for setting a nodal force vector component corresponding to a diagonal component that is 0 in the overall stiffness matrix to 0;
(b-3-1) and the second state variable updating means, and a second state variable updating step of updating the state variable vector based on the search direction vector and the residual vector,
(b-3-2) the residual update unit updates the residual vector based on the search direction vector and the overall stiffness matrix;
(b-3-3) a search direction update step in which the search direction update means updates the search direction vector based on the residual vector and the overall stiffness matrix, and a square of the norm of the residual vector The structural optimum design method is characterized by repeating until the value becomes smaller than a set value or exceeds a predetermined number of iterations. The earlier in the process, the structure optimum design method according to claim 1, characterized in that the component sections point force vector corresponding to a row or column of the diagonal components of the global stiffness matrix is 0 to 0.   The structural optimal design method according to claim 1, wherein the first solving step is realized by any one of a sequential linear programming method, an optimality criterion method, and a sequential convex function approximation method. Program for executing the steps of the structure optimum design method according to the computer in any one of claims 1 to 3. Computer read-readable storage medium storing a program of claim 4.

Images