JP2020016936A - Arithmetic unit for structural optimization with unnecessary supporting structure - Google Patents
Arithmetic unit for structural optimization with unnecessary supporting structure Download PDFInfo
- Publication number
- JP2020016936A JP2020016936A JP2018137736A JP2018137736A JP2020016936A JP 2020016936 A JP2020016936 A JP 2020016936A JP 2018137736 A JP2018137736 A JP 2018137736A JP 2018137736 A JP2018137736 A JP 2018137736A JP 2020016936 A JP2020016936 A JP 2020016936A
- Authority
- JP
- Japan
- Prior art keywords
- matrix
- eigenvalue
- load
- design
- calculating
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Abstract
Description
本発明は、コンピュータを用いた構造体の構造最適化演算を実行する装置に係り、より詳細には、コンピュータによる数値解析を通じて、任意の設計条件・境界条件下に於いて最適な構造（所与の力学的条件下にて剛性が最大化される構造）を見出す装置に係る。 The present invention relates to an apparatus for performing a structure optimization operation of a structure using a computer, and more particularly, to an optimum structure (given condition) under arbitrary design conditions and boundary conditions through numerical analysis by a computer. (A structure in which the rigidity is maximized under the mechanical conditions of the present invention).
コンピュータを用いた有限要素法などの数値計算技術の発展により、任意の構造体の、任意の設計条件・境界条件下に於ける最適な構造又は形状をコンピュータ上で求める「構造最適化」演算が実施されるようになっている。かかるコンピュータを用いた構造最適化の手法によれば、車両のボディーを初め、種々の機械器具、建築物などの種々の構造体を実際に製作する前に、構造体が要求される特性を備えているか否かの把握、或いは、要求される若しくは好ましい性能を備えた構造体の設計案の取得が可能となる点で有利である。 With the development of numerical calculation technology such as the finite element method using a computer, the "structure optimization" operation for finding the optimal structure or shape of an arbitrary structure under arbitrary design conditions and boundary conditions on a computer It is being implemented. According to the structure optimization method using such a computer, before actually manufacturing various structures such as a vehicle body, various machinery and equipment, and buildings, the structures have the required characteristics. This is advantageous in that it is possible to determine whether or not the design is performed or to obtain a design proposal for a structure having required or preferable performance.
そのような「構造最適化」の一つとしては、「トポロジー最適化（位相最適化）」と称される手法が知られている（特許文献１、２、非特許文献１、２）。かかる「トポロジー最適化」では、端的に述べれば、有限要素法等による数値解析とその結果から得られる平均コンプライアンス等の目的関数又は評価関数を最適化する演算処理とを通じて、或る与えられた領域（設計領域）内に於いて、任意の設計条件・境界条件の下での最適な構造、例えば、任意の力学的条件下での剛性が最大化される構造、を与える材料分布が生成される。より具体的には、例えば、特許文献１では、初めに、多数の、それぞれ微小な孔を有する有限要素から構成される任意の寸法の設計領域を設定した上で、任意の境界条件と体積制約条件の下、各有限要素の孔の寸法を変化させながら設計領域内にて有限要素により形成される構造体の平均コンプライアンスとその（各要素の孔の大きさに対する）感度を反復して演算し、最小の平均コンプライアンスを与える有限要素内の孔の寸法の分布を見出し、剛性の最大化された構造体の形状を生成する手法が開示されている。特許文献２は、材料の密度がそれぞれ割り当てられた多数の有限要素から構成される任意の寸法の設計領域に於いて、任意の境界条件と体積制約条件の下、各有限要素の材料密度を変化させながら有限要素により形成される構造体の平均コンプライアンスとその（各要素の密度に対する）感度を反復して演算し、平均コンプライアンスを最小化する材料の密度分布を生成し、剛性の最大化された構造体の形状を生成する手法を開示している。特許文献３は、任意の構造体に於いて外力を受ける部分の構造最適化を行うべく、構造最適化の対象となる部分と構造体のその他の部分との間で荷重が伝達されるように構造体の数理モデルを設定することを提案している。 As one of such “structural optimizations”, a technique called “topology optimization (phase optimization)” is known (Patent Documents 1 and 2, Non-Patent Documents 1 and 2). In short, in the “topology optimization”, in short, a given area is obtained through a numerical analysis by a finite element method or the like and an arithmetic processing for optimizing an objective function or an evaluation function such as an average compliance obtained from the result. Within the (design area), a material distribution is generated that gives an optimal structure under any design conditions and boundary conditions, for example, a structure that maximizes rigidity under any mechanical conditions. . More specifically, for example, in Patent Document 1, first, a design area of an arbitrary size composed of a large number of finite elements each having a minute hole is set, and then an arbitrary boundary condition and a volume constraint are set. Under the conditions, the average compliance of the structure formed by the finite elements and its sensitivity (with respect to the hole size of each element) are repeatedly calculated in the design area while changing the dimensions of the holes of each finite element. Disclosed are techniques for finding the distribution of pore sizes within a finite element that gives the least average compliance and generating a stiffened structure shape. Patent Document 2 discloses that the material density of each finite element is changed under an arbitrary boundary condition and a volume constraint condition in a design area of an arbitrary size composed of a large number of finite elements each assigned a material density. The mean compliance of the structure formed by the finite elements and its sensitivity (for the density of each element) are iteratively calculated, producing a density distribution of the material that minimizes the average compliance and maximizing the stiffness. A method for generating a shape of a structure is disclosed. Patent Document 3 discloses that, in order to perform structural optimization of a portion of an arbitrary structure receiving an external force, a load is transmitted between a portion to be subjected to the structure optimization and another portion of the structure. It is proposed to set up a mathematical model of the structure.
また、任意の構造体のコンプライアンスを最小化するように構造体の形状を変化させる形状最適化演算も種々知られている。例えば、特許文献４では、形状最適化演算に於ける計算の処理負担を低減させるべく、フェーズフィールド法を適用して任意の構造体の形状の最適化演算を実行する構成が提案されている。 Also, various shape optimization operations for changing the shape of a structure so as to minimize the compliance of an arbitrary structure are known. For example, Patent Literature 4 proposes a configuration in which a phase field method is applied to perform an optimization operation on a shape of an arbitrary structure in order to reduce a processing load of calculation in the shape optimization operation.
更に、構造体に対して任意の荷重が作用した場合、即ち、構造体に対して作用する荷重が不定である場合に剛性が最大化される最適構造を求める「ロバスト構造最適化」も提案されている（非特許文献３、４）。かかるロバスト構造最適化に於いては、構造体の或る領域に作用する荷重の値と分布が不定な場合の最適な構造として、想定される最も危険な荷重（最悪荷重）に対して剛性を最大化する構造体の構造が演算により求められる。より具体的には、最悪荷重が付与されたときのコンプライアンスと荷重の分布は、構造体に於いて荷重の作用する領域（荷重作用領域）内の荷重分布を表すベクトルＦｌと荷重作用領域の変位分布を表す変位ベクトルＵｌとを関係づける方程式：
Ｕｌ＝Ｃ・Ｆｌ …（１）
に於ける局所的なコンプライアンス行列Ｃの（荷重ベクトルＦｌのノルム||Fl||＝１のときの）最大固有値とこれに対応する固有ベクトルにより与えられることが理論的に見出されているので、ロバスト構造最適化では、実際には、式（１）の局所コンプライアンス行列Ｃの最大固有値を最小化するように、構造体に於ける材料分布や形状の最適化演算が実行される。このロバスト構造最適化演算によれば、任意の構造体に於いて、任意の領域に最も変形を大きくする最悪荷重が作用した場合に最大の剛性を呈し変形が最も小さくなる構造が算定できることとなるので、如何なる分布にて荷重が作用したしても高い剛性を呈する最適な構造が提供できることが期待される。
Furthermore, "robust structure optimization" for finding an optimal structure that maximizes rigidity when an arbitrary load acts on the structure, that is, when the load acting on the structure is undefined, has also been proposed. (Non-Patent Documents 3 and 4). In such a robust structure optimization, as an optimal structure in a case where the value and distribution of a load acting on a certain region of the structure are indefinite, rigidity is assumed with respect to an assumed most dangerous load (worst load). The structure of the structure to be maximized is calculated. More specifically, the compliance and the distribution of the load when the worst load is applied are determined by a vector Fl representing the load distribution in a region where the load is applied (load application region) in the structure and a vector Fl of the load application region. Equation relating the displacement vector U l representing the displacement distribution:
U 1 = C · F 1 (1)
Has been theoretically found to be given by the maximum eigenvalue (when the norm of the load vector F 1 || F l || = 1) and the corresponding eigenvector of the local compliance matrix C at Therefore, in the robust structure optimization, in practice, the material distribution and the shape optimization operation in the structure are executed so as to minimize the maximum eigenvalue of the local compliance matrix C in Expression (1). According to this robust structure optimization calculation, it is possible to calculate a structure that exhibits maximum rigidity and minimizes deformation when a worst load that maximizes deformation is applied to an arbitrary region in an arbitrary structure. Therefore, it is expected that an optimum structure exhibiting high rigidity can be provided even when a load is applied in any distribution.
ところで、現実の任意の構造体に荷重が作用する状態を考えたとき、構造体を支持する構造或いは構造体を固定する支点の位置などが不定であることが多い。例えば、ユーザーの手からの荷重を受ける携帯電話の外殻など、実際の構造体は、どのように保持或いは支持されるかは不定であり、そのような場合、構造体に於ける支点の位置は不定となる。 By the way, in consideration of a state in which a load is applied to an actual arbitrary structure, a structure for supporting the structure or a position of a fulcrum for fixing the structure is often indefinite. For example, the actual structure, such as the outer shell of a mobile phone, which is under the load of the user's hand, is undefined as to how it will be held or supported, in which case the position of the fulcrum in the structure Becomes indefinite.
しかしながら、上記の如き従前のロバスト構造最適化又はその他の構造最適化の方法に於いては、構造体に於いて支持構造或いは支点を明確に規定する拘束条件を与えなければ、演算を実行できないようになっている。上記の如き構造最適化演算では、端的に述べれば、最適構造を見出すために必要な最小化されるべきコンプライアンス若しくはその他の評価関数又はその感度（材料分布又は形状の変化に対する評価関数の変化、設計変数についての評価関数の導関数）の決定のために、構造体に於ける変位ベクトルＵと荷重ベクトルＦとを関連付ける剛性方程式：
ＫＵ＝Ｆ …（２）
に於ける剛性行列Ｋの逆行列Ｋ-1を用いて、上記の局所コンプライアンス行列Ｃの算出、並びに、変位ベクトルＵの算出が必要となる。ところが、構造体を支持する構造又は支点を与える拘束条件がない場合には、剛性行列Ｋの逆行列Ｋ-1が存在せず、局所コンプライアンス行列Ｃと変位ベクトルＵの演算が実行できないこととなっている。これは、物理的には、構造体を支持する構造がない状態では、構造体が剛体として自由に平行又は回転運動し得ることとなり（以下、「剛体運動」と称する。）、構造体内の各点の変位が式（２）の解として一意に決定されないことに相当する。
However, in the above-described conventional robust structure optimization or other structure optimization methods, the operation cannot be executed unless a constraint condition that clearly defines a supporting structure or a fulcrum in the structure is given. It has become. In the structure optimization operation as described above, in short, the compliance or other evaluation function to be minimized in order to find the optimum structure or its sensitivity (change of the evaluation function with respect to a change in material distribution or shape, design, etc.) Stiffness equation relating the displacement vector U and the load vector F in the structure to determine the derivative of the evaluation function for the variables):
KU = F (2)
It is necessary to calculate the local compliance matrix C and the displacement vector U using the inverse matrix K -1 of the stiffness matrix K in the above. However, when there is no constraint that gives a structure or a fulcrum that supports the structure, the inverse matrix K −1 of the rigidity matrix K does not exist, and the calculation of the local compliance matrix C and the displacement vector U cannot be performed. ing. This means that physically, in a state where there is no structure supporting the structure, the structure can freely move parallel or rotationally as a rigid body (hereinafter, referred to as "rigid body motion"), and each of the structures. This corresponds to that the displacement of the point is not uniquely determined as the solution of the equation (2).
この点に関し、本発明の発明者が研究したところ、支持構造又は支点のない構造体モデルの場合でも、局所コンプライアンス行列Ｃの算出並びに変位ベクトルＵの算出に於いて、剛性行列Ｋの逆行列Ｋ-1に代えて、剛性行列Ｋの「擬似逆行列（ムーア・ペンローズの逆行列または一般逆行列とも称される。）」Ｋ＋を用いると、構造体に於ける剛体運動を除いた変位ベクトル、かかる変位ベクトルにより生ずるコンプライアンス及びその感度を算出することができ、最小化されたコンプライアンスを与える最適な性能を有する構造体の構造又は形状が得られること、即ち、構造の最適化が可能であることが見出された。かかる知見は、本発明に於いて利用される。 In this regard, the inventor of the present invention has studied and found that even in the case of a support structure or a structure model having no fulcrum, the inverse matrix K of the stiffness matrix K is calculated in the calculation of the local compliance matrix C and the calculation of the displacement vector U. By using a “pseudo-inverse matrix (also called a Moore-Penrose inverse matrix or a general inverse matrix)” K + of the rigidity matrix K instead of −1 , a displacement vector excluding the rigid body motion in the structure is obtained. It is possible to calculate the compliance caused by such a displacement vector and its sensitivity, and to obtain the structure or shape of the structure having the optimal performance to give the minimized compliance, that is, to optimize the structure. Was found. Such knowledge is used in the present invention.
かくして、本発明の課題は、コンピュータを用いた数値解析技術により、任意の構造体に任意の荷重が付与される場合に構造体の剛性が最大化される最適な構造を見出す構造最適化演算を実行する場合に於いて、支持構造又は支点のない構造体の構造最適化が可能な手法を提供することである。 Thus, an object of the present invention is to use a computer-based numerical analysis technique to perform a structure optimization operation for finding an optimal structure that maximizes the rigidity of a structure when an arbitrary load is applied to the structure. It is an object of the present invention to provide a method capable of optimizing the structure of a structure having no supporting structure or fulcrum when executed.
本発明によれば、上記の課題は、コンピュータを用いた数値演算処理により設計領域内にて設定された条件の下で設計変数にて構造が決定される構造体にして支持構造のない構造体の最適な構造を見出す構造最適化演算装置であって、
前記設計領域の寸法、前記構造体に於いて荷重の作用する荷重作用領域の範囲、前記構造体を構成する材料の弾性特性を表す物性定数及び前記設計領域の体積に対する前記構造体の体積の割合である体積制約割合を含む条件を設定する条件設定手段と、
前記設計変数を用いて前記構造体の剛性行列を算出する剛性行列算出手段と、
前記剛性行列の擬似逆行列を算出する擬似逆行列算出手段と、
前記擬似逆行列を用いて前記構造体の前記荷重作用領域に於ける局所コンプライアンス行列を算出する局所コンプライアンス行列算出手段と、
前記局所コンプライアンス行列の少なくとも一つの固有値と対応する固有ベクトルとを算出する固有値解析手段と、
前記固有ベクトルと前記擬似逆行列とを用いて前記固有ベクトルを荷重ベクトルとしたときの前記構造体に於ける変位ベクトルを算出する変位ベクトル算出手段と、
前記変位ベクトルを用いて前記設計変数に対する前記少なくとも一つの固有値の感度を算出する固有値感度算出手段と、
前記体積制約割合を満たした状態で前記感度を用いて前記少なくとも一つの固有値が最小化されるよう前記設計変数を更新する設計変数更新手段と、
所定の終了条件が成立するまで、前記更新された設計変数を用いて、前記剛性行列算出手段、前記擬似逆行列算出手段、前記局所コンプライアンス行列算出手段、前記固有値解析手段、前記変位ベクトル算出手段、前記固有値感度算出手段及び前記設計変数更新手段の各処理を反復させる反復処理手段と、
前記設計変数により決定される構造体の最新の構造を出力又は表示する構造体出力手段と
を含む装置によって達成される。
According to the present invention, the above-described object is to provide a structure whose structure is determined by design variables under conditions set in a design area by a numerical operation process using a computer, and a structure having no support structure. A structure optimization arithmetic unit for finding an optimal structure of
The dimensions of the design area, the range of the load application area in which a load acts on the structure, physical property constants representing the elastic properties of the materials constituting the structure, and the ratio of the volume of the structure to the volume of the design area Condition setting means for setting a condition including a volume constraint ratio that is
Stiffness matrix calculation means for calculating a stiffness matrix of the structure using the design variables,
Pseudo-inverse matrix calculating means for calculating a pseudo-inverse of the rigidity matrix,
Local compliance matrix calculation means for calculating a local compliance matrix in the load application region of the structure using the pseudo inverse matrix,
Eigenvalue analysis means for calculating at least one eigenvalue and a corresponding eigenvector of the local compliance matrix,
Displacement vector calculating means for calculating a displacement vector in the structure when the eigenvector is a load vector using the eigenvector and the pseudo inverse matrix,
Eigenvalue sensitivity calculation means for calculating the sensitivity of the at least one eigenvalue to the design variable using the displacement vector,
Design variable updating means for updating the design variable so that the at least one eigenvalue is minimized using the sensitivity while the volume constraint ratio is satisfied,
Until a predetermined termination condition is satisfied, using the updated design variables, the stiffness matrix calculating means, the pseudo-inverse matrix calculating means, the local compliance matrix calculating means, the eigenvalue analyzing means, the displacement vector calculating means, Iterative processing means for repeating each processing of the eigenvalue sensitivity calculating means and the design variable updating means,
Structure output means for outputting or displaying the latest structure of the structure determined by the design variable.
上記の構成に於いて、「構造体」とは、車両のボディーを初め、種々の機械器具、建築物又はそれらを構成する任意の部材などとして使用される構造体であってよい。「構造体の最適な構造を見出す」とは、通常の構造最適化演算、例えば、トポロジー最適化演算、形状最適化演算などの場合と同様に、構造体に荷重が付与される場合の構造体の変形の程度が最小となる構造、より具体的には、構造体のコンプライアンスが最小となる構造を見出すことを意味している（なお、本明細書に於いて、「構造体の構造」というときには、構造体の形状、トポロジー或いは材料密度の分布をいうものとする。）。従って、上記の「設計領域」は、通常の構造最適化演算と同様に設定されてよく、構造体の構造を決定する「設計変数」は、構造最適化演算の態様によって適宜選択されてよい。構造最適化演算の手法としては、具体的には、例えば、所謂「密度法」を用いたトポロジー最適化演算（例えば、特許文献２、非特許文献１、２）が採用されてよく、その場合には、「設計変数」は、「設計領域」を任意に格子状に分割して得られる複数の節点を結んで構成される多面体（又は多角形）の要素に於いて割り当てられる０〜１の値を取る「材料密度」であり、演算過程に於いて更新されて、好ましくは、０又は１に収束していくことが期待される変数であってよい。そして、「設計変数」が「材料密度」である場合には、材料密度が実質的に非０である要素（例えば、小数点第一位を四捨五入して１になる要素）の集合が成す形状が、「設計領域」内にて生成される目的の「構造体」の形状となる。また、設計変数は、構造体の輪郭を画定する線又は面を決定する変数であってもよい。 In the above-described configuration, the “structure” may be a structure used as a vehicle body, various types of machinery, buildings, or any members constituting the same. "Finding the optimal structure of a structure" refers to a structure in which a load is applied to a structure in the same manner as in a normal structure optimization operation such as a topology optimization operation or a shape optimization operation. Means a structure that minimizes the degree of deformation of the structure, more specifically, a structure that minimizes the compliance of the structure (referred to as “structure of the structure” in this specification). Sometimes it refers to the distribution of the shape, topology or material density of the structure.) Therefore, the above-mentioned "design area" may be set in the same manner as in a normal structure optimization operation, and the "design variable" for determining the structure of the structure may be appropriately selected according to the mode of the structure optimization operation. As a method of the structure optimization calculation, specifically, for example, a topology optimization calculation using a so-called “density method” (for example, Patent Literature 2, Non-Patent Literatures 1 and 2) may be adopted. The “design variable” is a polyhedral (or polygon) element that is assigned by connecting a plurality of nodes obtained by arbitrarily dividing the “design area” into a lattice shape. This is a “material density” that takes a value, and may be a variable that is expected to be updated in the calculation process and preferably converge to 0 or 1. When the “design variable” is “material density”, the shape formed by a set of elements whose material density is substantially non-zero (for example, an element that becomes 1 by rounding off the first decimal place) is formed. , A target “structure” generated in the “design area”. Further, the design variable may be a variable that determines a line or a plane that defines the outline of the structure.
上記の本発明の装置に於ける「条件設定手段」では、設計領域の寸法、構造体に於いて荷重の作用する領域（荷重作用領域）の範囲、構造体を構成する材料の弾性特性を表す物性定数、設計領域の体積に対する構造体の体積の割合である体積制約割合を含む構造最適化演算に必要な初期条件が設定される。しかしながら、特に本発明の装置に於いては、最適化の対象となる構造体が、支持構造のない構造体、即ち、外部からの拘束を受けず、自由な剛体運動が許される状態にある構造体であるので、従前の構造最適化演算に際して通常設定される構造体に於ける変位に対する拘束条件又は境界条件は一切設定されなくてよいことは理解されるべきである。また、荷重条件について、構造体に於ける荷重作用領域の具体的な範囲は指定されるが、荷重作用領域内に於ける荷重の値或いは荷重の分布は、従前のロバスト構造最適化演算の場合と同様に不定のままであってよく、これらの条件は、後述の固有値解析演算の結果によって設定される。構造体を構成する材料の弾性特性を表す物性定数は、典型的には、ヤング率とポアソン比であってよいが、これに限定されない。「体積制約割合」は、最適化演算に於ける拘束条件となる（構造体の運動を拘束するものではない。）。 The "condition setting means" in the above-described apparatus of the present invention indicates the dimensions of the design region, the range of the region where the load acts on the structure (load application region), and the elastic characteristics of the material constituting the structure. Initial conditions necessary for a structure optimization operation including a physical property constant and a volume constraint ratio which is a ratio of the volume of the structure to the volume of the design area are set. However, in particular, in the device of the present invention, the structure to be optimized is a structure having no support structure, that is, a structure in which free rigid motion is allowed without being restrained from the outside. It is to be understood that, since it is a body, no constraint or boundary condition for the displacement in the structure normally set in the conventional structure optimization operation needs to be set. Although the specific range of the load application area in the structure is specified for the load condition, the value of the load or the distribution of the load in the load application area is determined by the conventional robust structure optimization calculation. The conditions may be set indefinitely in the same manner as described above, and these conditions are set according to the result of the eigenvalue analysis operation described later. The physical property constant representing the elastic property of the material constituting the structure may typically be Young's modulus and Poisson's ratio, but is not limited thereto. The “volume constraint ratio” is a constraint condition in the optimization calculation (it does not limit the motion of the structure).
上記の装置に於いて、
「剛性行列算出手段」にて算出される構造体の「剛性行列」は、上記の式（２）に於ける剛性行列Ｋと同様のものであってよく、具体的には、構造体に格子状に設定される複数の要素に於ける点（節点）の各々に作用する荷重を成分とする荷重ベクトルと各節点に於ける変位を成分とする変位ベクトルとを関係付ける剛性行列であってよい。各要素の弾性は、材料の弾性特性を表す物性定数と各要素の設計変数（例えば、材料密度）によって決定されるので、剛性行列は、各要素の設計変数の関数となり、本発明による最適化演算の過程で設計変数が更新されると共に剛性行列の各成分の値も変化することとなる。
「擬似逆行列算出手段」にて算出される構造体の剛性行列の「擬似逆行列」は、既に述べた如き、ムーア・ペンローズの逆行列または一般逆行列とも称される行列であり、剛性行列の成分を用いて所定のアルゴリズムにより一義的に算出される。具体的には、擬似逆行列は、プログラム言語に用意された関数コマンドを利用して算出されてよく、或いは、任意の近似方法を用いて算出することも可能である。
「局所コンプライアンス行列算出手段」にて算出される「局所コンプライアンス行列」とは、上記の式（１）に於ける局所コンプライアンス行列Ｃと同様のものであり、具体的には、荷重作用領域内に於ける各節点の荷重を成分とする局所的な荷重ベクトルと各節点の変位を成分とする局所的な変位ベクトルとを関係づける行列である。かかる局所コンプライアンス行列Ｃは、後の実施形態の欄にて詳細に説明される如く、本来は、構造体の全節点の荷重を表す荷重ベクトルと荷重作用領域内の局所的な荷重ベクトルとを関連付ける変換行列と構造体の剛性行列の逆行列とを用いて算出されるところ、本発明の装置に於いては、既に述べた如く、支持構造のない構造体が最適化の対象となっているので、構造体の剛性行列の逆行列が存在していない。そこで、「局所コンプライアンス行列算出手段」では、構造体の剛性行列の逆行列に代えて、構造体の剛性行列の擬似逆行列を用いて、局所コンプライアンス行列が算出される。
In the above device,
The “rigidity matrix” of the structure calculated by the “rigidity matrix calculating means” may be the same as the rigidity matrix K in the above equation (2). It may be a stiffness matrix that associates a load vector with a load acting on each of the points (nodes) of a plurality of elements set in a shape and a displacement vector with a displacement at each node as a component. . Since the elasticity of each element is determined by a physical property constant representing the elastic property of the material and a design variable (eg, material density) of each element, the stiffness matrix becomes a function of the design variable of each element, and is optimized by the present invention. In the course of the calculation, the design variables are updated, and the values of the components of the rigidity matrix change.
The "pseudo-inverse matrix" of the rigidity matrix of the structure calculated by the "pseudo-inverse matrix calculation means" is a matrix also referred to as a Moore-Penrose inverse matrix or a general inverse matrix, as described above. Is uniquely calculated by a predetermined algorithm using the components of Specifically, the pseudo-inverse matrix may be calculated using a function command prepared in a programming language, or may be calculated using an arbitrary approximation method.
The “local compliance matrix” calculated by the “local compliance matrix calculating means” is the same as the local compliance matrix C in the above equation (1). It is a matrix that associates a local load vector having a load at each node as a component and a local displacement vector having a displacement at each node as a component. The local compliance matrix C originally associates a load vector representing loads at all nodes of the structure with a local load vector in the load application region, as will be described in detail in the section of a later embodiment. Calculated using the transformation matrix and the inverse matrix of the rigidity matrix of the structure. In the apparatus of the present invention, as described above, since the structure without the support structure is the object of optimization, , There is no inverse matrix of the rigidity matrix of the structure. Therefore, the "local compliance matrix calculation means" calculates the local compliance matrix using a pseudo inverse matrix of the rigidity matrix of the structure instead of the inverse matrix of the rigidity matrix of the structure.
次いで、「固有値解析手段」に於いては、上記の局所コンプライアンス行列の固有値と固有ベクトルとが算出されるところ、かかる算出処理は、任意の態様にて実行されてよく、典型的には、プログラム言語に用意された関数コマンドが利用されてよい。なお、好適には、固有値解析手段に於いて得られる固有ベクトルはそのノルム＝１のものが算出される（以下、特に断らない限り、固有ベクトルは、ノルム＝１のベクトルであるものとする。）。局所コンプライアンス行列の固有値と固有ベクトルは、典型的には、複数組算出され、本発明の装置に於いては、最大固有値とその固有ベクトルだけではなく、値の大きいものから複数の組の固有値と固有ベクトルが、後の処理に用いられてよい。なお、ここで得られた固有値と固有ベクトルは、それぞれ、荷重作用領域内の変位ベクトルが荷重ベクトルに平行な場合のコンプライアンスと荷重ベクトルを表しており、固有値が大きいほど、対応する固有ベクトルが荷重ベクトルとして作用したときのコンプライアンスが大きく、変形の程度が大きくなることを示している。また、後の実施形態の欄に於いて説明される如く、支持構造のない構造体の剛性行列の擬似逆行列を用いて算出される局所コンプライアンス行列の固有ベクトルは、構造体に於いて荷重の合計とモーメントの合計とが０となり、剛体運動を発生しない荷重ベクトルを与え、局所コンプライアンス行列の固有値は、そのときのコンプライアンスとなる。 Next, in the “eigenvalue analysis means”, the eigenvalue and the eigenvector of the above local compliance matrix are calculated, and such calculation processing may be executed in any mode. May be used. Preferably, the eigenvector obtained by the eigenvalue analysis means has a norm of 1 (hereinafter, unless otherwise specified, the eigenvector is a vector of norm = 1). A plurality of sets of eigenvalues and eigenvectors of the local compliance matrix are typically calculated. In the apparatus of the present invention, not only the largest eigenvalue and its eigenvector but also a plurality of sets of eigenvalues and eigenvectors from the largest value are calculated. , May be used in subsequent processing. The eigenvalue and eigenvector obtained here represent the compliance and the load vector, respectively, when the displacement vector in the load application area is parallel to the load vector.As the eigenvalue increases, the corresponding eigenvector becomes the load vector. This shows that the compliance when acting is large and the degree of deformation is large. Further, as described in the section of the embodiment below, the eigenvector of the local compliance matrix calculated using the pseudo inverse matrix of the rigidity matrix of the structure without the support structure is the sum of the loads in the structure. And the sum of the moments become 0, giving a load vector that does not generate rigid body motion, and the eigenvalue of the local compliance matrix becomes the compliance at that time.
更に、上記の本発明の装置の「変位ベクトル算出手段」では、上記の固有ベクトルと擬似逆行列とを用いて構造体全体に於ける各節点の変位を表す「変位ベクトル」算出される。かかる「変位ベクトル」は、本来であれば、式（２）に於ける荷重ベクトルＦに固有ベクトルを代入し、式（２）の解Ｕを演算することにより、即ち、
Ｕ＝Ｋ−１Ｆ …（３）
を演算することにより算出されるが、既に述べた如く、本発明に於いては、式（２）の剛性行列Ｋの逆行列Ｋ−１が存在しない。そこで、本発明では、剛性行列Ｋの逆行列Ｋ−１に代えて、剛性行列Ｋの疑似逆行列Ｋ＋を用いて、
Ｕ＝Ｋ＋Ｆ …（４）
を演算することにより、Ｕとして「変位ベクトル」が算出される。なお、「変位ベクトル」は、局所コンプライアンス行列の複数得られた固有ベクトルのそれぞれに対して、算出されてよい。荷重ベクトルから変位ベクトルを算出する具体的な演算手法は、剛性行列の逆行列に代えて擬似逆行列を使用することの他は、任意の手法にて、例えば、この分野に於いて通常利用されている有限要素法の手法にて実行されてよい。上記の如く局所コンプライアンス行列Ｃの固有ベクトルを荷重ベクトルとして疑似逆行列Ｋ＋を用いて算出された変位ベクトルＵは、構造体に於ける剛体運動を除いた構造体内に於いて対称的に生ずる変位を表すこととなる。
Further, the "displacement vector calculating means" of the apparatus of the present invention calculates a "displacement vector" representing the displacement of each node in the entire structure using the eigenvector and the pseudo inverse matrix. The “displacement vector” is originally calculated by substituting the eigenvector for the load vector F in the equation (2) and calculating the solution U of the equation (2),
U = K −1 F (3)
As described above, in the present invention, there is no inverse matrix K −1 of the stiffness matrix K in Expression (2). Therefore, in the present invention, instead of the inverse matrix K −1 of the stiffness matrix K, a pseudo inverse matrix K + of the stiffness matrix K is used,
U = K + F (4)
Is calculated, the “displacement vector” is calculated as U. The “displacement vector” may be calculated for each of a plurality of obtained eigenvectors of the local compliance matrix. A specific calculation method for calculating a displacement vector from a load vector is to use a pseudo-inverse matrix instead of the inverse matrix of the stiffness matrix. It may be performed by a finite element method. As described above, the displacement vector U calculated using the pseudo-inverse matrix K + with the eigenvector of the local compliance matrix C as the load vector is the displacement symmetrically generated in the structure excluding the rigid body motion in the structure. Will be represented.
そして、上記の装置に於いて、変位ベクトルを用いて「固有値感度算出手段」にて設計変数に対する固有値、即ち、構造体のコンプライアンスの感度が算出され、その感度を用いて「設計変数更新手段」にて、構造体のコンプライアンスである固有値が最小化されるよう設計変数が更新される。即ち、本発明の装置に於いて、「固有値」（又はそれらの線形和）が構造最適化演算に於ける「評価関数」となる。また、ここにおいて、「感度」とは、設計変数の変化に対する固有値の変化量、即ち、設計変数についての固有値の微分に相当する量であり、変位ベクトルと、設計変数についての剛性行列の成分の微分値とから算出可能である。設計変数の更新は、従前の構造最適化演算に於いて用いられている手法（例えば、最適性基準法、逐次線形計画法、逐次凸関数近似法など）と同様であってよく、設計変数は、評価関数である固有値の感度を用いて、構造体が、初めに設定された体積制約割合（設計領域内にて最終的に形成される構造体の体積の、設計領域の体積に対する割合）を満たすことを数学的な拘束条件として、構造体の剛性が最大化されるよう更新されることとなる。 Then, in the above apparatus, the eigenvalue for the design variable, that is, the sensitivity of the compliance of the structure is calculated by the “eigenvalue sensitivity calculating means” using the displacement vector, and the “design variable updating means” is calculated using the sensitivity. In, the design variables are updated so that the eigenvalue that is the compliance of the structure is minimized. That is, in the apparatus of the present invention, the “eigenvalue” (or a linear sum thereof) becomes the “evaluation function” in the structure optimization operation. Here, the “sensitivity” is an amount of change of the eigenvalue with respect to the change of the design variable, that is, an amount corresponding to the differentiation of the eigenvalue of the design variable. It can be calculated from the differential value. The update of the design variables may be the same as the method used in the conventional structure optimization operation (eg, the optimality criterion method, the sequential linear programming method, the sequential convex function approximation method, etc.) Using the sensitivity of the eigenvalue, which is an evaluation function, the structure uses the initially set volume constraint ratio (the ratio of the volume of the structure finally formed in the design region to the volume of the design region). Satisfaction is a mathematical constraint, and the structure is updated so that the rigidity of the structure is maximized.
上記の本発明の装置に於いて、反復処理手段による更新された設計変数を用いた、剛性行列算出手段、擬似逆行列算出手段、局所コンプライアンス行列算出手段、固有値解析手段、変位ベクトル算出手段、固有値感度算出手段及び設計変数更新手段の各処理の反復は、所定の終了条件が成立するまで実行される。ここで、「所定の終了条件」は、目的の構造体の形状が得られるように任意に設定されてよいところ、具体的には、例えば、（ｉ）反復処理回数が任意に設定されてよい回数に到達した場合、（ii）構造体にて算出される評価関数（固有値又はそれらの線形和）の値が収束した場合（処理サイクル毎の評価関数の変化量が所定の小数以下となった場合）、（iii）更新される設計変数の値が収束した場合（処理サイクル毎の設計変数の変化量の最大値が所定の小数以下となった場合）などが選択されてよい。更に、本発明の装置に於いては、上記の如く「構造体出力手段」が設けられ、演算処理過程に於いて得られた設計変数により表される構造体の構造又は形状（例えば、材料密度の分布、有意な材料密度を有する有限要素の分布）がコンピュータのモニターの画面上や紙面にてグラフィックス表示されるようになっていてよい。 In the above-described apparatus of the present invention, the stiffness matrix calculating unit, the pseudo-inverse matrix calculating unit, the local compliance matrix calculating unit, the eigenvalue analyzing unit, the displacement vector calculating unit, the eigenvalue using the design variables updated by the iterative processing unit. The repetition of each process of the sensitivity calculating means and the design variable updating means is performed until a predetermined end condition is satisfied. Here, the “predetermined termination condition” may be arbitrarily set so as to obtain the shape of the target structure. Specifically, for example, (i) the number of repetition processes may be arbitrarily set. When the number of times is reached, (ii) when the value of the evaluation function (eigenvalue or linear sum thereof) calculated by the structure converges (the amount of change in the evaluation function for each processing cycle is less than or equal to a predetermined decimal number) Case), (iii) a case where the value of the design variable to be updated converges (a case where the maximum value of the change amount of the design variable for each processing cycle becomes equal to or smaller than a predetermined decimal number), or the like may be selected. Further, in the apparatus of the present invention, the "structure output means" is provided as described above, and the structure or shape (for example, material density) of the structure represented by the design variables obtained in the arithmetic processing process is provided. (A distribution of finite elements having a significant material density) may be displayed graphically on a screen of a computer monitor or on paper.
かくして、上記の本発明の装置では、基本的な構成としては、従前のロバスト構造最適化演算の場合と同様に、設計領域内に於ける構造体の構造を決定する設計変数を逐次的に更新しながら、局所コンプライアンス行列の固有値解析により構造体のコンプライアンスに相当する局所コンプライアンス行列の固有値とその固有値を与える固有ベクトルを決定し、その固有ベクトルを荷重ベクトルとして有限要素法等によって各要素の節点の変位ベクトルを算出し、かかる変位ベクトルを用いて評価関数である局所コンプライアンス行列の固有値が最小化されるように設計変数の値を探索し、そこで得られた設計変数にて表される形状が目的の構造体の形状として得られることとなる。かかる一連の演算処理を実行する構成に於いて、最適化の対象となる構造体が支持構造のない構造体である場合、構造体の剛性行列の逆行列がなく、局所コンプライアンス行列の算出及び変位ベクトルの算出ができないところ、本発明の装置では、後により詳細に説明される如く、構造体に支持構造のない場合でも算出することのできる構造体の剛性行列の擬似逆行列を用いると、支持構造のない構造体に於いて剛体運動を除いた変位ベクトル、コンプライアンス及びその感度を算出することができるという知見に基づき、局所コンプライアンス行列の算出及び変位ベクトルの算出に於いて、擬似逆行列を用い、これにより、支持構造のない構造体に於いて剛体運動を除いた構造体内で生ずる変位によるコンプライアンスを最小化する構造最適化演算の実行を可能にしている。 Thus, in the above-described apparatus of the present invention, as a basic configuration, the design variables for determining the structure of the structure in the design area are sequentially updated as in the case of the conventional robust structure optimization operation. The eigenvalue analysis of the local compliance matrix determines the eigenvalues of the local compliance matrix corresponding to the compliance of the structure and the eigenvectors that give the eigenvalues, and the eigenvectors are used as load vectors, and the displacement vectors of the nodes of the respective elements are determined by the finite element method or the like. Is calculated, and the value of the design variable is searched using the displacement vector so that the eigenvalue of the local compliance matrix, which is the evaluation function, is minimized. It will be obtained as a body shape. In the configuration for executing such a series of arithmetic processing, when the structure to be optimized is a structure without a support structure, there is no inverse matrix of the rigidity matrix of the structure, and calculation and displacement of the local compliance matrix are performed. Where the vector cannot be calculated, the apparatus of the present invention uses a pseudo-inverse matrix of the rigidity matrix of the structure that can be calculated even when the structure does not have a support structure, as will be described in more detail later. Based on the knowledge that displacement vectors excluding rigid body motion, compliance, and their sensitivities can be calculated in a structure without a structure, a pseudo-inverse matrix is used in the calculation of the local compliance matrix and the calculation of the displacement vector. Thus, in a structure without a support structure, a structure that minimizes compliance due to displacement occurring in the structure excluding the rigid body motion. And it enables the execution of the reduction operation.
ところで、上記の如く、局所コンプライアンス行列の固有値と固有ベクトルは複数組得られるところ、最適化の対象とされた構造体に支持構造があり、構造体の剛体運動が制限されている場合には、構造体の構造が最適化されていく過程で局所コンプライアンス行列の最大固有値は、常に或る一つの固有ベクトルにより与えられ、その一つの固有ベクトルが連続的に変化することとなる。即ち、支持構造のある構造体の場合には、構造最適化の過程に於いて、複数の固有値と固有ベクトルの組のうちで最大固有値を与える固有ベクトルが交代することはなく、最初に固有値を最大にした固有ベクトルが、その後の反復演算の間も、対応する固有値が最大となったまま連続的に変化することとなる。従って、変位ベクトルの算出が反復して実行される間に於いて、常に最大固有値の固有ベクトルを荷重ベクトルとして選択しておけば、算出される変位ベクトル及びそれを用いて算出される評価関数の感度も連続的に変化し、その感度を用いて設計変数の更新を実行することにより、構造体の構造は連続的に変化して最適な構造に収束することが可能となっていた。 By the way, as described above, a plurality of sets of eigenvalues and eigenvectors of the local compliance matrix can be obtained.If a structure targeted for optimization has a support structure and rigid motion of the structure is restricted, the structure In the process of optimizing the structure of the body, the maximum eigenvalue of the local compliance matrix is always given by a certain eigenvector, and the one eigenvector changes continuously. That is, in the case of a structure having a supporting structure, in the process of structural optimization, the eigenvector that gives the maximum eigenvalue among a plurality of eigenvalues and eigenvector sets does not alternate, and the eigenvalue is first maximized. The eigenvector changes continuously during the subsequent iterative operation while the corresponding eigenvalue is maximized. Therefore, if the eigenvector having the largest eigenvalue is always selected as the load vector while the calculation of the displacement vector is repeatedly executed, the sensitivity of the calculated displacement vector and the evaluation function calculated using the same is calculated. Also changes continuously, and the design variable is updated using the sensitivity, whereby the structure of the structure can change continuously and converge to an optimal structure.
一方、支持構造のない構造体の場合には、構造最適化の過程に於いて、局所コンプライアンス行列の複数の固有値と固有ベクトルの組のうちで最大固有値を与える固有ベクトルの交代が生じ得ることが明らかとなった。そして、かかる固有ベクトルの交代が生ずると、最大固有値を与える固有ベクトルが不連続に変化し、構造体の構造が最適な構造へ収束しづらくなる状況が発生した。そこで、本発明の装置に於いては、局所コンプライアンス行列の複数の固有値と固有ベクトルの組のうち、固有値の大きい方から複数個の組を選択し、それらの選択された固有値が全体的に最小化されるように設計変数の更新が実行されるようになっていてよい。そのために、本発明の装置に於いては、構造最適化の評価関数として、固有値の大きい方から選択された複数の固有値の線形和が採用されてよい。具体的な装置の構成に於いては、固有値の大きい方から複数個の組の固有ベクトルを選択した後、選択された組の固有ベクトルの各々に対して変位ベクトルの算出及び感度の算出が実行され、構造体全体の固有値の感度として、各固有値の感度の線形和が算出され、その構造体全体の固有値の感度を用いて、選択された固有値の線形和が評価関数として最小化されるように設計変数の更新が実行されてよい。なお、複数の固有値及びそれらの感度の線形和は、重み付き和であってよく、重み付き和に於ける各項の重みは、上記の如く、個々の固有値が大きいほど重み付き和に於ける対応する固有値に対する重みが大きくなるよう設定されてよい。このように重みを設定すると、より大きな固有値、即ち、より大きなコンプライアンスを与える荷重ベクトルに対する剛性を高くして、かかる荷重ベクトルに対する変形を最小化することが可能となる。 On the other hand, in the case of a structure without a support structure, it is clear that in the process of structural optimization, alternation of the eigenvector giving the largest eigenvalue among a plurality of pairs of eigenvalues and eigenvectors of the local compliance matrix may occur. became. Then, when such alternation of the eigenvectors occurs, the eigenvector giving the maximum eigenvalue changes discontinuously, and a situation has arisen in which the structure of the structure hardly converges to the optimal structure. Therefore, in the apparatus of the present invention, a plurality of sets of a plurality of eigenvalues and eigenvectors of the local compliance matrix are selected from the largest eigenvalue, and the selected eigenvalues are minimized as a whole. The update of the design variable may be performed so that the design variable is updated. For this purpose, in the apparatus of the present invention, a linear sum of a plurality of eigenvalues selected from the larger eigenvalues may be adopted as the evaluation function of the structure optimization. In the specific configuration of the device, after selecting a plurality of sets of eigenvectors from the largest eigenvalue, the calculation of the displacement vector and the calculation of the sensitivity is performed for each of the selected set of eigenvectors, The linear sum of the sensitivities of each eigenvalue is calculated as the sensitivity of the eigenvalues of the entire structure, and the design is made so that the linear sum of the selected eigenvalues is minimized as an evaluation function using the sensitivity of the eigenvalues of the entire structure. A variable update may be performed. Note that the linear sum of a plurality of eigenvalues and their sensitivities may be a weighted sum, and the weight of each term in the weighted sum is, as described above, the greater the individual eigenvalue, the more the weighted sum is The weight for the corresponding eigenvalue may be set to be large. By setting the weights in this way, it is possible to increase the rigidity of a load vector that gives a larger eigenvalue, that is, a larger compliance, and minimize the deformation of the load vector.
かくして、上記の本発明の装置の固有値感度算出手段に於いて、固有値を二つ以上用いるときには、固有値の感度が、局所コンプライアンス行列の複数の固有値のうちの値の大きい方から選択された複数の固有値の設計変数に対する感度の線形和であってよく、設計変数更新手段に於いて、固有値を二つ以上用いるときには、局所コンプライアンス行列の複数の固有値のうちの値の大きい方から選択された複数の固有値の線形和が最小化されるよう設計変数が更新されるようになっていてよい。また、上記の固有値及びその感度の線形和は、重み付き和であってよく、その場合の各項の重みは、個々の固有値が大きいほど重み付き和に於ける対応する固有値に対する重みが大きくなるよう設定されていてよい。上記の如く値の大きい複数の固有値が全体的に最小化されるように設計変数の更新を実行する場合には、構造体の構造がより安定的に最適な構造へ収束されることとなる。なお、設定された条件によって、最大固有値を与える固有値の交代がない場合或いはそのように想定される場合には、評価関数を最大固有値とし、設計変数に対する感度は、最大固有値の感度であってもよく、そのような場合も本発明の範囲に属することは理解されるべきである。（その場合には、局所コンプライアンス行列の固有値解析に於いて、最大固有値とそれに対応する固有ベクトルが選択される。） Thus, when two or more eigenvalues are used in the eigenvalue sensitivity calculation means of the above-described apparatus of the present invention, the sensitivity of the eigenvalues is determined based on the plurality of eigenvalues selected from the larger one of the eigenvalues of the local compliance matrix. A linear sum of the sensitivity of the eigenvalues to the design variables may be used. In the design variable updating means, when two or more eigenvalues are used, a plurality of eigenvalues selected from the larger one of the plurality of eigenvalues of the local compliance matrix are used. The design variables may be updated so that the linear sum of the eigenvalues is minimized. The linear sum of the eigenvalues and the sensitivities may be a weighted sum. In this case, the weight of each term is such that the larger the individual eigenvalue is, the larger the weight for the corresponding eigenvalue in the weighted sum is. May be set as follows. When the design variables are updated such that a plurality of large eigenvalues are minimized as a whole, the structure of the structure is more stably converged to the optimum structure. In addition, according to the set conditions, when there is no change of the eigenvalue that gives the maximum eigenvalue, or in such a case, the evaluation function is set to the maximum eigenvalue, and the sensitivity to the design variable is the sensitivity of the maximum eigenvalue. It should be understood that such cases are well within the scope of the present invention. (In that case, the maximum eigenvalue and the corresponding eigenvector are selected in the eigenvalue analysis of the local compliance matrix.)
上記の本発明によれば、支持構造又は支点のない構造体の場合、即ち、構造体の剛性行列に逆行列が存在しない場合でも構造最適化が可能な手法を提供される。かかる本発明の装置によれば、構造体に作用する荷重が不定である場合だけでなく、構造体を支持する構造が不定である場合でも、任意の荷重が付与される場合に構造体の剛性が最大化される最適な構造を見出すことが可能となると期待される。また、その際、支持構造のある構造体の場合と異なり、局所コンプライアンス行列の最大固有値を与える固有ベクトルが設計変数の更新に伴う構造体の構造の変化に伴って交代することがあり、かかる固有ベクトルの交代があると、構造体の構造の収束が困難となることがあるところ、本発明の装置は、局所コンプライアンス行列の複数の固有値が全体的に最小化されるように評価関数の感度及び設計変数の更新を実行するようになっていてよく、これにより、最大固有値を与える固有ベクトルの交代があっても、構造体の構造を収束させて、最適な構造が得られるようになっていてよい。 According to the present invention described above, there is provided a method capable of optimizing the structure even in the case of a structure having no support structure or fulcrum, that is, even when there is no inverse matrix in the rigidity matrix of the structure. According to such an apparatus of the present invention, not only when the load acting on the structure is indefinite, but also when the structure supporting the structure is indefinite, the rigidity of the structure when an arbitrary load is applied It is expected that it will be possible to find the optimal structure in which is maximized. Also, at this time, unlike the case of a structure having a support structure, the eigenvector giving the maximum eigenvalue of the local compliance matrix may alternate with the change in the structure of the structure due to the update of the design variable, and the eigenvector of the eigenvector may be changed. Where the alternation can make the convergence of the structure of the structure difficult, the apparatus of the present invention provides the sensitivity and design variables of the evaluation function such that the multiple eigenvalues of the local compliance matrix are globally minimized. May be executed, so that even if the eigenvector giving the maximum eigenvalue is changed, the structure of the structure may be made to converge to obtain an optimal structure.
本発明のその他の目的及び利点は、以下の本発明の好ましい実施形態の説明により明らかになるであろう。 Other objects and advantages of the present invention will become apparent from the following description of preferred embodiments of the present invention.
１…コンピュータ本体
２…コンピュータ端末
３…モニター
４…キーボード、マウス（入力装置）
1. Computer main body 2. Computer terminal 3. Monitor 4. Keyboard and mouse (input device)
コンピュータ装置の構成
本発明による構造最適化演算を実行する装置は、この分野で通常使われている形式の、図１に例示されている如き、コンピュータ装置１に於けるコンピュータ・プログラムの作動により実現されてよい。コンピュータ装置１には、通常の態様にて、双方向コモン・バスにより相互に連結されたＣＰＵ、記憶装置、入出力装置（Ｉ／Ｏ）が装備され、記憶装置は、本発明の演算で使用する演算処理を実行する各プログラムを記憶したメモリと、演算中に使用されるワークメモリ及びデータメモリを含んでいる。また、実施者によるコンピュータ装置１への指示及び計算結果その他の情報の表示及び出力は、コンピュータ装置１に接続されたコンピュータ端末装置２を通じて為される。コンピュータ端末装置２には、通常の態様にて、モニター３とキーボード並びにマウスといった入力装置４が設けられ、プログラムが起動されると、実施者は、プログラムの手順に従ってモニター３上の表示に従って、入力装置４を用いてコンピュータ装置１に各種の指示及び入力を行うとともに、モニター３上にてコンピュータ装置１からの演算状態及び演算結果等を視覚的に確認することが可能となる。なお、図示していないその他の周辺機器（結果を出力するプリンタ、計算条件及び演算結果情報等を入出力するための記憶装置など）がコンピュータ装置１及びコンピュータ端末装置２に装備されていてよいことは理解されるべきである。コンピュータ装置１を用いて、以下に述べる各種の処理又は演算を実行する際には、通常の態様にて、各種の処理又は演算に必要なプログラムが起動され、実施者は、コンピュータ端末装置２に於いて、プログラムに於いて準備された入力手順に従って、演算に必要なデータ、演算時の計算条件その他の各種設定を入力し、演算が開始される。そして、演算の実行中又は終了後に、演算結果が、適宜、コンピュータ端末装置２を通じて出力可能となる。
Configuration of Computer Apparatus An apparatus for performing a structure optimization operation according to the present invention is realized by the operation of a computer program in a computer apparatus 1 as shown in FIG. 1 of a type commonly used in this field. May be. The computer device 1 is equipped with a CPU, a storage device, and an input / output device (I / O) interconnected by a bidirectional common bus in a usual manner, and the storage device is used in the calculation of the present invention. And a work memory and a data memory used during the calculation. In addition, instructions by the practitioner to the computer device 1 and display and output of calculation results and other information are made through the computer terminal device 2 connected to the computer device 1. The computer terminal device 2 is provided with a monitor 3 and an input device 4 such as a keyboard and a mouse in a usual manner. When the program is started, the practitioner inputs data according to the display on the monitor 3 according to the procedure of the program. Using the device 4, various instructions and inputs can be made to the computer device 1, and it is possible to visually confirm the calculation state and the calculation result from the computer device 1 on the monitor 3. Note that other peripheral devices (not shown) such as a printer for outputting results, a storage device for inputting and outputting calculation conditions and calculation result information, and the like may be provided in the computer device 1 and the computer terminal device 2. Should be understood. When various processes or calculations described below are performed using the computer device 1, programs necessary for various processes or calculations are activated in a normal manner, and the practitioner operates the computer terminal device 2. According to the input procedure prepared in the program, data necessary for calculation, calculation conditions at the time of calculation and other various settings are input, and the calculation is started. Then, during or after the execution of the calculation, the calculation result can be output through the computer terminal device 2 as appropriate.
上記のコンピュータ装置１に於いて、本発明による構造最適化演算を行う場合、端的に述べれば、まず、入力装置４からの実施者の入力操作によって、最適化された構造体が生成される設計領域の寸法、構造体を形成する材料の特性値（ヤング率、ポアソン比）及び体積制約割合を含む設計条件や構造体に作用する荷重条件を含む境界条件などの構造最適化演算を実行するために必要な初期条件の設定が実行される。そうすると、コンピュータ装置１は、構造最適化演算処理の準備として、プログラムメモリに記憶されたプログラムを用いて、設計領域を、典型的には格子状に、複数の要素に分割する処理を実行すると共に、得られた要素の各々への設計変数の初期値の割り当てを実行する。かくして、かかる演算処理の準備が完了すると、コンピュータ装置１は、プログラムメモリに記憶されたプログラムに従い、（ａ）構造体の剛性行列Ｋ、擬似逆行列Ｋ＋、局所コンプライアンス行列Ｃの算出、（ｂ）局所コンプライアンス行列Ｃの固有値解析、（ｃ）各要素の各節点に生ずる変位の算出、（ｄ）設計領域内にて形成される構造体の評価関数の感度（各要素の設計変数の変化に対する評価関数の変化量）の算出、（ｅ）評価関数が最適化されるように各要素の設計変数を更新する処理を含む一連の演算処理を、後述の終了条件が成立するまで反復して実行し、これにより、設計変数により形成される構造体の構造又は形状が、逐次的に最適化された状態に近づいていくことが期待される。また、一連の演算処理の結果として得られた構造体の構造又は形状、設計変数の分布、局所コンプライアンス行列Ｃの固有値の推移、固有ベクトルの成分の分布などのグラフィック表示がモニター３上に表示されてよい。 In the case of performing the structure optimization operation according to the present invention in the computer apparatus 1 described above, in short, first, a design in which an optimized structure is generated by an input operation of a practitioner from the input device 4. To execute structural optimization calculations such as design conditions including the dimensions of the region, the characteristic values (Young's modulus and Poisson's ratio) of the material forming the structure and the volume constraint ratio, and the boundary conditions including the load conditions acting on the structure Of initial conditions necessary for the operation are executed. Then, the computer device 1 executes a process of dividing the design area into a plurality of elements, typically in a lattice pattern, using a program stored in the program memory, in preparation for the structure optimization arithmetic processing. And assigning initial values of design variables to each of the obtained elements. Thus, when the preparation for the arithmetic processing is completed, the computer apparatus 1 calculates (a) the stiffness matrix K, the pseudo inverse matrix K + , and the local compliance matrix C of the structure according to the program stored in the program memory; ) Eigenvalue analysis of the local compliance matrix C, (c) calculation of displacement occurring at each node of each element, (d) sensitivity of the evaluation function of the structure formed in the design area (with respect to changes in the design variables of each element) (E) calculating a change in the evaluation function, and (e) repeatedly executing a series of calculation processes including a process of updating design variables of each element so that the evaluation function is optimized until an end condition described later is satisfied. Accordingly, it is expected that the structure or shape of the structure formed by the design variables will gradually approach an optimized state. Further, graphic displays such as the structure or shape of the structure obtained as a result of a series of arithmetic processing, distribution of design variables, transition of the eigenvalue of the local compliance matrix C, and distribution of the components of the eigenvector are displayed on the monitor 3. Good.
ロバスト構造最適化演算の概要
本発明の装置の実施形態に於ける構造最適化演算では、基本的には、ロバスト構造最適化演算の手法（非特許文献３、４）に従って、不定の荷重の存在下に於いて、剛性を最大にする、即ち、変形を最小にする構造体の構造を最適な構造として求めるために、構造体の変形程度を表す構造体のコンプライアンスを最大にする荷重（最悪荷重）の存在下に於いて、構造体のコンプライアンスが最小化されるように、設計領域内の構造体の体積が所与の体積制約割合の条件を満たすことを条件として、構造体の構造を決定する設計変数の更新する演算処理が実行され、更新された設計変数により表される構造体の構造が目的の構造体の構造として見出されるようになっていてよい。
Outline of Robust Structure Optimization Operation In the structure optimization operation in the embodiment of the apparatus of the present invention, basically, the existence of an indefinite load is determined according to the robust structure optimization operation method (Non-Patent Documents 3 and 4). Below, in order to maximize the rigidity, that is, to determine the structure of the structure that minimizes the deformation as the optimal structure, a load that maximizes the compliance of the structure indicating the degree of deformation of the structure (the worst load) ), The structure of the structure is determined, provided that the volume of the structure in the design domain satisfies the given volume constraint ratio so that the compliance of the structure is minimized The operation of updating the design variable to be executed may be executed, and the structure of the structure represented by the updated design variable may be found as the structure of the target structure.
より具体的には、まず、ロバスト構造最適化演算に於ける最小化の対象となる評価関数となるコンプライアンスｌは、
ｌ＝ＦＴ・Ｕ …（５）
により与えられる。ここに於いて、ＦとＵは、構造体に於いて格子状の要素を設定した場合の要素の節点に於ける荷重を成分とする荷重ベクトル（上付きの添え字Ｔは、転置を表す。）と変位を成分とする変位ベクトルであり、既に述べた如く、構造体に於ける剛性行列Ｋを用いた式（２）［ＫＵ＝Ｆ］にて互いに関係付けられる。上記の式（５）のコンプライアンスｌに於いて、荷重ベクトルＦが０以外の値を持つのは、構造体に於いて荷重の作用する「荷重作用領域」内の節点のみとなるので、構造体のコンプライアンスｌは、荷重作用領域内の荷重分布を表す荷重ベクトルＦｌと荷重作用領域の変位分布を表す変位ベクトルＵｌを用いて
ｌ＝Ｆｌ Ｔ・Ｕｌ …（６）
により表される。そして、荷重ベクトルＦｌと変位ベクトルＵｌは、式（１）により局所コンプライアンス行列Ｃを用いて互いに関係付けられるので、コンプライアンスｌは、
ｌ＝Ｆｌ Ｔ・Ｃ・Ｆｌ …（７）
により与えられる。このコンプライアンスｌに関し、その値が最大となるとき、即ち、構造体に任意の荷重が作用する場合にその変形が最大となるときは、荷重ベクトルＦｌが局所コンプライアンス行列Ｃの最大固有値を与える固有ベクトルとなっているときであり、荷重ベクトルＦｌのノルム||Fl||＝１としたとき、コンプライアンスｌは、局所コンプライアンス行列Ｃの最大固有値となることが、Ｒａｙｌｅｉｇｈ−Ｒｉｔｚの定理により理論的に示されている（非特許文献３、４参照）。従って、不定の荷重に対して構造体の剛性を最大化する構造最適化演算に於いては、荷重ベクトルが局所コンプライアンス行列Ｃの最大固有値を与える固有ベクトルとなっているときのコンプライアンスｌを最小化するように構造体の構造を決定する設計変数が更新されることとなる。
More specifically, first, the compliance l which is an evaluation function to be minimized in the robust structure optimization operation is
l = F TU (5)
Given by Here, F and U are load vectors (the superscript suffix T indicates transposition) that has a load at a node of the element when a lattice-like element is set in the structure. ) And a displacement vector having displacement as components, and as described above, are related to each other by equation (2) [KU = F] using the rigidity matrix K in the structure. In the compliance 1 of the above equation (5), the load vector F has a value other than 0 only at the node in the “load application area” where the load acts on the structure. Is calculated using a load vector F 1 representing a load distribution in the load application area and a displacement vector U 1 representing a displacement distribution in the load application area. 1 = F 1 T · U 1 (6)
Is represented by Then, the load vector F 1 and the displacement vector U 1 are related to each other using the local compliance matrix C according to the equation (1).
l = F l T · C · F l (7)
Given by When the value of the compliance 1 is maximized, that is, when the deformation is maximized when an arbitrary load is applied to the structure, the load vector F 1 is an eigenvector that gives the maximum eigenvalue of the local compliance matrix C. is when that is the, when a norm || F l || = 1 the load vector F l, compliance l may be the maximum eigenvalue of the local compliance matrix C, the theoretical theorem of Rayleigh-Ritz (See Non-Patent Documents 3 and 4). Therefore, in the structural optimization operation for maximizing the rigidity of the structure with respect to an indefinite load, the compliance l when the load vector is the eigenvector that gives the maximum eigenvalue of the local compliance matrix C is minimized. Thus, the design variables for determining the structure of the structure are updated.
ところで、式（１）に於ける局所コンプライアンス行列Ｃに関して、荷重作用領域内の荷重ベクトルＦｌと構造体全体の荷重ベクトルＦとを関連づける変換行列Ｈを
Ｆ＝Ｈ・Ｆｌ …（８ａ）
Ｆ＝ＨＴ・Ｆ …（８ｂ）
と置くと、荷重作用領域内の変位ベクトルＵｌは、変換行列Ｈと構造体の剛性行列Ｋの逆行列Ｋ−１を用いて、
Ｕｌ＝ＨＴ・Ｕ＝ＨＴ・Ｋ−１・Ｆ＝ＨＴ・Ｋ−１・Ｈ・Ｆｌ …（９）
と表されるので、式（１）から、局所コンプライアンス行列Ｃは、
Ｃ＝ＨＴ・Ｋ−１・Ｈ …（１０）
により与えられる。なお、変換行列Ｈは、設計領域に対して指定される荷重作用領域の範囲に基づいて決定される。
Incidentally, wherein with respect to at local compliance matrix C in (1), the transformation matrix H relating the load vector F l and the entire structure load vector F of the load acting in the area F = H · F l ... (8a)
F = H T · F (8b)
, The displacement vector U 1 in the load application area is calculated using the transformation matrix H and the inverse matrix K −1 of the rigidity matrix K of the structure.
U l = H T U = H T K- 1 · F = H T K- 1 · H · F 1 (9)
From Equation (1), the local compliance matrix C is
C = H T · K −1 · H (10)
Given by Note that the conversion matrix H is determined based on the range of the load application area specified for the design area.
従って、ロバスト構造最適化演算では、設定された初期条件を用いて式（２）の剛性行列Ｋを決定した後、その逆行列Ｋ−１を用いて式（１０）により局所コンプライアンス行列Ｃを算出し、局所コンプライアンス行列Ｃの固有値解析を実行し、そこで得られた最大固有値を与える固有ベクトルを荷重ベクトルＦとして用いて
Ｕ＝Ｋ−１・Ｆ …（１１）
により、式（２）を解いて構造体の変位ベクトルＵを算出し、しかる後に、算出された変位ベクトルＵを用いて、コンプライアンスｌを最小化するように設計変数の更新が実行される。
Therefore, in the robust structure optimization calculation, after determining the stiffness matrix K of Expression (2) using the set initial conditions, the local compliance matrix C is calculated by Expression (10) using the inverse matrix K −1. Then, the eigenvalue analysis of the local compliance matrix C is performed, and the obtained eigenvector giving the maximum eigenvalue is used as the load vector F, and U = K −1 · F (11)
Then, the displacement vector U of the structure is calculated by solving the equation (2), and thereafter, the design variables are updated using the calculated displacement vector U so as to minimize the compliance l.
支持構造のない構造体の構造最適化演算
（ａ）支持構造のない構造体に於ける変位ベクトルの算出-擬似逆行列Ｋ＋を用いた演算
「発明の概要」の欄に於いて既に述べた如く、現実の任意の構造体に於いては、ユーザーの手からの荷重を受ける携帯電話の外殻など、構造体を支持し或いは固定する構造や支点の位置が不定な場合がある。従って、構造最適化演算に於いては、荷重だけでなく、支持構造も不定である場合、即ち、支持構造のない状態の構造体についても、構造最適化演算が実行できることが好ましい。
Structural optimization calculation of a structure without a support structure (a) Calculation of a displacement vector in a structure without a support structure-calculation using a pseudo-inverse matrix K + As already described in the “Summary of the Invention” section As described above, in any actual structure, the structure for supporting or fixing the structure or the position of a fulcrum, such as the outer shell of a mobile phone, which receives a load from the user's hand, may be indeterminate. Therefore, in the structure optimization calculation, it is preferable that not only the load but also the support structure is undetermined, that is, it is preferable that the structure optimization calculation can be performed even for a structure having no support structure.
しかしながら、構造体に於いて支持構造がない場合、構造体は自由な平行移動及び回転（剛体運動又は剛体変位）が許されるので、或る荷重ベクトルＦが与えられた場合の式（２）［ＫＵ＝Ｆ］の解である変位ベクトルＵは、一意に定まらなくないこととなる。即ち、支持構造がない構造体に於いて、数学的には、剛性行列Ｋの逆行列Ｋ−１が存在しないため、上記のロバスト構造最適化演算の実行に際して、逆行列Ｋ−１を用いる式（１０）による局所コンプライアンス行列Ｃの算出及び式（１１）による変位ベクトルＵの算出を実行することができないこととなっている。 However, when there is no support structure in the structure, the structure is allowed to freely translate and rotate (rigid body movement or rigid body displacement). Therefore, when a certain load vector F is given, equation (2) [ KU = F], the displacement vector U is not uniquely determined. That is, in a structure having no support structure, mathematically, there is no inverse matrix K- 1 of the stiffness matrix K. Therefore, when the above robust structure optimization operation is performed, an expression using the inverse matrix K- 1 is used. The calculation of the local compliance matrix C by (10) and the calculation of the displacement vector U by Expression (11) cannot be performed.
ところで、式（２）［ＫＵ＝Ｆ］の如き連立一次方程式に於いて、方程式の係数を成分とする行列Ｋに逆行列Ｋ−１が存在しない場合に擬似逆行列Ｋ＋を用いることにより方程式の解を得ることが可能である。擬似逆行列Ｋ＋の性質から、式（２）［ＫＵ＝Ｆ］が解を有するとき、解Ｕは、任意のベクトルＵｗを用いて、
Ｕ＝Ｋ＋・Ｆ＋（Ｉ−Ｋ＋・Ｋ）Ｕｗ …（１２）
により表される（Ｉは、単位行列である。）。ここで、構造体の剛性方程式である式（２）［ＫＵ＝Ｆ］が静的な状態を表し、荷重ベクトルＦによって平行及び回転について加速運動が生じていないとすると、式（１２）の右辺の第２項（Ｉ−Ｋ＋・Ｋ）・Ｕｗは、剛性行列Ｋを乗算すると、擬似逆行列Ｋ＋の定義から
Ｋ・（Ｉ−Ｋ＋・Ｋ）・Ｕｗ＝０ …（１３）
となり、構造体に対して力が作用していない状態でも生じ得る構造体の剛体運動を表す変位ベクトルであると考えられる。一方、式（１２）の右辺の第１項Ｋ＋・Ｆは、この場合、剛体運動以外の構造体内に生ずる変位であると考えられ、従って、荷重ベクトルＦは、剛体運動を生じさせない構造体内に於いて対称的に変位を生ずる荷重ベクトルとなっていると考えられる（もし対照的でない変位があったとすると、支持構造がないので、剛体運動となっているはずである。）。即ち、荷重ベクトルＦが、図２に模式的に例示されている如く、構造体に対して荷重の合計とモーメントの合計とがそれぞれ０となるように作用している荷重ベクトルであるときには、式（１２）の右辺の第１項Ｋ＋・Ｆが構造体内部の対称的な変位を表していることとなる。従って、式（９）に於いて、逆行列Ｋ−１に代えて擬似逆行列Ｋ＋を用いた場合、即ち、変位ベクトルＵｌが
Ｕｌ＝ＨＴ・Ｋ＋・Ｆ＝ＨＴ・Ｋ＋・Ｈ・Ｆｌ …（１４）
により表された場合、荷重ベクトルＦが剛体運動を生じない荷重ベクトルであるとき、Ｋ＋・Ｆが剛体運動による変位を含まないので、変位ベクトルＵｌは、構造体内部の対称的な変位となっていると考えられる。
By the way, in a simultaneous linear equation such as the equation (2) [KU = F], when the inverse matrix K −1 does not exist in the matrix K having the coefficients of the equation as components, the pseudo inverse matrix K + is used. It is possible to obtain the solution of From the property of the pseudo inverse matrix K + , when Equation (2) [KU = F] has a solution, the solution U is obtained by using an arbitrary vector Uw.
U = K + .F + (I-K + .K) Uw (12)
(I is a unit matrix). Here, assuming that the equation (2) [KU = F], which is a rigidity equation of the structure, represents a static state, and that no acceleration motion occurs in parallel and rotation by the load vector F, the right side of the equation (12) the second term of (I-K + · K) · Uw , when multiplied by the stiffness matrix K, K · from the pseudo inverse matrix K + definition (I-K + · K) · Uw = 0 ... (13)
It can be considered that the displacement vector represents a rigid body motion of the structure that can occur even when no force acts on the structure. On the other hand, the first term K + · F on the right side of the equation (12) is considered to be a displacement occurring in the structure other than the rigid body motion in this case. Is considered to be a load vector that produces a symmetrical displacement at (2) (if there is an unsymmetrical displacement, there should be a rigid body motion because there is no support structure). That is, when the load vector F is a load vector acting so that the total of the load and the total of the moment on the structure become 0, as schematically illustrated in FIG. The first term K + · F on the right side of (12) indicates a symmetric displacement inside the structure. Therefore, in the equation (9), when using a pseudo-inverse matrix K + in place of the inverse matrix K -1, i.e., the displacement vector U l is U l = H T · K + · F = H T · K + · H · F 1 (14)
If it represented by, when the load vector F is a load vector that does not cause the rigid body motion, since K + · F does not include the displacement caused by rigid motion, the displacement vector U l is the symmetrical displacement of the structure section It is thought that it has become.
更に、上記の式（１４）を参照して、式（１０）に代えて、局所コンプライアンス行列Ｃを、擬似逆行列Ｋ＋を用いて
Ｃ＝ＨＴ・Ｋ＋・Ｈ …（１５）
により決定した場合に、或る荷重ベクトルＦｌｉ（ただし、||Flｉ||＝１）が局所コンプライアンス行列Ｃの一つの固有値λｉを与える固有ベクトルφｉであるとき、即ち、
Ｕｌｉ＝Ｃ・Ｆｌｉ＝λｉＦｌｉ（Ｃ・φｉ＝λｉφｉ） …（１６）
であるとき、荷重ベクトルＦｌｉは、構造体内部の対称的な変位ベクトルＵｌｉに平行となるので、固有ベクトルφｉは、構造体に対して荷重の合計とモーメントの合計とがそれぞれ０となるように作用している荷重ベクトルを与えることとなる。実際、図３に模式的に描かれている如き、直方体に於いて、その底面（ｙ＝０）を荷重作用領域として設定して（ｙ＝０の面に不定荷重が作用すると設定したことに相当する。）、式（１５）による局所コンプライアンス行列Ｃの固有ベクトルを算出すると、図４（Ａ）の（ａ）〜（ｆ）に描かれている如く、固有ベクトルの成分により表されるｙ＝０の面内に於ける荷重分布は、荷重の合計とモーメントの合計とがそれぞれ０となる滑らかな分布となることが見出された。即ち、局所コンプライアンス行列Ｃの固有ベクトルを荷重ベクトルとした場合、荷重ベクトルと変位ベクトルとは平行となるので、構造体には対称的な変位が生ずることとなる。
Furthermore, with reference to the above equation (14), wherein in place of (10), the local compliance matrix C, and using a pseudo-inverse matrix K + C = H T · K + · H ... (15)
When a certain load vector F li (|| F li || = 1) is an eigen vector φ i that gives one eigen value λi of the local compliance matrix C,
U li = C · F li = λ i F li (C · φ i = λ i φ i ) (16)
Since the load vector F li is parallel to the symmetric displacement vector U li inside the structure, the eigenvector φ i is such that the sum of the load and the sum of the moment with respect to the structure are each 0. Is given. Actually, as shown schematically in FIG. 3, in the rectangular parallelepiped, the bottom surface (y = 0) is set as a load application area (an indeterminate load is set to act on the y = 0 surface). When the eigenvector of the local compliance matrix C is calculated by Expression (15), y = 0 represented by the eigenvector component is drawn as shown in (a) to (f) of FIG. It has been found that the load distribution in the plane is a smooth distribution in which the total of the load and the total of the moment are each 0. That is, when the eigenvector of the local compliance matrix C is a load vector, the load vector and the displacement vector are parallel to each other, so that a symmetric displacement occurs in the structure.
かくして、本実施形態に於ける支持構造のない構造体の構造最適化演算では、剛性行列Ｋの逆行列Ｋ−１ではなく、擬似逆行列Ｋ＋を用いることにより、局所コンプライアンス行列Ｃは、式（１５）により算出され、局所コンプライアンス行列Ｃの固有ベクトルは、構造体内に於いて対称的な分布の荷重ベクトルＦを与え、変位ベクトルＵが、
Ｕ＝Ｋ＋・Ｆ …（４）
により算出される（Ｕは、構造体内に対称的な変位分布を与えるベクトルとなる。）。
Thus, in the structure optimization operation of the structure having no support structure in the present embodiment, the local compliance matrix C is obtained by using the pseudo inverse matrix K + instead of the inverse matrix K −1 of the stiffness matrix K. The eigenvector of the local compliance matrix C calculated by (15) gives a load vector F having a symmetric distribution in the structure, and the displacement vector U
U = K + · F (4)
(U is a vector that gives a symmetrical displacement distribution within the structure.)
（ｂ）支持構造のない構造体に於けるコンプライアンスの変化
局所コンプライアンス行列Ｃに対しては複数組の固有値と固有ベクトルが得られるところ、支持構造のある構造体の場合、図５（Ａ）に模式的に描かれている如く、構造最適化演算の過程に於いて、最初の段階から、一つの固有ベクトル（φ０）が、徐々に変化しながら、最大のコンプライアンス、即ち、最大固有値（λ０）を与えるようになっている。構造体が支持構造を持つ場合、固有ベクトルにより与えられる荷重分布が対称的となる条件がなく、固有ベクトルに一方向に向かう傾向のあるベクトル（φ０）が存在し、そのような固有ベクトル（φ０）が最大固有値（λ０）を与える荷重ベクトルとなると、支点などの支持構造とその付近に大きな変形を生み、コンプライアンスを著しく増大させ、かくして、構造最適化演算の反復が進んでも、常に最大固有値を与える固有ベクトルであり続けることとなる。即ち、構造最適化演算の過程に於いて、複数存在する固有ベクトルのうち、最大固有値を与える固有ベクトルは通常交代しないので、最大固有値を与える荷重分布及び変位分布は徐々に変化することとなる。従って、従前のロバスト構造最適化演算では、構造体のコンプライアンスを最小化する設計変数の更新に於いて使用する変位ベクトルＵの算出のために、常に、局所コンプライアンス行列Ｃの固有値解析により得られた固有値のうちの最大のものを与える固有ベクトルを荷重ベクトルとして選択していれば、構造体の構造は、安定的に最適な構造へ収束することとなっていた。
(B) Change in compliance in a structure without a support structure A plurality of sets of eigenvalues and eigenvectors are obtained for a local compliance matrix C. In the case of a structure with a support structure, FIG. As depicted in the drawing, in the course of the structure optimization operation, from the first stage, one eigenvector (φ 0 ) gradually changes to the maximum compliance, that is, the maximum eigenvalue (λ 0 ). Is to give. When the structure has a support structure, there is no condition that the load distribution given by the eigenvector is symmetric, and there is a vector (φ 0 ) that tends to one direction in the eigenvector, and such an eigenvector (φ 0 ) Becomes a load vector that gives the maximum eigenvalue (λ 0 ), which causes a large deformation in the support structure such as a fulcrum and its vicinity, greatly increases the compliance, and thus the maximum eigenvalue is always increased even if the repetition of the structure optimization operation proceeds. It will continue to be the given eigenvector. That is, in the process of the structure optimization operation, among the plurality of eigenvectors, the eigenvector that gives the maximum eigenvalue does not normally change, so that the load distribution and the displacement distribution that give the maximum eigenvalue gradually change. Therefore, in the conventional robust structure optimization operation, the eigenvalue analysis of the local compliance matrix C was always used for the calculation of the displacement vector U used in updating the design variables for minimizing the compliance of the structure. If the eigenvector giving the largest one of the eigenvalues is selected as the load vector, the structure of the structure will stably converge to the optimal structure.
一方、支持構造のない構造体の場合、（擬似逆行列Ｋ＋を用いて得られた）局所コンプライアンス行列Ｃの固有ベクトルに於いては、図２或いは図４の例示からも理解される如く、複数得られる固有ベクトルは、それぞれの態様にて、構造体に於いて荷重の合計及びモーメントの合計が共に０となる対称的な荷重分布を保存しながら、構造体の構造最適化演算の過程に於いて徐々に変化することとなり、いずれも、一方向に向かうベクトルになる傾向がなく、変形が局所的に一方向に生ずることがない（図５（Ｂ）に示されているように、固有値λ０を与える一方向に向かう固有ベクトルφ０が存在しない。）。こういった支持構造のない構造体に於ける性質により、構造最適化演算の過程では、図５（Ｂ）に描かれているように、最大固有値を与える固有ベクトルが分布の態様の異なるベクトルへ代わってしまう現象、即ち、最大固有値を与える固有ベクトルの交代（最大固有値がλ１からλ２へ代わるため、最大固有値を与える固有ベクトルがφ１からφ２に交代する。）が生じ得ることが見出された。かかる最大固有値を与える固有ベクトルの交代が生ずると、最大固有値を与える固有ベクトルを用いて算出される変位ベクトルの分布が不連続的に変化し、これにより、構造体のコンプライアンスを最小化するべく更新される設計変数も不連続的に変化することとなり、構造体の構造が安定的に最適な構造へ収束し難くなる場合が生じた。 On the other hand, in the case of a structure having no support structure, the eigenvector of the local compliance matrix C (obtained using the pseudo inverse matrix K + ) has a plurality of eigenvectors as can be understood from the example of FIG. 2 or FIG. In each aspect, the obtained eigenvectors preserve the symmetric load distribution in which the sum of the loads and the sum of the moments in the structure are both 0, and in the course of the structure optimization calculation of the structure. In each case, there is no tendency to become a vector in one direction, and no deformation occurs locally in one direction (eigenvalue λ 0 as shown in FIG. 5B). There is no eigenvector φ 0 going in one direction that gives Due to such a property of the structure having no supporting structure, in the process of the structure optimization operation, as shown in FIG. 5B, the eigenvector giving the maximum eigenvalue is replaced by a vector having a different distribution mode. and thus a phenomenon, i.e., (since the maximum eigenvalue is replaced to lambda 2 from lambda 1, eigenvector gives the maximum eigenvalue alternates from phi 1 to phi 2.) Substitution of eigenvectors giving the maximum eigenvalue been found that may occur Was. When such an alternation of the eigenvector giving the maximum eigenvalue occurs, the distribution of the displacement vector calculated using the eigenvector giving the maximum eigenvalue changes discontinuously, thereby being updated to minimize the compliance of the structure. The design variables also change discontinuously, which sometimes makes it difficult for the structure of the structure to stably converge to an optimal structure.
そこで、本実施形態では、局所コンプライアンス行列Ｃから得られる複数の固有ベクトルにより与えられる複数の固有値が全体的に最小化されるように設計変数の更新が実行されてよい。具体的には、局所コンプライアンス行列Ｃの固有値解析を実行した後、まず、固有値の大きい方から所定の数Ｎの固有値λｉと固有ベクトルφｉの組が選択され、
λｔ＝ｗ１λ１＋ｗ２λ２＋…＋ｗＮλＮ …（１７）
により与えられる固有値の線形和（重み付き和）λｔが評価関数として最小化されるように設計変数の更新が実行されてよい。ここにおいて、ｗｉは、線形和λｔに於ける各項λｉについての重みである。かかる重みｗｉは、具体的には、例えば、下記のように設定されてよい。
ｗｉ＝（λｉ−λＮ＋１）μ／Σ（λｉ−λＮ＋１）μ …（１８）
（Σは、ｉ＝１〜Ｎまでの総和であり、μはべき指数である。）
式（１８）の重みｗｉによれば、固有値λｉが大きいほど、線形和λｔに於ける寄与が大きくなり、最小化に於ける寄与が増大することとなるので、構造体の構造がより速やかに最適な構造へ収束することが期待される。かくして、設計変数に対する最小化される評価関数の感度は、各固有値の感度∂λｉ／∂ρｅを用いて（ρｅは要素ｅに割り当てられた設計変数）
ｓｅ＝−Σｗｉ（∂λｉ／∂ρｅ） …（１９）
（Σは、ｉ＝１〜Ｎまでの総和）
により与えられる。
Thus, in the present embodiment, design variables may be updated such that a plurality of eigenvalues given by a plurality of eigenvectors obtained from the local compliance matrix C are entirely minimized. Specifically, after performing the eigenvalue analysis of the local compliance matrix C, first, a set of a predetermined number N of eigenvalues λ i and eigenvectors φ i is selected from the largest eigenvalue,
λ t = w 1 λ 1 + w 2 λ 2 + ... + w N λ N (17)
Updating the design variables may be performed as is minimized as a linear sum (weighted sum) lambda t evaluation function given eigenvalues by. Here, w i is the weight for in terms lambda i in linear sum lambda t. Specifically, the weight wi may be set as follows, for example.
w i = (λ i −λ N + 1 ) μ / Σ (λ i −λ N + 1 ) μ (18)
(Σ is the sum of i = 1 to N, and μ is a power exponent.)
According to the weight w i in the equation (18), the larger the eigenvalue λ i , the larger the contribution in the linear sum λ t and the larger the contribution in the minimization. It is expected to converge to the optimal structure more quickly. Thus, the sensitivity of the evaluation function to be minimized to the design variable is calculated using the sensitivity of each eigenvalue ∂λ i / ∂ρ e (where ρ e is the design variable assigned to the element e).
s e = -Σw i (∂λ i / ∂ρ e) ... (19)
(Σ is the sum of i = 1 to N)
Given by
構造最適化演算手法の概要（密度法を用いたトポロジー最適化演算の概要）
本発明の装置の実施形態に於ける具体的な構造最適化演算は、上記の如く、擬似逆行列を利用し、評価関数として複数の固有値の線形和を採用しながら、任意のトポロジー構造最適化や形状最適化などの任意のアルゴリズムに従って実行されてよいことは理解されるべきである。本実施形態に於いては、例えば、以下に概説される密度法（Density-based Approach）を用いたトポロジー最適化演算が採用されてよい。
Outline of structure optimization calculation method (Overview of topology optimization calculation using density method)
As described above, a specific structure optimization operation in the embodiment of the device of the present invention uses a pseudo-inverse matrix and employs a linear sum of a plurality of eigenvalues as an evaluation function, and performs arbitrary topology structure optimization. It should be understood that the algorithm may be performed according to any algorithm, such as shape and shape optimization. In the present embodiment, for example, a topology optimization calculation using a density method (Density-based Approach) outlined below may be employed.
密度法を用いたトポロジー最適化演算に於いては、具体的には、まず、三次元空間に於いて任意の寸法の設計領域（例えば、図３の如き直方体の領域など）が設定され、設計領域を格子状に複数の要素に分割し、各要素に対する設計変数として各要素内材料の密度（材料密度）ρｅ（ｅは、要素番号）が
０≦ρｅ≦１ …（２０）
の条件を満たすよう設定され、材料密度ρｎのときの要素の材料のヤング率Ｅｎが、
Ｅｅ（ρｅ）＝ρｅ ＰＥｏ …（２１）
と設定される。ここで、Ｐは、ペナルティパラメータと称する値（＞１）であり、Ｅｏは、材料密度ρｅ＝１のときの材料のヤング率である。なお、後に説明される変位ベクトルＵの演算に於ける特異性を回避するために、ヤング率Ｅｎは、
Ｅｅ（ρｅ）＝Ｅmin＋ρｅ Ｐ(Ｅｏ−Ｅmin) …（２２）
と修正されて用いられてよい（非特許文献１、２）。ここで、Ｅminは、材料密度ρｅ＝０のときの仮想的な材料のヤング率である。そして、評価関数として構造体のコンプライアンスが採用され、有意な材料密度ρｅを有する要素から構成される構造体のコンプライアンスが最小化されるように要素の材料密度ρｅを反復して更新する演算処理が実行され、設計領域内に於ける有意な材料密度を有する要素の分布の形状が、目的の構造体の形状、即ち、剛性ができるだけ高くなる構造体の形状、として見出される。本実施形態に於いては、最適化の評価関数である構造体のコンプライアンスは、構造体の局所コンプライアンス行列Ｃの固有値又は固有値の線形和（式（１７））となる。
In the topology optimization calculation using the density method, first, a design area of an arbitrary size (for example, a rectangular parallelepiped area as shown in FIG. 3) is set in a three-dimensional space. The region is divided into a plurality of elements in a grid pattern, and the density (material density) ρ e (e is an element number) of the material in each element as a design variable for each element is 0 ≦ ρ e ≦ 1 (20)
The condition is satisfied setting, the Young's modulus E n of the material of the element at the time of material density [rho n is,
E e (ρ e ) = ρ e P Eo (21)
Is set. Here, P is a value (> 1) called a penalty parameter, and Eo is the Young's modulus of the material when the material density ρ e = 1. In order to avoid in specificity to the calculation of the displacement vector U that is described later, the Young's modulus E n is
E e (ρ e) = Emin + ρ e P (Eo-Emin) ... (22)
(Non-Patent Documents 1 and 2). Here, Emin is the virtual Young's modulus of the material when the material density ρ e = 0. Then, operation compliance structure as an evaluation function is employed to update iteratively the material density [rho e elements as compliance structures comprised elements with significant material density [rho e is minimized The process is performed and the shape of the distribution of elements having significant material density within the design domain is found as the shape of the structure of interest, that is, the structure with the highest possible stiffness. In the present embodiment, the compliance of the structure, which is the evaluation function of the optimization, is the eigenvalue of the local compliance matrix C of the structure or a linear sum of the eigenvalues (Equation (17)).
なお、材料密度ρｎの更新は、設計領域の体積に対する構造体の体積の割合Ｖflacが所与の体積制約割合の条件を満たすように実行される。設計領域内の構造体の体積は、要素ｅの体積ｖｅを用いて、Σρｅ・ｖｅ（ここで、Σは、全要素についての総和である。）により与えられるので、材料密度ρｅは、体積制約割合Ｖflacの条件として、
Σρｅ・ｖｅ≦Ｖflac・Ｖｔ …（２３）
（ここで、Ｖｔは、設計領域の体積）
が満たされるように更新されることとなる。また、演算結果に於いて、各要素の材料密度の値ができるだけ０又は１に収束するように材料密度ρｅには、任意のフィルター処理を施した値が用いられてよい。
The update of the material density ρ n is performed such that the ratio Vflac of the volume of the structure to the volume of the design region satisfies the condition of the given volume constraint ratio. The volume of the structure design area, using the volume v e elements e, Σρ e · v e (where, sigma is. The sum of all elements) so provided by, the material density [rho e Is the condition of the volume constraint ratio Vflac,
Σρ e · v e ≦ Vflac · Vt ... (23)
(Where Vt is the volume of the design area)
Will be updated so that is satisfied. Further, in the calculation result, the material density [rho e to converge only 0 or 1 can the value of the material density of each element may be used a value subjected to any filtering.
トポロジー最適化演算の具体的な処理手順に於いては、まず、初期状態の設定として、設計領域内にて生成される構造体に対する拘束条件、構造体に対して作用する荷重条件、体積制約割合（設計領域の体積に対する生成されるべき構造体の体積の割合）、材料密度ρｅの初期値、構造体を構成する材料の特性値（ヤング率、ポアソン比）を含む構造最適化演算のための初期条件が実施者により設定される。そして、設定された設計領域に対して、任意のプログラムを用いて要素の設定が為される。材料密度ρｅの初期値は、設計領域内全域に於いて均一な値であってよい。また、本実施形態に於いては、最適化の対象は、支持構造のない構造体であるので、構造体に於ける自由な剛体運動を制限する拘束条件は設定されない。更に、本実施形態は、ロバスト構造最適化の手法に従って荷重条件を設定するので、構造体に作用する荷重条件については、荷重作用領域は設定されるが、かかる領域内での荷重分布は局所コンプライアンス行列Ｃの固有値解析により得られた固有ベクトルによって決定されるので、荷重の具体的な値及びその分布は初期条件としては設定されない。 In the specific processing procedure of the topology optimization calculation, first, as the initial state setting, the constraint condition for the structure generated in the design area, the load condition acting on the structure, the volume constraint ratio (ratio of the volume of the structure to be produced to the volume of the design domain), an initial value of the material density [rho e, characteristic values of the material constituting the structure (Young's modulus, Poisson's ratio) for structural optimization operation comprising Are set by the practitioner. Then, elements are set in the set design area using an arbitrary program. The initial value of the material density [rho e may be a uniform value at the design area throughout. Further, in the present embodiment, since the optimization target is a structure having no support structure, a constraint condition for restricting a free rigid motion in the structure is not set. Further, in the present embodiment, since the load condition is set according to the robust structure optimization method, the load application region is set for the load condition acting on the structure, but the load distribution in this region is the local compliance. Since the load is determined by the eigenvector obtained by the eigenvalue analysis of the matrix C, the specific value of the load and its distribution are not set as initial conditions.
次いで、下記の（ａ）、（ｂ）、（ｃ）の演算処理が所定の終了条件が成立するまで反復して実行される。 Next, the following calculation processes (a), (b), and (c) are repeatedly executed until a predetermined termination condition is satisfied.
（ａ）荷重条件の設定−構造体に最悪荷重が作用している条件を見出すために、まず、設計領域内の各要素の材料密度ρｅから各要素のヤング率Ｅｅ（ρｅ）を計算し、かかるヤング率を用いて設計領域内に形成されている構造体の剛性行列Ｋが決定される。剛性行列Ｋに於いて使用される各要素のヤング率Ｅｅには、式（２１）又は（２２）が用いられる。そして、剛性行列Ｋを用いて、擬似逆行列Ｋ＋、局所コンプライアンス行列Ｃが算出され、局所コンプライアンス行列Ｃの固有値解析を実行し、複数組の固有値λｉ及び固有ベクトルφｉのうち、固有値の大きい方から所定の数Ｎの組が選択され、選択された固有ベクトルφｉがそれぞれ荷重ベクトルＦｉに設定される。なお、擬似逆行列Ｋ＋及び局所コンプライアンス行列Ｃは、任意のプラグラムに従って計算されてよい（プログラム言語に用意されている関数を利用してもよく、任意のアルゴリズムによって算出されてもよい。）。また、局所コンプライアンス行列Ｃの固有値解析も任意のプラグラムに従って計算されてよい（プログラム言語に用意されている関数を利用してもよく、任意のアルゴリズムによって算出されてもよい。） (A) Setting of load condition—In order to find out the condition in which the worst load is acting on the structure, first, the Young's modulus E e (ρ e ) of each element is determined from the material density ρ e of each element in the design area. The calculated rigidity matrix K of the structure formed in the design area is determined using the Young's modulus. The Young's modulus E e of each element to be employed in the stiffness matrix K, the formula (21) or (22) is used. Then, a pseudo inverse matrix K + and a local compliance matrix C are calculated using the stiffness matrix K, and an eigenvalue analysis of the local compliance matrix C is performed, and among a plurality of sets of eigenvalues λ i and eigenvectors φ i , the eigenvalues having large eigenvalues are calculated. A set of a predetermined number N is selected from the two, and the selected eigenvectors φ i are respectively set as the load vectors F i . The pseudo inverse matrix K + and the local compliance matrix C may be calculated according to an arbitrary program (a function prepared in a programming language may be used, or may be calculated by an arbitrary algorithm). Also, the eigenvalue analysis of the local compliance matrix C may be calculated according to an arbitrary program (a function prepared in a programming language may be used, or may be calculated by an arbitrary algorithm).
（ｂ）変位ベクトルの算出−上記の荷重条件の下、設計領域内に於ける各要素の各節点の変位が算出される。具体的な算出方法は、有限要素法等に従った任意のプログラムを用いたものであってよい。既に述べた如く、本実施形態に於いては、最適化の対象となる構造体に支持構造がなく、構造体の全要素の剛性行列Ｋの逆行列Ｋ−１が存在しないので、各要素の各節点の変位は、剛性行列Ｋの擬似逆行列Ｋ＋を用いて、式（４）によって算出される。また、最適化の評価関数として、局所コンプライアンス行列Ｃの複数の選択された固有値λｉの線形和が用いられる場合には、各節点の変位は、複数の選択された固有ベクトルの各々を荷重ベクトルとして用いた場合のそれぞれについて算出される（即ち、評価関数がＮ個の固有値の線形和（式（１７））であるときには、Ｎ個の対応する固有ベクトルを荷重ベクトルとした場合に対応して、Ｎ個の変位ベクトルが算出されることとなる。）。 (B) Calculation of displacement vector-The displacement of each node of each element in the design area is calculated under the above load conditions. A specific calculation method may use an arbitrary program according to the finite element method or the like. As described above, in the present embodiment, there is no support structure in the structure to be optimized and there is no inverse matrix K- 1 of the rigidity matrix K of all elements of the structure. The displacement of each node is calculated by Expression (4) using a pseudo inverse matrix K + of the rigidity matrix K. Further, when a linear sum of a plurality of selected eigenvalues λ i of the local compliance matrix C is used as the optimization evaluation function, displacement of each node is obtained by using each of the plurality of selected eigenvectors as a load vector. (Ie, when the evaluation function is a linear sum of N eigenvalues (Equation (17)), corresponding to the case where N corresponding eigenvectors are used as load vectors, Will be calculated.)
（ｃ）材料密度の更新−上記の演算により得られた設計領域内に全要素の節点変位ベクトルＵと全要素の剛性行列Ｋとを用いて、式（１９）により構造体の材料密度ρｅに対する評価関数の感度（材料密度ρｅについての評価関数の微分）が算出され、更に、材料密度ρｅに対する構造体の評価関数の感度を用いて、体積制約割合条件（式（２３））が満たされることを拘束条件として構造体の評価関数が最小化されるように新たな材料密度ρｅを算出する演算処理、即ち、材料密度ρｎの更新処理が実行される。材料密度ρｅの更新処理は、例えば、最適性基準法、逐次線形計画法、逐次凸関数近似法などの手法に従って実行されてよい。具体的には、材料密度ρｅについての各固有値の感度∂λｉ／∂ρｅは、式（２２）を参照して、
∂λｉ／∂ρｅ＝−Ｕｉ Ｔ[pρｅ p-1(Ｅｏ−Ｅmin)Ｋ０ ｅ]Ｕｉ …（２４）
により与えられる。ここで、Ｕｉは、固有値λｉを与える固有ベクトルφｉを荷重ベクトルとした場合の変位ベクトルであり、Ｋ０ ｅは、要素ｅの単位剛性行列（ヤング率＝１と設定した剛性行列）である。従って、固有値λ１〜λＮを選択したときには、式（２４）、式（１８）及び式（１９）を用いて、各要素ｅの設計変数に対する評価関数の感度ｓｅが算出される。そして、例えば、新しい材料密度ρｅ (k+1)は、それまでの材料密度ρｅ (k)を用いて（kは、反復処理の回数）、
ρｅ (k+1)＝Ｂｅ (k)η・ρｅ (k) …（２５ａ）
或いは、
ρｅ (k+1)＝max(ρｅ (k)-m, min(ρｅ (k)+m, Ｂｅ (k) η・ρｅ (k)))…（２５ｂ）
の形式にて与えられる。ここで、Ｂｎ (k)は、感度ｓｅを用いて、
Ｂｅ (k)＝Ｈｂ(ｓｅ)[Λ（∂ｖ（Ρ(k)）／∂ρｅ）]−１ …（２６）
により与えられる。ここで、Ｈｂ(ｓｅ)は、ρｅを滑らかにするフィルター関数であり、∂ｖ（Ρ(k)）／∂ρｅは、式（２３）の体積制約条件
ｖ（Ρ(k)）＝Σρｅ (k)・ｖｅ−Ｖflac・Ｖｔ≦０
のρｅについての微分であり（∂ｖ（Ρ(k)）／∂ρｅ＝ｖｅ）、Λは、ラグランジュの未定係数であり、二分法などによって体積制約条件（式（２３））を満たす値が決定されてよい。ｍは、1回の更新に於ける材料密度の変化限度であり、ηは、Ｂｅ (k)のべき指数であり、＜１である。なお、Ｂｅ (k)ηの演算負荷を低減するために、
Ｂｅ (k)η〜１＋η（Ｂｅ (k)−１） …（２７）
による近似演算が用いられてよい。具体的な演算処理のアルゴリズムは、例えば、非特許文献１、２などに記載されているものが使用されてよい。
(C) Updating the material density—using the nodal displacement vectors U of all elements and the stiffness matrix K of all elements in the design area obtained by the above calculation, the material density ρ e of the structural body is obtained by equation (19). the sensitivity of the evaluation function (differential of the evaluation function for the material density [rho e) is calculated for further using the sensitivity of the evaluation function of the structure with respect to the material density [rho e, volumetric constraint ratio condition (equation (23)) is An arithmetic process of calculating a new material density ρ e , that is, an update process of the material density ρ n is executed so that the evaluation function of the structure is minimized with the satisfaction of the constraint as a constraint condition. The update processing of the material density ρ e may be executed according to, for example, a method such as an optimality criterion method, a sequential linear programming method, and a sequential convex function approximation method. Specifically, the sensitivity ∂λ i / ∂ρ e of each eigenvalue of material density [rho e, with reference to equation (22),
∂λ i / ∂ρ e = −U i T [pρ ep -1 (Eo−Emin) K 0 e ] U i (24)
Given by Here, U i is a displacement vector when the eigen vector φ i giving the eigen value λ i is a load vector, and K 0 e is a unit stiffness matrix of element e (a stiffness matrix set to Young's modulus = 1). is there. Therefore, when the eigenvalues λ 1 to λ N are selected, the sensitivity se of the evaluation function to the design variable of each element e is calculated using Expressions (24), (18), and (19). Then, for example, the new material density ρ e (k + 1) is calculated using the material density ρ e (k) up to that point (k is the number of iterations),
ρ e (k + 1) = B e (k) η · ρ e (k) (25a)
Or,
ρ e (k + 1) = max (ρ e (k) -m, min (ρ e (k) + m, Be (k) η · ρ e (k) )) ... (25b)
Is given in the form Here, B n (k), using the sensitivity s e,
B e (k) = Hb ( s e) [Λ (∂v (Ρ (k)) / ∂ρ e)] -1 ... (26)
Given by Here, Hb (s e ) is a filter function for smoothing ρ e , and ∂v (Ρ (k) ) / ∂ρ e is a volume constraint condition v (Ρ (k) ) in Expression (23). = Σρ e (k) · v e -Vflac · Vt ≦ 0
It is a differential of about ρ e a (∂v (Ρ (k)) / ∂ρ e = v e), Λ is the undetermined coefficients of Lagrange, volume constraints, such as by dichotomy (equation (23)) The value to be satisfied may be determined. m is the change limit once in material density to update, eta is B e (k) is an index powers of a <1. In order to reduce the calculation load of B e (k) η,
B e (k) η １1 + η (B e (k) −1) (27)
May be used. As a specific algorithm of the arithmetic processing, for example, those described in Non-Patent Documents 1 and 2 and the like may be used.
上記の（ａ）〜（ｃ）の演算処理は、（ｃ）の処理にて新たな材料密度ρｎ (k+1)が得られる度に、それを材料密度ρｎ (k)として用いて、繰り返し実行される。かかる反復処理を実行されるにつれて、材料密度ρｎの分布は、体積制約割合条件（２３）を満たしつつ、構造体の評価関数を最小化する方向、即ち、構造体の剛性が高くなる方向に変化していくことが期待される。そして、かかる処理が反復されていく過程で、好ましくは、各要素の材料密度ρｎの値は、０又は１に収束され、これにより、設計領域内に於いて有意な材料密度（ρｎ≒１）を有する要素の分布の形状が体積制約割合を満たす条件下で剛性が高くなった構造体の形状として見出されることとなる。なお、（ａ）〜（ｃ）の反復処理を終了するための所定の終了条件としては、（ｉ）反復処理回数が任意に設定されてよい回数に到達した場合、（ii）構造体の評価関数の値が収束した場合（処理サイクル毎の評価関数の変化量が所定の小数以下となった場合）、（iii）更新される材料密度ρｎの値が収束した場合（処理サイクル毎の材料密度の変化量の最大値が所定の小数以下となった場合）などが選択されてよい。（ii）又は（iii）を終了条件として採用する場合には、結果として得られる構造体の形状は、体積制約割合条件の下で、最適化されていることが期待される。 In the above-described arithmetic processing (a) to (c), each time a new material density ρ n (k + 1) is obtained in the processing of (c), the new material density ρ n (k + 1) is used as the material density ρ n (k). Is executed repeatedly. As such an iterative process is performed, the distribution of the material density ρ n is directed in a direction in which the evaluation function of the structure is minimized, that is, in a direction in which the rigidity of the structure increases, while satisfying the volume constraint ratio condition (23). It is expected to change. Then, in the process of repeating such processing, preferably, the value of the material density ρ n of each element is converged to 0 or 1, whereby the significant material density (ρ nに) in the design region is obtained. Under the condition that the shape of the distribution of the element having 1) satisfies the volume restriction ratio, it is found as the shape of the structure having increased rigidity. In addition, as predetermined termination conditions for ending the repetition processing of (a) to (c), (i) when the number of repetition processing reaches a number that can be arbitrarily set, (ii) evaluation of the structure When the value of the function converges (when the amount of change of the evaluation function in each processing cycle becomes a predetermined decimal or less), (iii) When the value of the updated material density ρ n converges (the material in each processing cycle) (A case where the maximum value of the change amount of the density is equal to or less than a predetermined decimal number) may be selected. When (ii) or (iii) is adopted as the termination condition, it is expected that the shape of the resulting structure is optimized under the volume constraint ratio condition.
装置の作動
図６を参照して、本実施形態の装置に於ける構造最適化演算の処理手順は、次のように実行されてよい。ロバスト最適化演算を実行するために必要な初期条件の設定が実行される（ステップ１）。既に触れた如く、これらの初期条件の設定は、演算の実施者により入力装置４を通じて実行される。しかる後、コンピュータの作動により、まず、剛性行列Ｋの演算（ステップ２）、擬似逆行列Ｋ＋の演算（ステップ３）、擬似逆行列Ｋ＋を用いた局所コンプライアンス行列Ｃの演算（ステップ４：式（１５））が順に実行され、局所コンプライアンス行列Ｃの固有値解析が実行されて、局所コンプライアンス行列Ｃの所定数の組の固有値λｉと固有ベクトルφｉが選択される（ステップ５）。局所コンプライアンス行列Ｃの固有値と固有ベクトルが選択されると、選択された固有ベクトルφｉをそれぞれ荷重ベクトルとして、疑似逆行列Ｋ＋を用いて、有限要素法による設計領域内の各要素の各節点の変位（変位ベクトルＵｉ）が演算される一方（ステップ６：式（４））、選択された固有値λｉを用いて、評価関数感度に於ける各固有値の重みｗｉが演算される（ステップ７：式（１８））。そして、固有ベクトル毎の変位ベクトルＵｉと剛性行列Ｋの成分を用いて式（２４）により、各固有値λｉの材料変数に対する感度（∂λｉ／∂ρｅ）が算出され、それらの感度と重みｗｉを用いて、各要素ｅの評価関数（固有値の線形和）の感度ｓｅが算出される（ステップ８：式（１９）、固有値感度算出手段）。そうすると、構造体の評価関数の感度ｓｅを用いて、上記の態様にて、トポロジー最適化の手法に従った材料密度（密度関数）の更新演算が実行され、設計領域内の新たな材料密度の分布が得られることとなる（ステップ９：式（２５ａ）又は式（２５ｂ））。かくして、一連の演算処理が完了すると、設計領域内の材料密度の分布、有意な材料密度を有する要素の分布、又は、構造体の形状がコンピュータのモニター３上に表示されてよい（ステップ１０）。
Operation of Apparatus With reference to FIG. 6, the processing procedure of the structure optimization operation in the apparatus of the present embodiment may be executed as follows. Initial conditions necessary for executing the robust optimization operation are set (step 1). As mentioned above, the setting of these initial conditions is executed by the operator of the calculation through the input device 4. Thereafter, the computer operates to calculate the stiffness matrix K (step 2), calculate the pseudo inverse matrix K + (step 3), and calculate the local compliance matrix C using the pseudo inverse matrix K + (step 4: Equation (15)) is sequentially executed, eigenvalue analysis of the local compliance matrix C is performed, and a predetermined number of sets of eigenvalues λ i and eigenvectors φ i of the local compliance matrix C are selected (step 5). When the eigenvalues and eigenvectors of the local compliance matrix C are selected, the displacement of each node of each element in the design area by the finite element method using the pseudo-inverse matrix K + with the selected eigenvectors φ i as load vectors, respectively. While (displacement vector U i ) is calculated (step 6: equation (4)), weight wi of each eigenvalue in the evaluation function sensitivity is calculated using the selected eigenvalue λ i (step 7: Equation (18). Then, the sensitivity (∂λ i / ∂ρ e ) of each eigenvalue λ i to the material variable is calculated by Expression (24) using the displacement vector U i for each eigen vector and the components of the stiffness matrix K. using the weight wi, sensitivity s e of the evaluation function of each element e (linear sum of the eigenvalues) is calculated (step 8: equation (19), the eigenvalue sensitivity calculating means). Then, by using the sensitivity s e of the evaluation function of the structure, in the above embodiment, the update operation of the material density in accordance with the method of the topology optimization (density function) is executed, a new material density within the design area (Step 9: Equation (25a) or Equation (25b)). Thus, when a series of arithmetic processing is completed, the distribution of the material density in the design area, the distribution of the elements having a significant material density, or the shape of the structure may be displayed on the monitor 3 of the computer (Step 10). .
上記の一連の処理が実行されると、既に述べた如き所定の終了条件が成立しているか否かが判定され（ステップ１１）、所定の終了条件が成立していなければ、ステップ２〜ステップ１０の実行が繰り返され、所定の終了条件が成立すると、演算処理が終了する。 When the above series of processing is executed, it is determined whether or not the predetermined end condition as described above is satisfied (step 11). If the predetermined end condition is not satisfied, steps 2 to 10 are performed. Is repeated, and when a predetermined end condition is satisfied, the arithmetic processing ends.
計算実験例
上記に説明した本発明の装置の有効性を検証するために、図３に例示する直方体を設計領域１０とし、荷重作用領域として、直方体の底面（ｙ＝０の面）に上下方向に不定荷重が作用する領域とし、上記の処理手順に従って、計算実験を実行し、構造体が最適化されることを確認した。なお、設計領域の寸法は、ｘ×ｙ×ｚ＝１４×７×１８と設定し、一つの要素の体積は、ｖｅ＝１と設定し、体積制約割合Ｖflacは、０．３と設定した。また、構造体を構成する材料のヤング率は、Ｅｏ＝１、Ｅmin＝１０−９とし、ポアソン比は、０．３とし、ペナルティパラメータは、Ｐ＝３とした。また、式（１７）の評価関数λｔに於いて選択する固有値と固有ベクトルの組数は、Ｎ＝４とし、式（１８）に於いて、μ＝１とした。更に、式（２５ａ）、（２５ｂ）に於いて、η＝０．１５、ｍ＝０．０５と設定した。
Calculation Experimental Example In order to verify the effectiveness of the above-described apparatus of the present invention, the rectangular parallelepiped illustrated in FIG. , A calculation experiment was performed according to the above-described processing procedure, and it was confirmed that the structure was optimized. The dimensions of the design area were set as x × y × z = 14 × 7 × 18, the volume of one element was set as v e = 1, and the volume constraint ratio Vflac was set as 0.3. . The Young's modulus of the material constituting the structure was Eo = 1, Emin = 10 −9 , the Poisson's ratio was 0.3, and the penalty parameter was P = 3. Further, the number of sets of eigenvalues and eigenvectors selected in the evaluation function λt in Expression (17) is set to N = 4, and in Expression (18), μ = 1. Further, in equations (25a) and (25b), η = 0.15 and m = 0.05.
まず、図４（Ａ）は、既に説明されている如く、構造最適化演算開始前の図３に示す直方体に於ける局所コンプライアンス行列Ｃの固有値解析によって得られた固有ベクトルの成分のｙ＝０の面内に於ける分布を示している（初期条件に於いて、ｙ＝０の面にのみ荷重が生ずるように設定しているので、ｙ＝０以外の領域の荷重は全て０である。）。図４（Ａ）の（ａ）〜（ｆ）の分布は、値が大きい方から順に六つの固有値を与える固有ベクトルの分布である。上記の如く、構造最適化演算開始前の段階で、支持構造のない構造体に於いて算出される局所コンプライアンス行列Ｃの固有ベクトルは、構造体に於いて一方向の運動を発生させない対称的な荷重分布、即ち、荷重の合計とモーメントの合計との両方が０となる荷重分布を与えることが示された。なお、構造体の荷重作用領域に於ける変位ベクトルは、局所コンプライアンス行列Ｃの固有ベクトルと平行であるので、構造体内に生ずる変位も対称的になることが理解される。 First, as described above, FIG. 4A shows the case where y = 0 of the component of the eigenvector obtained by the eigenvalue analysis of the local compliance matrix C in the rectangular parallelepiped shown in FIG. 3 before the start of the structure optimization operation. The distribution in the plane is shown (in the initial condition, the load is set to occur only on the plane of y = 0, so the load in the area other than y = 0 is all zero). . The distributions (a) to (f) in FIG. 4A are distributions of eigenvectors that give six eigenvalues in descending order of value. As described above, at the stage before the start of the structural optimization operation, the eigenvector of the local compliance matrix C calculated for the structure having no support structure is a symmetrical load that does not generate a unidirectional motion in the structure. It was shown that the distribution, that is, a load distribution in which both the total load and the total moment were zero. Since the displacement vector in the load application region of the structure is parallel to the eigenvector of the local compliance matrix C, it is understood that the displacement generated in the structure is also symmetric.
一方、図４（Ｂ）は、図３に示す直方体に於いて本実施形態による構造最適化演算に於いて設計変数の更新を１００回実行した後の構造体に於ける局所コンプライアンス行列Ｃの固有値解析によって得られた、値が大きい方から六つの固有値を与える固有ベクトルを示している。図４（Ｂ）に於ける（ａ）〜（ｆ）は、図４（Ａ）の（ａ）〜（ｆ）の固有ベクトルに対応する固有ベクトルの分布である。図４（Ｂ）に示した固有ベクトルの分布に於いても局所コンプライアンス行列Ｃの固有ベクトルは荷重の合計とモーメントの合計との両方が０となる対称的な分布となった。また、図４（Ｂ）と図４（Ａ）とを比較して、構造最適化後の固有ベクトルは、構造最適化前の固有ベクトルに比して振幅又は起伏が緩やかであることが把握される。これは、構造最適化により、構造体の剛性が高まり、コンプライアンスが低減された、即ち、変形の程度が小さくなったことを示している。 On the other hand, FIG. 4B shows the eigenvalues of the local compliance matrix C in the rectangular parallelepiped shown in FIG. 3 after the design variables have been updated 100 times in the structural optimization operation according to the present embodiment. It shows the eigenvectors obtained by analysis and giving the six eigenvalues from the larger value. 4A to 4F show the distribution of the eigenvectors corresponding to the eigenvectors of FIGS. 4A to 4F. Also in the distribution of the eigenvectors shown in FIG. 4B, the eigenvectors of the local compliance matrix C have a symmetric distribution in which both the total of the load and the total of the moment are zero. Further, by comparing FIG. 4B and FIG. 4A, it is understood that the eigenvector after the structure optimization has a smaller amplitude or undulation than the eigenvector before the structure optimization. This indicates that the structural optimization has increased the rigidity of the structure and reduced compliance, that is, reduced the degree of deformation.
図７は、図３に示す直方体から上記の手順に従って設計変数の更新を１００回実行した後に得られた構造体の形状を示している。図に於いては、材料密度が０．５以上になった要素が表示されている。色の濃さは、材料密度の高さに対応する。図示の如く、不定荷重の条件下にて支持構造のない構造体に於いて剛性がより高く且つ体積制約条件を満たす形状が得られ、本発明の装置により支持構造のない構造体の構造最適化が可能であることが示された。 FIG. 7 shows the shape of the structure obtained after the design variables are updated 100 times from the rectangular parallelepiped shown in FIG. 3 in accordance with the above procedure. In the figure, elements whose material density is 0.5 or more are displayed. The color intensity corresponds to the height of the material density. As shown in the figure, a structure having a higher rigidity and satisfying the volume constraint condition is obtained in a structure without a support structure under an indefinite load condition, and the apparatus of the present invention optimizes the structure of the structure without a support structure. Has been shown to be possible.
以上の説明は、本発明の実施の形態に関連してなされているが、当業者にとつて多くの修正及び変更が容易に可能であり、本発明は、上記に例示された実施形態のみに限定されるものではなく、本発明の概念から逸脱することなく種々の装置に適用されることは明らかであろう。 Although the above description has been made in connection with the embodiments of the present invention, many modifications and changes are easily possible for those skilled in the art, and the present invention is limited to only the above-exemplified embodiments. It will be apparent that the invention is not limited and can be applied to various devices without departing from the concept of the invention.
Claims (1)
前記設計領域の寸法、前記構造体に於いて荷重の作用する荷重作用領域の範囲、前記構造体を構成する材料の弾性特性を表す物性定数及び前記設計領域の体積に対する前記構造体の体積の割合である体積制約割合を含む条件を設定する条件設定手段と、
前記設計変数を用いて前記構造体の剛性行列を算出する剛性行列算出手段と、
前記剛性行列の擬似逆行列を算出する擬似逆行列算出手段と、
前記擬似逆行列を用いて前記構造体の前記荷重作用領域に於ける局所コンプライアンス行列を算出する局所コンプライアンス行列算出手段と、
前記局所コンプライアンス行列の少なくとも一つの固有値と対応する固有ベクトルとを算出する固有値解析手段と、
前記固有ベクトルと前記擬似逆行列とを用いて前記固有ベクトルを荷重ベクトルとしたときの前記構造体に於ける変位ベクトルを算出する変位ベクトル算出手段と、
前記変位ベクトルを用いて前記設計変数に対する前記少なくとも一つの固有値の感度を算出する固有値感度算出手段と、
前記体積制約割合を満たした状態で前記感度を用いて前記少なくとも一つの固有値が最小化されるよう前記設計変数を更新する設計変数更新手段と、
所定の終了条件が成立するまで、前記更新された設計変数を用いて、前記剛性行列算出手段、前記擬似逆行列算出手段、前記局所コンプライアンス行列算出手段、前記固有値解析手段、前記変位ベクトル算出手段、前記固有値感度算出手段及び前記設計変数更新手段の各処理を反復させる反復処理手段と、
前記設計変数により決定される構造体の最新の構造を出力又は表示する構造体出力手段と
を含む装置。 Structural optimization arithmetic unit that finds the optimal structure of a structure that has no support structure as a structure whose structure is determined by design variables under conditions set in the design area by numerical processing using a computer And
The dimensions of the design area, the range of the load application area in which a load acts on the structure, physical property constants representing the elastic properties of the materials constituting the structure, and the ratio of the volume of the structure to the volume of the design area Condition setting means for setting a condition including a volume constraint ratio that is
Stiffness matrix calculation means for calculating a stiffness matrix of the structure using the design variables,
Pseudo-inverse matrix calculating means for calculating a pseudo-inverse of the rigidity matrix,
Local compliance matrix calculation means for calculating a local compliance matrix in the load application region of the structure using the pseudo inverse matrix,
Eigenvalue analysis means for calculating at least one eigenvalue and a corresponding eigenvector of the local compliance matrix,
Displacement vector calculating means for calculating a displacement vector in the structure when the eigenvector is a load vector using the eigenvector and the pseudo inverse matrix,
Eigenvalue sensitivity calculation means for calculating the sensitivity of the at least one eigenvalue to the design variable using the displacement vector,
Design variable updating means for updating the design variable so that the at least one eigenvalue is minimized using the sensitivity while the volume constraint ratio is satisfied,
Until a predetermined termination condition is satisfied, using the updated design variables, the stiffness matrix calculating means, the pseudo-inverse matrix calculating means, the local compliance matrix calculating means, the eigenvalue analyzing means, the displacement vector calculating means, Iterative processing means for repeating each processing of the eigenvalue sensitivity calculating means and the design variable updating means,
A structure output unit that outputs or displays the latest structure of the structure determined by the design variable.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2018137736A JP2020016936A (en) | 2018-07-23 | 2018-07-23 | Arithmetic unit for structural optimization with unnecessary supporting structure |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2018137736A JP2020016936A (en) | 2018-07-23 | 2018-07-23 | Arithmetic unit for structural optimization with unnecessary supporting structure |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JP2020016936A true JP2020016936A (en) | 2020-01-30 |
Family
ID=69581804
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2018137736A Pending JP2020016936A (en) | 2018-07-23 | 2018-07-23 | Arithmetic unit for structural optimization with unnecessary supporting structure |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP2020016936A (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111444579A (en) * | 2020-03-11 | 2020-07-24 | 华中科技大学 | Composite material structure optimization design method considering manufacturability |
| CN114896728A (en) * | 2022-05-06 | 2022-08-12 | 大连理工大学 | Method and device for identifying external loads of structure, computer equipment and storage medium |
| CN117709164A (en) * | 2023-12-22 | 2024-03-15 | 河海大学 | Self-supporting structure topology optimization design method based on phase field method |
-
2018
- 2018-07-23 JP JP2018137736A patent/JP2020016936A/en active Pending
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111444579A (en) * | 2020-03-11 | 2020-07-24 | 华中科技大学 | Composite material structure optimization design method considering manufacturability |
| CN111444579B (en) * | 2020-03-11 | 2022-04-12 | 华中科技大学 | Composite material structure optimization design method considering manufacturability |
| CN114896728A (en) * | 2022-05-06 | 2022-08-12 | 大连理工大学 | Method and device for identifying external loads of structure, computer equipment and storage medium |
| CN117709164A (en) * | 2023-12-22 | 2024-03-15 | 河海大学 | Self-supporting structure topology optimization design method based on phase field method |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Jahangiry et al. | An isogeometrical approach to structural level set topology optimization | |
| Takezawa et al. | Topology optimization for worst load conditions based on the eigenvalue analysis of an aggregated linear system | |
| JP5377501B2 (en) | Structure optimization device, structure optimization method, and structure optimization program | |
| Wang et al. | A level set method for structural topology optimization | |
| US10850495B2 (en) | Topology optimization with microstructures | |
| Wang et al. | Isogeometric analysis for parameterized LSM-based structural topology optimization | |
| Luo et al. | Shape and topology optimization of compliant mechanisms using a parameterization level set method | |
| Guest et al. | Reducing dimensionality in topology optimization using adaptive design variable fields | |
| Zuo et al. | Evolutionary topology optimization of structures with multiple displacement and frequency constraints | |
| Rong et al. | Topological optimization design of structures under random excitations using SQP method | |
| JP2020016936A (en) | Arithmetic unit for structural optimization with unnecessary supporting structure | |
| Zuo et al. | Evolutionary topology optimization of continuum structures with a global displacement control | |
| Zuo et al. | An improved bi-directional evolutionary topology optimization method for frequencies | |
| Changizi et al. | Stress-based topology optimization of steel-frame structures using members with standard cross sections: Gradient-based approach | |
| Munk et al. | A simple alternative formulation for structural optimisation with dynamic and buckling objectives | |
| CN109002614B (en) | Improved level set topology optimization method for stable pore forming | |
| Modak | Model updating using uncorrelated modes | |
| JP2021101330A (en) | Geometrical dimension control in optimization | |
| Zhao et al. | A PEM-based topology optimization for structures subjected to stationary random excitations | |
| Rad et al. | Reliability based bi-directional evolutionary topology optimization of geometric and material nonlinear analysis with imperfections | |
| Gunwant et al. | Topology optimization of continuum structures using optimality criterion approach in ANSYS | |
| Soares | Nonlinear dynamic analysis considering explicit and implicit time marching techniques with adaptive time integration parameters | |
| Groh et al. | Localised post-buckling states of axially compressed cylinders and their energy barriers | |
| Zhang et al. | Finite strain topology optimization with nonlinear stability constraints | |
| Jang et al. | Design space optimization using design space adjustment and refinement |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20180816 |