WO2013065764A1 - 解析装置、解析方法及び解析プログラム - Google Patents
解析装置、解析方法及び解析プログラム Download PDFInfo
- Publication number
- WO2013065764A1 WO2013065764A1 PCT/JP2012/078262 JP2012078262W WO2013065764A1 WO 2013065764 A1 WO2013065764 A1 WO 2013065764A1 JP 2012078262 W JP2012078262 W JP 2012078262W WO 2013065764 A1 WO2013065764 A1 WO 2013065764A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- equation
- boundary condition
- variable
- equations
- dual
- Prior art date
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Definitions
- the present invention relates to a technique for analyzing a motion or a state of a continuum by representing a motion of an analysis target such as a continuum by a differential equation and calculating a solution according to a boundary condition.
- an object of the present invention is to provide an analysis apparatus, an analysis method, and an analysis program that can calculate a solution of a non-self-adjoint problem represented by simultaneous differential equations.
- the flowchart which shows the operation example of the analyzer in 4th Embodiment The figure which shows the example of the element of 16 node rectangle Diagram showing the analysis results with a circular condition under uniform gravity and with “zero displacement and zero surface force” as the boundary condition Figure showing the analysis result when the upper boundary condition is “zero displacement and zero surface force” for a square under uniform gravity
- the problem you want to solve is called the “Primal Problem”.
- the homogeneous boundary condition [Homogeneous Boundary Condition] is determined from the inhomogeneous boundary condition [Inhomogeneous Boundary Condition].
- the solution function of the main problem is called the main solution [Primal Solution], and this variation is called the main variation [Primal Variation].
- the adjoint differential operator [Adjoint Differential Operator] is determined from the original differential operator [Original Differential Operator], and the homogeneous adjoint boundary condition [Homogeneous Adjoint Boundary Condition] ] Is determined.
- the homogeneous adjoint boundary condition [Homogeneous Adjoint Boundary Condition] ]
- One of the elements is the inhomogeneous adjoint boundary condition [Inhomogeneous Adjoint Boundary
- the problem that gives this is called the “Dual Problem”.
- the solution function of a dual problem is called a dual solution [Dual Solution], and this variation is called a dual variation.
- both groups are tied together.
- the stress component is The components of the outward unit normal vector of the boundary surface are n x and n y , and the surface forces p x and p y are obtained from Cauchy's formula,
- the function groups ⁇ j and ⁇ j * should be eigenfunctions. If attention is paid to the equation (68), the external force term f Hj is composed of the function group ⁇ j * and the external force term f Hj * is composed of the function group ⁇ j , so that the coefficients c k and c k are obtained through the inner product. * It can be seen that there is a relationship between them.
- equation (70) is substituted into equations (45), (48), and (51) representing stress, It becomes. Therefore, if the stress generated by the function group ⁇ j * is expressed as ⁇ Ex * , ⁇ Ey * , ⁇ Exy * , It becomes. Surface forces caused by the function group ⁇ j * p Ex *, p Ey *, It becomes. This is, on the surface forces p Ei * by the set of functions ⁇ j *, is a preparation for to satisfy the boundary conditions.
- Equation (98) is called the primary simultaneous eigenvalue problem [Primal Simultaneous Eigenvalue Problem]
- equation (99) is called the dual simultaneous eigenvalue problem [Dual Simultaneous Eigenvalue Problem].
- the function ⁇ i is called a main eigenfunction [Primal Eigenfunction]
- ⁇ i * is called a dual eigenfunction [Dual Eigenfunction].
- the operators on the left are called the main differential operator [Primal Differential Operator] and the dual differential operator [Dual Differential Operator], respectively. Both (main, dual) differential operators are self-adjoint differential operators.
- a homogeneous boundary condition is imposed on the function ⁇ i
- a homogeneous adjoint boundary condition is imposed on the function ⁇ i * . Since there is a function ⁇ i * on the right side of the main simultaneous differential equation (87), the function ⁇ i satisfies the homogeneous boundary condition and the sum ⁇ j L ij ⁇ j of the differential coefficients satisfies the homogeneous adjoint boundary condition There is a need. These are collectively referred to as “Primal Boundary Condition”.
- ⁇ i in Eq. (98) and ⁇ i * in Eq. (99) imposed by these boundary conditions have orthogonality, and constitute basis functions of Hilbert space. For example, if m-th order eigenvalues and main eigenfunctions are set as ⁇ m and ⁇ im and n-th order eigenvalues and main eigenfunctions are set as ⁇ n and ⁇ in, Therefore, it is known that the main eigenfunction ⁇ i is orthogonal.
- the main simultaneous differential equation (87) is obtained from the following equations (93) and (97): It becomes.
- the dual simultaneous differential equation (88) has no deformation, It is.
- the functions ⁇ i and ⁇ i * used in the eigenfunction method in the next chapter are combinations that satisfy the (main, dual) simultaneous differential equations (105) and (106). However, ⁇ i remains as it is, and a combination in which the signs of ⁇ and ⁇ i * are reversed satisfies the requirement, but this combination is excluded.
- equation (111) is an eigenfunction method in the Hilbert space in which the displacement u Hj is expressed by the main eigenfunction ⁇ i and the external force term f Hi is expressed by the dual eigenfunction ⁇ i * .
- equation (115) is an eigenfunction method in the Hilbert space in which the displacement u Hj * is expressed by a dual eigenfunction ⁇ i * and the external force term f Hi * is expressed by a main eigenfunction ⁇ i .
- FIG. 1 shows the dual principle of simultaneous partial differential equations.
- the basis for the name “dual” is the above-mentioned feature “front and back integration / reversal use”.
- Clapeyron's theorem representing the law of conservation of energy If we take the variation in It becomes. From the above equation and (132) It becomes. This is the same form as the principle of complementary virtual work, but is obtained as a result of applying the energy conservation law to (123), and should be recognized as another form of the variational direct method. If the area obtained by removing the surface force specified area from the boundary surface is Cu, the displacement on Cu is the specified value. This is equivalent to the variation being zero. This is the principle of minimum work.
- linear elastic finite elements are formulated by Eq. (132), but no variational calculation is performed for h-method elements. If the variation symbol ⁇ is removed from Eq. (132) and multiplied by 1 ⁇ 2, the energy conservation law of Eq. (134) is matched, but the h-method element follows this. If a variational operation is performed using equation (123), an element equivalent to the p method can be obtained, and a new finite element can be constructed from equations (122) and (127).
- the inside of the boundary is expressed in the range of the inequality sign, and the boundary is expressed by the equal sign.
- In-element displacement Displacement that is the main solution is divided into a term that expresses boundary displacement and a term that has zero displacement at the boundary. It expresses.
- a subscript A represents a term representing boundary displacement
- a subscript o represents a term where displacement is zero at the boundary.
- the displacement u j as the main solution is represented by the sum of the term u Bj that satisfies the inhomogeneous boundary condition and the term u Hj that satisfies the homogeneous boundary condition. Expressed as a matrix, It becomes. The details are determined when the boundary conditions are given.
- it may be a Legendre ⁇ function or a trigonometric function.
- the function system is Becomes the trial function ⁇ ok .
- the differential coefficient of displacement is It becomes.
- the matrix on the right side is It becomes.
- the strain component can be calculated from Eqs. (4) to (6). It expresses.
- the stress component is calculated from equations (7) to (9). It becomes. This represents the stress component of equations (10) to (12).
- [I] is a unit matrix.
- the above equation represents nodal displacements U, nodal force F, the relationship must be established between the coefficient C o. In combination with the in-element displacement u of (173), It becomes.
- the matrix [K b ] is obtained by extracting a corresponding column from the matrices [K U ], [ ⁇ I] according to the arrangement of the known parts U b , F b , and similarly, the unknown parts U v .
- the matrix [K according to the order of F v v ], [-I] to extract the corresponding column, the matrix [K v ] is obtained.
- the matrix [Y b ], [Y v ] is obtained as follows.
- the operation of extracting a column from the matrix [Y v ] is performed by the matrix [ It can be expressed by extracting a column of c v ].
- Extract the relevant column from [c v ] and compose [c ve ], [c vh ] is defined by extracting a column of c v ].
- the matrix [ ⁇ 0 ] of the equation (263) is The operation of extracting a column from the matrix [ ⁇ 0 ] can be expressed by extracting the column of the matrix [I].
- the corresponding column is extracted from the matrix [I] to form [I oe ] and [I oh ], Therefore, the matrices [Y e ] and [Y h ] are It becomes.
- the function matrix And coefficient matrix If the matrix [Y e ] and [Y h ] are defined, It becomes.
- the first term is a function with a known part s b , which is a displacement [u B ] that satisfies the inhomogeneous boundary condition.
- the second term is a function having an unknown part s h, the displacement [u H] which satisfies the homogeneous boundary conditions. That is, Defined as And It becomes.
- the above equation shows the details of equation (174), and the displacement [u] is expressed by the displacement [u B ] that satisfies the inhomogeneous boundary condition and the displacement [u H ] that satisfies the homogeneous boundary condition. I understand that it was done.
- the matrix [K v * ] is obtained.
- the matrix [ Y b * ] and [Y v * ] are obtained as follows.
- the matrices [ ⁇ A ] and [0] of the equation (318) are The operation of extracting a column from the matrix [ ⁇ A ] is to extract the column of [I], and similarly, the operation of extracting a column from the matrix [0] is to extract the column of [0]. Can be expressed.
- a corresponding column is extracted from [I] to form [I Ub * ] and [I Uv * ].
- a corresponding column is extracted from “0” to form [0 Fb * ] and [0 Fv * ].
- the unknown parts s v * and c o * are rearranged to be unknown parts s e * and s h *, and the columns of the matrices [K v * ] and [K 0 ] are changed accordingly, and [K e * ] , [K h * ].
- the columns of the matrices [Yv * ] and [ ⁇ 0 ] are exchanged to be [Y e * ] and [Y h * ].
- equation (324) becomes It becomes.
- the purpose of this procedure is to erase the unknown part s e * so that [K e * ] has an inverse matrix. Therefore, [K e * ] is a square matrix.
- the unknown parts s e * and s h * are constructed by the following procedure. (1) n ve * pieces to fall to an unknown part s v * unknown part s e of the *, collected in the first half of the unknown part s e *. (2) the n oe * pieces to fall to an unknown part s e * of the unknown part c o *, gather in the second half of the unknown part s e *. (3) the n vh * pieces to fall to an unknown part s h * of the unknown part s v *, collected in the first half of the unknown part s h *.
- the operation of extracting a column from the matrix [Y v * ] can be expressed by extracting the column of the matrix [c v * ]. Extract the corresponding column from [c v * ] and [c ve * ], Configure [c vh * ] Define Using the unit matrix [I], the matrix [ ⁇ 0 ] of the equation (324) is The operation of extracting a column from the matrix [ ⁇ 0 ] can be expressed by extracting the column of the matrix [I].
- the corresponding column is extracted from the matrix [I] to construct [I oe * ], [I oh * ], Therefore, the matrices [Y e * ] and [Y h * ] are It becomes. In addition, the coefficient matrix , The matrices [Y e * ] and [Y h * ] are It becomes.
- the first term is a function with a known part s b * and a displacement that satisfies the inhomogeneous adjoint boundary condition [ u B * ].
- the second term is a function having an unknown part s h * , which is a displacement [u H * ] that satisfies the homogeneous adjoint boundary condition. That is, Defined as And It becomes.
- the above equation shows the details of the equation (191).
- the displacement [u * ] is a displacement [u B * ] that satisfies the inhomogeneous adjoint boundary condition and a displacement [u H that satisfies the homogeneous adjoint boundary condition. It can be seen that it is expressed as * ].
- equations (390) to (393) Since the main trial function [ ⁇ h ] satisfies the homogeneous boundary condition and the dual trial function [ ⁇ h * ] satisfies the homogeneous adjoint boundary condition, It becomes. That is, It is. Using equations (390) to (393), equations (388) and (389) are It becomes. In order for this to hold for any variation, It must be. By substituting equations (398) and (399) into each other, It becomes. By solving this, eigenvalues ⁇ and eigenvectors ⁇ e h ⁇ , ⁇ e h * ⁇ can be obtained. On the other hand, if the equations (398) and (399) are combined into one, It is. The equal sign in equation (403) is according to equation (395).
- eigenvalues ⁇ and eigenvectors ⁇ e h ⁇ , ⁇ e h * ⁇ can be obtained. Since [A s ] and [B s ] are both symmetric matrices, the eigenvalue ⁇ is a real number, and the eigenvector ⁇ e s ⁇ is orthogonal with [A s ] and [B s ] interposed therebetween.
- the eigenfunctions ⁇ and ⁇ * ⁇ are obtained by returning the components ⁇ e h ⁇ and ⁇ e h * ⁇ of the eigen vector ⁇ e s ⁇ belonging to the eigen value ⁇ to the expressions (315) and (375).
- the unknown part s h is determined in this way, the unknown part s e is determined by the equation (300). Therefore, the displacement ⁇ u ⁇ is completely determined from the equation (299) or (301). If the displacement is determined, the stress distribution, strain energy distribution, etc. inside the element can be obtained.
- the equation (432) is established. Since this means a homogeneous differential equation, the coefficient ⁇ a 0 ⁇ associated with the solution is arbitrary. In the previous section, we showed how to determine the arbitraryness by adding boundary conditions on the edges, but there is also a method in which the designer who performs the calculation decides himself. A system that displays in a diagram what deformation state or stress state is to be obtained according to the coefficient ⁇ a 0 ⁇ , and the designer himself / herself so that the boundary conditions can be realized while viewing this diagram. The coefficient ⁇ a 0 ⁇ is determined.
- the known quantities are ⁇ U 1 ⁇ and ⁇ F 2 ⁇ , and the unknown quantities are ⁇ U 2 ⁇ and ⁇ F 1 ⁇ . Disassembling the above equation and, It expresses. (438) Thus, an unknown quantity ⁇ U 2 ⁇ is determined. Substituting this into equation (437) Thus, an unknown quantity ⁇ F 1 ⁇ is determined. That is, all unknown quantities were obtained. Such a procedure It is established as a conventional means of the finite element method.
- equation (434) is obtained as an equivalent equation from the energy conservation law without including the calculation of variation, nothing needs to be bound by self-adjoint boundary conditions, and non-self-adjoint It should be good to solve the above equation for various boundary conditions.
- the matrix [K b ] is determined according to the arrangement of the known parts U b and F b .
- the matrix [K v ] is obtained.
- the deformation of the whole system is governed by the mode coefficient ⁇ a 0 ⁇
- the variational principle equivalent to the differential equation of the element that is, the least square method of the equations (122) to (128) is used. Use it to determine the appropriate ⁇ a 0 ⁇ . If the element is not bent at an appropriate value around ⁇ a 0 ⁇ , it can be considered that an accurate analysis result is obtained.
- the mode coefficient ⁇ a 0 ⁇ can be determined by substituting ⁇ u A ⁇ for ⁇ u ⁇ in equation (421).
- the functional of equation (128) is created for the entire system, and the equation (127) can be constructed by differentiating with the mode coefficient ⁇ a 0 ⁇ . As a result, the mode coefficient ⁇ a 0 ⁇ can be determined.
- the configuration coefficient is displayed so that the mode coefficient ⁇ a 0 ⁇ is linked to the slider and the designer adjusts the mode coefficient by the slider so that the deformation diagram is displayed. This makes it possible to bring out the designer's image. By refining the boundary conditions by narrowing down the possible deformation states, it is possible to study to obtain a better structure.
- FIG. 2 is a functional block diagram illustrating an example of the configuration of the analysis apparatus of the present embodiment.
- the analysis device 10 analysis system 10) illustrated in FIG. 2 can be configured by a computer including the calculation unit 1.
- the calculation unit 1 inputs design data of the analysis object, boundary condition data, and differential equation data including differential operators for the analysis object, calculates the solution of the differential equation to be analyzed, and outputs it as analysis result data.
- the calculation unit 1 calculates a solution u j according to the following equation.
- the following formula is the above formula (12 Same as 2).
- the calculation unit 1 inputs data indicating the shape and material of the analysis object as design data, and reads the data of the original differential operator L ij used for calculation based on the design data.
- the solution of the differential equation represented by the above equation is calculated using the primitive differential operator L ij . This makes it possible to calculate solutions for non-self-adjoint problems.
- the force f i act to be analyzed, but u i can be a displacement of the analyzed, f i and u i is not limited thereto.
- f i and u i may indicate other than force and displacement.
- an appropriate value corresponding to the analysis target is f i and u i can be set.
- the computing unit 1 may calculate the solution u j by the following equation.
- the following formula is the above formula (123) It is the same. In this case, the solution of the self-adjoint problem can be calculated.
- the calculation unit 1 can also calculate the solution u j to be analyzed in the following equation by the variational direct method.
- the following formula is the same as the above (127).
- the calculation of the solution u j by the variation direct method can calculate a solution in which the variation is zero in the functional function ⁇ in the following equation.
- FIG. 3 is a functional block diagram illustrating an example of the configuration of the analysis apparatus according to the second embodiment. 3 includes a setting unit 11, an adjoining boundary condition calculation unit 12, and a non-self-adjoint calculation unit 13.
- the setting unit 11 sets the primordial differential operator to be analyzed and the boundary condition of displacement.
- the setting unit 11 can receive an input of data representing a configuration or material to be analyzed as design data from a user, and can determine an original differential operator based on the configuration or material to be analyzed.
- the setting unit 11 can accept input of boundary condition data from the user.
- the accompanying boundary condition calculation unit 12 calculates the accompanying boundary condition from the boundary condition received by the setting unit 11. For example, the adjoint boundary condition calculation unit 12 multiplies the differential equation including the primitive differential operator included in the differential equation data read by the setting unit 11 by the dual eigenfunction ⁇ i * or the dual displacement u i *, that is, the inner product. Can be calculated from the boundary term obtained by partial integration.
- the accompanying boundary condition can be calculated using the above equation (91).
- u Hi when the homogeneous boundary conditions imposed on the p Hi imposed phi i, p Ei, such that the boundary section becomes zero ⁇ i *, p Ei * Can be obtained as an adjoining boundary condition.
- the non-self-adjoint calculation unit 13 calculates the main simultaneous differential equation and the dual simultaneous differential equation from the primitive differential operator, and uses the main simultaneous differential equation and the dual simultaneous differential equation, the boundary condition and the adjoining boundary condition, The solution of simultaneous differential equations is calculated by obtaining eigenfunctions and dual eigenfunctions.
- the main simultaneous differential equation can be defined as the above equation (105), and the dual simultaneous differential equation can be defined as the above equation (106).
- the adjoint differential operator L ij * can be calculated by the above equation (90).
- the non-self-adjoint calculation unit 13 uses the main simultaneous eigenvalue problem (the above formula (105)), the dual simultaneous differential equation (the above formula (106)), the boundary condition, and the adjoint boundary condition. 98)) and the dual simultaneous eigenvalue problem (the above equation (99)) can be defined to calculate the main eigenfunction ⁇ and the dual eigenfunction ⁇ * .
- the solution of the non-self-adjoint calculation unit 13 is not limited to the one using the main simultaneous eigenvalue problem (the above equation (98)) and the dual simultaneous eigenvalue problem (the above equation (99)).
- the main eigenfunction ⁇ and the dual eigenfunction ⁇ * can be calculated using the above equations (122) and (123).
- the main eigenfunction ⁇ and the dual eigenfunction ⁇ * can be calculated both analytically and approximately.
- the approximate eigenfunction can be calculated by the above formulas (122) and (123).
- the non-self-adjoint calculation unit 13 calculates a solution of the simultaneous equations of the main problem using the main eigenfunction ⁇ and the dual eigenfunction ⁇ * .
- the non-self-adjoint calculation unit 13 can calculate the coefficient c k * of the dual problem by the above equation (115), and calculates the solution function u Hj * of the dual problem by substituting this into the above (70). Can do.
- FIG. 4 is a flowchart showing an operation example of the analysis apparatus 10a shown in FIG.
- the setting unit 11 receives an input of design data from the user (S1).
- the design data is not limited to a specific object or format.
- the analysis target is a solid
- input of data representing the shape or material of the analysis target can be accepted.
- the analysis target is a fluid
- input of data representing the physical properties of the fluid can be accepted.
- input of data indicating the configuration and physical properties of each element may be accepted.
- the setting unit 11 can acquire design data not only by a user input but also by reading design data from an accessible recording medium or by downloading via a network.
- the setting unit 11 determines the original differential operator based on the design data acquired in S1, and determines the differential equation of the main problem (S2). For example, the setting unit 11 extracts a differential operator corresponding to the shape of the analysis target indicated by the design data from pre-recorded data, and sets the differential equation including the extracted differential operator as the differential equation of the main problem. Can do. Thus, the differential equation of the main problem can be determined according to the design data.
- the original differential operator is L ij will be described.
- the setting unit 11 further acquires boundary condition data (S3).
- the boundary condition data is data related to the boundary condition to be analyzed, and is not limited to a specific target or format. For example, data indicating at least one condition of displacement and force on the boundary surface of the analysis object can be input as the boundary condition. In the invention disclosed in the present application, it is possible to set both the displacement and the force on the boundary surface, which has been impossible in the past.
- the adjoining boundary condition calculation unit 12 generates adjoining boundary condition data using the differential equation set in S2 and the boundary condition data acquired in S3.
- the adjoining boundary condition can be calculated using the above equation (91).
- the boundary terms represented by the above formula (60) are expressed by the above formulas (208) and (229) indicating the boundary displacement, dual displacement, surface force, and dual surface force, respectively. ), Expression (207), and expression (228).
- the non-self-adjoint calculation unit 13 generates data indicating the main simultaneous differential equation (for example, the above equation (105)) and the dual simultaneous differential equation (the above equation (106)) from the primary differential operator L ij (S5).
- the non-self-adjoint calculation unit 13 uses the main eigenfunction ⁇ , the main simultaneous differential equation and the dual simultaneous differential equation calculated in S5, the boundary condition acquired in S2, and the adjoint boundary condition calculated in S4. And the dual eigenfunction ⁇ * is calculated (S6).
- the non-self-adjoint calculation unit 13 uses the main eigenfunction ⁇ and the dual eigenfunction ⁇ * calculated in S6 to calculate the solution function u Hj of the simultaneous equations of the main problem by, for example, the above equation (111) ) And equation (69) (S7). Further, the non-self-adjoint calculation unit 13 uses the main eigenfunction ⁇ and the dual eigenfunction ⁇ * to obtain the solution function u Hj * of the simultaneous equations of the dual problem, for example, the above formulas (115) and (70). Can be calculated. In this way, the calculated solution is output as an analysis result to a computer display (not shown) or the like (S8). A specific example of the analysis result will be described later.
- the operation example of the analysis device 10a has been described above, but the processing of the analysis device 10a is not limited to the above example.
- This embodiment can be said to be one specific example of the first embodiment.
- FIG. 5 is a functional block diagram illustrating an example of the configuration of the analysis apparatus according to the third embodiment.
- the analysis apparatus 10b illustrated in FIG. 5 includes a setting unit 11, an adjoint boundary condition calculation unit 12, a self-adjoint determination unit 15, a non-self-adjoint calculation unit 13, and a self-adjoint calculation unit 14.
- the setting unit 11, the accompanying boundary condition calculating unit 12, and the non-self accompanying calculating unit 13 can be configured in the same manner as in the first embodiment.
- the analysis apparatus 10b includes a self-adjoint determination unit 15 that determines whether or not the boundary condition matches the adjoint boundary condition calculated by the adjoint boundary condition calculation unit.
- the self-adjoint determination unit 15 determines that the boundary condition matches the adjoint boundary condition
- the self-adjoint calculation unit 14 solves the self-adjoint problem by obtaining a self-adjoint eigenfunction of the self-adjoint problem from the original differential operator.
- the non-self-adjoint operation unit 13 uses the main simultaneous differential equation and the dual simultaneous differential equation, and the boundary condition and the adjoint boundary condition. Calculate the solution of simultaneous differential equations by finding the function and dual eigenfunction.
- the adjoint differential operator L ij * can be calculated by partial integration of an integral sum or inner product obtained by multiplying a differential equation including a primitive differential operator by a dual eigenfunction ⁇ i * or a dual displacement u i * .
- the adjoint differential operator L ij * can be calculated by the above formula (62) or formula (90).
- the self-adjoint calculation unit 14 generates a simultaneous differential equation of the self-adjoint problem from the original differential operator and calculates a solution of the simultaneous differential equation.
- the self-adjoint eigenfunction ⁇ i can be calculated using the above equation (107) or equation (123) as a simultaneous differential equation of the self-adjoint problem.
- the self-adjoint calculation unit 14 calculates c k by, for example, an expression in which the self-adjoint eigenfunction ⁇ ik is replaced with ⁇ ik * in the above formula (111), and substitutes this c k into the above formula (69).
- the solution function u Hj of the self-adjoint problem can be obtained.
- FIG. 6 is a flowchart showing an operation example of the analysis apparatus 10b shown in FIG.
- the analysis apparatus 10b can execute the processing of S1 to S4 in the same manner as the processing of S1 to 4 in the second embodiment.
- the self-accompanying determination unit 15 compares the data indicating the boundary condition acquired in S1 with the accompanying boundary condition data generated in S4 to determine whether or not the boundary condition matches the accompanying boundary condition. Determine. If the self-accompaniment determination unit 15 determines in S21 that they do not match (that is, it is non-self-accompaniment), the non-self-accompanying operation unit 13 executes the processes of S5 to S7.
- the processing of S5 to S7 may be the same as S4 to S7 in the second embodiment.
- the self-adjoint calculation unit 14 If the self-adjoint determination unit 15 determines in S21 that the boundary condition and the non-boundary condition match (that is, is self-adjoint), the self-adjoint calculation unit 14 generates self-adjoint differential equation data from the primitive differential operator ( S22). For example, the self-adjoint calculation unit 14 can set the above formula (107) as a self-adjoint differential equation. In this case, the self-adjoint calculation unit 14 calculates the self-adjoint eigenfunction ⁇ i using the original differential operator L ij acquired in S3 as the self-adjoint differential operator and the boundary condition set in S2 as the self-adjoint boundary condition.
- the self-adjoint calculation unit 14 calculates the solution function uj of the self-adjoint simultaneous differential equation using the self-adjoint eigenfunction ⁇ i (S24). .
- the result output of S8 may be the same as S8 in the second embodiment.
- This embodiment can be said to be one specific example of the first or second embodiment.
- This embodiment is an example when the analysis apparatus is applied to an analysis using a finite element method (FEM).
- FEM finite element method
- the original differential operator and the boundary condition are set for a finite element.
- the configuration and analysis processing of the analysis apparatus described in the first to third embodiments can be similarly applied to finite elements.
- the present embodiment which uses the analysis processing shown in the first to third embodiments as a finite element, is different from the conventional finite element method will be described.
- a solution is calculated by a variational direct method (for example, the first implementation described above).
- Form a variational direct method
- the p method is a technique for approximating intra-element displacement by a higher-order function, and a function that appropriately approximates intra-element displacement is found and used.
- the main eigenfunction and the dual eigenfunction are calculated, and with these, the displacement in the element is intended to be explained appropriately.
- FIG. 7 is a flowchart illustrating an operation example of the analysis apparatus according to the present embodiment.
- the configuration of the analysis apparatus in the present embodiment can be the same as in FIG. 2 or FIG. Note that the flowchart shown in FIG. 7 is an example, and the present invention is not limited to this.
- the setting unit 11 first receives input of data indicating the shape to be analyzed (S31). For example, data indicating the shape of a finite element is input.
- the setting unit 11 creates a weight constant (S32).
- the weight constant w is preferably used so that the eigenvalue ⁇ is determined as a dimensionless quantity.
- the object to be analyzed is a ring
- the rank of the differential operator is second
- it is preferable that w 1 / r 2 if the representative length is r.
- the length of some characteristic part of the object is adopted as the representative length.
- the setting unit 11 inputs the nodal displacement and the nodal force (S33).
- the nodal displacement u and the dual displacement u * can be represented by a boundary function u A representing the displacement of the element boundary and a correction function u 0
- the correction function u 0 can be represented by the sum of the trial function ⁇ 0 (for example, The above formula (173) and formula (190)).
- the setting unit 11 sets a trial function (S34).
- the trial function can be set, for example, as in the above equation (179).
- the displacement at the boundary surface is the above equation (208)
- the dual displacement is the above equation (229)
- the surface force is the above equation (207)
- the dual surface force is It can be represented by the above formula (228).
- the accompanying boundary condition calculation unit 12 calculates the accompanying boundary condition from the boundary condition. For example, an integral sum or inner product obtained by multiplying a differential equation including a primitive differential operator by a dual displacement u i * is substituted into a boundary term R obtained by partial integration, and an adjoint boundary condition Can be calculated.
- the boundary term R can be expressed by the above formula (250).
- the formula (250) by dividing the matrix elements of the nodal force F to a known portion F b and the unknown portion F v sorting sorts separately nodal displacements U in known portion U b and unknown portion U v.
- the dual nodal force F * is also known part F b * and unknown part F v *. If the dual displacement U * is also divided into the known part U b * and the unknown part U v * , the following equation (251) is obtained.
- U b * and F b * corresponding to the unknown parts F v and U v must be zero, so that the adjoint boundary condition of Expression (254) is obtained.
- the coefficient vector ⁇ e h ⁇ is expressed by the above equation (122) and the equation (123).
- ⁇ e h * ⁇ the approximate eigenfunctions ⁇ and ⁇ * ⁇ can be obtained (S39).
- the main trial function [ ⁇ h ] can be expressed as the above equation (313).
- the function matrix ⁇ is defined by the above equation (293)
- the coefficient matrix [c e ] [c h ] is defined by the above equations (294) and (295).
- the dual trial function [ ⁇ h * ] can be expressed as, for example, the above equation (373).
- the coefficient matrix [c e * ] [c h * ] is defined by the above equations (354) and (355).
- the non-self-adjoint calculation unit 13 calculates the solution function of the simultaneous equations of the main problem by the eigenfunction method using the main eigenfunction and dual eigenfunction calculated in S39 (S40).
- the calculated solution outputs the analysis result to a computer display (not shown) or the like (S41).
- a computer display not shown
- S41 the like
- Equation (1) approximate eigenfunctions variation of ⁇ ⁇ calculated in S37, the approximation dual eigenfunction ⁇ phi * ⁇ variation of calculated at S38 and ⁇ .delta..phi * ⁇ Eigenvalues and eigenvectors can be obtained by solving the equations (eg, the above equations (400) and (401) or the above equation (402)) obtained by multiplying and integrating the main simultaneous differential equation and the dual simultaneous differential equation. Using these, the displacement ⁇ u ⁇ of the main problem is obtained.
- equations eg, the above equations (400) and (401) or the above equation (402)
- [Fifth embodiment] 67 includes an initial equation determination unit 21, a boundary condition determination unit 22, and a calculation unit 23.
- the initial equation determination unit 21 reads data indicating the structure of the system to be processed and the properties of the components, and determines an equation representing the system.
- This equation can be, for example, n equations including a variable representing a physical quantity to be obtained.
- the equation can be represented by, for example, matrix data or differential equation data, and the data representation format is not particularly limited.
- the variable can be represented by data such as a vector or a function, for example, but is not particularly limited to an expression format.
- the initial equation determination unit 21 can accept, for example, designation of a system to be processed from the user, and can select a system corresponding to the designated system from among the simultaneous differential equations recorded in advance. Or the structure which receives simultaneous simultaneous differential equations themselves from a user may be sufficient.
- the boundary condition determination unit 22 reads a value representing a physical quantity as data indicating the boundary condition and determines the boundary condition.
- the boundary condition determination unit 22 also receives a boundary condition that specifies an external force and a displacement for one direction of one point, that is, one degree of freedom. Even if a plurality of physical quantities are defined as boundary conditions for one degree of freedom by the function of the calculation unit 23 described later, it can be handled. Further, as a boundary condition, not only a known part but also an unknown part can be specified. For example, it is possible to accept a setting in which neither displacement nor external force is unknown for a certain node.
- the computing unit 23 transforms the n initial equations into 2n variable equations or equations including 2n equations, and in the modified 2n variable equations or equations including 2n equations, A known part including a variable that becomes known according to a condition and an unknown part including an unknown variable are determined, and a solution of the equation is calculated for the unknown part.
- the equation can be transformed into a format including twice as many equations as the initial equation, as shown in Equation (963).
- an equation represented by a matrix and a variable vector is determined as an initial equation, for example, as shown in Expression (442), the matrix can be transformed into a matrix including a variable vector having twice the degree of freedom of the variable vector. it can.
- the calculation unit 23 can output the calculated solution or information obtained from the solution. For example, when an object is displaced by receiving a force in a system to be processed, it is possible to output a deformation state, a stress distribution, an error distribution thereof, and the like as information obtained from the solution.
- the initial equation determination unit 21 can determine, for example, a simultaneous differential equation that represents a relationship between a force acting on an object that is a component of the system and a variable that represents a physical quantity of the object.
- the physical quantity of the object can be, for example, displacement or speed, but is not limited to a specific one.
- the calculation unit 23 generates data indicating the 2n equations using the differential operators of the simultaneous differential equations determined by the initial equation determination unit 21 and the adjoint differential operators determined from the differential operators, and generates 2n equations. Calculate the solution.
- the calculation unit 23 newly introduces dual simultaneous equations (including n equations) including dual variables to the differential operators and variables of the simultaneous differential equations of the main problem representing the target system. As a result, the solution of 2n equations is calculated. In this way, by using the adjoint differential operator determined from the differential operator, the 2n equations are generated and the solutions of the 2n equations are calculated, thereby setting a plurality of physical quantities as boundary conditions for one degree of freedom. If so, the solution can be calculated.
- the inner product of a variable and a differential variable of a simultaneous differential equation and a dual variable causes the variable and the differential operator of the dual simultaneous differential equation to act on the dual variable. It is equal to the inner product with the thing.
- examples of obtaining dual differential equations for various simultaneous differential equations are disclosed.
- Information regarding the dual differential equation can be recorded in advance in association with the differential equation.
- the calculation unit 23 can read out information on the dual differential equation corresponding to the differential equation determined by the initial equation determination unit 21 and use it for the calculation.
- the calculation unit 23 can calculate the solution of the equation by transforming the initial equation into an equation including 2n variables without calculating the adjoint operator.
- the initial equation determination unit 21 can determine an initial equation indicating a physical quantity at a node of a system component by determining two n-dimensional variable vectors and an n-row matrix.
- an n-row matrix representing a relationship between an n-dimensional vector of a variable representing the displacement of the node and an n-dimensional vector of a variable representing the force of the node can be generated. This can be said to represent the differential equation to be satisfied by the system using an n-dimensional vector and an n-row matrix.
- the initial equation can be determined in the form as shown in Expression (441), Expression (449), and the like.
- n equations are represented by an n-by-n matrix K and n-dimensional vector variable vectors ⁇ U ⁇ , ⁇ F ⁇ .
- the boundary condition determination unit 22 is configured to be able to determine boundary conditions in which the number of degrees of freedom of a variable whose value is known is different from the number of degrees of freedom of a variable whose value is unknown. For example, even when a user specifies that a plurality of physical quantities are unknown or known for one degree of freedom as a boundary condition, an error such as improper setting by the user is not generated. Can do.
- the calculation unit 23 generates a 2n-dimensional vector based on the two variable vectors. Furthermore, after transforming the matrix into a 2n-column matrix based on the 2n-dimensional vector variable, among the 2n-dimensional vector variable, a known part including a variable that is known by a boundary condition, and an unknown variable And determine the unknown part. At this time, the degree of freedom of the unknown part and the degree of freedom of the known part may not be the same. For example, the calculation unit 23 can determine the unknown part and the known part without adding a restriction that the degrees of freedom of the unknown part and the known part are the same.
- the calculation unit 23 modifies the 2n-column matrix and the 2n-dimensional vector so that the unknown part variable is represented by the known part variable, and uses the deformed matrix to change the unknown part variable. calculate. For example, the calculation unit 23 can rearrange the elements of the variable vector so that the known part and the unknown part of the variable are grouped, and can rearrange the elements of the matrix accordingly. This specific example is shown in equations (441) to (447) in section 9.2.
- the computing unit 23 can calculate solutions for both self-associated problems and non-self-associated problems by executing the above processing.
- the boundary condition that the displacement and force are zero for the x direction of a certain node or the boundary condition that the displacement and force are unknown for the y direction of a certain node. be able to. It is possible to calculate a solution even under a non-self-adjoint boundary condition by an operation that does not limit the number of degrees of freedom of the unknown part and the known part with the equation in which the variables are multiplied by n as described above.
- the mode coefficient for the homogeneous solution is determined so that the variation is zero, and the solution is obtained using the determined mode coefficient. it can. Further, it is also possible to adopt a configuration in which an input of a value near the determined mode coefficient is received from the user, a solution is calculated using the mode coefficient of the input value, and the solution or information obtained from the solution is output. Note that these configurations are not necessarily provided in the information processing apparatus provided with the calculation unit 23. For a plurality of solutions calculated by the calculation unit 23, a process for obtaining a solution using the above-described mode coefficient is performed in other configurations. The information processing apparatus may execute. In addition, another information processing apparatus may execute a process of receiving an input of a value near the mode coefficient from the user and calculating and outputting a solution using the mode coefficient of the input value.
- the solution is expressed by multiplying the special solution and the homogeneous solution by an arbitrary mode coefficient. be able to.
- the differential equation can be approximated as a whole structure system by the mode coefficient obtained when the solution expressed using this mode coefficient is applied to the equation in which the variation is zero in the functional ⁇ of the equation (128). A satisfactory solution is obtained. Examples of processing when there are a plurality of such solutions are described in, for example, sections 9.4, 9.5, and 11.3.10.
- the calculation unit 23 calculates the solution of the equation with the degree of freedom doubled in a situation where the degrees of freedom of the unknown part and the known part are not necessarily the same.
- a plurality of solutions may be calculated. Even in such a case, an appropriate result can be output by calculating a representative solution from a plurality of solutions and outputting the solution within a predetermined range including this solution, as in the above example.
- a user interface is provided that can convey the results of a plurality of solutions in a recognizable form by displaying the solutions varied by a user operation.
- An information processing apparatus 10d illustrated in FIG. 68 includes an initial equation determination unit 21, a boundary condition determination unit 22, and a determination unit 33.
- the initial equation determination unit 21 and the boundary condition determination unit 22 can be configured similarly to the initial equation determination unit 21 and the boundary condition determination unit 22 of FIG.
- the determination unit 33 determines whether or not the boundary condition of the simultaneous differential equations determined by the initial equation determination unit 21 is a self-adjoint condition.
- the determination unit 33 is obtained by partial integration of the integral sum, that is, the inner product, of the simultaneous differential equations and the dual variables.
- the condition of the dual variable that makes the boundary term zero can be determined as the adjoining boundary condition.
- the determining unit 33 can determine that the adjoining condition is a self-adjoining condition.
- an expression for calculating the boundary term in association with the differential equation can be recorded in advance in a recording unit accessible by the information processing apparatus 10d.
- the determination unit 33 can read the boundary term expression corresponding to the differential equation determined by the initial equation determination unit 21 from the recording unit and calculate the boundary term.
- the boundary term R is changed to the transposed matrix of the nodal external force F and the displacement U, the dual nodal external force F *, and the dual displacement U *. It can be expressed by the product of the matrix.
- the nodal external force F, the displacement U, the dual nodal external force F *, and the dual displacement U * are divided into the unknown part and the known part, respectively. , An adjoint boundary condition is obtained.
- the adjoint boundary condition is obtained by sorting the unknown part and the known part of each variable separately in the boundary term expressed by the product of the transpose matrix of the variable of the differential equation and the matrix of the dual variable.
- the expression (251) or the expression (1814) can be used, for example, when the operator and the boundary condition of the initial differential equation are set to finite elements and the analysis is executed.
- the determination unit can output, for example, information based on the determination result to the user as to whether or not it is a self-accompanying condition. For example, when it is determined that the boundary condition does not coincide with the accompanying boundary condition and is a non-self accompanying condition, a message can be output indicating that the analysis target is non-self accompanying in the boundary condition.
- the determination unit 33 needs to have a function capable of handling a non-self-adjoint problem for analysis, or a guide for analyzing a non-self-adjoint problem such as where to obtain or obtain software that can analyze a non-self-adjoint problem. Can be output.
- the determination unit 33 can also control subsequent calculations based on information on whether the adjoining boundary condition and the boundary condition match. For example, when the adjoining boundary condition and the boundary condition match, the analysis process using the conventional finite element method is executed, and when the adjoining boundary condition and the boundary condition do not match, as described above, It is possible to cause the information processing apparatus to execute processing capable of handling non-self-association.
- a process that can handle non-self-adjoint is, for example, transforming an initial equation into an equation including 2n equations, including a modified 2n equation, and a known part including a known variable and an unknown variable. The unknown part to be included is determined, and the solution of the equation is calculated for the unknown part.
- the self-adjoint condition may be determined without calculating the adjoint boundary condition.
- the number n1 of the known portions Ub of the nodal displacement U is set to the unknown portion Fv of the nodal external force.
- the number n2 of known parts Fb of the nodal external force F is equal to the number of unknown parts Uv of the nodal displacement. Therefore, the determination unit 33 can determine whether or not it is a self-adjoining condition based on whether or not the degree of freedom of the unknown part of one variable is equal to the degree of freedom of the known part of another variable.
- a configuration in which the determination unit 33 in FIG. 68 is added to the information processing apparatus in FIG. 68 is also possible. For example, based on the determination result of the determination unit 33, it is possible to control whether or not the calculation unit 23 performs a calculation using a 2n-dimensional equation. As a result, the calculation unit 23 can appropriately switch between a non-self-adjoint operation using a 2n-dimensional equation and a self-adjoint operation using an n-dimensional equation to improve processing efficiency.
- the boundary condition determination unit 22 can accept designation of known and unknown variables from the user.
- the determination unit 33 determines the unknown part and the known part of the variable based on the designation received from the user, and whether the degree of freedom of the unknown part of one variable is equal to the degree of freedom of the known part of the other variable. Outputs whether or not.
- the calculation unit 23 uses the n-row matrix representing the differential operator. Calculating the unknown part of the variable.
- the differential equation includes 2n variables and 2n equations.
- the model is transformed into a form, and the solution of the equation is calculated for the unknown part among the transformed 2n equations.
- FIG. 69 is a flowchart showing an example of processing for setting a nodal displacement vector and a nodal external force vector on the premise that the boundary condition satisfies the self-adjoint condition in the analysis using the conventional finite element.
- the information processing apparatus first receives an input of nodal displacement (S101).
- S101 nodal displacement
- the degree of freedom in which the displacement is known and its numerical value are input from the user.
- data indicating the degree of freedom and numerical value of the nodal displacement may be read.
- the degree of freedom is specified by, for example, the position and direction of the node.
- the degree of freedom of the known portion Ub of the nodal displacement and its numerical value are determined (S102).
- the degree of freedom and the numerical value of the known part Ub of the nodal displacement are determined, the degree of freedom of the unknown part Fv of the nodal external force can be automatically determined (S103). This is because, when the self-adjoint condition is satisfied, the degree of freedom of the unknown part Fv is equal to the known part Ub.
- the degree of freedom of the unknown part Uv of the nodal displacement is obtained by subtracting the degree of freedom of the known part Ub from the total degree of freedom (S104).
- the degree of freedom of the known part Fb is equal to the degree of freedom of the unknown part Uv. Therefore, the degree of freedom of the known part Fb is determined, and the numerical value of the known part Fb is assumed to be zero (S105). .
- the information processing apparatus corrects the value of the known part Fb after confirming that the known part Fb of the nodal external force does not contradict the degree of freedom.
- the degree of freedom of the external force where no value is specified.
- the degree of freedom satisfies the self-adjoint boundary condition, it can be analyzed by the conventional finite element method based on the set nodal displacement and nodal external force (S108).
- FIG. 70 is a flowchart showing an example of processing for setting a nodal displacement vector and a nodal external force vector that can be handled even when the boundary condition is a non-self-adjoint condition in an analysis using a finite element.
- the boundary condition determination unit 22 of the information processing apparatus accepts input of degrees of freedom of node displacement and numerical values (S201), and determines and records the degrees of freedom and numerical values of the known part Ub of node displacement ( S202).
- the unknown part Fv (S203) of the nodal external force, the unknown part Uv (S204) of the nodal displacement, and the known part Fb (S205) are automatically calculated.
- zero is set to the degree of freedom of the nodal external force for which no value is specified.
- the boundary condition determination unit 22 receives the degree of freedom of the nodal external force and numerical values (S206), and corrects the degree of freedom and the numerical values of the nodal external force Fb to the input degrees of freedom and numerical values (S207). Further, the boundary condition determination unit 22 accepts an input of an unknown part of the nodal displacement (S208).
- the degree of freedom of a node whose displacement is unknown is input. For example, a configuration in which the user can input by designating a node and a direction on the screen may be used, or a configuration in which data indicating the degree of freedom of the unknown part may be read.
- the boundary condition determination unit 22 corrects the degree of freedom of the unknown part Uv based on the input information (S209).
- the boundary condition determination unit 22 accepts an input of an unknown part of the nodal external force (S210), and corrects the degree of freedom of the unknown part Fv based on the input information (S211). In S208 to S211, it is also possible to accept the designation that both the displacement and the external force are unknown parts for the same degree of freedom. That is, the boundary condition determination unit 22 can set a condition in which different variables are unknown for the same degree of freedom.
- the determination unit 33 determines whether or not the degree of freedom satisfies the self-adjoint boundary condition (S212), and displays the determination result (S213).
- the determination unit 33 can determine, for example, whether or not the number of degrees of freedom of the known part Ub of the nodal displacement is equal to the number of unknown parts Fv of the nodal external force. If they are equal, the self-adjoining boundary condition is satisfied, and if they are not equal, it can be determined that the self-adjoining boundary condition is not satisfied.
- Non-self-adjoint boundary conditions can be analyzed (S215).
- the self-adjoint boundary condition can be analyzed by a conventional finite element method (S214).
- the analysis processing of S214 or S215 may be executed by another device.
- the determination unit 33 determines that the self-adjoint boundary condition is satisfied in S212
- the determination unit 33 performs the analysis process using the conventional finite element method in S214 and determines that the self-adjoint boundary condition is not satisfied.
- information can be output to enable the analysis processing in step S215.
- it is possible to output information such as how to obtain software and information for executing the processing of S215.
- FIG. 71 is a flowchart showing another example of processing for setting a nodal displacement vector and a nodal external force vector that can be handled even when the boundary condition is a non-self-adjoint condition.
- input of degrees of freedom and numerical values of nodal displacement and nodal external force is accepted, and the degree of freedom without input is set as an unknown part.
- the boundary condition determination unit 22 receives the input of the degree of freedom of the node displacement and the numerical value (S301), and based on the input, determines the degree of freedom and the numerical value of the known part Ub of the node displacement (S302), and the unknown The degree of freedom of the part Uv is also determined (S303). Similarly, the input of the external force of the nodal point and the input of a numerical value are accepted (S304), and based on the input, the degree of freedom and the numerical value of the known part Fb of the nodal external force are determined (S305). The degree is also determined (S306).
- the determination process (S307), the display process (S308), and the analysis process (S309, 310) can be executed in the same manner as in FIG.
- the nodal displacement and the degree of freedom of the nodal external force that are not particularly input are set as unknown parts.
- the user does not need to explicitly specify the unknown part, but for example, it is necessary to input to explicitly set a value of zero for a node where the displacement is desired to be zero.
- the variable of the degree of freedom in which the user did not designate the unknown part or the known part is set to some value as the known part. This eliminates the need for the user to set whether it is unknown or known for all degrees of freedom.
- the information processing apparatus can be configured such that, for example, the user can select which of the processes in FIG. 70 or 71 is used for input.
- the processing examples in FIGS. 70 and 71 correspond to the two processing methods (1) and (2) described in 11.7.1.
- FIG. 72 is a flowchart illustrating an execution example of analysis of a non-self-adjoint problem using a finite element.
- the information processing apparatus 10c inputs a calculation model (S401). For example, a finite element model, shape, material, element, node, etc. of the structure to be processed are input.
- the initial equation determination unit 21 creates an element stiffness matrix using the data input in S401, and creates an overall stiffness matrix by superimposing the element stiffness matrices (S405).
- the overall stiffness matrix created here is an example of an initial equation expressed by a matrix.
- S401 for example, two n-dimensional variable vectors indicating physical quantities at nodes of a system having n degrees of freedom and n rows of element stiffness matrices are determined.
- the matrix can be represented by the formula (441) or the formula (449).
- the boundary condition determination unit 22 inputs gravity load data (S403) and creates a gravity load (S405).
- S403 for example, data indicating the value of gravity or the value of displacement applied to each degree of freedom specified by the node and direction is input.
- the boundary condition determination unit 22 also inputs data specifying an unknown part of the nodal displacement and an unknown part of the nodal external force (S404). Based on the data input in S403 and S404, a known part and an unknown part of the nodal displacement and a known part and an unknown part of the nodal external force are created (S405). Note that the processes of S403, S404, and S405 can be executed in the same manner as the process shown in FIG. 70, for example.
- the calculation unit 23 deforms the rigid element matrix to create a matrix used for calculation (S406). For example, based on two variable vectors with n degrees of freedom, the calculation unit 23 generates a 2n-dimensional vector, and transforms the n-row matrix into a 2n-column matrix in accordance with the 2n-dimensional vector variable. . For example, it can be transformed into a form as shown in Formula (450).
- the computing unit 23 determines a known part including a variable that is known by a boundary condition and an unknown part including an unknown variable among the 2n-dimensional vector variables.
- the 2n-column matrix and the 2n-dimensional vector are transformed so that the unknown part variable is represented by the known part variable.
- the variable of the unknown part can be calculated using a matrix transformed as in Expression (452).
- the computing unit 23 calculates the solution of the matrix equation created in S406 (S407).
- the calculation unit 23 creates a special solution portion for the nodal displacement and a special solution portion for the nodal external force (S408).
- a vector representing the unknown part and the known part of the special solution can be created as in Expression (936).
- the calculation unit 23 defines a homogeneous mode coefficient (S409). For example, if the number of homogeneous solutions is n o , an n-dimensional vector ⁇ a o ⁇ can be defined as a mode variable.
- the calculation unit 23 creates a homogeneous solution portion of the nodal displacement and a homogeneous solution portion of the nodal external force (S410). For example, as shown in Expression (937), a vector representing the unknown part and the known part of the homogeneous solution is created.
- the calculation unit 23 determines the unknown part of the nodal displacement solution using, for example, Expression (938), and determines the unknown part of the nodal external force solution, for example, using Expression (939). If the solution can be determined uniquely as a result of the process of S411, the result is output and the analysis process is terminated (S412).
- the optimal value of ⁇ a o ⁇ is calculated by the method of least squares.
- the calculation unit 23 creates a matrix such as Expression (960) and Expression (961) (S413).
- the calculation unit 23 determines the mode coefficient by solving the equation (962) for the mode coefficient ⁇ a o ⁇ (S414).
- the computing unit 23 outputs a solution corresponding to a predetermined range of mode coefficients centered on the determined mode coefficient (S415). For example, an input of a value near the determined mode coefficient is received from the user, a solution is calculated using the mode coefficient of the input value, and the calculated solution or information obtained from the solution is output.
- a deformation state, a stress distribution, an error distribution, and the like are output.
- a user interface can be provided so that the user can designate a value within a predetermined range including the mode coefficient determined in S414 by operating ⁇ a o ⁇ with a slider.
- the program for executing the processing in S408 to S410 when there are a plurality of solutions and the program for executing the processing of S401 to S407 may be provided separately. That is, a program for executing at least one step of S408 to S410 or a recording medium on which the program is recorded is also included in the embodiment of the present invention.
- FIG. 9 is a diagram showing a result of determining how to support the outer side in order to make the inner boundary condition “zero displacement and zero surface force” in a ring under uniform gravity. .
- the left diagram shows the stress ⁇ r
- the center diagram shows the stress ⁇
- the right diagram shows the shear stress ⁇ r ⁇ .
- the result shown in FIG. 9 is a result calculated analytically without making the ring a finite element.
- An example of a specific main simultaneous differential equation used for the analysis is shown below.
- FIG. 10 is a diagram showing an analysis result in a case of a square under uniform gravity and the boundary condition of the upper side is “displacement is zero and surface force is zero”. That is, it is the result of the problem of how to support the other side in order to make both the displacement of the upper side of the square and the surface force zero under uniform gravity.
- the left figure shows the stress in the horizontal direction (x direction)
- the middle figure shows the stress in the vertical direction (y direction)
- the right figure shows the shear stress.
- the displacement of the upper side is kept at zero by the large stresses ⁇ x and ⁇ y acting on the lower part of the square.
- the square is formulated as a finite element, and there is only one element.
- P method finite element As a method that can be analyzed with one element.
- the analysis shown in FIG. 10 is an example of a non-self-associated boundary condition, and the P-method finite element method cannot be analyzed under such a boundary condition.
- Whether a boundary condition is given to an edge or a node depends on the formulation, but here it is formulated to give a boundary condition to a node.
- FIG. 11 is a diagram showing an analysis result of in-plane deformation of a square flat plate with a fixed periphery.
- the left figure in FIG. 11 shows the stress in the horizontal direction (x direction), the middle figure shows the stress in the vertical direction (y direction), and the right figure shows the shear stress. This is a self-associated boundary condition.
- the analysis result shown in FIG. 11 is obtained by performing an analysis equivalent to the P method finite element.
- the analysis method of the present invention is considered to include the P method. That is, if the finite element method is coded based on the present invention, the P method can be used seamlessly.
- FIG. 12 is a diagram showing a result of analyzing deformation due to its own weight when the periphery of a square out-of-plane deformed flat plate element is fixed.
- the square is formulated as a finite element, and there is only one element.
- FIG. 13A is a diagram of a model of a left-end free-fixed beam
- FIG. 13B is a diagram illustrating main eigenfunctions ⁇ 1 , ⁇ 2 , ⁇ 3 and dual eigenfunctions ⁇ 1 * , ⁇ 2 * , The calculation result of ⁇ 3 * is shown.
- FIG. 3C is a diagram showing a solution (1 mode, 2 mode and 3 mode) and an analytical solution obtained by the eigenfunction method. The three-mode solution almost overlaps with the analytical solution.
- FIG. 14 and FIG. 15 are diagrams showing analysis results when a square virtual boundary is defined in an infinite region and a spring (FIG. 14) or vorticity (FIG. 15) is arranged at the center.
- the analysis results shown in FIGS. 14 and 15 show the state of the potential flow calculated by the eigenfunction method with the Dirac ⁇ function given to the right side of the above equation (12) or (13). From these figures, it can be confirmed that the internal flow is close to the analytical solution.
- FIG. 17 shows the respective eigenfunctions normalized so that the maximum value is 1.
- FIG. 18 shows a comparison between the analytical solution of Eq. (525) and the result obtained by the eigenfunction method of Eq. (494).
- the broken line indicates the analytical solution
- the solid line indicates the eigenfunction method.
- 50M indicates that 50 modes have been calculated. It can be seen that the greater the number of modes, the better the accuracy of the eigenfunction method, and the displacement u ⁇ agrees well with about 10 modes.
- FIG. 19 shows a comparison between the analytical solution of equation (528) and the result obtained by the eigenfunction method of equation (494).
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 50M indicates that 50 modes have been calculated. It can be seen that the greater the number of modes, the better the eigenfunction method accuracy, and the displacement u ⁇ agrees well with about 10 modes.
- FIG. 20 shows a comparison between the analytical solution of equation (531) and the result obtained by the eigenfunction method of equation (494).
- the red broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 50M indicates that 50 modes have been calculated. It can be seen that the greater the number of modes, the better the eigenfunction method accuracy, and the displacement u ⁇ agrees well in about 10 modes.
- FIG. 24 shows a comparison between the analytical solution of equation (625) and the result obtained by the eigenfunction method of equation (578).
- the broken line indicates the analytical solution
- the solid line indicates the eigenfunction method.
- 56M indicates that up to 56 modes have been calculated. It can be seen that the greater the number of modes, the better the eigenfunction method accuracy, and the displacement u ⁇ agrees well with about 10 modes.
- FIG. 25 shows a comparison between the analytical solution of equation (632) and the result obtained by the eigenfunction method of equation (584).
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 50M indicates that 50 modes have been calculated. It can be seen that the greater the number of modes, the better the eigenfunction method accuracy, and the displacement u ⁇ agrees well with about 10 modes.
- FIG. 27 shows the respective eigenfunctions normalized so that the maximum value is 1.
- FIG. 28 shows a comparison between the analytical solution of equation (706) and the result obtained by the eigenfunction method of equation (702).
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 30M indicates that up to 30 modes have been calculated. It can be seen that the greater the number of modes, the better the accuracy of the eigenfunction method, and the displacement u ⁇ agrees well in about 3 modes.
- FIG. 29 shows a comparison between the analytical solution of equation (709) and the result obtained by the eigenfunction method of equation (702).
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 30M indicates that up to 30 modes have been calculated. It can be seen that the greater the number of modes, the better the eigenfunction method accuracy, and the displacement u ⁇ agrees well in about 3 modes.
- FIG. 30 shows a comparison between the analytical solution of equation (713) and the result obtained by the eigenfunction method of equation (704).
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 30M indicates that up to 30 modes have been calculated. It can be seen that the greater the number of modes, the better the eigenfunction method accuracy, and the displacement u ⁇ agrees well in about 3 modes.
- FIG. 31 shows a comparison between the analytical solution of equation (715) and the result obtained by the eigenfunction method of equation (704).
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 30M indicates that up to 30 modes have been calculated. It can be seen that the greater the number of modes, the better the accuracy of the eigenfunction method, and the displacement u ⁇ agrees well in about 3 modes.
- FIG. 35 shows a comparison between the analytical solution of equation (784) and the result obtained by the eigenfunction method of equation (780).
- each function value is plotted by a negative sign. Only the coordinate value on the vertical axis was read with the sign reversed.
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 32M indicates that 32 modes have been calculated. The more modes, the more The accuracy of the eigenfunction method is improved, and the displacement u ⁇ agrees well in about 3 modes.
- FIG. 36 shows a comparison between the analytical solution of equation (787) and the result obtained by the eigenfunction method of equation (780).
- each function value is plotted by a negative sign. Only the coordinate value on the vertical axis is read in reverse.
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 32M indicates that 32 modes have been calculated. The more modes, the more The accuracy of the eigenfunction method is improved, and the displacement u ⁇ agrees well in about 3 modes.
- FIG. 37 shows a comparison between the analytical solution of Equation (790) and the result obtained by the eigenfunction method of Equation (782).
- each function value is plotted with a negative sign. Only the coordinate value on the vertical axis is read in reverse.
- the red broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 32M indicates that 32 modes have been calculated. It can be seen that the greater the number of modes, the better the accuracy of the eigenfunction method, and the displacement u ⁇ agrees well in about 3 modes.
- FIG. 38 shows a comparison between the analytical solution of Equation (793) and the result obtained by the eigenfunction method of Equation (782).
- each function value is plotted by a negative sign. Only the coordinate value on the vertical axis is read in reverse.
- the broken line indicates the analytical solution
- the black solid line indicates the eigenfunction method.
- 32M indicates that 32 modes have been calculated. The more modes, the more The accuracy of the eigenfunction method is improved, and the displacement u ⁇ agrees well in about 3 modes.
- [Other analysis results] 39 is a diagram illustrating one of solutions other than the solution illustrated in FIG. 10 among the plurality of solutions when the plurality of solutions are generated due to the eigenvalue being zero in the analysis illustrated in FIG. 10. .
- the designer can study a method that can be realized in design.
- 40A and 40B are examples of a plurality of solutions calculated when the displacement of the upper surface of the square plate under uniform gravity and the displacement that makes the surface force zero are analyzed using a finite element.
- a square is formed by 6 ⁇ 6 ⁇ finite elements in the vertical and horizontal directions.
- a circle means a node.
- a total of 10 modes are shown. Any combination of these is the solution.
- a plurality of solutions can be obtained by the calculation shown in the section 9.2 above.
- the mode coefficient is linked to the slider so that a deformation diagram is displayed as the designer adjusts the mode coefficient by the slider, the designer can have any support method. It will be possible to examine the designer's image. In addition, the designer can narrow down the feasible support methods and then calculate again.
- the calculation unit 1, the setting unit 11, the adjoint boundary condition calculation unit 12, the non-self-adjoint calculation unit 13, the self-adjoint determination unit 15, and the self-adjoint calculation unit 14 are executed by the computer processor from the memory.
- the program and data are read and data processing is executed according to the program.
- the calculation unit 1 may acquire design data and boundary condition data by receiving input from a user via an input device and a user interface, or may be acquired by reading from a memory.
- the differential equation data is preferably recorded in advance in a memory accessible by the calculation unit 1.
- the function of the calculation unit 1 can be realized by a CPU of a single computer, or can be realized by a plurality of computers capable of data communication with each other.
- a program that causes a computer to execute the processing of the computing unit 1 and a non-transitrory that records the program ) is also included in the embodiment of the present invention.
- present invention is not limited to the first to fourth embodiments.
- the present invention can be applied to all types of analysis using differential equations such as structural analysis, vibration analysis, electromagnetic field analysis, material analysis, fluid analysis, temperature analysis, sound field analysis, electromagnetic analysis, circuit simulator, and the like.
- the differential equation of elasticity is expressed by equations (1) and (2). It becomes. Using these equations, the left side of (123) is It is. Using these stress components, if equation (819) is expressed, It becomes. From the partial integration formula (817), It becomes. The variational direct method of (123) means that this is zero. Get Both Eqs. (132) and (829) are none other than the variational direct method, but in the past they were called the virtual work principle.
- Shape Function Figure 8 in Section 8.1 shows the node arrangement of a 16-node quadrilateral element. Although this can be used for the h method, the node arrangement of a 4-node quadrilateral element, which is a more typical element, is shown in FIG.
- strain component is calculated from equations (4) to (6).
- Strain energy Strain energy U is expressed by equation (130). It is. This is used on the left side of the energy conservation law (837) and doubles to represent the left side of (838).
- the strain energy eng inside the element is Can be obtained at In the above equation, the energy U and the nodal displacement U are the same symbol, which is confusing, but the U on the left side represents the strain energy, and the two ⁇ U ⁇ appearing on the right side represent the nodal displacement.
- the integral term in the above equation It can be calculated by.
- Element stiffness matrix [K elm ] per unit thickness This represents the energy equivalence between the nodal displacement U and the equivalent nodal force F.
- Equation (865) After dividing the object to be analyzed into an enormous number of elements, formula (865) is created for each element, and the overall system equation is constructed while paying attention to the shared nodes. It becomes.
- N is the total number of degrees of freedom of the entire system. Equation (434) expresses the total degree of freedom N as n.
- the method of recognizing this as a balance of forces is the linear elastic finite element method by the h method, that is, the conventional method.
- Equation (865) If you take a close look at the procedure for obtaining Equation (865), you can see that the variational calculation is not used anywhere and no virtual displacement appears. Many textbooks explain that Eq. (829) is calculated and Eq. (865) is obtained, but in reality it is not, and the laws of energy conservation (837) and (838) are used. There is. Therefore, it makes more sense to understand that (866) and (434) represent energy equivalence rather than understanding force balance.
- the equivalent nodal force F per unit thickness is given by It is. Each is obtained by integrating the surface force on the element boundary.
- the equivalent nodal force F per unit thickness is given by equation (861). It is. This is obtained by placing the strain energy inside the element equal to the work by the equivalent nodal force.
- the stiffness matrix for the new finite element consists of two [K U ] and [K O ], whereas the stiffness matrix for the h-method element consists of one [K eng ]. From the details, that is, the calculation methods of (212), (213) and (853), the two new and old stiffness matrices are completely different in nature, and the effect of the correction function, ie, [K O ] can be ignored. For example, you can see that they are not the same.
- In-element displacement u j is expressed by equation (845). It is. As isoparametric element, it causes the interpolation function number chi Ai matches the shape functions N i. U ij is a nodal displacement, and the arrangement of various matrices is determined according to the nodal displacement. From equations (192) and (193), It is. Therefore, It becomes.
- equation (870) can be calculated, It becomes. here, obviously, It is. Also, It is.
- Equation (870) For the element stiffness matrix [K elm ] per unit thickness, substituting (876) into (864), It is common to define Using the above equation is a bit cumbersome because it is necessary to rewrite the matrix elements in the plane stress state and the plane strain state. However, it is necessary to show the formulas often seen in textbooks and to know the correspondence between them and formula (893), so the matrix elements of both are described below. As a result, you will find that (893), which does not need to distinguish between the two, is easier to handle.
- the nodal displacement U is calculated from the equations (934) and (938). If we define and, It is. Here, n is the number of degrees of freedom of one element. Since it is confusing when the same symbol n is used, replace the symbol as follows. If the total number of degrees of freedom when the entire structure is divided into meshes and expressed in terms of h-method elements is n T , Eq. (943) is It becomes. The above equation represents the displacement distribution within the element, and it can be seen that it is a function of the mode coefficient ⁇ a O ⁇ .
- formulation using conventional elements (442)
- the solution is obtained as follows. It is a solution that can handle both unique and multiple solutions, which is a structural feature. These items are common features consistent with both solutions.
- Figure 42 shows the deformation and stress distribution when a uniform gravity load fy is applied.
- the load coefficients c x and cy are The displacement is made dimensionless with the outer diameter R.
- Subscript B is added to a term that satisfies the inhomogeneous boundary condition, and u Bj is added, and a term that satisfies the homogeneous boundary condition is added to subscript H and u Hj is set.
- the main displacement u j is the sum, Get The inner product of this and the function u Hi * that satisfies the homogeneous adjoint boundary condition is It becomes.
- This consists of the functional system of (1110), and the combination of ⁇ x and ⁇ y that satisfies (1128) satisfies the above simultaneous differential equation.
- (1163) and (1164) are as follows: It becomes. This consists of the functional system of (1112), and the combination of ⁇ x and ⁇ y that satisfies (1128) satisfies the simultaneous differential equation of the above equation.
- the equations (1089) and (1090) of the dual simultaneous differential equations are Can be transformed.
- ⁇ 2 in (1100) Get Substituting this into the dual simultaneous differential equations (1089) and (1090), It becomes.
- This is composed of the functional system of (1110), and the combination of ⁇ x * and ⁇ y * that satisfies (1144) satisfies the above simultaneous differential equation.
- the solution ⁇ x * , ⁇ y * is expressed by equation (1111) (1195) and (1196) are as follows: It becomes. This indicates that the combination of ⁇ x * and ⁇ y * , which is composed of the function system of the formula (1112) and satisfies the formula (1144), satisfies the simultaneous differential equation of the above formula.
- the function ⁇ oe is a function that is vertically symmetric (symmetric with respect to the x axis) and antisymmetric with respect to the left and right (antisymmetric with respect to the y axis).
- the odd function (odd) is in the x direction
- the even function (even) is in the y direction.
- the functions are antisymmetric vertically (antisymmetric to the x axis) and antisymmetric to the left and right (antisymmetric to the y axis). That is, the odd function (odd) is in the x direction and the odd function (odd) is in the y direction. Take the initial letter and add the subscript oo.
- the function ⁇ eo is a function that is antisymmetric vertically (antisymmetric to the x axis) and symmetrical (symmetric to the y axis). That is, an even function (even) in the x direction and an odd function (odd) in the y direction. Take the initial letter and add the subscript eo.
- [AS] asymmetric form of the function phi xoo the x-direction displacements (1214) equation, when the displacement in the y direction (1210) represented in the form of a function phi yee of formula, variations in x-axis as 44 in FIG. (asymmetry), symmetry about the y-axis. Taking the initial letter, such a transformation is called mode AS. A function combination AS representing this deformation is expressed as Define with subscript AS.
- mode SS Taking the initial letter, such a transformation is called mode SS.
- the function combination ⁇ ss representing this deformation is expressed as As shown, it is defined with a subscript SS.
- Eigenvalues omega 1, represented as omega are collectively omega 2, (1108), a solution with a first kind modified Bessel function I m of (1112) equation, in accordance with (1209) equation, It expresses.
- the main variables of these functions are r and ⁇ , and n and ⁇ are auxiliary variables. If all these variables are written, the formula becomes long, so only the variables of interest are written below. For example, if r changes on both sides of the equation but r, ⁇ , and ⁇ do not change, the argument is expressed as (n). Then the derivative is Is obtained.
- the eigenvalues ⁇ 1 and ⁇ 2 are generically expressed as ⁇ , and a solution using the second type modified Bessel function K m of the equations (1108) and (1112) is expressed by the equation (1209), It expresses.
- the main variables of these functions are r and ⁇ , and n and ⁇ are auxiliary variables. If all these variables are written, the formula becomes long, so only the variables of interest are written below. For example, if r changes on both sides of the equation but r, ⁇ , and ⁇ do not change, the argument is expressed as (n). Then the derivative is Is obtained. Furthermore, from equations (1358) to (1361), Is obtained.
- the mode SA includes J 1SA in the expression (1250), J 2SA in the expression (1262), Y 1SA in the expression (1294), Y 2SA in the expression (1306), I 1SA in the expression (1338), and I 1SA in the expression (1350).
- the mode SS includes J 1SS in the expression (1252), J 2SS in the expression (1264), Y 1SS in the expression (1296), Y 2SS in the expression (1308), I 1SS in the expression (1340), and I 1SS in the expression (1352).
- a total of 8 sets belong to I 2SS , K 1SS of equation (1384), and K 2SS of equation (1396). Since it is a solution function for any integer n, the total number is 8n.
- (1) For J 1SS It becomes. Convert (1462) by (969) and (970), Get (2) For J2SS , It becomes. (1466) is converted by (969) and (970), Get (3) For Y 1SS , It becomes.
- the mode AA includes J 1AA in the formula (1253), J 2AA in the formula (1265), Y 1AA in the formula (1297), Y 2AA in the formula (1309), I 1AA in the formula (1341), I 1AA in the formula (1353) A total of 8 sets of I 2AA , K 1AA of the formula (1385), and K 2AA of the formula (1397) belong. Since it is a solution function for any integer n, the total number is 8n. (1) For J 1AA , It becomes. (1494) is transformed by (969) and (970), Get (2) For J2AA , It becomes. (1498) is transformed by (969) and (970), Get (3) For Y 1AA , Get (4) For Y 2AA , It becomes.
- the main eigenfunction, dual eigenfunction, and stress function based on them are collectively denoted as F.
- the coefficient c J1 is added to J 1AS in the formula (1251)
- the coefficient c J2 is set to J 2AS in the formula (1263)
- the coefficient c Y1 is set to the Y 1AS in the formula (1295)
- the coefficient is set to Y 2AS in the formula (1307).
- the c Y2 the coefficients c I1 to (1339) formula I 1AS, the coefficients c I1 to (1351) formula I 2AS, a coefficient c K1 to (1383) equation K 1AS, (1395) equation K 2AS Are multiplied by the coefficient c K2 and added together to make it equal to F.
- J 1AS is expressed by (1432) and (1433)
- J 2AS is expressed by (1436) and (1437)
- Y 1AS (1440) and (1441) Y 2AS in (1444) and (1445)
- I 1AS in (1448) and (1449) I 2AS in (1452) and (1453)
- K 1AS is expressed by (1456) and (1457)
- K 2AS is expressed by (1460) and (1461), respectively.
- the first and second columns of A are obviously, It is.
- the third and fourth columns of A are obviously, It is.
- the fifth and sixth columns of A are obviously, It is.
- the seventh and eighth columns of A are obviously, It is.
- the determinant of A is a function of ⁇ 1 and ⁇ 2 , but from the equations (1533) and (1534), it is eventually a function of the eigenvalue ⁇ .
- Eigenvalue ⁇ is obtained by solving the eigen equation with determinant zero.
- an eigenvector c F corresponding to the eigenvalue ⁇ is obtained from the equation (1531).
- the eigen equations and eigenvectors are very long and unfortunately cannot be described here. However, the solved results can be described numerically.
- the four eigenvalues ⁇ from the lower order are: In the vicinity. If the eigenvector c F corresponding to each eigenvalue ⁇ is obtained and substituted into the equation (1530), F is determined. F includes a main eigenfunction, a dual eigenfunction, and a stress function based on them.
- the upper row shows the main eigenfunction, and the surface force is zero and the displacement is zero on the inner side.
- the bottom is a dual eigenfunction, with zero surface force and zero displacement at the outer edge.
- the upper deformation is very similar to that of the analytical solution in Fig. 42, but the stress distribution is quite different. For example, while the distribution of the stresses sigma tau in the analysis solution has been seen, it is seen the distribution of similar stresses sigma theta is the main eigenfunction.
- the upper row shows the main eigenfunction, and the surface force is zero and the displacement is zero on the inner side.
- the bottom is a dual eigenfunction, with zero surface force and zero displacement at the outer edge. Paying attention to ⁇ ⁇ by the main eigenfunction, it can be seen that the distribution is very similar to the state of ⁇ ⁇ in the first mode. In addition, distribution of ⁇ ⁇ is it can also be seen, which is similar to the analytical solution. Therefore, it can be expected that the difference between the first mode and the second mode will be similar to the analytical solution.
- the state of deformation is somewhat complicated.
- the solution is obtained by the eigenfunction method in section 5.1. Like the analytical solution of 11.4.10 Section, weighting factors c x, a c y, The displacement is made dimensionless with the outer diameter R.
- Figure 52 shows the calculation results up to the third mode. It can be seen that the distribution of stress ⁇ ⁇ is approaching the analytical solution.
- Fig. 53 shows the calculation results using up to the 10th mode. It can be seen that the stress distribution is approaching the analytical solution.
- Figure 54 shows the calculation results using up to the 30th mode. The results are almost the same as the results up to the 10th mode, and the stress distribution and deformation are generally close to the analytical solution.
- Figure 55 shows the dimensionless velocity with the observation area as a square.
- Fig. 56 shows the dimensionless speed of the observation area as a square.
- FIGS. 57 to 64 show the first mode of the main eigenfunction and dual eigenfunction (1) to (8).
- residual flow is expressed by u Bj of the homogeneous solution part.
- the special solution part is expressed only by the eigenfunction F 1 11SS of the equation (1704), and the eigenvalue is
Abstract
Description
ここで,連立偏微分方程式で表された非自己随伴問題を解く為の,新しい固有関数法を示す。次に,本手法と最小2乗法,変分直接法,弾性学のエネルギ原理との関係について論じる。その後,仮想仕事の原理及び有限要素定式化の深部について論じる。さらに、2次元弾性体の有限要素の定式化について述べる。
2. 連立偏微分方程式
2.1 用語と造語
固有関数と微分作用素や境界条件との関係について,必要な用語や造語を図1に示す。
として,表の群に含まれる。主問題の解関数を主解[Primal Solution]と呼び,この変分
を主変分[Primal Variation]と呼ぶ。
する条件からは同次随伴境界条件[Homogeneous Adjoint Boundary Condition]が定まる。同じ同次随伴境界条件を持つ問題は多数あり,同様に加群を成す。これを裏の群[Tails Group]と呼ぶ。その要素のひとつを非同次随伴境界条件[Inhomogeneous Adjoint Boundary
Condition]と呼び,これを与える問題を双対問題[Dual Problem]と呼ぶ。双対問題の解
関数を双対解[Dual Solution]と呼び,この変分を双対変分[Dual Variation]と呼ぶ。内
積を通じて,両群は表裏一体に結ばれる。
2次元弾性体の静的釣合方程式は,直交座標系で,
平面応力,平面ひずみ状態共に,変位ひずみ関係は,
態を統一的に取り扱えるようになる。(1),(2)式の微分作用素を原初微分作用素Lijと
呼び,
Lij *=Lij
となる場合を自己随伴な微分作用素と呼ぶ。
2次元ポテンシャル流れの問題を,速度ポテンシャルを用いずに,速度だけで解こうとすれば,直交座標系で,
Lij * ≠Lij (22)
となる場合を非自己随伴な微分作用素と呼ぶ。
3.1 境界条件の同次化
以降は,前章の弾性体の例に題材を絞って記す。主解となる変位ux,uyを主変位と呼びu1,u2と書く。方程式の右辺をf1,f2と書けば,連立偏微分方程式(1),(2)
は,
双対解となる変位を双対変位uj *と呼ぶ。前節同様に,非同次随伴境界条件を満たす項uBj *と同次随伴境界条件を満たす項uHj *の和で,双対変位uj *を,
uHj *で表現されるから,最終的な未知関数はuHj *のみである。同様に,uBj *,pBj *は非同次随伴境界条件を満たす既知関数であるが,pBj *はuBj *で表現されるから,最終的な既知関数はuBj *のみである。
境界表面をc,内部領域をsとする。方程式(23)に双対変位ui *を乗じた積分和,即ち,内積は,
(59),(64)式右辺の作用素に着目して,双対問題の連立偏微分方程式
4.1関数群
同次境界条件を満足する関数群をφjとおき,この群のk番目の関数をφjkと書いて,解関数uHjをその和で,
応力を表す(27),(30),(33)式に(69)式を代入すれば,
る。
4.1節の知見により,同次化された微分方程式(40)に注目して,主連立微分方程式[Primal Simultaneous Differential Equations]
*,pEi *に課される為,境界項REはゼロとなる。よって,(88)式により,(90)式は,
λ*=λ (93)
となるようにλ,λ*を定めれば,
分係数の和ΣjLij *φj *が同次境界条件を満たす必要がある。これらを合わせて双対境界条件[Dual Boundary Condition]と呼ぶ。(主,双対)境界条件は共に自己随伴な境界条件である。
主連立微分方程式 (87)は,(93),(97)式より,
である。次章の固有関数法に用いる関数φi,φi *は,(主,双対)連立微分方程式(105),(106)を満たす組合せとする。ただし,φiはそのままで,λとφi *の符号を反転させた組合せも要件を満たすが,この組合せは排除する。
5.1固有関数法
5.1.1 主問題
同次化された微分方程式(40)に(69)式を代入して,
間での固有関数法である。
同次化された微分方程式(66)に(70)式を代入して,
となる。上式は,双対連立微分方程式(106)により,
となる。関数φiと内積をとって直交性を利用すれば,
方, (113)式に戻せば外力項fHi *が再構成され検証に役立つ。即ち,(115)式の解法
は,変位uHj *を双対固有関数φi *で,外力項fHi *を主固有関数φiで表現する,ヒルベルト空間での固有関数法である。
ここで改めて図1を熟視すれば,あたかもコインの如く表裏一体の構造は,原初微分作用素がその材質に,同次境界条件が表裏紋様の礎に当り,これに内在する固有値と(主,双対)固有関数が解法の主役を演ずる事が判る。
している。「双対」の名を冠した根拠は,上記の特徴「表裏一体・反転利用」にある。
主変分[Primal Variation]δuiは(24)式,双対変分[Dual Variation]δui *は(42)式より,それぞれ,
内積をとる手順を複数組み合わせる事は,(121)式の双対変分δui *と内積をとる事と
等価である。即ち,
致するから,(120),(121)式も同等になり,(122)式は,
一方,(105),(120)式により,主変分δuiの微分係数の和について,
物体力が指定される問題を扱う。(128)式より,基本的には全ての問題に対して,最小
2乗法が変分原理となる。
PEi=0 or φi=0 (129)
となる条件は自己随伴な境界条件である事が判る。即ち,表面力および変位のx,y方向
の成分について,各成分ごとにどちらか一方だけを与える問題は自己随伴問題となる。従って,この場合には,(123)式を用いても良い。単位厚さ当りの全ひずみエネルギUは,
分直接法として認識すべきである。境界表面から変位指定領域を取り除いた領域をcpとすれば,cp上で表面力は指定値となるから,上式は,汎関数Πを,
用した結果得られたものであるから,変分直接法の別形式として認識すべきである。境界表面から表面力指定領域を取り除いた領域をCuとすれば,Cu上で変位は指定値となるから,上式は,汎関数Πを,
(1) 連立偏微分方程式(23)は,(111)式の固有関数法によって解ける。(128)式より,最
小2乗法は変分原理である。非自己随伴問題では(122),(127)式が変分直接法となり,自
己随伴問題では(123),(127)式が変分直接法になる。
一形式の仮想仕事の原理では,δuiは「幾何学的境界条件を満たす,すべての関数」で
仮想変位と称される。力学的境界条件を満たすとは限らない関数に対して変分記号δを使う事は不適切であるから,仮想変位をδ~uiと書くと,仮想仕事の原理とは,
)式による恩恵である。
えばp法と同等の要素が得られ,(122),(127)式からは新たな有限要素を構成できる。
8.1 形状関数
例えば,16節点四角形要素の節点番号を図8のように表すと,その形状関数Ni(i=1~16)は,
内挿関数χAi(i=1~n~)を用いて要素内部の変位uAjを
試行関数ψok(k=1~l~)は,境界での変位がゼロとなる関数系を用いる。例え
ば,ルジャンドル陪関数でも良いし,三角関数でも良い。三角関数の場合には,関数系を
双対解となる双対変位uj *を,境界変位を表現する項と,境界で変位がゼロとなる項とに分けて,
ひずみ成分は,(4)~(6)式から,
境界表面での変位は(208),双対変位は(229)式であり,表面力は(207),双対表面力は(228)式である。これらを用いて,(60)式の境界項Rは,
のとき,境界項Rをゼロとする条件より,同次随伴境界条件が得られ,
界条件であり,一致しない場合には非自己随伴境界条件である。
等価節点力の定義式(216)を変形して,方程式
、
入れ替わって、
v]、[-I]から該当列を抽出すれば,行列[Kv]が得られる.また,行列[Yb]、
[Yv]は、次の様にして得られる.単位行列[I]とゼロ行列[0]を用いて,(257)式の行列[χA]、[0]を
に,行列[0]から列を抽出する操作は[0]の列を抜き出すことで表現できる.既知部Ubと未知部Uvの並びに従って、[I]から該当列を抽出して[IUb]、[IUv]を構成する。同様に、既知部Fbと未知部Fvの並びに従って,「0」から該当列を抽出して[0Fb]、[0Fv]を構成する。
である。そして、未知部sv、coと既知部sbの間には、
、[K0]の列を入れ替えて[Ke]、[Kh]とする。同様に,行列[Yv]、[ψ0]の列を入れ替えて[Ye]、[Yh]とする.その結果,(263)式は,
(1)未知部svのうち未知部seに入るものnve個を、未知部seの前半に集める。
(2)未知部coのうち未知部seに入るものnoe個を、未知部seの後半に集める。
(3)未知部svのうち未知部shに入るものnvh個を、未知部shの前半に集める。
(4)未知部coのうち未知部shに入るものnoh個を、未知部shの後半に集める。
よって、未知部の個数は、
cv]の列を抜き出すことで表現できる。[cv]から該当列を抽出して[cve]、[cvh]を構成して、
表現できる.行列[I]から該当列を抽出して[Ioe]、[Ioh]を構成して、
位[uB]と、同次境界条件を満足する変位[uH]とで表現された事が判る。
(310)、(313)式により、行列[χB]、[ψh]が関数部[Γ]と係数部[Hb]、[Hh]とに分離された事が判る。また、変位{u}が、
され,任意の未知部shによって同次境界条件を満たす関数が構成されることを知る。未
知部shに対して変位[uH]が同次境界条件を満たすのは、行列[ψh]による効果である
から,係数ベクトルを{eh}として,主固有関数{φ}を、
等価節点力の定義式(235)を変形して,方程式
入れ替わって、
I]から該当列を抽出して得られるもので,同様に、未知部Uv *,Fv *の並びに従って行列[Kv]、[-I]から該当列を抽出すれば,行列[Kv *]が得られる.また,行列[
Yb *]、[Yv *]は、次の様にして得られる。単位行列[I]とゼロ行列[0]を用いて,(318)式の行列[χA]、[0]を
に,行列[0]から列を抽出する操作は[0]の列を抜き出すことで表現できる.既知部Ub *と未知部Uv *の並びに従って、[I]から該当列を抽出して[IUb*]、[IUv*]を構成する。同様に、既知部Fb*と未知部Fv*の並びに従って,「0」から該当列を抽出し
て[0Fb*]、[0Fv*]を構成する。
れる状態である。そして、未知部sv *、co *と既知部sb *の間には、
(1)未知部sv *のうち未知部se *に入るものnve *個を、未知部se *の前半に集める。
(2)未知部co *のうち未知部se *に入るものnoe *個を、未知部se *の後半に集める。
(3)未知部sv *のうち未知部sh *に入るものnvh *個を、未知部sh *の前半に集める。
(4)未知部co *のうち未知部sh *に入るものnoh *個を、未知部sh *の後半に集める。
よって、未知部の個数は、
[cv *]の列を抜き出すことで表現できる。[cv *]から該当列を抽出して[cve *]、
[cvh *]を構成して、
表現できる.行列[I]から該当列を抽出して[Ioe *]、[Ioh *]を構成して、
uB *]となる.第2項は未知部sh *を備えた関数であり,同次随伴境界条件を満足する変
位[uH *]となる。即ち、
(370)、(373)式により、行列[χB *]、[ψh *]が関数部[Γ]と係数部[Hb *]、[Hh *]とに分離された事が判る。また、変位{u*}が、
主連立微分方程式(105)、双対連立微分方程式(106)を行列で表現して、それぞれ、
となる。(315),(375)式の固有関数{φ}、{φ*}及び(316)、(376)式の変分{δφ}、{δφ*}を代入して、
ができる。逆に,(398),(399)式をひとつにまとめると、
る。[As]、[Bs]は共に対称行列であるから、固有値λは実数となり,固有ベクトル{es}は[As]、[Bs]を挟んで直交する。
クトル{eh}、{eh *}の組み合わせについて,注意深く取り扱う必要が生じる.一方
,(402)式で計算すれば,固有値問題のサイズは大きくなるが,固有ベクトル{eh}、{eh *}の組み合わせも同時に得られるので、組み合わせの確認は不要となる。どちらの方法を用いるかは、メモリサイズや計算速度に応じて決めれば良いが、定式化の都合からは、(402)式の方が扱い易いので、以降は、これを考え方の基準にする。
す事によるもので、4.4 節で述べた事でもある。[As]がゼロ空間を持つ場合には,固
有値がゼロとなる。[Ah]のゼロ空間をm0 個とする。
戻すことにより,固有関数{φ}、{φ*}が得られる。
(390)式の[Ah]を、(293),(313),(373)式により変形すれば、
形すれば,
、(373)式により変形すれば、
主問題の連立偏微分方程式(23)を行列で表現して、
{ap}、{a0}と定義して、(406),(415)式を利用すれば、
している。
(422)、(423)式より、連立偏微分方程式は、
{a0}の関数となるので,節点情報の他に,辺上の境界条件を加えることで{a0}を定めることができる。
応力分布、ひずみエネルギ分布なども得られる。
式を意味しているので,その解に付随の係数{a0}は任意である。前節では、その任意
性を辺上の境界条件を付加することで決定する方法を示したが,計算を行う設計者が自ら決る方法もある。係数{a0}に応じてどのような変形状態,あるいは応力状態となるか
を図で表示するシステムを具備して,この図を見ながら実現可能な境界条件となるように,設計者が自ら係数{a0}を決めるのである。
9.1 全体系の方程式
7 章の(4)項にて,従来の線形弾性有限要素,特にh法要素では変分の演算は行われて
おらず,(134)式のエネルギ保存則と一致することを述べた。h法要素を複数個組み合わ
せて,全体系の運動方程式は、
が変位指定領域、n2個が外力指定領域で,
有限要素法の常套手段として確立されている。
れていることに着目すれば、何も自己随伴な境界条件に縛られる必要はなく、非自己随伴な境界条件に対して上式を解いても、良いはずである。
全体系の運動方程式(434)を,単位行列[I]を用いて,
を
から該当列を抽出して得られるもので、同様に、未知部Uv、Fvの並びに従って行列[K]、[-I]から該当列を抽出すれば,行列[Kv]が得られる。上式からは、
できる。即ち、行列[Kv]が、解の複数存在する劣決定系になるのか、解が存在しない
可能性のある優決定系になるのかに応じて、線形代数学の知見に基づき、未知部svを得
ることができる。解が複数存在する場合には、8.18節のゼロ固有値の場合に相当するので、設計者自身によって解を決定することができる。
全体系の運動方程式(434)の節点外力Fを、重力荷重など物体力を表す項Fgと表面力を表す項Fsとに分けて、
解が複数存在する場合には,(452)式を満足する1つの特解{sp}と、(452)式の右辺
をゼロとおいた同次解{s0}が得られる。同次解の個数n0は、
がった状態となって,不正な解析結果と見なされることが良くある。これは、全体系として,エネルギ保存則と等価な方程式(434)を満たすものの、要素内部で微分方程式を満足
しないことに起因する。
の微分方程式に等価な変分原理、即ち、(122)式~(128)式の最小自乗法を用いることで、適切な{a0}を決める.適切な{a0}の周辺の値で要素が折れ曲がった状態にならない場合には、まずまず正確な解析結果が得られていると考えることができる。
従来型の有限要素では,要素内部の変位は(166)式の{uA}で表現される。{uA}は
節点変位Uで表され,Uは(455)式のようにモード係数{a0}によって表される。
数を調節するのに合わせて、変形図が表示されるように構成すれば、どのような支持方法があり得るのか検討できて、設計者のイメージを湧き立たせることが可能となる。実現可能な変形状態を絞り込んで境界条件を再定義することで、より良い構造を得るための検討が可能となる。
図2は、本実施形態の解析装置の構成の一例を示す機能ブロック図である。図2に示す解析装置10(解析システム10)は、演算部1を備えるコンピュータにより構成することができる。演算部1は、解析対象物の設計データ、境界条件データ及び解析対象に対する微分作用素を含む微分方程式データを入力し、解析対象の微分方程式の解を計算して、解析結果データとして出力です。ここで、演算部1は、双対変位ui *の変分のことを双対変分δui *とするとき,下記式によって解ujを計算する。下記式は、上述した式(12
2)と同じである。
とすることができるが、fi及びuiは、これに限られない。解析対象及び微分作用素に応じて、fi及びuiは力及び変位以外のものを示す場合がありうる。例えば、構造解析、振動解析、電磁場解析、材料解析、流体解析、温度解析、音場解析、電磁気学解析、又は回路シミュレータ等の様々な解析において、解析対象に応じた適切な値がfi及びuiに設定することができる。
演算部1は、下記式によって解ujを計算してもよい。下記式は上述した式(123)
と同様である。この場合、自己随伴問題の解を計算することができる。
演算部1は、変分直接法により下記式において解析対象の解ujを計算することもでき
る。下記式は、上記(127)と同様である。
解を計算することができる。
図3は、第2の実施形態における解析装置の構成の一例を示す機能ブロック図である。図3に示す解析装置10aは、設定部11、随伴境界条件演算部12、及び非自己随伴演算部13を備える。
位もゼロ(φi=0)の境界条件が課されていた場合、pEi *は任意となる。この場合、随伴境界条件は、境界表面の変位、力とも任意となる。このような場合、境界条件と随伴境界条件が一致しないため、非自己随伴問題となるただし、上記式(106)の作用により,Lji *φj *がφiに課された境界条件を満たさねばならない。これは前記双対境界条件に含まれる。非自己随伴問題の場合には、前記主境界条件と前記双対境界条件とを利用することにより、適切に問題が解ける。
方程式(上記式(105))、双対連立微分方程式(上記式(106))、境界条件、随伴境界条件を用いて、主連立固有値問題(上記式(98))と双対連立固有値問題(上記式(99))を定義して、主固有関数φ、及び双対固有関数φ*を計算することができる
。
る。その場合、主固有関数φ、及び双対固有関数φ*は、解析的にも近似的にも計算する
ことができる。例えば、後述するように、上記式(122)及び式(123)によって近似固有関数を計算することができる。
立方程式の解を計算する。例えば、上記式(111)で示される係数ckを算出し、これ
を上記式(69)に代入して、解関数uHjを計算することができる。また、非自己随伴演算部13は、上記式(115)により双対問題の係数ck *を計算でき、これを上記(70)に代入することで双対問題の解関数uHj *を計算することができる。
図4は、図3に示す解析装置10aの動作例を示すフローチャートである。図4に示す例では、まず、設定部11は、ユーザからの設計データの入力を受け付ける(S1)。設計データは特定の対象や形式に限定されない。例えば、解析対象が固体の場合は、解析対象の形状や材料を表すデータの入力を受け付けることができる。解析対象が流体の場合は
、流体の物性を表すデータの入力を受け付けることができる。また、有限要素法による解析の場合は、要素ごとの構成や物性を示すデータの入力を受け付けてもよい。また、設定部11は、ユーザからの入力のみならず、例えば、アクセス可能な記録媒体から設計データを読み込む、あるいはネットワークを介するダウンロード等により、設計データを取得することができる。
て、固有関数法により、主問題の連立方程式の解関数uHjを、例えば、上記式(111)及び式(69)により、計算する(S7)。また、非自己随伴演算部13は、主固有関数φ、及び双対固有関数φ*を用いて、双対問題の連立方程式の解関数uHj *を、例えば、上記式(115)及び式(70)により、計算することができる。このようにして、計算された解は、解析結果をコンピュータのディプレイ(図示せず)などに出力される(S8)。解析結果の具体例については、後述する。
図5は、第3の実施形態における解析装置の構成の一例を示す機能ブロック図である。図5に示す解析装置10bは、設定部11、随伴境界条件演算部12、自己随伴判定部15、非自己随伴演算部13、及び、自己随伴演算部14を備える。設定部11、随伴境界条件演算部12、及び非自己随伴演算部13は、上記第1の実施形態の同様に構成することができる。
対固有関数φi *または双対変位ui *を乗じた積分和すなわち内積を、部分積分することにより計算することができる。例えば、上記式(62)又は式(90)により、随伴微分作用素Lij *を計算することができる。
随伴演算部14は、例えば、自己随伴固有関数φikを上記式(111)のφik *と置き換
えた式により、ckを計算し、このckを上記式(69)に代入することにより、自己随伴問題の解関数uHjを求めることができる。
図6は、図5に示す解析装置10bの動作例を示すフローチャートである。図6において、解析装置10bは、S1~S4の処理を、上記第2の実施形態におけるS1~4の処理と同様に実行することができる。S21において、自己随伴判定部15は、S1で取得された境界条件を示すデータと、S4で生成された随伴境界条件データとを比較することにより、境界条件が随伴境界条件と一致するか否かを判定する。S21で自己随伴判定部15が、一致しないと(すなわち、非自己随伴である)と判定すると、非自己随伴演算部13は、S5~S7の処理を実行する。S5~S7の処理は、上記第2の実施形態におけるS4~S7と同様にしてもよい。
随伴固有関数φiを用いて、自己随伴連立微分方程式の解関数ujを計算する(S24)
。S8の結果出力は、上記第2の実施形態におけるS8と同様にしてもよい。
本実施の形態は、解析装置を、有限要素法(FEM)を用いた解析に適用した場合の例である。本実施形態では、原初微分作用素及び境界条件は、有限要素に対して設定される。例えば、上記の第1~第3の実施形態で説明した解析装置の構成及び解析処理を、有限要素に対しても、同様に、適用することができる。
ば、w=1/r2とするのが好ましい。ここで、代表長さには、物体のどこか特徴的な部分の長さが採用される。円環の場合には、半径を代表長さにとるのが好ましい。この場合、相似図形の固有値λが一致するので、相似図形の物体については、計算を節約することができる。
u0は、試行関数ψ0の和で表すことができる(例えば、上記式(173)及び式(190))。設定部11は試行関数を設定する(S34)。試行関数は、例えば、上記式(179)のように設定することができる。例えば、図8に示すような16節点四角形の要素の場合、境界表面での変位は上記式(208)、双対変位は上記式(229)、表面力は上記式(207)、双対表面力は上記式(228)で表すことができる。
に、双対変位U*も既知部Ub *と未知部Uv *に分けてならべかえると式(251)のよう
になる。Rをゼロとするには、未知部Fv、Uvに対応するUb *、Fb *がゼロでなければならないので、式(254)の随伴境界条件が得られる。
(S37、S38)、近似主固有関数を{φ}、近似双対固有関数を{φ*}、それぞれ
の係数ベクトルを{eh}、{eh *}とおいて、
{φ}=[ψh]{eh}
{φ*}=[ψh *]{eh *}
と表し、これらを主連立微分方程式(上記式(105))及び双対連立微分方程式(上記式(106))へ代入して、上記式(122)及び式(123)によって係数ベクトル{eh}および{eh *}を得た後、近似固有関数{φ}および{φ*}を得ることができる(S39)。
場合、関数行列Γを上記式(293)で定義し、係数行列[ce][ch]を上記式(294)(295)で定義している。双対試行関数[ψh *]は、例えば、上記式(373)のように表すことができる。この場合、係数行列[ce *][ch *]を上記式(354)(355)で定義している。
図67に示す情報処理装置10cは、初期方程式決定部21、境界条件決定部22、及び演算部23を備える。初期方程式決定部21は、処理対象となる系の構造及び構成要素の性質を示すデータを読み込んで、系を表す方程式を決定する。この方程式は、例えば、求めるべき物理量を表す変数を含むn個の方程式とすることができる。方程式は、例えば、行列のデータ又は微分方程式のデータで表すことができ、データの表現形式は特に限定されない。変数は、例えば、ベクトル又は関数等のデータで表すことができるが、表現形式に特に限定されない。
初期方程式決定部21は、例えば、前記系の構成要素である物体に作用する力と、前記物体の物理量を表す変数との関係を表す連立微分方程式を決定することができる。物体の物理量は、例えば、変位又は速度等とすることができるが、特定のものに限定されない。
また、演算部23は、随伴作用素を計算せずに、初期方程式を2n個の変数を含む方程式に変形し、方程式の解を計算することができる。例えば、初期方程式決定部21が、2つのn次元の変数ベクトルとn行の行列を決定することで、系の構成要素の節点における物理量を示す初期方程式を決定することができる。例えば、節点の変位を表す変数のn次元ベクトルと、節点の力を表す変数のn次元ベクトルとの関係を表すn行の行列を生成することができる。これは、系が満たすべき微分方程式を、n次元ベクトルとn行の行列を用いて表したものと言える。例えば、式(441)や式(449)等に示すような形で初期方程式を決定することができる。式(441)の例では、n行n列の行列Kと、n次元ベクトルの変数ベクトル{U},{F}でn個の方程式が表される。
境界条件決定部22は、方程式の変数において、値が既知となる変数の自由度の数と、値が未知となる変数の自由度の数が異なる境界条件も決定可能に構成される。例えば、ユーザから、境界条件として、1自由度に対して複数の物理量が未知または既知であると指定された場合も、ユーザの設定が不適切である旨等のエラーは出さないよう構成することができる。
上記演算部23は、上記処理を実行することにより自己随伴な問題と、非自己随伴な問題の両方について解を計算することができる。例えば、1つの自由度、例えば、ある節点のx方向については変位も力もゼロであるという境界条件が設定された場合や、ある節点のy方向については変位も力も未知であるといった境界条件も受け入れることができる。上記のような、変数をn倍にした方程式で、未知部と既知部の自由度の数に制限をかけない演算により、非自己随伴境界条件の下でも解を計算することが可能になる。
本実施形態では、上記演算部23では、自由度を2倍にした方程式を、未知部と既知部の自由度が必ずしも同じでない状況で方程式の解を計算する。この演算部23では、複数の解が算出される場合がある。その場合にも、上記例のように、複数解のうち代表的な解を計算し、この解を含む所定の範囲で解を出力することで適切な結果を出力することができる。また、ユーザの操作により解を変動させて表示することで、複数解の結果を認識可能な形態で伝えることができるユーザインタフェースが提供される。
図68に示す情報処理装置10dは、初期方程式決定部21、境界条件決定部22、及び判定部33を備える。初期方程式決定部21及び境界条件決定部22は、図67の初期方程式決定部21及び境界条件決定部22と同様に構成することができる。判定部33は、初期方程式決定部21で決定された連立微分方程式の境界条件が、自己随伴条件か否かを判定する。
判定部33は、例えば、連立微分方程式の変数と同じ数の双対変数と、双対連立微分方程式とを定義した場合に、連立微分方程式と双対変数との積分和すなわち内積を部分積分して得られる境界項を計算し、境界条件決定部22で決定された境界条件の下で、当該境界項をゼロにする双対変数の条件を、随伴境界条件として決定することができる。判定部33は、随伴境界条件が境界条件と一致する場合は、自己随伴条件であると判断できる。
なお、上記例のように、随伴境界条件を計算しなくても、自己随伴条件の判定ができる場合がある。例えば、11.7.1に記載されているように、上記の式(251)において、自己随伴境界条件となる場合、節点変位Uの既知部Ubの個数n1は、節点外力の未知部Fvに等しく、節点外力Fの既知部Fbの個数n2は、節点変位の未知部Uvの個数に等しくなる。したがって、判定部33は、1つの変数の未知部の自由度が、他の変数の既知部の自由度と等しいか否かによって、自己随伴条件か否かを判定することができる。
なお、上記図68の判定部33を図68の情報処理装置に追加した構成も可能である。例えば、判定部33の判定結果に基づいて、演算部23が、2n次元の方程式を用いた演算を行うか否かを制御することができる。これにより、演算部23は、2n次元の方程式を用いた非自己随伴の演算と、n次元の方程式を用いた自己随伴の演算とを適切に切り替えて、処理効率を向上させることができる。
(自己随伴条件に対応する設定例)
図69は、従来の有限要素を用いた解析において、境界条件が自己随伴条件を満たすことを前提として、節点変位ベクトル、節点外力ベクトルを設定する処理例を示すフローチャートである。図69に示す例では、情報処理装置は、まず、節点変位の入力を受け付ける(S101)。例えば、変位がわかっている自由度とその数値がユーザから入力される。あるいは、節点変位の自由度と数値を示すデータが読み込まれてもよい。自由度は、例えば、節点の位置と方向により特定される。節点変位が入力されると、節点変位の既知部Ubの自由度と、その数値が確定する(S102)。節点変位の既知部Ubの自由度及び数値が確定すると自動的に節点外力の未知部Fvの自由度を確定することができる(S103)。自己随伴条件を満たす場合、未知部Fvの自由度は、既知部Ubに等しいからである。
図70は、有限要素を用いた解析において、境界条件が非自己随伴条件である場合にも対応可能な節点変位ベクトル、節点外力ベクトルを設定する処理例を示すフローチャートである。図70に示す例では、情報処理装置の境界条件決定部22は、節点変位の自由度及び数値の入力を受け付け(S201)、節点変位の既知部Ubの自由度と数値を決定し記録する(S202)。既知部Ubが決まると節点外力の未知部Fv(S203)、節点変位の未知部Uv(S204)、及び既知部Fb(S205)が自動的に計算される。S205では、値が指定されていない節点外力の自由度にはゼロが設定される。
図71は、境界条件が非自己随伴条件である場合にも対応可能な節点変位ベクトル、節点外力ベクトルを設定する処理の他の例を示すフローチャートである。図71に示す例では、節点変位および節点外力の自由度と数値の入力を受け付け、入力がなかった自由度は未知部として設定する。境界条件決定部22は、節点変位の自由度と数値の入力を受け付け(S301)、入力に基づいて、節点変位の既知部Ubの自由度と数値を確定し(S302)、これに伴って未知部Uvの自由度も確定する(S303)。同様にして、節点外力の自由と数値の入力を受け付け(S304)、入力に基づいて、節点外力の既知部Fbの自由度と数値を確定し(S305)、これに伴って未知部Fvの自由度も確定する(S306)。判定処理(S307)、表示処理(S308)、解析処理(S309、310)は、図70の場合と同様に実行することができる。
図72は、有限要素を用いた、非自己随伴問題の解析の実行例を示すフローチャートである。図72に示す例では、まず、情報処理装置10cは、計算モデルを入力する(S401)。例えば、処理対象となる構造の有限要素モデル、形状、材料、要素、節点等が入力される。
図9は、一様重力下の円環で,内辺の境界条件を「変位がゼロかつ表面力もゼロ」とするには,外辺を如何に支持すべきかを求めた結果を示す図である。従来は求解不能とされてきた問題であるが,新たに解析解を得て,これと固有関数法による計算結果とを比較し,解析解に近い事を確認した。円環の上部と下部に大きな垂直応力σrが作用し,左右部
には剪断応力τrθが作用する事で,重力荷重を受け持っている様子が判る。内辺は真円
に保たれ,表面力も作用していない事が確認できる。なお、図9において、左図は応力σr、中央図は応力σθ、右図は剪断応力τrθを示す。図9に示す結果は、円環を有限要素とせず、解析的に計算した結果である。以下に、解析に用いた具体的な主連立微分方程式の一例を示す。
σyが作用することによって上辺の変位がゼロに保たれている。ここでは、正方形を有限
要素として定式化しており、要素はひとつだけである。ひとつの要素で解析できる方法としてP法有限要素がある。図10に示す解析は、非自己随伴な境界条件の例であり、P法有限要素法ではこのような境界条件での解析はできない。また、境界条件を辺に与えるか、節点に与えるかは定式化の都合によるが、ここでは、節点に境界条件を与えるよう定式化している。
3Cは、固有関数法により求めた解(1モード、2モード及び3モード)および解析解(Analytical Solution)を示す図である。3モードの解は、解析解とほぼ重なっている。
とを示している.モード数が多いほど,固有関数法の精度は良くなり,10 モード程度で
変位 u^は良く一致していることがわかる。
たことを示している。モード数が多いほど、固有関数法の精度は良くなり、10 モード程
度で変位 u^は良く一致していることがわかる。
算したことを示している。モード数が多いほど,固有関数法の精度は良くなり、10 モー
ド程度で変位 u^は良く一致していることがわかる.
していない。
図中、破線が解析解を、実線が固有関数法を示している。56M は56 モードまで計算した
ことを示している。モード数が多いほど、固有関数法の精度は良くなり、10 モード程度
で変位 u^は良く一致していることがわかる。
たことを示している。モード数が多いほど、固有関数法の精度は良くなり、10 モード程
度で変位 u^は良く一致していることがわかる。
たことを示している。モード数が多いほど、固有関数法の精度は良くなり、3 モード程度で変位 u^は良く一致していることがわかる。
たことを示している。モード数が多いほど、固有関数法の精度は良くなり、3 モード程度で変位 u^は良く一致していることがわかる。
たことを示している。モード数が多いほど、固有関数法の精度は良くなり、3 モード程度で変位 u^は良く一致していることがわかる。
たことを示している。モード数が多いほど、固有関数法の精度は良くなり、3 モード程度で変位 u^は良く一致していることがわかる。
10.4.1 微分方程式
図32の梁の静たわみにて,梁の長さをl、曲げ剛性をEIとする。剪断力F(x)を、曲げモーメントをM(x)とする。変位をu(x)、分布荷重を f(x)とすると、
作用素Lを
固有関数法の精度は良くなり、3 モード程度で変位 u^ は良く一致していることがわか
る。
固有関数法の精度は良くなり、3 モード程度で変位 u^は良く一致していることがわかる。
線の向きを図32の方向に合わせるため,各関数値には負号を乗じて作図している。縦軸の座標値のみ、正負反転して読まれたし。図中、赤色破線が解析解を、黒色実線が固有関数法を示している。32M は32 モードまで計算したことを示している。モード数が多いほ
ど、固有関数法の精度は良くなり、3 モード程度で変位 u^は良く一致していることがわかる。
固有関数法の精度は良くなり、3 モード程度で変位 u^は良く一致していることがわかる。
図39は、図10に示した解析において、固有値がゼロになることにより複数の解が発生する場合の、複数の解のうちの図10に示す解以外の解の一つを示す図である。このような複数の解を設計者に提示することにより、設計者は、設計的に実現可能な方法を検討することができる。
上記第1~4の実施形態における演算部1、設定部11、随伴境界条件演算部12、非自己随伴演算部13、自己随伴判定部15、自己随伴演算部14は、コンピュータのプロセッサがメモリから、プログラム及びデータを読み込んでプログラムに従ってデータの処理を実行する。演算部1は、設計データ、境界条件データを、入力デバイス及びユーザインターフェースを介して、ユーザからの入力を受け付けることにより取得してもよいし、メモリから読み込むことで取得してもよい。微分方程式データは、例えば、演算部1がアクセス可能なメモリに予め記録しておくことが好ましい。
)な記録媒体も本願発明の実施形態に含まれる。
11.1 節5.2の「変分直接法と最小2乗法」における式変形の詳細
11.1.1 主問題の変分直接法
同次化された微分方程式(40)式
以上のことから,(795)式を変形した式
主連立微分方程式(105)
同次化された微分方程式(66)式
双対連立微分方程式(106)
11.2.1 部分積分公式
2次元領域に直交座標x,yをとり,物体の境界表面をC,内部領域をSとする.境界表面の外向き単位法線ベクトルの成分をnx,nyとする.任意の関数f,f*をx,yの関数として部分積分すると,
ひずみエネルギUは(130)式で表され,
有限要素法の一般的教科書では,(829)式を仮想仕事の原理と認識して要素の剛性行列が構成されるように記されているが,7章(4)で述べたように,実際のところは,エネルギ保存則を表す(838)式が使われている.次章にて,詳しく述べる.
11.3.1 形状関数
8.1節の図8では,16節点四角形要素の節点配置を示した.h法にこれを用いることもできるが,より代表的な要素である4節点四角形要素の節点配置を図41に示す.その形状関数Ni(i=1~4)は,
8.3節の(163)式では,変位ujを行列で,
外力のなす仕事wrkは,
エネルギ保存則(837)式に従い,(854)式と(859)式を等置して,
新しい有限要素では,単位厚さあたりの等価節点力Fは,(214)式で,
これに対して,従来の有限要素法即ちh法の要素では,
単位厚さあたりの等価節点力Fは,(861)式で,
図41の節点配置で,節点1,2,3,4の座標をそれぞれ(-a,-b),(a,-b),(a,b),(-a,b)として,アスペクト比κが,
9.5節で示した変分原理の計算方法を補足説明する.まずは,(251)式と同様に,節点変位Uを既知部Ubと未知部Uvに分けて,
同次解に対するモード係数を{aO}とおいて,解{sv}は(455)式で,
連立偏微分方程式(23)は,
本節は,10章に含める方が良い.
11.4.1 微分方程式
2次元弾性体の静的釣合方程式は,直交座標系で,
微分方程式(1),(2)の連立をはずして,重調和方程式
極座標系の微分方程式(996),(997)の一般解をψ(r,θ)とおいて,
(1000)式を(999)式に代入して,
なお,上式ではm=0の場合に不都合が生じるので,
さらに,(1006)式ではm=1の場合にも不都合が生じるので,
以上のことから,(1006),(1011),(1016)式を微分方程式(1001)の解の形式とすれば良い.
2次元円環の外径をR ,内径をγR とする.内径比率γの範囲は,
内辺の変位をゼロとする境界条件は,
微分方程式(1),(2)と境界条件式(1018),(1019)を満たす解では,(1016)式の形式は不要となるので,微分方程式(996),(997)の一般解ux,uy を,それぞれ(1006),(1011)式で表す.
計算の結果,関数を,
・|i=1,j=1| のとき,(1029),(1030)式を利用して,
・|i=1,j=2| のとき,(1031)式を利用して,
・ |i=2,j=1|のとき,(1032)式を利用して,
・ |i=2,j=2|のとき,(1029)式,(1030)式を利用して,
境界項Rをゼロとする条件から,随伴境界条件が定まる.
以下にその例を示すが,円環の外辺の半径をr=R ,内辺の半径をr=γRとしており,境界項Rと同じ記号が現れるのでまぎらわしいが,両者の意味合いが大きく異なるので判別は可能である.
円環にて,内辺(r=γR )で表面力がゼロかつ変位もゼロとする境界条件は,(1018),(1019)式であり,
円環にて,内辺(r=γR)で変位がゼロ,外辺(r=R )でも変位がゼロとする境界条件は,
非同次境界条件を満たす項に添字Bを付けてuBjとおき,同次境界条件を満たす項に添字Hを付けてuHjとおく.主変位ujをその和で,
・|i=1,j=1|のとき,(1029),(1030)式を利用して,
・ |i=1,j=2|のとき,(1031)式を利用して,
・ |i=2,j=1|のとき,(1032)式を利用して,
・ |i=2,j=2|のとき,(1029),(1030)式を利用して,
この例のように,
代表長さを外径Rとして,重み定数を
主連立固有値問題は,(98)式
(1094)式と(1098)式は同じ微分方程式であるから,結局,φx,φyとφx *,φy *は,同じ関数セットから構成されることがわかる.よって,(1094)式を解く事に照準を合わせて,
さて,θの関数には 2πの周期性を期待してmを整数とおき,変数分離解を求めることにして,微分方程式
同様にして,微分方程式
主連立固有値問題の(1092),(1093)式を変形して,
双対連立固有値問題の(1096),(1097)式を変形して,
一方,双対連立固有値問題の(1096),(1097)式は,
主連立微分方程式の(1087),(1088)式を変形して,
解 φx *,φy *を,(1105)式
一方,主連立微分方程式の(1087),(1088)式は,
双対連立微分方程式の(1089),(1090)式を変形して,
を満たさねばならないから,(1179),(1180)は,
解φx *,φy *を,(1107)式
一方,双対連立微分方程式の(1089),(1090)式は,
解φx *,φy *を,(1111)式
(1106),(1108),(1110),(1112)式の解を用いるにあたって,整数mを偶数と奇数の場合に分けるのが便利で、nを整数として、
x方向の変位が(1210)式の形の関数φxeeで,y方向の変位が(1214)式の形の関数φyooで表されるとき,変形は図43の様にx軸に対称(symmetry),y軸に非対称(asymmetry)になる.頭文字をとって,このような変形をモードSAと呼ぶ.この変形を表す関数の組み合わせφSAを、
x方向の変位が(1214)式の形の関数φxooで、y方向の変位が(1210)式の形の関数φyeeで表されるとき,変形は図の44の様にx軸に非対称(asymmetry),y軸に対称(symmetry)になる.頭文字をとって,このような変形をモードASと呼ぶ.この変形を表す関数の組み合わせφASを,
x方向の変位が(1212)式の形の関数φxoeで,y方向の変位が(1216)式の形の関数φyeoで表されるとき,変形は図45の様にx軸に対称(symmetry),y軸に対称(symmetry)になる.頭文字をとって,このような変形をモードSSと呼ぶ.この変形を表す関数の組み合わせφssを,
x方向の変位が(1216)式の形の関数φxeoで,y方向の変位が(1212)式の形の関数φyoeで表されるとき,変形は図46の様にx軸に非対称(asymmetry),y軸に非対称(asymmetry)になる.頭文字をとって,このような変形をモードAAと呼ぶ.この変形を表す関数の組み合わせφAAを,
固有値ω1、ω2を総称してωと表し,(1106),(1110)式の第1種ベッセル関数Jmを用いた解を,(1209)式に従い,
固有値ω1、ω2を総称してωと表し、(1106),(1110)式の第2種ベッセル関数Ymを用いた解を、(1209)式に従い、
モードSAには,(1250)式のJ1SA,(1262)式のJ2SA,(1294)式のY1SA,(1306)式のY2SA,(1338)式のI1SA,(1350)式のI2SA,(1382)式のK1SA,(1394)式のK2SA,の合計8セットが属する。任意の整数nに対しても解関数となるから,総数は8nとなる。
(4)Y2SAについては
(5) I1SAについては,
(6) I2SAについては,
(7) I1SAについては,
(8) I2SAについては,
モードASには,(1251)式のJ1AS,(1263)式のJ2AS ,(1295)式のY1AS ,(1307)式のY2AS ,(1339)式のI1AS,(1351)式のI2AS,(1383)式のK1AS,(1395)式のK2AS ,の合計8セットが属する.任意の整数nに対しても解関数となるから,総数は 8nとなる.
(1) J1ASについては,
(2) J2ASについては,
(3) Y1ASについては,
(4) Y2ASについては,
(5) I1ASについては,
(6) I2ASについては,
(7) K1ASについては,
(8) K2ASについては,
モードSSには,(1252)式のJ1SS,(1264)式のJ2SS ,(1296)式のY1SS,(1308)式のY2SS,(1340)式のI1SS,(1352)式のI2SS,(1384)式のK1SS,(1396)式の K2SS,の合計8セットが属する.任意の整数nに対しても解関数となるから,総数は 8nとなる.
(1) J1SSについては,
(2) J2SSについては,
(3) Y1SSについては,
(4) Y2SSについては,
(5) I1SSについては,
(6) I2SSについては,
(7) K1SSについては,
(8) K2SSについては,
モードAAには,(1253)式のJ1AA,(1265)式のJ2AA,(1297)式のY1AA,(1309)式のY2AA,(1341)式のI1AA,(1353)式のI2AA,(1385)式のK1AA,(1397)式のK2AA,の合計8セットが属する.任意の整数nに対しても解関数となるから,総数は8nとなる.
(1) J1AAについては,
(2) J2AAについては,
(3) Y1AAについては,
(4) Y2AAについては,
(5) I1AAについては,
(6) I2AAについては,
(7) K1AAについては,
(8) K2AAについては,
11.4.10節の解析解と比較するために,11.4.9節で与えた境界条件で問題を解く.即ち,y方向の一様重力下の円環で,内辺の境界条件を,表面力がゼロ,変位もゼロとするには,外辺をどのように支持すれば良いか求める.従来は求解不能とされた問題である.
固有関数に課される境界条件は,
第3モードまで用いた計算結果を図52に示す.応力στの分布が解析解に近づいていることがわかる.
本節は,10章に含める方が良い.
11.5.1 微分方程式
2次元渦なし流れの問題は,直交座標系で,連続の式
湧き出しをディラックの 関数で表現した場合,(18)式および(1545)式の右辺は 関数となる.また,渦度をディラックの 関数で表現した場合,(19)式および(1547)式の右辺は 関数となる.湧き出しや渦度の特性によって,方程式を切り替える必要を無くすためには,(18)式を
主解となる速度ux,uyを主速度と呼び u1,u2と書く.方程式の右辺をf1,f2と書けば,連立偏微分方程式(1548),(1549)式は,(23)式の形となり,
11.4.8節では重調和方程式の解を得たが,その幾つかは調和方程式(ラプラス方程式)を満足する.微分作用素の極座標形式は(992)式
(1552)式から得られる微分方程式
特殊な場合として,(1552)式でν=0とすれば,微分方程式
調和方程式(ラプラス方程式)
(1566)式から得られる微分方程式
特殊な場合として,(1572)式でν=0とすれば,微分方程式
座標原点にある湧き出しの強さQ[m2/s]を として,速度ポテンシャルφは,(1564)式より,
座標原点にある渦度の強さ(循環)をΓ[m2/s]として,速度ポテンシャルφは,(1564)式より,
非同次境界条件を満たす項に添字Bを付けてuBjとおき,同次境界条件を満たす項に添字Hを付けてuHjとおく.主速度 をその和で,
(1669)式と(1673)式は同じ微分方程式であるから,結局,φx,φyと、φx *,φy * は,同じ関数セットから構成されることがわかる.よって,(1669)式を解く事に照準を合わせて,
となる.これら,主・双対連立微分方程式を同時に満たす解は,(1593)式の組み合せから,以下の8つ得られる.
5.1節の固有関数法により解が求められる.得られた解に,11.5.7節と同様の無次元化を施して,流れの様子を示す.
観測領域を正方形として,固有関数のモード番号m×nを最大30×30とした計算を行った.外力項fiが(1698)式であることから,(1704)式の固有関数fl 11SSのみで特解部分uHjを表現できる.無次元速度の様子を図65に示す.この結果は,(1599),(1600)式の特解である.
図65は、解析解による様子(図55)と良く似ており,原点から流体が湧き出す様子が表現できている.しかし、特解部分のみ表現しているため,解析解との差がどのようなものか調べておく必要がある.(1595)式の解析解から,固有関数解を差し引いた様子を図66に示す.
図66より,境界周辺部での差が大きいことがわかる.これを残差流と呼ぶ.残差流を表現するのは,同次解部分のuBjである.特解部は(1704)式の固有関数F1 11SSのみで表現されており、固有値は、
Claims (16)
- 処理対象となる系の構造及び系の構成要素の性質を示すデータを読み込んで、読み込んだデータに基づいて、前記系を表す方程式であって、求めるべき物理量を表す変数を含むn個の初期方程式を決定する初期方程式決定部と、
境界条件を示すデータとして、前記物理量を表す値を読み込んで、境界条件を決定する境界条件決定部と、
前記n個の初期方程式を、2n個の変数の方程式又は2n個の方程式を含む方程式に変形し、前記変形した2n個の変数の方程式又は2n個の方程式を含む方程式において、前記境界条件によって既知となる変数を含む既知部と、未知の変数を含む未知部とを決定し、当該未知部について方程式の解を計算する演算部を備える、情報処理装置。 - 初期方程式決定部は、前記物理量を表す変数を含むn個微分方程式を決定し、
前記演算部は、前記初期方程式決定部が決定した前記微分方程式の微分作用素と、当該微分作用素から決まる随伴微分作用素を用いて、前記2n個の方程式を示すデータを生成し、2n個の方程式の解を計算することにより、少なくとも1つの前記物理量を出力する、請求項1に記載の情報処理装置。 - 前記初期方程式決定部は、前記系の構成要素の節点における物理量を示すn次元の変数ベクトルを2つと、n行の行列を含むn個の方程式を前記初期方程式として決定し、
前記境界条件決定部は、前記変数ベクトルにおいて、値が既知となる変数の自由度の数と、値が未知となる変数の自由度の数が異なる境界条件も決定可能とし、
前記演算部は、前記2つの変数ベクトルに基づいて2n次元のベクトルを生成し、当該2n次元のベクトルの変数に基づいて前記n行の行列を2n列の行列に変形した上で、当該2n次元のベクトルの変数のうち、前記境界条件によって既知となる変数を含む既知部分と、未知の変数を含む未知部分であって、必ずしも自由度が同じでない既知部分と未知部分を決定し、前記未知部分の変数を、前記既知部分の変数で表す形式になるよう前記2n列の行列及び前記2n次元のベクトルを変形して、変形した行列を用いて、前記未知部分の変数を計算する、請求項1に記載の情報処理装置。 - 前記決定したモード係数近傍の値の入力をユーザから受け付け、入力された値のモード係数を用いて解を計算し、前記解または前記解から得られる情報を出力する、請求項4に記載の情報処理装置。
- 処理対象となる系の構造及び系の構成要素の性質を示すデータを読み込んで、読み込んだデータに基づいて、前記系を表す微分方程式であって、求めるべき物理量を表す変数を含む微分方程式を初期方程式として決定する初期方程式決定部と、
境界条件を示すデータとして、前記物理量を表す値を読み込んで、境界条件を決定する境界条件決定部と、
前記微分方程式の前記変数と同じ数の双対変数と、双対微分方程式を定義した場合に、前記微分方程式の微分作用素を前記変数へ作用させたものと前記双対変数との積分和すなわち内積を部分積分して得られる境界項を計算し、前記境界条件の下で、当該境界項をゼロにする双対変数の条件を、随伴境界条件として決定する随伴境界条件決定部と、
前記随伴境界条件と前記境界条件とを比較した結果を出力する判定部とを備え、
前記微分方程式の微分作用素を前記変数へ作用させたものと前記双対変数との内積は、前記変数と前記双対微分方程式の微分作用素を前記双対変数へ作用させたものとの内積に等しい関係にある、情報処理装置。 - 処理対象となる系の複数の要素における物理量を表す2つのn次元の変数と、n行の行列で表されるn個の初期方程式を決定する設定部と、
境界条件を示すデータとして、前記物理量を表す値を読み込んで、境界条件を決定する境界条件決定部と、
前記2つのn次元の変数のうち、前記境界条件によって既知となる変数を含む既知部と、未知の変数を含む未知部を決定し、1つの変数の未知部の自由度が、他の変数の既知部の自由度と等しくない場合に、前記境界条件が非自己随伴である旨の情報を出力する判定部を備える、情報処理装置。 - 解析対象の原初微分作用素及び変数の境界条件を設定する設定部と、 前記境界条件から随伴境界条件を計算する随伴境界条件演算部と、 原初微分作用素から主微分作用素及び双対微分作用素を計算し、主連立微分方程式及び双対連立微分方程式、並びに、前記境界条件及び前記随伴境界条件を用いて、主固有関数及び双対固有関数を求めることにより、連立微分方程式の解を計算する非自己随伴演算部と、を備える情報処理装置。
- 前記境界条件と前記随伴境界条件が一致するか否かを判定する自己随伴判定部と、 前記境界条件と前記随伴境界条件が一致すると判定された場合、原初微分作用素から自己随伴問題の自己随伴固有関数を求めることにより、前記自己随伴問題の解を計算する自己随伴演算部をさらに備え、 前記境界条件と前記随伴境界条件が一致しないと判定された場合、前記非自己随伴演算部が、前記主連立微分方程式及び前記双対連立微分方程式、並びに、前記境界条件及び前記随伴境界条件を用いて、主固有関数及び双対固有関数を求めることにより、連立微分方程式の解を計算する、請求項9に記載の情報処理装置。
- 前記境界条件決定部又は前記設定部は、ユーザから前記物理量を表す変数で値が未知の部分を示す情報の入力を受け付け、当該情報を用いて前記境界条件を決定する、請求項1~3、6~10のいずれか1項に記載の情報処理装置。
- コンピュータが、処理対象となる系の構造及び系の構成要素の性質を示すデータを読み込んで、読み込んだデータに基づいて、前記系を表す方程式であって、求めるべき物理量を表す変数を含むn個の初期方程式を決定する初期方程式決定工程と、
コンピュータが、境界条件を示すデータとして、前記物理量を表す値を読み込んで、境界条件を決定する境界条件決定工程と、
コンピュータが、前記n個の方程式を、2n個の変数の又は2n個の方程式を含む方程式に変形し、前記変形した2n個の変数の方程式又は2n個の方程式を含む方程式において、前記境界条件によって既知となる変数を含む既知部と、未知の変数を含む未知部とを決定し、当該未知部について方程式の解を計算する演算工程とを有する、情報処理方法。 - 処理対象となる系の構造及び系の構成要素の性質を示すデータを読み込んで、読み込んだデータに基づいて、前記系を表す方程式であって、求めるべき物理量を表す変数を含むn個の初期方程式を決定する初期方程式決定処理と、
境界条件を示すデータとして、前記物理量を表す値を読み込んで、境界条件を決定する境界条件決定処理と、
前記n個の方程式を、2n個の変数の方程式又は2n個の方程式を含む方程式に変形し、前記変形した2n個の変数の方程式又は2n個の方程式を含む方程式において、前記境界条件によって既知となる変数を含む既知部と、未知の変数を含む未知部とを決定し、当該未知部について方程式の解を計算する演算処理とをコンピュータに実行させる、情報処理プログラム。 - コンピュータが、処理対象となる系の構造及び系の構成要素の性質を示すデータを読み込んで、読み込んだデータに基づいて、前記系を表す微分方程式であって、求めるべき物理量を表す変数を含むn個の微分方程式を初期方程式として決定する初期方程式決定工程と、
コンピュータが、境界条件を示すデータとして、前記物理量を表す値を読み込んで、境界条件を決定する境界条件決定工程と、
コンピュータが、前記微分方程式の前記変数と同じ数の双対変数と、双対微分方程式を定義した場合に、前記微分方程式の微分作用素を前記変数に作用させたものと前記双対変数との積分和すなわち内積を部分積分して得られる境界項を計算し、前記境界条件の下で、当該境界項をゼロにする双対変数の条件を、随伴境界条件として決定する随伴境界条件決定工程と、
コンピュータが、前記随伴境界条件と前記境界条件とを比較した結果を出力する出力工程を有し、
前記微分方程式の微分作用素を前記変数に作用させたものと前記双対変数との内積は、前記変数と前記双対微分方程式の微分作用素を前記双対変数に作用させたものとの内積に等しい関係にある、情報処理方法。 - 処理対象となる系の構造及び系の構成要素の性質を示すデータを読み込んで、読み込んだデータに基づいて、前記系を表す微分方程式であって、求めるべき物理量を表す変数を含むn個の微分方程式を初期方程式として決定する初期方程式決定処理と、
境界条件を示すデータとして、前記物理量を表す値を読み込んで、境界条件を決定する境界条件決定処理と、
前記微分方程式の前記変数と同じ数の双対変数と双対連立微分方程式を定義した場合に、前記微分方程式の微分作用素を前記変数に作用させたものと前記双対変数との積分和すなわち内積を部分積分して得られる境界項を計算し、前記境界条件の下で、当該境界項をゼロにする双対変数の条件を、随伴境界条件として決定する随伴境界条件決定処理と、
前記随伴境界条件と前記境界条件とを比較した結果を出力する出力処理をコンピュータに実行させ、
前記微分方程式の微分作用素を前記変数に作用させたものと前記双対変数との内積は、前記変数と前記双対微分方程式の微分作用素を前記双対変数に作用させたものとの内積に等しい関係にある、情報処理プログラム。
Priority Applications (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US14/354,632 US20140379314A1 (en) | 2011-10-31 | 2012-10-31 | Analyzer, analysis method, and analysis program |
EP12846477.3A EP2838040A4 (en) | 2011-10-31 | 2012-10-31 | DEVICE, METHOD, AND ANALYSIS SOFTWARE |
JP2013541831A JP5886309B2 (ja) | 2011-10-31 | 2012-10-31 | 解析装置、解析方法及び解析プログラム |
Applications Claiming Priority (4)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2011-239812 | 2011-10-31 | ||
JP2011239812 | 2011-10-31 | ||
JP2012085968 | 2012-04-04 | ||
JP2012-085968 | 2014-04-04 |
Publications (1)
Publication Number | Publication Date |
---|---|
WO2013065764A1 true WO2013065764A1 (ja) | 2013-05-10 |
Family
ID=48192104
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
PCT/JP2012/078262 WO2013065764A1 (ja) | 2011-10-31 | 2012-10-31 | 解析装置、解析方法及び解析プログラム |
Country Status (4)
Country | Link |
---|---|
US (1) | US20140379314A1 (ja) |
EP (1) | EP2838040A4 (ja) |
JP (1) | JP5886309B2 (ja) |
WO (1) | WO2013065764A1 (ja) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US10488553B2 (en) * | 2016-04-01 | 2019-11-26 | Baker Hughes, A Ge Company, Llc | Stress tensor computation using Mindlin formulation |
CN106777825B (zh) * | 2017-01-24 | 2020-03-27 | 湖南科技大学 | 一种基于谱有限元的矩形板振动模态计算方法 |
US11113436B1 (en) * | 2019-11-21 | 2021-09-07 | United States Of America As Represented By The Administrator Of Nasa | Method for simulation of flow in fluid flow network having one-dimensional and multi-dimensional flow components |
CN115081148B (zh) * | 2022-07-20 | 2022-11-15 | 上海索辰信息科技股份有限公司 | 一种基于势能理论的加筋板等效参数确定方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6718291B1 (en) * | 1999-07-02 | 2004-04-06 | Vadim Shapiro | Mesh-free method and system for modeling and analysis |
US7519518B2 (en) * | 2000-08-02 | 2009-04-14 | Comsol Ab | Method for assembling the finite element discretization of arbitrary weak equations, involving local or non-local multiphysics couplings |
-
2012
- 2012-10-31 EP EP12846477.3A patent/EP2838040A4/en not_active Withdrawn
- 2012-10-31 WO PCT/JP2012/078262 patent/WO2013065764A1/ja active Application Filing
- 2012-10-31 US US14/354,632 patent/US20140379314A1/en not_active Abandoned
- 2012-10-31 JP JP2013541831A patent/JP5886309B2/ja active Active
Non-Patent Citations (8)
Title |
---|
CHEN, CHUN HSIUNG ET AL.: "The Variational Principle for Non-Self-Adjoint Electromagnetic Problems", IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, vol. MTT-28, no. 8, August 1980 (1980-08-01), pages 878 - 886, XP055067259 * |
ERHARD SCHMIDT: "Zur Theorie der linearen und nichtlinearen Integralgleichungen", vol. 63, 1907, MATH. ANN., pages: 433 - 476 |
FUJIKAZU SAKAI: "Yugen Yosoho to Sabunho no Tokasei Oyobi Aru Risanka Shuho", JAPAN SOCIETY OF CIVIL ENGINEERS RONBUN HOKOKUSHU, December 1973 (1973-12-01), pages 39 - 52, XP055153927, Retrieved from the Internet <URL:http://library.jsce.or.jp/jsce/open/00037/220/220-126063.pdf> [retrieved on 20121116] * |
KIYOSHIGE MATSUMURA ET AL.: "2 Jigen Mitei Shinsuicho Mondai no Jiko Zuihan Henbun Genri ni Tsuite", 2004 NENDO SHUKI ZOSEN SANGAKKAI RENGO TAIKAI, THE KANSAI SOCIETY OF NAVAL ARCHITECTS, November 2004 (2004-11-01), pages 95 - 98, XP008173810 * |
KYUICHIRO WASHIZU; YOSHIYUKI YAMAMOTO; TADAHIKO KAWAI: "The Method of Weighted Residuals and Variational Principles", 1974, BAIFUKAN CO., LTD. |
See also references of EP2838040A4 |
T. BODDINGTON: "Variational Methods in Mathematical Physics", 1964, PERGAMON PRESS |
TOSHIYA SAITO; SHIGEYA MARUYAMA: "Methods of Mathematical Physics", vol. 1, 1959, TOKYOTOSHO CO., LTD. |
Also Published As
Publication number | Publication date |
---|---|
EP2838040A4 (en) | 2015-06-03 |
JP5886309B2 (ja) | 2016-03-16 |
JPWO2013065764A1 (ja) | 2015-04-02 |
EP2838040A1 (en) | 2015-02-18 |
US20140379314A1 (en) | 2014-12-25 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Bacciocchi et al. | Vibration analysis of variable thickness plates and shells by the generalized differential quadrature method | |
Liu et al. | Smoothed finite element methods | |
Ferreira et al. | Analysis of laminated shells by a sinusoidal shear deformation theory and radial basis functions collocation, accounting for through-the-thickness deformations | |
Liu et al. | The finite element method: a practical course | |
Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods | |
Reddy | Theory and analysis of elastic plates and shells | |
Kanno | Nonsmooth mechanics and convex optimization | |
Zhao et al. | Modal analysis of a rotating thin plate via absolute nodal coordinate formulation | |
Li et al. | Hamiltonian system-based new analytic free vibration solutions of cylindrical shell panels | |
JP5886309B2 (ja) | 解析装置、解析方法及び解析プログラム | |
Millithaler et al. | Structural dynamics of electric machine stators: Modelling guidelines and identification of three-dimensional equivalent material properties for multi-layered orthotropic laminates | |
Avery et al. | A computationally tractable framework for nonlinear dynamic multiscale modeling of membrane woven fabrics | |
Nguyen-Minh et al. | Static analysis of corrugated panels using homogenization models and a cell-based smoothed mindlin plate element (CS-MIN3) | |
Yin et al. | Degenerate grazing bifurcations in a three-degree-of-freedom impact oscillator | |
Cinefra | Formulation of 3D finite elements using curvilinear coordinates | |
Shamloofard et al. | Development of a shell superelement for large deformation and free vibration analysis of composite spherical shells | |
Mohazzab et al. | A spectral collocation solution for in-plane eigenvalue analysis of skew plates | |
Reille et al. | Learning data-driven reduced elastic and inelastic models of spot-welded patches | |
Pinho et al. | Nonlinear static analysis of thin shallow and non-shallow shells using tensor formulation | |
Jekel et al. | Obtaining a hyperelastic non-linear orthotropic material model via inverse bubble inflation analysis | |
Ansari et al. | Nonlinear bending analysis of first-order shear deformable microscale plates using a strain gradient quadrilateral element | |
Wang et al. | A physics-informed machine learning approach of improving RANS predicted reynolds stresses | |
Buljak et al. | Proper orthogonal decomposition and radial basis functions for fast simulations | |
Kim et al. | Nondimensional analysis of a two-dimensional shear deformable beam in absolute nodal coordinate formulation | |
Valle et al. | An efficient general curvilinear coordinates finite element method for the linear dynamic study of thickness-independent shells |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 12846477 Country of ref document: EP Kind code of ref document: A1 |
|
DPE1 | Request for preliminary examination filed after expiration of 19th month from priority date (pct application filed from 20040101) | ||
ENP | Entry into the national phase |
Ref document number: 2013541831 Country of ref document: JP Kind code of ref document: A |
|
NENP | Non-entry into the national phase |
Ref country code: DE |
|
WWE | Wipo information: entry into national phase |
Ref document number: 2012846477 Country of ref document: EP |
|
WWE | Wipo information: entry into national phase |
Ref document number: 14354632 Country of ref document: US |