JP3978377B2 - Molding simulation analysis method - Google Patents

Description

[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a technique for analyzing deformation due to press forming of a plate material by computer simulation.
[0002]
[Prior art]
There is a field in which products are manufactured by forming a plate material by plastic working using a mold. One example is in the field of manufacturing automobiles.
[0003]
In the molding of plate materials, mold correction is performed as necessary to improve various factors such as product formability, product productivity, product shape accuracy, and product surface quality (surface accuracy). . There is a strong demand for reduction in man-hours and time required for mold correction. For example, for mold correction for the purpose of improving formability and productivity, reduction of man-hours and time has been steadily achieved in recent years. is there.
[0004]
However, with regard to mold correction for the purpose of improving shape accuracy and surface quality, man-hours and time have not been reduced so much, especially the problem of poor shape accuracy is the man-hour for mold correction. This hinders the reduction of production preparation lead time.
[0005]
The shape accuracy of the product depends on the material of the product. For example, in the field of manufacturing automobiles, recently, the use of high-tensile steel plates and aluminum alloys as product materials is becoming popular for reducing the weight of automobiles and improving collision safety.
[0006]
On the other hand, there is a spring back of the product after molding as a cause of the poor shape accuracy of the product. The phenomenon of springback is a phenomenon in which when the material is unloaded after being loaded in plastic working such as press molding, the material elastically recovers and its shape changes. This amount of springback depends not only on the shape of the product but also on the type of material. For example, in the case of a high-strength steel plate and an aluminum alloy, there is a tendency that the amount of springback is larger than that of a material generally used conventionally.
[0007]
The following methods have already been proposed in order to reduce the man-hours and time required for mold correction for the purpose of improving the shape accuracy. This is because the state of deformation of the plate material by press forming is analyzed by numerical simulation, and the result of the analysis is used in the mold design stage to predict in advance the shape accuracy of the product. Thus, the mold is corrected on the computer so that the shape accuracy of the product is within the allowable range.
[0008]
Even if this method is adopted, the mold is corrected by machining or the like so as to improve the shape accuracy after measuring the shape accuracy of the product by molding using the actually manufactured mold. The process is not completely omitted. However, if that method is adopted, it is possible to virtually correct the mold, while the man-hours and time for the correction are compared with the man-hours and time for actually correcting the mold. Can be easily reduced. Therefore, if this method is adopted, it becomes easy to reduce the man-hours and time for designing and correcting the mold as a whole.
[0009]
In order to embody this technique, the following springback amount prediction method has already been proposed. In order to improve the accuracy of prediction of poor shape accuracy, we studied both aspects of material modeling and calculation technology, and modeled the nonlinearity of material characteristics at the time of material unloading and introduced it into FEM calculation as numerical analysis calculation It is a thing. This type of springback amount prediction method includes a molding simulation step and a springback analysis step, and one conventional example is disclosed in Japanese Patent Application Laid-Open No. 2000-312933.
[0010]
As a result of carrying out an experiment by embodying this kind of springback amount prediction method, the present inventors have confirmed that the springback amount can be predicted with a precision that can be used in mold design for a thin steel plate.
[0011]
The inventors further conducted an experiment for predicting the springback amount of a thick plate by the same springback amount prediction method. As an example of a product manufactured by molding a thick plate, there are functional parts such as a body frame and a suspension arm.
[0012]
[Problems to be solved by the invention]
However, the present inventors have found that it is difficult to perform the simulation of deformation due to press molding of a plate material with higher accuracy than the case of a thin plate when the plate material is a thick plate. I noticed.
[0013]
[Means for Solving the Problems and Effects of the Invention]
The following aspects are obtained by the present invention. Each mode is divided into sections, each section is numbered, and is described in a form that cites other section numbers as necessary. This is to facilitate understanding of some of the technical features described herein and some of the combinations thereof. The technical features and combinations of the technical features described herein are It should not be construed as limited.
[0014]
(1) The punch is brought into contact with the bent inner surface which is one of the two surfaces of the plate material to be bent and deformed while supporting the outer surface with at least one of the die and the pad in a pressurized state. A method of analyzing the deformation of a plate material when press-molding the plate material by computer simulation,
The shell element model modeled on the computer is used as an aggregate of a plurality of shell elements in which the plate members are arranged in the plane direction thereof, and the thickness of each shell element is considered in consideration of the thickness direction stress of each shell element. A molding simulation analysis method including a molding simulation step of analyzing in-plane stress that is stress in a plane direction.
[0015]
There is already a technique for predicting deformation due to press forming of a plate material by simulation.
[0016]
In this forming simulation, a plate material is virtually divided into a plurality of solid elements along both the surface direction and the plate thickness direction and modeled, and numerical analysis such as FEM is performed based on the solid element model. It is classified into a method to perform and a method in which a plate material is virtually divided into a plurality of shell elements in the plane direction and modeled, and a numerical analysis such as FEM is performed based on the shell element model.
[0017]
In general, the previous method is superior to the latter method in that it can easily improve the prediction accuracy of plate deformation, while the latter method can easily reduce the calculation time required for numerical analysis. This is superior to the previous method. Therefore, there is a demand for development of a method that is practically excellent in both prediction accuracy and calculation time.
[0018]
When the latter method is adopted, the conventional shell element is designed under the following analysis conditions (a) to (d) based on the Reissner-Mindlin plate theory.
[0019]
Hereinafter, the analysis conditions will be described. For convenience of explanation, the direction in which the neutral surface of the plate material extends, that is, the bending direction is referred to as the x direction, the plate thickness direction is referred to as the z direction, and the direction perpendicular to both directions is the y direction (cross section). Also called direction). Among them, the x direction and the y direction are both in-plane directions for the plate material, and the remaining z directions can be referred to as out-of-plane directions in the sense of intersecting the surface of the plate material.
[0020]
(A) Assuming a plane stress state for the shell element, the vertical stress σz in the thickness direction in the shell element is ignored and regarded as 0.
[0021]
(B) For each shell element, the plane retainability, that is, the cross section perpendicular to the shell element axis before bending deformation retains the plane even after the bending change. Therefore, the cross-sectional rotation angle θy (rotation angle about the y axis) of the shell element does not depend on the position in the plate thickness direction, that is, the z coordinate value.
[0022]
Here, the section rotation angle θy means the angle of the side of the shell element with respect to the normal of the neutral plane of the shell element. More precisely, if the plane cross-section perpendicular to the longitudinal axis of the shell element before bending deformation is the pre-deformation cross-section, and the cross-section to which the pre-deformation cross-section moved after the bending change is the post-deformation cross-section, It means the angle formed between the tangent line at an arbitrary point on the cross section and the pre-deformation cross section.
[0023]
(C) Since the cross-sectional rotation angle θy does not depend on the z coordinate value, the vertical strain εx in the bending direction in the shell element is linearly distributed according to the position in the plate thickness direction, and the shear strain γxz in the plate thickness direction is It is a constant value that is not zero.
[0024]
(D) Due to the planar retainability of the shell element, its neutral plane always coincides with the plate thickness center plane.
[0025]
In this conventional simulation simulation using shell elements, if the bending of the plate material is relatively weak like thin plate molding, the above analysis conditions are appropriate, so the deformation can be calculated accurately. it can.
[0026]
However, as a result of research by the present inventors, it has been found that it is difficult to accurately predict the deformation of a thick plate by this conventional forming simulation. In the case of thick plates, it is not appropriate to adopt the above analysis conditions, and the stress distribution of each shell element (described by normal stress σ and shear stress τ) and deformation state (normal strain ε and shear strain) It is necessary to establish a new analysis model and analysis technology that can reproduce (with γ) with higher accuracy.
[0027]
For that purpose, paying attention to the fact that the calculated values of the bending direction stress σx (similar to the cross-sectional direction stress σy) do not agree with each other when the plate thickness direction stress σz is considered and ignored, In consideration of this, it is important to calculate the bending direction stress σx (the same applies to the cross-sectional direction stress σy).
[0028]
Based on this knowledge, according to the method according to this section, a plate thickness of each shell element is obtained by using a shell element model modeled on a computer as an aggregate of a plurality of shell elements in which a plate is arranged in the plane direction. The in-plane stress that is the stress in the plane direction of each shell element is analyzed in consideration of the directional stress, and the deformation of the plate material is analyzed using the analysis result of the in-plane stress.
[0029]
Thus, according to this method, the deformation of the plate material is analyzed using the shell element model, so that the calculation time for simulating the deformation of the plate material can be shortened compared to the case of using the solid element model. Becomes easy.
[0030]
Moreover, according to this method, unlike the case where the deformation of the plate material is analyzed using the above-described conventional shell element model, the analysis is performed in consideration of the stress in the plate thickness direction. It becomes easy to approach when using a solid element model.
[0031]
Thus, according to this method, it is easy to achieve both reduction in calculation time and improvement in analysis accuracy.
[0032]
An example of the use of this method is to analyze the springback of the plate material based on the analysis result of the forming simulation, but this method can be applied to other uses.
[0033]
The method according to this section and each of the following items is of course effective when the plate is a thick plate, but can be applied not only to a thick plate forming simulation but also to a thin plate forming simulation. is there.
[0034]
In this section and the following sections, the term “stress” includes at least one of normal stress and shear stress, and the term “in-plane stress” includes at least one of bending direction stress and cross-sectional direction stress. The term “strain” includes at least one of normal strain and shear strain.
[0035]
The “molding simulation process” in this section analyzes the in-plane stress considering the thickness direction stress, and then uses the in-plane stress of the in-plane stress and the thickness direction stress, but uses the thickness direction stress. However, the present invention can be implemented in a mode in which the deformation of the plate material is analyzed, or in a mode in which the deformation of the plate material is analyzed using both the in-plane stress and the thickness direction stress. In the former mode, it is easier to analyze the deformation of the plate material by a simpler algorithm than the latter mode.
[0036]
(2) The molding simulation step includes:
For each shell element, a provisional value analysis step of analyzing a provisional value of the in-plane stress assuming a plane stress state in which the plate thickness direction stress is zero;
A final value analysis step for analyzing the final value of the in-plane stress for each shell element based on the analyzed provisional value and considering the balance of forces in the plate thickness direction;
The molding simulation analysis method according to item (1), including:
[0037]
When the shell element model is used, it is difficult to directly calculate the thickness direction stress σz. However, for example, assuming a minute element forming a partial cylinder, the plate thickness in the plate material to be bent and deformed is assumed. If the balance of force in the direction (radial direction) is expressed by an equation, the relationship between the plate thickness direction stress σz and the bending direction stress σx (similar to the cross-sectional direction stress σy) can be derived. Therefore, if this relationship is used, the plate thickness direction stress σz can be calculated from the bending direction stress σx (similar to the cross-sectional direction stress σy) calculated using the shell element model.
[0038]
Based on such knowledge, according to the method according to this section, the provisional value of the in-plane stress is analyzed for each shell element assuming a plane stress state in which the thickness direction stress is zero. Further, for each shell element, the final value of the in-plane stress is analyzed based on the analyzed provisional value and considering the balance of forces in the thickness direction.
[0039]
Therefore, according to this method, the in-plane stress calculated using the shell element model is effectively used even though the shell element model is difficult to directly calculate the thickness direction stress. By using it, it is possible to calculate the in-plane stress in consideration of the thickness direction stress.
[0040]
(3) The final value analysis step is such that the cross-sectional force, which is the sum of the in-plane stresses on the cross-section of each shell element, substantially matches the cross-sectional force of the shell element adjacent to the cross-section. The molding simulation analysis method according to item (2), further including a final value determination step of determining the final value.
[0041]
The method according to the item (2) is implemented in such a manner that the final value of the in-plane stress is determined as a value different from the provisional value by considering the balance of forces in the thickness direction for each shell element. Is possible.
[0042]
For each shell element, the provisional value of the in-plane stress, that is, the final value of the in-plane stress so as to result in a value different from the in-plane stress calculated by ignoring the thickness direction stress in the plane stress state. The determination is that the cross-sectional force, which is the sum of the in-plane stresses on the cross-section of each shell element, may differ from the cross-sectional force for the shell element adjacent to that cross-section. Such a situation is contrary to the in-plane balancing at the nodes shared by adjacent shell elements.
[0043]
Based on such knowledge, according to the method according to this section, the cross-sectional force, which is the sum of the in-plane stresses on the cross-section of each shell element, is substantially equal to the cross-sectional force of the shell element adjacent to the cross-section. The final value of the in-plane stress is determined so as to match.
[0044]
Therefore, according to this method, the continuity of the cross-sectional force is maintained between the adjacent shell elements even though the in-plane stress calculated in the plane stress state is changed in consideration of the thickness direction stress. It becomes possible.
[0045]
(4) The cross section in which the final value determining step is the sum of the in-plane stresses in the plane stress state, the cross-sectional force being the sum of the in-plane stresses considering the plate thickness direction stress on the cross section of each shell element (3) The step of determining the position of the neutral plane of each shell element so as to be substantially maintained by force, and determining the final value based on the determined neutral plane position. Molding simulation analysis method.
[0046]
According to this method, the position of the neutral plane is determined so that the cross-sectional force, which is the sum of the in-plane stresses acting on the cross-section of each shell element, is substantially maintained, and as a result, between the adjacent shell elements. It becomes possible to maintain the continuity of the cross-sectional force in step.
[0047]
(5) The final value analyzing step includes a step of analyzing the final value based on the analyzed provisional value and the contact pressure that the plate material receives from the punch on the bent inner surface during the press molding (2 The molding simulation analysis method according to any one of items (4) to (4).
[0048]
In order to accurately simulate the deformation due to press forming of a plate material, the inventors need to analyze the in-plane stress in consideration of the contact force with which the punch contacts the bent inner surface of the plate material during forming. I noticed.
[0049]
Based on this knowledge, according to the method according to this section, the final value of the in-plane stress is analyzed based on the provisional value of the in-plane stress and the contact pressure that the plate material receives from the punch on the bent inner surface during press forming.
[0050]
(6) The molding simulation step includes:
For each shell element, a provisional value analysis step of analyzing a provisional value of the in-plane stress assuming a plane stress state in which the plate thickness direction stress is zero;
For each shell element, based on the analyzed provisional value, identify a stress approximation function that approximately defines the distribution of the thickness direction stress in the thickness direction, and use the identified stress approximation function A final value analysis step of analyzing the final value of the in-plane stress;
The molding simulation analysis method according to item (1), including:
[0051]
This method can be considered to have a relationship constituting one embodiment of the method according to the item (2). Specifically, the implementation of the method according to this section specifies a stress approximation function that approximately defines the distribution of stress in the thickness direction in the thickness direction, and uses the identified stress approximation function. In the aspect of analyzing the final value of the in-plane stress, this is equivalent to the execution of the method according to the item (2).
[0052]
(7) The final value analysis step
A stress approximation function specifying step for specifying the stress approximation function based on the analyzed provisional value;
Using the identified stress approximation function, an in-plane stress determination step for determining an in-plane stress established together with the plate thickness direction stress as an in-plane stress,
Based on the determined in-plane stress, the cross-sectional force, which is the sum of the in-plane stresses on the cross section of each shell element, substantially matches the cross-sectional force of the shell element adjacent to the cross section. A final value determining step for determining the final value;
The molding simulation analysis method according to item (6) including:
[0053]
According to this method, by using the stress approximation function, the in-plane stress not considering the plate thickness direction stress is changed to the in-plane stress considering the plate thickness direction stress, that is, the in-plane stress. Furthermore, the final value of the in-plane stress is determined by considering the balance of forces between the shell elements based on the in-plane stress.
[0054]
Therefore, according to this method, the in-plane stress not considering the thickness direction stress is assumed as a provisional value, and the final in-plane stress is determined by considering both the thickness direction stress and the balance of forces between the shell elements. The value is determined.
[0055]
(8) In the cross-sectional force in which the correction step is a total of the in-plane stresses in the plane stress state, the cross-sectional force that is the sum of the in-plane stresses in consideration of the plate thickness direction stress on the cross section of each shell element The molding according to (7), including a step of determining a position of a neutral surface of each shell element so as to be substantially maintained, and determining the final value based on the determined neutral surface position. Simulation analysis method.
[0056]
This method has the same relationship with the method according to item (7) as the method according to item (4) has the same relationship as the method according to item (3).
[0057]
(9) The final value analysis step includes a step of specifying the stress approximation function based on the analyzed provisional value and a contact pressure that the plate member receives from the punch on the inner surface of the bending during the press molding ( The molding simulation analysis method according to any one of 2) to (8).
[0058]
The present inventors can specify the stress approximation function that approximately defines the distribution of the stress in the thickness direction in the thickness direction by considering only the in-plane stress. It has been found that the analysis accuracy of deformation due to press forming of the plate can be easily improved if the punch is specified in consideration of the contact force with which the punch contacts the bent inner surface of the plate.
[0059]
Based on such knowledge, according to the method according to this section, the stress approximation function of the stress in the plate thickness direction is based on the provisional value of the in-plane stress and the contact pressure that the plate material receives from the punch on the bent inner surface during press molding. Identified.
[0060]
(10) The molding simulation analysis method according to any one of (6) to (9), wherein the stress approximation function formulates the distribution on the premise of force balance in the plate thickness direction for each shell element. .
[0061]
(11) Each of the shell elements is defined by a plurality of nodes, and the plurality of nodes do not include nodes that are not essential for geometrically identifying the corresponding shell elements (1) to (10) The molding simulation analysis method according to any one of the items.
[0062]
According to this method, the plurality of nodes defining each shell element are constituted by the minimum number of nodes necessary for geometrically specifying the same shell element. Therefore, according to this method, it is possible to easily reduce the number of nodes that must be set for the entire plate, and it is possible to easily reduce the calculation time by reducing the number of nodes.
[0063]
(12) In the forming simulation step, the computer as an aggregate of a plurality of solid elements in which equivalent plate materials equivalent to the plate materials are arranged in both the surface direction and the plate thickness direction with respect to the severity of bending by the press forming. Using the solid element model modeled above (or the calculation result using the solid element model and the measurement result by the model experiment), on the plane cross section perpendicular to the longitudinal axis of the equivalent plate before the bending deformation Analyze the section rotation angle, which is the angle at which each point is rotated by the bending deformation, depending on the position in the thickness direction of the equivalent plate, and analyze the distortion of the equivalent plate using the analysis result of the section rotation angle. The molding simulation analysis method according to any one of (1) to (11), including a strain analysis step.
[0064]
As a result of various studies, the present inventors have found that, in the thick plate molding, unlike the thin plate molding, the plane retaining property is not established for each shell element. I also realized that it was important to make it dependent.
[0065]
Furthermore, the present inventors have also realized that it is important to consider the nonlinear distribution in the thickness direction of the strain by making the cross-sectional rotation angle θy dependent on the z coordinate value.
[0066]
Furthermore, the present inventors have the same bending severity by press molding (for example, it can be expressed by a bending index R / t, which is a value obtained by dividing the bending radius R of the plate material by the plate thickness t). As far as the bending radius R or the thickness t is different, the fact that the cross-sectional rotation angle and the distribution of strain in the thickness direction are common is also noticed.
[0067]
Therefore, when it is necessary to perform a forming simulation for a plurality of plate materials, it is not indispensable to analyze the section rotation angle (or model experiment) using the above-described solid element model for each plate material.
[0068]
Based on these findings, according to the method according to this section, regarding the severity of bending by press forming, an equivalent plate equivalent to the plate to be subjected to forming simulation is obtained in both the surface direction and the plate thickness direction. Using a solid element model (or a calculation result using a solid element model and a measurement result obtained by a model experiment) modeled on a computer as an assembly of a plurality of solid elements arranged side by side, the longitudinal direction of the equivalent plate before bending deformation Analysis is made so that the cross-section rotation angle, which is the angle at which each point on the plane cross-section perpendicular to the axis rotates by bending deformation, depends on the thickness direction position of the equivalent plate. Further, the distortion of the equivalent plate material is analyzed using the analysis result of the section rotation angle.
[0069]
In this section, “equivalent plate” is used as a term that does not exclude the plate itself for which in-plane stress is to be calculated. Further, in this section, “strain” includes at least one of bending direction strain and plate thickness direction strain, and the strain in each direction includes at least one of vertical strain and shear strain.
[0070]
(13) The severeness of the bending is defined to be stronger as the bending radius of the plate after bending deformation is smaller and to be stronger as the plate thickness of the plate is thicker. Molding simulation analysis method.
[0071]
(14) The molding simulation analysis method according to (13), wherein the severity of bending is expressed by a bending index obtained by dividing a bending radius of the plate material after the bending deformation by a plate thickness of the plate material.
[0072]
(15) The forming simulation analysis method according to the item (14), wherein the method is performed on a plate material having the bending index of 3 or less.
[0073]
(16) A strain approximation function specifying step for specifying a strain approximation function that approximately defines a distribution of the strain in the plate thickness direction based on the analysis result of the section rotation angle for the equivalent plate. Including
In-plane stress analysis in which the forming simulation step uses the shell element model with the specified strain approximation function for each shell element and analyzes the in-plane stress in consideration of the plate thickness direction stress. The molding simulation analysis method according to item (12), including a step.
[0074]
According to this method, the strain approximation function is based on the section rotation angle calculated in association with the bending severity using the solid element model (or the calculation result using the solid element model and the measurement result by the model experiment). And the in-plane stress is analyzed by using the shell element model together with the identified strain approximation function.
[0075]
Therefore, according to this method, although the shell element model is used to analyze the in-plane stress, it is easy to analyze the in-plane stress in consideration of the nonlinear distribution of strain in the plate thickness direction.
[0076]
(17) Further, the method includes a springback analysis step of analyzing, by a computer, a springback that is a phenomenon in which the plate material is elastically recovered as the plate material is unloaded after the plate material is plastically processed by press molding.
The molding simulation step analyzes residual stress remaining in the plate material at the time of unloading and residual strain remaining in the plate material at the time of unloading,
The forming simulation analysis method according to any one of (1) to (16), wherein the springback analysis step analyzes the springback using analysis values of the residual stress and residual strain.
[0077]
According to this method, it is possible to analyze the springback of the plate material with high accuracy by using the in-plane stress analyzed with high accuracy in consideration of the plate thickness direction stress.
[0078]
(18) Furthermore,
A specimen having substantially the same material as the plate material is loaded in one direction of the tensile direction and the compression direction and plastically deformed, then unloaded, and further loaded in the opposite direction to be plastically deformed. An experimental value acquisition step for repeating the acquisition of the experimental value for the stress-strain relationship of the test piece every time the prestrain, which is a plastic strain at the start of unloading of the test piece, is changed,
A material for obtaining the material characteristic value necessary for analyzing the springback in association with the prestrain by approximating the stress-strain relationship represented by the obtained experimental value with a linear material hardening model. Characteristic value acquisition process and
Including, and
When the analysis value of the residual strain is viewed as the current value of the pre-strain, the springback analysis step uses the material characteristic value obtained in association with the current value of the pre-strain, and the residual strain and the residual The forming simulation analysis method according to item (17), wherein the springback is analyzed according to the linear material hardening model using each analysis value of stress.
[0079]
(19) A program executed by a computer to implement the molding simulation analysis method according to any one of (1) to (18).
[0080]
If this program is executed by a computer, the same operational effects as the method according to any one of the above items (1) to (18) can be realized.
[0081]
The “program” in this section can be interpreted to include not only a combination of instructions executed by a computer to fulfill its function, but also files and data processed by each instruction.
[0082]
In addition, it is not indispensable to create the “program” in this section as a dedicated program created in order to fulfill the functions of each section or as a combination of dedicated subprograms (modules, etc.). It is possible to create a combination of a subprogram and a general-purpose subprogram. In some cases, a combination of a plurality of general-purpose subprograms without a dedicated subprogram may be generated.
[0083]
(20) A recording medium on which the program according to item (19) is recorded in a computer-readable manner.
[0084]
If the program recorded on the recording medium is executed by a computer, the same operation and effect as the method according to any one of the items (1) to (18) can be realized.
[0085]
The “recording medium” in this section can adopt various formats, such as a magnetic recording medium such as a flexible disk, an optical recording medium such as a CD and a CD-ROM, a magneto-optical recording medium such as an MO, and a ROM. At least one such as an unremovable storage can be employed.
[0086]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, one of more specific embodiments of the present invention will be described in detail with reference to the drawings.
[0087]
One embodiment of the present invention is a springback analysis method, and FIG. 1 shows a springback analysis device (hereinafter simply referred to as “analysis device”) suitable for implementing the springback analysis method. .
[0088]
This analysis device includes an input device 10, a computer 12, an output device 14, and an external storage device 16. The input device 10 is configured to include a mouse, a keyboard, and the like. The computer 12 includes a processor 20 such as a CPU, a memory 22 such as a ROM, a RAM, and a hard disk, and a bus 24 that connects the processor 20 and the memory 22 to each other. The output device 14 is configured as a display, a printer, a plotter, or the like.
[0089]
The external storage device 16 can be loaded with a recording medium such as a CD-ROM 30 and a writable MO (magneto-optical disk) 32. In the loaded state, data can be read from and written to the recording medium as necessary. Done.
[0090]
In the present embodiment, a program necessary for the analysis of the springback is stored in the CD-ROM 30, and data necessary for the analysis is stored in the MO 32. Programs necessary for the analysis of the springback include a material characteristic value acquisition program, a molding simulation program, a material characteristic value determination program, a backstress calculation program, and a springback analysis program.
[0091]
At the time of analysis of springback, necessary programs and data are read from the CD-ROM 30 and MO 32 and transferred to the RAM or hard disk of the computer 12, and then the processor 20 executes the programs.
[0092]
FIG. 2 is a process diagram illustrating the springback analysis method according to the present embodiment.
[0093]
If it demonstrates roughly, the springback analysis method will analyze the springback which arises in the board | plate material, after pressing the thick steel plate as a board | plate material. The press molding is performed such that there is a portion where a relatively high plastic strain occurs, such as a plastic strain of 0.15 and 0.20.
[0094]
In addition, hereinafter, the term “stress” means both normal stress and shear stress unless otherwise specified. Similarly, the term “strain” also means both normal strain and shear strain. Means.
[0095]
FIG. 3 illustrates a press molding machine 100. This press molding machine 100 is capable of bending a thick plate as a workpiece W into a U shape without using a wrinkle presser. In the press molding machine 100, an upper ram 102 and a lower ram 104 are opposed to each other with a gap therebetween, and a punch 110, a die 112, and a pad 114 are placed between the upper ram 102 and the lower ram 104. Arranged to sandwich W.
[0096]
The punch 110 receives a force in a direction approaching the die 112 from the upper ram 102 and is pressed against the workpiece W at a preset speed. The die 112 is fixed to the stage 120. The pad 114 is supported by the lower ram 104 via the hydraulic cylinder 122, and receives a force in a direction approaching the punch 110 from the hydraulic cylinder 122. The die 112 is attached to the stage 120 by a die set 124.
[0097]
In the springback analysis method of the present embodiment, by using a test piece that is substantially the same as the material to be press-molded, experimental values of the stress-strain relationship of the material are acquired as a plurality of discrete values, Based on the obtained experimental values, material characteristic values necessary for analyzing the springback of the material are obtained. This material property value can be said to be information necessary to define the stress-strain relationship of the material so that the Bauschinger effect of the material can be accurately expressed.
[0098]
In addition, in this springback analysis method, in order to analyze the stress-strain relationship remaining in the material immediately after press forming, a forming simulation is performed to simulate deformation of the material by press forming by using the finite element method. Is called.
[0099]
In this springback analysis method, the finite element method is then used,
(1) The material property value acquired by executing the material property value acquisition program,
(2) Residual stress and residual strain calculated by execution of the forming simulation program, neutral surface position for each element of the plate, and plate thickness for each element of the plate,
(3) By using the back stress calculated by executing the back stress calculation program, the spring back of the material is analyzed.
[0100]
More specifically, in the present embodiment, as shown in FIG. 2, first, in step S1, a test piece having substantially the same material as the material to be press-molded is prepared by the operator. Thereafter, by using a tester (not shown), a tensile-compression test is performed on the test piece in a uniaxial stress state.
[0101]
Specifically, this test involves unloading the specimen after it is plastically deformed by loading it in the tensile direction, and then plastically deforming by loading it in the opposite direction, that is, in the compression direction. Experimental values of the stress-strain relationship are obtained as a plurality of discrete values.
[0102]
The plurality of experimental values are acquired every time the pre-strain, which is a plastic strain at the start of unloading of the test piece, is changed. FIG. 4 shows graphs representing stress-strain relationships during unloading and compression for five prestrains, respectively.
[0103]
Next, in step S2, the material characteristic value acquisition program is executed by the computer 12, whereby the material characteristic value is determined for each pre-strain based on the plurality of acquired experimental values. The contents of the material property value acquisition program will be described in detail later.
[0104]
Thereafter, in step S3, data necessary for executing the molding simulation, that is, data for numerical analysis is created by the operator. Data for numerical analysis includes data for virtually dividing a plate material to be press-molded into a plurality of shell elements, and a load applied to the plate material during press molding (including loads applied to the plate material from punches and pads) And a boundary condition for a plurality of shell elements.
[0105]
In the present embodiment, the plurality of nodes defining each shell element is configured not to include nodes that are not essential for geometrically specifying each shell element. An example of such a shell element is a quadrilateral element that forms a quadrilateral and has four vertices as four nodes. By defining each shell element as described above, an increase in the number of nodes set in the plate material is avoided, and as a result, an increase in calculation time accompanying an increase in the number of nodes is also avoided.
[0106]
Subsequently, in step S4, the computer 12 executes the above-described forming simulation program, thereby analyzing the stress-strain relationship remaining in the plate material immediately after press forming. Furthermore, the shape of the material immediately after press molding is also analyzed.
[0107]
This forming simulation program is a form that uses an isotropic hardening model as a material hardening model that approximates the stress-strain relationship of materials. This forming simulation program stores an initial strain, which is a plastic strain existing in each shell element immediately after press forming, in the memory 22 as a residual strain in association with the number j of each shell element of the plate material. Details of this molding simulation program will be described later.
[0108]
This molding simulation program is manufactured by LSTC and sold by the Japan Research Institute under the name “LS-DYNA3D” and by ESI and sold under the name “PAM-STAMP”. And a new step is added to them.
[0109]
Thereafter, in step S5, the material characteristic value determination program is executed by the computer 12, and by using the plurality of material characteristic values acquired in step S2, the plasticity associated with the molding of each shell element of the plate material is obtained. A material property value corresponding to the strain is determined. The contents of the material property value determination program will be described in detail later.
[0110]
Subsequently, in step S6, the computer 12 executes the back stress calculation program, whereby the back stress α is calculated for each shell element of the plate material. Details of the backstress calculation program will be described later.
[0111]
Thereafter, in step S7, the computer 12 executes the springback analysis program, thereby analyzing the springback. This springback analysis program is clear from the above explanation,
(1) The material property value determined for each shell element of the plate material by executing the material property value determination program,
(2) Residual stress and residual strain calculated for each shell element by execution of the molding simulation program, neutral surface position r₀And thickness t₀When,
(3) Back stress calculated for each shell element by executing the back stress calculation program
Is a program for analyzing the springback of a plate material.
[0112]
An example of this springback analysis program is manufactured by the Japan Research Institute and sold under the name “JOH-NIKE3D”.
[0113]
This completes one execution of the springback analysis method.
[0114]
The springback analysis program is designed so that the stress-strain relationship of the material whose springback is to be analyzed can be defined in various ways. One method is to approximate the stress-strain relationship with a linear material hardening model using a parallelogram. In FIG. 5, when the test piece was loaded in the tensile direction (positive load) and plastically deformed, then unloaded, and subsequently loaded in the reverse direction, that is, in the compression direction (reverse load), and plastically deformed again, Plastic strain ε^PAnd the state where the stress σ changes are shown by parallelograms.
[0115]
Here, various symbols in FIG.₀"Indicates the initial yield stress, which is the yield stress when the material begins its first plastic deformation. "Y₁"Indicates the re-yield stress, which is the stress when the material re-yields. “H” indicates a plastic hardening coefficient. The plastic hardening coefficient H is expressed by tan φ when the slope of a straight line graph (third linear portion 54 described later) representing the stress-strain relationship in the plastic region at the time of positive load and reverse load is expressed by angle φ. “Β” is the re-yield stress Y₁Is the initial yield stress Y₀ The material softening coefficient that defines the degree of decrease is shown. "Ε^P ₀ Is the plastic strain ε at the start of unloading^PThe initial strain is shown.
[0116]
The springback analysis program can select an appropriate one from an isotropic hardening model, a kinematic hardening model, and a composite hardening model as a material hardening model that approximates the stress-strain relationship.
[0117]
This selection is made by determining the value of the material softening coefficient β. Specifically, as shown in FIG. 5, when β = 0, the kinematic hardening model is selected, and at this time, the re-yield stress Y₁Is
Y₁= -Y₀+ H · ε^P ₀
Is defined by the expression Further, when β = 1, an isotropic hardening model is selected, and at this time, the re-yield stress Y₁ Is
Y₁= -Y₀-H ・ ε^P ₀
Is defined by the expression When 0 <β <1, a composite hardening model is selected, and at this time, the re-yield stress Y₁ Is
Y₁= -Y₀+ H · (1-2β) ε^P ₀
Is defined by the expression
[0118]
The graph of FIG. 6 relates to the material to be pressed, but the same graph will be drawn for the specimen. When the graph representing the stress-strain relationship of the test piece is approximated by a parallelogram, the first linear portion 50 representing the stress-strain relationship in the elastic region of the test piece at the time of positive load, and at the time of unloading Alternatively, the second linear portion 52 representing the stress-strain relationship in the elastic region at the time of reverse loading is parallel to each other, and the third linear portion 54 representing the stress-strain relationship in the plastic region at the positive load of the test piece The stress-strain relationship is approximated so that the fourth linear portion 56 representing the stress-strain relationship in the plastic region during reverse loading is parallel to each other.
[0119]
On the other hand, in the graph representing the stress-strain relationship, the region where the Bauschinger effect appears is the unloading start point P._S To yield point P_Y Elastic range up to and its re-yield point P_Y And the plastic region from the end point of reverse load.
[0120]
Therefore, in the present embodiment, a polygonal line that approximately represents the stress-strain relationship in these two regions, that is, a geometric feature of the second straight line portion 52 and the fourth straight line portion 56 connected to each other is acquired. And based on the obtained geometric characteristics, the initial yield stress Y is directly applied to the specimen and indirectly to the material.₀And re-yield stress Y₁And the plastic hardening coefficient H and the material softening coefficient β are obtained.
[0121]
Specifically, the re-yield stress Y₁Is obtained as the stress σ at the break point of the broken line. The plastic hardening coefficient H is obtained by using the fact that the third straight line portion 54 is parallel to the fourth straight line portion 56 to obtain the gradient of the fourth straight line portion 56 at an angle φ and calculating its tan φ. Obtained by. The material softening coefficient β is the initial yield stress Y₀And re-yield stress Y₁, Plastic hardening coefficient H and pre-strain ε^P _gAnd a formula similar to the previous formula, i.e.
Y₁= -Y₀+ H · (1-2β) ε^P _g
It is obtained by substituting into the following equation and solving for β. Initial yield stress Y₀Is obtained based on experimental values at the time of positive load.
[0122]
This springback analysis program uses these initial yield stresses Y₀, Re-yield stress Y₁, Designed to analyze springback by utilizing the elastic modulus E of material (selectively representing Young's modulus E and shear elastic modulus G) in addition to plastic hardening coefficient H and material softening coefficient β. ing. That is, in this springback analysis program, the initial yield stress Y₀, Re-yield stress Y₁The plastic hardening coefficient H, the material softening coefficient β, and the elastic coefficient E constitute “material characteristic values”.
[0123]
By the way, what is taken on the horizontal axis of the graph of FIG. 5 is not the total strain ε, which is the sum of elastic strain and plastic strain, but plastic strain ε.^PIt is. Therefore, for example, the graph of FIG. 2 extends in parallel with the axis representing the stress σ in the elastic region at the time of positive load. In contrast, the horizontal axis of the graph of FIG. 6 is the total strain ε. Therefore, for example, the graph of FIG. 2 is inclined with respect to the axis representing the stress σ in the elastic region at the time of positive load. Therefore, in these two graphs, the correspondence relationship between the first to fourth linear portions 50, 52, 54, and 56 is as shown in these two figures.
[0124]
Therefore, when the gradient of the first straight line portion 50 in FIG. 6 is represented by an angle γ, the value obtained as tan γ represents the Young's modulus E of the material. On the other hand, it is assumed that the first straight part 50 and the second straight part 52 are parallel to each other. Therefore, in the present embodiment, the Young's modulus E is acquired from the gradient of the second straight part 52.
[0125]
From the acquired Young's modulus E, a shear elastic modulus G for the same material is derived according to a predetermined rule. Thereby, the Young's modulus E and the shear elastic modulus G, which are collectively referred to as the elastic modulus E, are acquired for the same material.
[0126]
The stress-strain relationship is the initial strain ε of the material^P ₀(Prestrain of specimen ε^P _gCorresponding to). The geometrical feature of the parallelogram that approximately represents the stress-strain relationship is the initial strain ε^P ₀Therefore, the material characteristic value is also the initial strain ε.^P ₀It depends on. Therefore, in the present embodiment, as described above, the pre-strain ε of the test piece^P _gThe tension-compression test is repeated each time the is changed.
[0127]
In FIG. 7, the linear material hardening model used in this embodiment is shown by a solid line graph, the isotropic hardening model is shown by a two-dot chain line graph, and the actual stress-strain relationship is shown by a broken line graph. Yes. The actual stress-strain relationship is consistent with the relationship represented by the isotropic hardening model (two-dot chain line) in the elastic region and plastic region under positive load, and in the elastic region and reverse plastic region during reverse load. Corresponds to the relationship (solid line) represented by the material hardening model of this embodiment. The material hardening model used in this embodiment is linear, but the pre-strain ε^P _gIt is defined that the shape changes according to. Note that the parallelogram representing the linear kinematic hardening model used in the present embodiment is omitted in FIG.
[0128]
In this embodiment, a kinematic hardening model or a composite hardening model that can express the Bauschinger effect can be used as the material hardening model, and the expected initial strain ε for each shell element of the material.^P ₀Depending on the material curing model with different shapes is used. However, the expected initial strain ε for each shell element^P ₀Even if they are different from each other, it is maintained that the material hardening model is expressed by a parallelogram.
[0129]
Therefore, as shown in FIG. 7, the linear material hardening model used in the present embodiment is not entirely consistent with the graph representing the actual stress-strain relationship, but the region where the Bauschinger effect appears, that is, The region that undergoes elastic deformation after unloading coincides with the region that undergoes plastic deformation with sufficiently high accuracy.
[0130]
For this reason, in the present embodiment, the stress-strain relationship of the material is approximated by a linear material hardening model using a parallelogram. The Bauschinger effect can be accurately simulated not only in the strain region but also in the high strain region, and as a result, the springback can be analyzed with high accuracy.
[0131]
As shown in FIG. 7, the broken line is acquired so as to match the graph representing the actual stress-strain relationship during unloading and reverse loading as much as possible. Specifically, on the graph, the provisional setting of three points that are not on the same straight line and that are separated from each other by a certain distance or more may be extracted as a line in the graph. It is performed sequentially from one end point of a part to the other end point, and when the angle of the temporary polygonal line composed of these temporary three points becomes maximum, these three points are changed to the final three points. At the same time, the final broken line composed of these three points is determined as the final broken line.
[0132]
These three points are composed of a first actual measurement point 60 that is closest to the unloading start point, a second actual measurement point 62 that is next closest, and a third actual measurement point 64 that is farthest away. Each stress-strain relationship at these three actual measurement points 60, 62, and 64 is located on a graph representing the actual stress-strain relationship, and therefore coincides with the actual measurement value.
[0133]
In the example shown in FIG. 7, of the two straight portions constituting the broken line, the one near the unloading start point (not shown) is not completely the same as the second straight portion 52 of the parallelogram. The other straight line partially coincides with the parallel straight fourth straight part 56, but not completely.
[0134]
Here, the contents of the material property value acquisition program will be described in detail.
[0135]
FIG. 8 shows a flowchart of this material property value acquisition program. First, in step S101 (hereinafter simply expressed as “S101”, the same applies to other steps), data representing experimental values is taken from the MO 32 into the memory 22 of the computer 12.
[0136]
Next, in S102, the pre-strain ε^P _gThe number i to be assigned to is set to 1. Thereafter, in S103, the current pre-strain ε^P _gData representing the stress-strain relationship related to (i) is taken from the memory 22.
[0137]
Subsequently, in S104, based on the captured data, the current pre-strain ε^P _gA graph representing the stress-strain relationship related to (i) is assumed, and the graph is approximated by a broken line. Thereby, the three measured points 60, 62, 64 on the parallelogram are extracted as described above. Thereafter, in S105, the gradient of the second straight line 52 connecting the second actual measurement point 62 and the third actual measurement point 64 is calculated by the angle φ, whereby the present pre-strain ε^P _gA plastic hardening coefficient H (i) (= tan φ) corresponding to (i) is calculated.
[0138]
Subsequently, in S106, the pre-strain ε of this time is obtained by using the above-described equation shown by the equation (1) in FIG.^P _gA material softening coefficient β (i) corresponding to (i) is calculated. In the present embodiment, the calculation is performed on the assumption that the material softening coefficient β is in the range of 0 <β <1, that is, the material hardening model is a composite hardening model. Hereinafter, the calculation method will be specifically described.
[0139]
When the material softening coefficient β is calculated using the above equation (1), the first condition shown by the equation (2) in FIG.^P _g The second condition that the condition represented by the expression (3) is satisfied is adopted as the constraint condition. Further, “Y” in the above formula (1)₁“And“ H ”are obtained by experimental values.
[0140]
Re-yield stress Y₁For various types of experiments, prestrain ε^P _gHas been confirmed to show a monotonic increase with respect to^P _gWhen = 0
Y₁= -Y₀
And all the prestrains ε^P _gabout,
Y₁> -Y₀
The relationship is established.
[0141]
When the first condition of the above formula (2) is considered and the above formula (1) is used, the formula (4) is obtained. This equation (4) can be transformed into equation (5). When this equation (5) is used and the second condition shown by the above equation (3) is taken into consideration, equation (6) is obtained, and equation (7) is obtained from this equation (6).
[0142]
Further, when the second condition shown in the above equation (3) is considered and the above equation (1) is solved for “β”, equation (8) is obtained.
[0143]
In this embodiment, it is determined whether or not the condition expressed by the above formula (7) is satisfied. If not, the present pre-strain ε is set so as to be satisfied.^P _g(I) is modified. For example, this pre-strain ε^P _g(I)
1.01 × (−Y₁/ H)
Or by selecting the smallest of the experimental values that satisfy the condition shown by the above equation (7).
[0144]
If the condition shown in the above equation (7) is satisfied, the present pre-strain ε modified as necessary^P _g(I) and initial yield stress Y obtained from experimental values₀, Re-yield stress Y₁Then, the material softening coefficient β (i) of this time is calculated by substituting the plastic hardening coefficient H (i) into the above formula (1).
[0145]
Further, it is determined whether or not the calculated material softening coefficient β (i) satisfies the first condition expressed by the above formula (2), that is, whether or not it is within the range of 0 <β <1. If so, the calculated material softening coefficient β (i) is the final value. If the first condition is not satisfied, the material softening coefficient β (i) is corrected so that the first condition is satisfied, and the initial yield stress Y is adjusted accordingly.₀Will also be fixed.
[0146]
Material softening coefficient β (i) and initial yield stress Y₀For example, when the original material softening coefficient β (i) is 0 or less, the material softening coefficient β (i) is changed to a set value larger than 0, for example, 0.1, and the material softening is performed. Initial yield stress Y when coefficient β (i) is 0.1₀Is calculated using the above equation (1). The calculated initial yield stress Y₀, The latest value β (i) of the material softening coefficient β and the initial yield stress Y₀Are adopted as final values as they are.
[0147]
In contrast, the initial yield stress Y under the modified material softening coefficient β (i)₀Does not satisfy the condition expressed by the above equation (3), the material softening coefficient β (i) is increased by, for example, 0.1, and similarly, a new initial yield stress Y₀Is calculated. Such a process is repeated within a range where the material softening coefficient β does not exceed 1, so that the material softening coefficient β (i) is a re-yield stress Y based on experimental values.₁And the plastic hardening coefficient H (i) is corrected so as to satisfy all of the above formulas (1) to (3).
[0148]
As described above, the correction method when the original material softening coefficient β (i) is 0 or less has been described.
[0149]
FIG. 10 is a flowchart showing details of S106 as a material softening coefficient calculation routine.
[0150]
First, in S251, by using the extracted second actual measurement point 62, the re-yield stress Y₁Is determined. Next, in S252, the determined re-yield stress Y₁And the calculated plastic hardening coefficient H (i) and the current pre-strain ε^P _gIt is determined whether or not (i) satisfies the condition shown by the equation (7) in FIG.
[0151]
If it is assumed that it is satisfied this time, the determination is YES, and in S253, the re-yield stress Y determined above is determined.₁And the calculated plastic hardening coefficient H (i) and the initial yield stress Y obtained from experimental values.₀Is substituted into the equation (1) to calculate the current material softening coefficient β (i). Thereafter, the process proceeds to S254.
[0152]
On the other hand, this time, if it is assumed that the condition expressed by the equation (7) is not satisfied, the determination in S252 is NO, and in S255, the current pre-strain ε^P _gAs described above, (i) is corrected so as to satisfy the condition shown in the equation (7).
[0153]
Thereafter, in S253, the corrected pre-strain ε^P _gThe material softening coefficient β (i) of this time is calculated by using the equation (1) under (i). Subsequently, the process proceeds to S254.
[0154]
In any case, after that, in S254, it is determined whether or not the calculated material softening coefficient β (i) satisfies the condition shown by the equation (2) in FIG. This time, if it is assumed that the condition is satisfied, the determination is YES, and the calculated material softening coefficient β (i) and the initial yield stress Y based on the experimental value are obtained.₀And are adopted as final values as they are. This completes one execution of this routine.
[0155]
On the other hand, if it is assumed that it is not satisfied this time, the determination is NO, and in S256, the current material softening coefficient β (i) is corrected as described above. Thereafter, in S257, the initial yield stress Y is obtained by using the above equation (1) under the corrected material softening coefficient β (i).₀Is calculated.
[0156]
Subsequently, in S258, the calculated initial yield stress Y₀Whether or not is greater than 0 is determined. This time, if it is assumed that it is greater than 0, the determination is YES, the corrected material softening coefficient β (i), and the calculated initial yield stress Y₀Are adopted as final values. Thus, one execution of this routine is completed.
[0157]
In addition, this time, the calculated initial yield stress Y₀If NO is less than or equal to 0, the determination in S258 is NO, and the process returns to S256. In S256, the material softening coefficient β (i) is corrected again as described above. Thereafter, in S257, a new initial yield stress Y is added below the corrected material softening coefficient β (i).₀Is calculated.
[0158]
Subsequently, in S258, the calculated initial yield stress Y in the same manner as described above.₀Whether or not is greater than 0 is determined. As a result of repeating S256 to S258 several times, if the determination in S258 is YES, the material softening coefficient β (i) and the initial yield stress Y that are finally corrected₀Are adopted as final values. This completes one execution of this routine.
[0159]
If the execution of S106 is completed as described above, then, in S107 of FIG. 8, the gradient of the second straight line 52 connecting the first measured point 60 and the second measured point 62 is calculated by the angle γ. tanγ is the current prestrain ε^P _gYoung E (i) corresponding to (i) is calculated. From the calculated Young's modulus E (i), a shear elastic modulus G (i) is derived according to a predetermined rule. In this way, Young's modulus E (i) and shear modulus G (i) are obtained.
[0160]
Subsequently, in S108, it is determined whether or not the number i is greater than or equal to the maximum value iMAX. It is determined whether all the experimental values have been used to determine the material property value. If it is assumed that the number i is not greater than or equal to the maximum value iMAX this time, the determination is no, the number i is incremented by 1 in S109, and then the process proceeds to S103. Hereinafter, S103 to S108 are the new pre-strain ε.^P _gIt will be executed in connection with (i).
[0161]
As a result of repeating the execution of S103 to S109 several times, if the determination in S108 is YES, in S110, the pre-strain ε^P _gThe change area is divided into a plurality of sections at regular intervals between the set minimum value and maximum value.
[0162]
At this time, pre-strain ε corresponding to each category^P _gIf there is no material characteristic value corresponding to^P _gMultiple prestrains ε near^P _gThen, necessary material property values are generated by interpolating those having material property values. As a result, pre-strain ε^P _gAs shown in the tabular form in FIG.^P _gThe plurality of representative values are obtained discretely.
[0163]
This completes one execution of the material property value acquisition program.
[0164]
Next, the contents of the molding simulation program will be described in detail.
[0165]
FIG. 12 conceptually shows the contents of the molding simulation program in a flowchart, and details of step S403 are shown in flowcharts as element stress calculation routines in FIGS. 13 and 14. The contents of the molding simulation program will be described below. Prior to that, the premise theory and principle will be described.
[0166]
In this forming simulation program, deformation due to press forming of a plate material is analyzed using a finite element method of time increment type.
[0167]
Here, as is well known, the finite element method uses a stiffness matrix that divides a material into a plurality of finite elements and correlates stress and displacement at each node of each element, and solves static problems. This is a numerical analysis method in which iterative calculations of stiffness equations, which are simultaneous equations, are repeated until the solution converges to a true solution, subject to the balance between external force and internal force at each node of each element. In the present embodiment, each element of the overall stiffness matrix formed by a collection of a plurality of element stiffness matrices is an unknown number, which is a solution to be obtained, and repeated calculation is performed until the solution converges to a true solution. .
[0168]
In the present embodiment, the finite element method is executed in combination with the incremental calculation method. Incremental calculation methods, as is well known, divide time into small intervals to analyze nonlinear behavior and material behavior during dynamic processes, and then calculate each physical quantity (displacement, This is a numerical analysis method that advances the analysis up to a predetermined time by sequentially accumulating increments of rotation angle, stress, and strain.
[0169]
In the present embodiment, the material behavior accompanying a single operation of the punch during press forming is divided into a plurality of calculation steps each having a minute time, and a series of forming simulations is performed by the completion of the plurality of calculation steps. Ends. Furthermore, in this embodiment, a dynamic explicit method is employed as the incremental calculation method.
[0170]
By this dynamic explicit method, a dynamic equation of motion equivalent to the equation of motion of the spring is solved. Specifically, the acceleration of each node of each shell element is calculated, the velocity of each node is calculated by time integration of the acceleration, and the displacement of each node is calculated by time integration of the velocity.
[0171]
In this forming simulation program, stress and strain are further analyzed based on a stress update algorithm for an elastic-plastic body. This stress update algorithm analyzes the stress increments for the strain increments that occurred during each calculation step from time t to time t + Δt for each shell element under the Von Mises yield condition. This stress update algorithm is embodied as an element stress calculation routine which will be described later with reference to FIGS. 13 and 14.
[0172]
In this stress update algorithm, the displacement increment u and the rotation angle increment θ in the current calculation step are calculated. The displacement increment u should originally be expressed as “du” or “Δu”, but is simply expressed as “u” here. The displacement increment u is a general term for the displacement increments u, v, and w of each component. Similarly, the rotation angle increment θ should originally be expressed as “dθ” or “Δθ”, but is simply expressed as “θ” here. The rotation angle increment θ is a general term for the rotation angle increments θy, θy, θz of each component.
[0173]
Specifically, the calculation of the displacement increment u is performed by using an incremental type overall stiffness equation considering force balance. On the other hand, the calculation of the rotation angle increment θ is performed by using an incremental type overall stiffness equation considering moment balance. In any case, the total stiffness equation includes a total stiffness matrix having an unknown component, and the unknown component of the total stiffness matrix is equal to the calculated value up to the previous calculation step in the current calculation step. . In this embodiment, an unknown component of the entire stiffness matrix is obtained by the finite element method.
[0174]
The strain increment is the isotropic strain increment ε_mAnd deviation strain increment e_ijAnd separated. Its isotropic strain increment ε_mIs the strain increment for each component ε_ijAnd its strain increment ε_ijIs calculated using a strain approximation function described later.
[0175]
Isotropic strain increment ε_mFrom hydrostatic pressure increment 3Kε_m(K: bulk modulus) is determined. This hydrostatic pressure increment 3Kε_mThe hydrostatic pressure component σ up to the previous calculation step_ij ⁰By adding hydrostatic pressure component σ_mIs updated. That is, the updated hydrostatic pressure component σ_m ^{t + Δt}Is calculated, which is shown in FIG. 15 by equation (11).
[0176]
Strain increment for each component ε_ijTo isotropic strain increment ε_mBy subtracting, deviation strain increment e for each component_ijIs required.
[0177]
This deviation strain increment e_ijDeviation stress component 2Ge_ij(G: shear modulus) is determined. Therefore, the deviation stress component 2Ge_ijIs obtained in the elastic region in a plane stress state. This deviation stress component 2Ge_ijThe deviation stress component S up to the previous calculation step_ij ⁰Is added to the trial stress (trial value of the deviation stress S) S_ij ^TrialIs required. That is, the updated trial stress S_ij ^TrialIs calculated, which is shown in equation (12).
[0178]
This trial stress S_ij ^TrialIs substituted into the equation (13) to obtain a trial solution σ of equivalent stress^TrialIs required.
[0179]
Trial solution σ of this equivalent stress^TrialIs within the elastic region and the material of each shell element has not yet yielded, the trial stress S_ij ^TrialIs the updated deviating stress component (current deviating stress component) S_ij ^FinalIt is said.
[0180]
On the other hand, as shown in FIG. 16A representing the stress space, the trial solution σ of the equivalent stress^TrialIs reached, it is determined that the material of each shell element has yielded. In this case, as shown in FIG. 15 (b), the trial stress S is obtained by the radial return algorithm (radius pullback method) shown in FIG._ij ^TrialIs pulled back onto the yield surface.
[0181]
In the radial return algorithm, the equivalent plastic strain increment ε^PFrom the equivalent stress σ representing the radius of the yield surface and the final solution S_ij ^FinalThe provisional value of is calculated. Further, the two calculated values are compared, and the final solution S when the equivalent stress σ is matched._ij ^FinalIs the final solution S_ij ^FinalThe final value is determined.
[0182]
In addition, in the case of adopting the above-mentioned von Mises yield condition, for example, by using the equation (16) in FIG.^PAnd then the final solution S by using equation (14)_ij ^FinalCan be calculated.
[0183]
If the material of each shell element yields, the final solution S pulled back_ij ^FinalIs obtained as the updated deviation stress component.
[0184]
Then, as shown in FIG. 15 by the equation (15), the updated hydrostatic pressure component σ_m ^{t + Δt}And the updated deviation stress component S_ij ^FinalAnd stress σ updated for each component_ij ^{t + Δt}Is calculated.
[0185]
One calculation step of the stress update described above is sequentially executed for all shell elements constituting the plate material.
[0186]
In the conventional shell element, it is assumed that the planarity of the cross section is maintained before and after the bending deformation. Therefore, the neutral surface of the shell element always coincides with the center plane of the plate thickness, and the cross section rotation angle θy is The value depends on the x-coordinate value representing the bending direction position, but does not depend on the z-coordinate value representing the plate thickness direction position.
[0187]
Accordingly, the cross-sectional rotation angle θy is expressed by equation (21) in FIG. 17, the vertical strain εx in the bending direction is expressed by equation (22), and the shear strain γxz in the thickness direction is expressed by equation (23). . In these equations, “u” represents a displacement in the x direction, and “w” represents a displacement in the z direction.
[0188]
Therefore, when a conventional shell element is used, as shown in FIG. 18, the vertical strain εx is distributed linearly in the plate thickness direction, while the shear strain γxz is distributed so as not to change along the plate thickness direction.
[0189]
However, according to the study by the present inventors, such an assumption is appropriate when the plate material to be bent is a thin plate, but the degree of bending of the plate material is severe as in the case of a thick plate. In some cases it turned out to lose validity.
[0190]
FIG. 19 is a graph showing how the relationship between the cross-sectional rotation angle θy and the distance ds from the inner surface of the plate changes depending on the bending severity Lb. This graph is created based on the experimental results of the present inventors. Here, the bending severity Lb is expressed by using a bending index R / t which is a value obtained by dividing the bending radius R by the plate thickness t, and the bending severity becomes stronger as the value is smaller.
[0191]
According to the graph of FIG. 19, the region where the bending severity Lb is 3 or more (the region where the bending severity is weak and in the same figure, the lower side of the horizontal line associated with the notation “Lb = 3” ), The section rotation angle θy hardly changes regardless of the position in the plate thickness direction, and it can be seen that the assumption about the conventional shell element is valid. However, in the region where the bending severity Lb is smaller than 3 (the region where bending is severe and the region above the horizontal line), the cross-sectional rotation angle θy changes significantly depending on the position in the plate thickness direction. It can be seen that the assumptions about conventional shell elements are not valid.
[0192]
Therefore, the inventors calculated the cross-sectional rotation angle θy of an arbitrary part in the thick plate by simulation. In this simulation, the plank is virtually divided into a plurality of solid elements. Simulation is performed using a solid element model.
[0193]
A program to be executed for this simulation is incorporated in the molding simulation program. In FIG. 20, the program is conceptually represented as a strain approximation function specifying routine.
[0194]
In this strain approximation function specifying routine, first, in S601, for each of a plurality of representative plate materials (potential equivalent plate materials) having different bending indices R / t, the displacement u (each of each solid element is determined by the finite element method). Component displacements u, v, w) and rotation angles θ (rotation angles θy, θy, θz of each component) are calculated. Next, in S602, the strain is calculated in association with the z coordinate value based on the calculated displacement u and the rotation angle θ and using the equations (32) and (33) in FIG. .
[0195]
Subsequently, in S603, based on the relationship between the calculated strain and the z-coordinate value, a strain approximation function that approximately defines the distribution of strain in the plate thickness direction is specified. Thereafter, in S604, the specified strain approximation function is associated with the bending index R / t of each plate material. This completes one execution of the strain approximation function specifying routine.
[0196]
It is not essential to execute the strain approximation function specifying routine for each plate material to be analyzed for springback. When a strain approximation function has already been specified for substantially the same bending severity Lb by past execution of this strain approximation function specifying routine, the strain approximation function is used by using the specified strain approximation function. This is because new execution of the function specifying routine can be omitted.
[0197]
In order to verify the accuracy of the simulation using the solid element model, the inventors measured the relationship between the cross-section rotation angle θy and the distance ds from the inner surface of the bend by experiment for a thick plate. Furthermore, the present inventors compared the experimental value representing the measured relationship with the calculated value by the above simulation.
[0198]
FIG. 21 shows the relationship between the cross-sectional rotation angle θy and the distance ds for both thick plates by both experimental values and calculated values. As is clear from the figure, there is a good agreement between the experimental values and the calculated values. Therefore, it is possible to use a calculated value instead of the experimental value for the thick plate. Therefore, even without an experiment, the relationship between the cross-sectional rotation angle θy and the distance ds, that is, the z-coordinate value of the cross-sectional rotation angle θy Dependencies can be acquired.
[0199]
FIG. 22 is a graph showing the relationship between the shear strain γxz (an example of the strain) in the plate thickness direction calculated from the calculated values by the simulation and the distance ds from the inner surface of the bending in association with each value of the bending index R / t. It is shown in The graph further shows a strain approximation function that approximates the relationship between the shear strain γxz and the distance ds when the values of the bending index R / t are 1.83, 1.61, and 2.05, respectively. Has been. This strain approximation function is an expression that defines the shear strain γxz by a function having a variable z representing the distance ds as a variable. This distortion approximation function can be configured as a polynomial which is, for example, a quadratic function, a cubic function, a quartic function, or a quintic function.
[0200]
By using this strain approximation function, it is possible to calculate the distribution of strain in the thickness direction with almost the same accuracy as in the case of using a solid element even in a molding simulation using a shell element. FIG. 23 is a graph showing a calculated value by a forming simulation using a solid element and a calculated value by a forming simulation using a shell element of the present embodiment regarding the shear strain γxz. Both are in good agreement.
[0201]
FIG. 23 also shows a calculated value by a forming simulation using a shell element of a conventional example for comparison. In the shell element of the conventional example, it is assumed that the shear strain γxz is a constant value that does not depend on the z coordinate value.
[0202]
FIG. 24 is a graph showing the distribution of the shear strain γxz calculated using the shell element of the present embodiment and the distribution of the shear strain γxz calculated using the shell element of the conventional example. The distribution of the shear strain γxz calculated using the shell element of the present embodiment shows an example when the strain approximation function is a quadratic function. As described above, when the shell element of the present embodiment is used, the actual material characteristics of the thick plate, that is, the actual distribution of the shear strain γxz, can be accurately reproduced as compared with the case of using the conventional shell element. The
[0203]
Unlike the shell element of the conventional example, for the shell element of the present embodiment, the cross-sectional rotation angle θy depends not only on the x coordinate value but also on the z coordinate value, which is expressed by equation (31) in FIG. Has been.
[0204]
Due to such dependency, for the shell element of the present embodiment, the vertical strain εx in the bending direction is expressed by the equation (32) in FIG. 25, and the shear strain γxz in the thickness direction is expressed by the equation (33). .
[0205]
In FIG. 26, the calculated value using the shell element of this embodiment is shown regarding the distribution in the plate thickness direction of each of the shear strain γxz and the normal strain εx. The distribution of the shear strain γxz shows an example when the strain approximation function is a higher order function than the second order. According to the shell element of the present embodiment, the cross-sectional rotation angle θy can be made to depend on the z coordinate value, and as a result, the actual distribution is accurately reproduced not only for the shear strain γxz but also for the vertical strain εx.
[0206]
As described above, the shell element of the present embodiment and the shell element of the conventional example are compared with each other with respect to the strain. Next, the vertical stress σz in the thickness direction is compared with each other.
[0207]
As described above, in the conventional shell element, as shown in FIG. 27, the normal stress σz is ignored and assumed to be zero. On the other hand, the normal stress σz is considered in the shell element of the present embodiment.
[0208]
The vertical stress σz in the plate thickness direction can be considered as the sum of the contribution σzs due to bending deformation of the thick plate and the contribution σzc due to the contact pressure p received by the thick plate from the punch. The contribution σzc is caused by shear deformation generated in the thick plate by the contact pressure p.
[0209]
The vertical stress σz can be analyzed from the vertical stress σx by paying attention to the balance of forces in the z direction (radial direction of the thick cylinder) of the microelements in the thick plate (thick cylinder). If the balance of the force is described in polar coordinates and discretized in terms of time, it is expressed by equation (41) in FIG. “R” in the equation (41) is a variable on the z-coordinate axis and can be replaced with “z”.
[0210]
In the present embodiment, first, the vertical stress σx is tentatively analyzed while ignoring the vertical stress σz, and based on the provisional value and using the above formula (41), the vertical stress σz in the thickness direction is calculated. A stress approximation function that approximately represents the distribution is identified. The provisional value of the vertical stress σx coincides with the vertical stress σx in the plane stress state of the thick plate.
[0211]
Specifically, the stress approximation function can be specified as follows.
[0212]
First, the basic form of the stress approximation function is selected. This basic form can be constituted by an n-order function or an exponential function. This basic shape differs depending on whether or not the contact pressure that the bending inner surface of the thick plate receives from the punch is taken into consideration.
[0213]
Equation (42) in FIG. 28 is an example of a stress approximation function that does not consider the contact pressure from the punch. In the stress approximation function expressed by the equation (42), “a” represents the z-direction position of the bent inner surface of the thick plate, and “b” represents the z-direction position of the bent outer surface of the thick plate. . The z-direction position of the neutral plane is r₀, T₀Respectively, “a” is expressed by equation (43), and “b” is expressed by equation (44).
[0214]
Therefore, the normal stress σz expressed by the above formula (42) is 0 on both surfaces of the thick plate, which indicates that the contact pressure that the bent inner surface of the thick plate receives from the punch is ignored. . A method for specifying the stress approximation function in consideration of the contact pressure will be described later.
[0215]
If the vertical stress σz is expressed by the stress approximation function f (z) expressed by the equation (42), the equation (41) is transformed into the equation (45).
[0216]
The above equation (42) can be expanded as represented by equation (46). Further, the first derivative f ′ (z) of the stress approximation function f (z) is expressed by Expression (47). By substituting these equations (46) and (47) into the above equation (45), equation (48) is obtained.
[0217]
When the equations (43) and (44) are substituted for the respective coefficients in the equation (48) and the variable substitution using the variables a0, a1, l, m is performed, the respective coefficients in the equation (48) are illustrated. They are arranged as shown in Equation (49) of 28. If the result of this arrangement is used, equation (48) in FIG. 28 is transformed into equation (50) in FIG.
[0218]
When the variable substitution shown in Expression (51) is performed, Expression (50) is transformed into Expression (52). In this equation (52), “X”, “Y”, “W”, and “Z” are known, while “l”, “m”, and “A ′” are unknown.
[0219]
By the way, in the present embodiment, each shell element is associated with a plurality of layers stacked in the thickness direction thereof, and a function (each of each shell element) representing the distribution of the vertical stress σx in the thickness direction. The position of the layer, that is, the z-coordinate value is defined as a variable). Accordingly, in this embodiment, the vertical stress σx is obtained for each z coordinate value for each shell element, and the least squares method is expressed by the equation (52) on the assumption of the relationship between the z coordinate value and the vertical stress σx. Is applied to the unknown l, m, A ′.
[0220]
Equation (53) in FIG. 30 represents the sum of squares for equation (52) in FIG. If this equation (53) is partially differentiated with respect to the unknowns l, m, A ′, equation (54) is obtained, which can be transformed into equation (55). Solving the three simultaneous equations represented by the equation (55), the unknowns l, m, and A 'are obtained, and thereby the stress approximation function f (z) is specified.
[0221]
In addition, it is possible to select the normal stress σz, the shear stress τxz, or both as variables for which the thickness direction distribution is defined by the function.
[0222]
The case where the stress approximation function f (z) is a quartic function has been described above. Next, the case of a cubic function will be described.
[0223]
If the normal stress σz is expressed by a stress approximation function f (z) expressed by the equation (56) in FIG. 31, the first derivative f ′ (z) of the stress approximation function f (z) is expressed by the equation (57). ). By substituting the equations (56) and (57) into the equation (45), the equation (58) is obtained. If this equation (58) is arranged with respect to the variable c, the equation (59) is obtained.
[0224]
When the variable substitution shown in the equation (60) is performed, the above equation (59) is transformed into the equation (61). In this equation (61), “X”, “Y”, and “Z” are known, while “c” and “A ′” are unknown. For each shell element, the unknowns c and A ′ are calculated by applying the least squares method to Equation (61) on the assumption of the relationship between the z coordinate value of each layer and the vertical stress σx of each layer. The sum of squares of equation (61) is represented by equation (62) in FIG.
[0225]
If the equation (62) is partially differentiated with respect to the unknown c, the equation (63) is obtained, which can be transformed into the equation (64). On the other hand, if the equation (62) is partially differentiated with respect to the unknown A ′, the equation (65) is obtained, which can be transformed into the equation (66). By solving the simultaneous equations of the equations (64) and (66), the unknowns c and A ′ are obtained, and thereby the stress approximation function f (z) is specified.
[0226]
The method for specifying the approximate function of the vertical stress σz without considering the contact pressure from the punch has been described above. However, when the contact pressure is considered, the basic form of the approximate function is, for example, the expression (67 in FIG. ). In that case, the normal stress σz is equal to the contact pressure p on the bent inner surface (z = a) of the thick plate, while the normal stress σz is zero on the bent outer surface (z = b). This is expressed by equation (68).
[0227]
Another example for obtaining the stress approximation function f (z) in consideration of the contact pressure p will be described. In this example, the vertical stress σz is expressed as the sum of the contribution σzc of the contact pressure p and the contribution σzs of bending deformation as shown in the equation (69) of FIG. The contribution σzc is formulated by a cubic function shown in Expression (70), for example. This contribution σzc is equal to the contact pressure p on the inner surface of the bend and 0 on the outer surface of the bend, and this is expressed by Expression (71). On the other hand, the contribution σzs is expressed by Expression (72).
[0228]
If the equation (70) is partially differentiated with respect to the variable r (equal to the variable z), the equation (73) is obtained. Therefore, in this case, the balance of force in the plate thickness direction (radial direction) in the microcylindrical element is expressed by Expression (74). By solving the equation (74) by the same method as described above, the stress approximation function f (z) that approximates the normal stress σz is obtained.
[0229]
Even if the contact pressure p from the punch is not taken into consideration in order to obtain the vertical stress σz or not, the provisional value of the vertical stress σx based on the vertical stress σz (as described above) The value in the plane stress state) is corrected. Hereinafter, this correction method will be described.
[0230]
According to the Prandtl-Reuss flow law, equation (81) in FIG. 35 is established for the deviation stress S of each component. When the normal stress σz is ignored, the equation (82) is established for the deviation stress s of each component. On the other hand, when the vertical stress σz is taken into account, the equation (83) is established for the deviation stress s of each component. “Σx” and “σy” in the equation (82) mean the x-direction normal stress and the y-direction normal stress when the normal stress σz is not considered, respectively, while the “σx ′” and “σy” in the equation (83) “” Means the x-direction normal stress and the y-direction normal stress when the normal stress σz is considered.
[0231]
Since the value on the right side of the first expression in the expression (82) matches the value on the right side of the first expression in the expression (83), the expression (84) is derived. Further, since the value on the right side of the second equation in the equation (82) and the value on the right side of the second equation in the equation (83) match each other, the equation (85) is derived.
[0232]
From these equations (84) and (85), equation (86) is derived for the normal stress σx ′, and equation (87) is derived for the normal stress σy ′. According to the equations (86) and (87), it can be seen that the normal stresses σx ′ and σy ′ considering the normal stress σz increase by the amount of the normal stress σz with respect to the normal stresses σx and σy not considered. . This increase is due to the fact that the cross-sectional force that is the sum of the vertical stresses σx ′ and σy ′ acting on the entire cross section of the thick plate (the cross-sectional force that takes into account the vertical stress σz) is the sum of the vertical stresses σx and σy. This means that the force (cross-sectional force in a plane stress state) increases by the sum of the vertical stress σz.
[0233]
Therefore, in the present embodiment, the position r in the plate thickness direction of the neutral surface is maintained so that the cross-sectional force considering the normal stress σz is maintained at the cross-sectional force in the plane stress state.₀Is changed. Specifically, the amount of movement of the neutral plane from the center position of the plate thickness is determined so that the cross-sectional force considering the vertical stress σz decreases by an increase in the vertical stress σ.
[0234]
FIG. 36 is a graph showing an example of a method for determining the movement amount of the neutral plane, taking as an example the case where the determination is performed only with respect to the normal stress σx. In the figure, the vertical stress σx corresponding to the neutral plane before movement is indicated by a solid line graph, and the vertical stress σx corresponding to the neutral plane after movement is indicated by a two-dot chain line graph.
[0235]
In FIG. 36, when the solid line graph is moved by the distance L1 in the direction in which the neutral plane approaches the bending inner surface, the value on the bending outer surface of the vertical stress σz and the distance L1 are equal on the tension side (+ side). The sectional force, which is the sum of the vertical stress σx, increases by the amount, while on the compression side (− side), the sum of the vertical stress σx is equal to the product of the value of the vertical stress σz on the inner surface of the bending and the distance L1. Some cross-sectional force is reduced. It can be considered that the increase on the tension side corresponds to the decrease on the compression side.
[0236]
Therefore, if the function of the stress on the outermost surface of the plate material is obtained and the increase amount of the normal stress σx calculated in advance is used, the distance L1 is uniquely determined, whereby the final neutral surface The position is determined.
[0237]
The example shown in FIG. 36 is an example in which the sign of the vertical stress σz is negative. Therefore, the position of the neutral surface is moved in the direction in which the cross-sectional force that is the sum of the vertical stress σx increases. That is, in this example, the decrease in the cross-sectional force, which is the sum of the vertical stress σx by considering the normal stress σz, is compensated by the movement of the neutral plane, and as a result, it is not considered when considering the normal stress σz. In some cases, the sectional force, which is the sum of the vertical stresses σx, does not change. This avoids an unnatural situation where the cross-sectional force between two shell elements sharing the same node becomes discontinuous due to the consideration of the vertical stress σz.
[0238]
In FIG. 37, the calculated value for the bending direction stress σx is represented by a graph. These calculated values are
(1) A calculated value using a shell element model, which represents a bending direction stress σx that does not consider the thickness direction stress σz;
(2) A calculated value using a shell element model, a bending direction stress σx taking into account the thickness direction stress σz;
(3) A calculated value using a shell element model, which is a bending direction stress σx that takes into account the thickness balance stress σz and force balance
Is included.
[0239]
In FIG. 38, the value of the thickness direction stress σz analyzed using the solid element model is also represented by a graph.
[0240]
FIG. 38 is a graph showing some calculated values for the bending direction stress σx. These calculated values include a calculated value using the shell element of the present embodiment and a calculated value using the shell element of the conventional example. In the same figure, for comparison, the calculated values are also shown in a graph using solid elements. This calculated value can be considered to be equivalent to the actually measured value of the bending direction stress σ. Therefore, according to the figure, the bending direction stress σx calculated using the shell element of the present embodiment reflects the measured value more accurately than the bending direction stress σx calculated using the shell element of the conventional example. I understand.
[0241]
The contents of the molding simulation program in FIG. 12 have been schematically described above, but will be specifically described below.
[0242]
In this molding simulation program, first, in S401, data necessary for molding analysis is read from the MO 32 and the analysis is initialized.
[0243]
Next, in S402, boundary conditions regarding the load acting on the plate material are applied. In this embodiment, the boundary condition regarding the load which acts on a board | plate material from a pad is applied.
[0244]
Subsequently, in S403, the stress acting on the target shell element that is the current execution target among the plurality of shell elements constituting the plate material is calculated. Details of S403 are shown in a flowchart as an element stress calculation routine in FIGS. This element stress calculation routine is designed to embody the stress update algorithm.
[0245]
In this element stress calculation routine, first, in S501, as described above, the total stiffness matrix is used in the incremental stiffness equation in a state where each component is equal to the previous value. A displacement increment u and a rotation angle increment θ are calculated. Based on these values, the strain at the thickness center plane is calculated. Next, in S502, the strain increment ε for each component based on the strain approximation function._ijIs calculated. Where strain increment ε_ijIncludes both normal and shear strains.
[0246]
Thereafter, in S503, the calculated strain increment ε_ijIsotropic strain increment ε_mIs calculated. Subsequently, in S504, the calculated isotropic strain increment ε_mBased on the hydrostatic pressure increment 3Kε_mIs calculated. Thereafter, in S505, the stress σ for each component in the previous calculation step_ij ⁰And hydrostatic pressure increment 3Kε_mAnd updated hydrostatic pressure σ as described above._m ^{t + ΔT}Is calculated.
[0247]
Subsequently, in S506, the calculated strain increment ε for each component._ijAnd isotropic strain increment ε_mAnd deviation strain increment e for each component_ijIs calculated. Subsequently, in S507, the deviation stress increment 2Ge in the plane stress state and in the elastic region._ijIs calculated.
[0248]
Thereafter, in S508, the deviation stress S for each component in the previous calculation step._ij ⁰And deviation stress increment 2Ge_ijAnd updated trial stress S as described above._ij ^{t + ΔT}Is calculated. Subsequently, in S509, the trial stress S_ij ^{t + ΔT}Based on the calculated value of^TrialIs calculated.
[0249]
Thereafter, in S510 of FIG. 14, it is determined whether or not the plate material has yielded on at least one of both surfaces thereof. If it is assumed that it has surrendered this time, the determination is yes, and the process proceeds to step S511 and subsequent steps.
[0250]
In S511, the equivalent plastic strain increment ε is obtained by using the equation (16) in FIG.^PIncrement Δε^PIs calculated. Then, in S512, the calculated increment Δε^PIs equivalent stress σ, and corresponding stress σ and trial stress S_ij ^TrialFinal solution S corresponding to_ij ^FinalIs calculated using equation (14) in FIG. In the present embodiment, the final solution S_ij ^FinalIs an example of the “provisional value” in item (6).
[0251]
In addition, when the von Mises yield condition is not adopted, the radial return algorithm can be realized by repetitive execution of S511 and S512.
[0252]
For example, at each execution of S511, the equivalent plastic strain increment ε^PIs a small increment Δε^PIs increased by On the other hand, at each execution of S512, the equivalent plastic strain increment ε^PEquivalent stress σ corresponding to the provisional value (calculated by execution of S510) is calculated as the corresponding equivalent stress σ. Furthermore, the corresponding equivalent stress σ and trial stress S_ij ^TrialFinal solution S corresponding to_ij ^FinalThe provisional value of is calculated.
[0253]
In S512, the corresponding equivalent stress σ and the final solution S_ij ^FinalIt is determined whether or not the provisional values are sufficiently consistent with each other. If there is a good agreement, the equivalent plastic strain increment ε^PProvisional value and final solution S_ij ^FinalThe provisional values are the final values. Final solution S_ij ^FinalThe final value of means the updated stress component. In contrast, the corresponding equivalent stress σ and the final solution S_ij ^FinalIf the provisional values of the two do not sufficiently match each other, the process returns to S511, and the equivalent plastic strain increment ε^PA small incremental Δε with a constant provisional value of^PIs increased by
[0254]
In any case, the equivalent plastic strain increment ε by the radial return algorithm^PAnd the updated deviation stress component S_ij ^FinalThen, in S513, the updated deviation stress component S is obtained._ij ^FinalAnd the updated hydrostatic pressure component σ_m ^{t + Δt}Based on the above, the stress σ updated for each component is obtained by using the equation (15) in FIG._ij ^{t + Δt}Is calculated.
[0255]
As described above, the case where the plate material yields on at least one of both surfaces thereof has been described. However, if the plate material does not yield on any surface, the determination in S510 is NO, and the process proceeds to S514. In S514, the trial stress S updated up to the previous calculation step._ij ^TrialAnd the updated hydrostatic pressure component σ_m ^{t + Δt}Based on the above, the stress σ updated for each component is obtained by using the equation (17) in FIG._ij ^{t + Δt}Is calculated.
[0256]
  If the plate yields, execute S511 to S513Thereafter, in S515, the sum of the bending direction stress σx (value in the plane stress state) acting on the cross section of the shell element of interest this time and the sum of the cross direction stress σy (value in the plane stress state) are Calculated as force. These cross-sectional forces have values calculated assuming that the thickness direction stress σz is zero.
[0257]
Subsequently, in S516, a stress approximation function of the thickness direction stress σz is specified by the above-mentioned least square method. Thereafter, in S517, a bending direction stress σx and a cross-sectional direction stress σy satisfying the specified stress approximation function are calculated. In the present embodiment, the calculated bending direction stress σx and cross-sectional direction stress σy are examples of the “temporary in-plane stress” in the item (6).
[0258]
In addition, this element stress calculation routine has a function of performing analysis in consideration of the plate thickness direction stress σz in advance although analysis using a shell element is performed. In this embodiment, first, assuming that the plate thickness direction stress σz is 0, the bending direction stress σx and the sectional direction stress σy are obtained, and then, from the bending direction stress σx and the sectional direction stress σy, By using the stress approximation function, the plate thickness direction stress σz is obtained. Thereafter, by using the above function, the bending direction stress σx and the cross-sectional direction stress σy are obtained again on the assumption of the obtained thickness direction stress σz, that is, taking the plate thickness direction stress σz into consideration.
[0259]
Subsequently, in S518, as described above, the neutral plane is set so that the total sum of the calculated stresses σx and σy (values in consideration of the normal stress σz) matches the calculated sectional force. Is moved. Specifically, the neutral plane position r₀Is corrected. This correction value affects the calculation of strain using the equations (32) and (33) in FIG. 25, the specification of the stress approximation function, the specification of the strain approximation function, the analysis of the springback, and the like.
[0260]
That is, in the present embodiment, the corrected neutral plane position r₀The bending direction stress σx and the cross-sectional direction stress σy defined by the above respectively constitute an example of the “final value” in the item (6). Further, the values of the bending direction stress σx and the cross-sectional direction stress σy when it is determined in S409 in FIG. 12 that the series of analysis has been completed mean the residual stress. The residual stress is stored in the memory 22 in association with each shell element.
[0261]
Thereafter, in S519, the updated plastic strain ε is obtained by using the equation (19) of FIG. 15 based on the strain increment of each component up to the previous calculation step and the current deviation strain increment.^P _ij ^{t + Δt}Is calculated. “Δλ” in the equation (19) is obtained from the equations (18) and (19). The plastic strain ε when it is determined in S409 in FIG.^P _ij ^{t + Δt}Means the residual strain. The residual strain is stored in the memory 22 in association with each shell element.
[0262]
Subsequently, in S520, the updated plastic strain ε^P _ij ^{t + Δt}Based on the thickness t of the shell element of interest this time₀Is updated. This updated value affects the specification of the stress approximation function, the specification of the strain approximation function, the analysis of the springback, and the like.
[0263]
Thereafter, in S521, it is determined whether or not the execution of S501 to S520 has been performed for all shell elements of the plate material. If it is assumed that the execution of S501 to S520 is not performed for all the shell elements, the determination is NO, and in S522, the next shell element is set as a new focused shell element, and then the process returns to S501. On the other hand, if it is assumed that the execution of S501 to S520 has been performed for all the shell elements, the determination in S521 is YES, and in S523, the finite element method is executed for the shell element model that approximates the plate material.
[0264]
Specifically, an overall stiffness matrix is constructed from the element stiffness matrix for each shell element, and simultaneous stiffness equations are solved while considering the balance between external force and internal force at each node.
[0265]
In S524, it is determined whether or not the solution of the stiffness equation has converged. If it does not converge, the determination is no, the process returns to S501 in FIG. 13, and S501 to S520 are executed using the new overall stiffness matrix. On the other hand, if the solution of the stiffness equation has converged, the determination in S524 is YES, and one execution of this element stress calculation routine is completed.
[0266]
Thereafter, in S404 of FIG. 12, the contact force acting on the plate material from the punch is calculated under consideration of the deformation of the plate material and the shape of the punch under the boundary condition.
[0267]
Subsequently, in S405, the acceleration at each node of each shell element is updated by taking into account the boundary condition regarding the machining speed of the punch. Thereafter, in S406, the acceleration of each node of each shell element is updated by taking into account the analysis conditions relating to the contact between the punch and the plate material. Subsequently, in S407, the speed of each node is updated by performing time integration on the node acceleration, and further, the displacement of each node is updated by performing time integration on the calculated speed. Is done.
[0268]
Thereafter, in S408, the analysis result is output and necessary post processing is performed. Subsequently, in S409, it is determined whether or not the current analysis is completed. One execution of the loop composed of S402 to S410 corresponds to one time increment (one calculation step) in the incremental calculation method, and in terms of phenomenon, the punch enters a small distance into the die. Corresponding to
[0269]
Therefore, in S409, it is determined whether or not the elapsed time from the first execution start time of the loop has reached a length corresponding to the predetermined number of executions of the loop, and the current analysis is completed. It is determined whether or not. If it is assumed that the analysis has not ended yet at the current time t, the determination is NO, and the updated current time t + Δt (Δt: time increment) is calculated in S410. Thereafter, the process returns to S402.
[0270]
As a result of repeating the loop consisting of S402 to S410 several times, if the current analysis is completed, the determination in S409 is YES, and one execution of this molding simulation program is completed.
[0271]
Next, the contents of the material property value determination program will be described in detail.
[0272]
FIG. 39 is a flowchart showing this material property value determination program. First, in S201, a number j assigned in order to a plurality of shell elements in which the material is virtually divided is set to 1. Next, in S202, the initial strain ε corresponding to the current shell element j is stored from the memory 22.^P ₀(J) is read out.
[0273]
Thereafter, in S203, the initial strain ε of this time^P ₀It is determined to which of the plurality of pre-strain categories (j) belongs, and the material property value corresponding to the pre-strain category determined to belong to the number j of the current shell element is associated with each other. Data is stored in the memory 22. Therefore, the spring back analysis program can analyze the spring back of the material by referring to the data and using the material characteristic value selected for each shell element.
[0274]
Subsequently, in S204, the current initial strain ε^P ₀Necessary modifications are made to (j). Specifically, the pre-strain ε strictly corresponding to the material property value selected in S203 above.^P _gIs strictly the initial strain ε^P ₀Considering the fact that it does not always match (j), if it does not match, this initial strain ε^P ₀(J) is modified. In addition, this shell element number j and the corrected initial strain ε^P ₀Data that associates (j) with each other is stored in the memory 22.
[0275]
Therefore, the spring back analysis program can analyze the spring back of the material by referring to the data and using the material characteristic value corrected for each shell element.
[0276]
Thereafter, in S205, it is determined whether or not the number j is greater than or equal to the maximum value jMAX, that is, whether or not all the material characteristic values have been selected for a plurality of shell elements of the material. If it is assumed that the number j is not greater than or equal to the maximum value jMAX this time, the determination is no, the number j is incremented by 1 in S206, and then the process proceeds to S202. Thereafter, S202 to S206 are executed for the new shell element.
[0277]
As a result of repeating the execution of S202 to S206 several times, if the determination in S205 is YES, one execution of this material property value determination program is completed.
[0278]
Next, the contents of the back stress calculation program will be described in detail.
[0279]
FIG. 40 is a flowchart showing the back stress calculation program. First, in S301, the number k given to the plurality of shell elements of the material in order is set to 1. Next, in S302, the initial stress-strain relationship corresponding to the current shell element k is read from the memory 22. This initial stress-strain relationship means the relationship between the stress and strain remaining in the current shell element k of the material immediately after press forming, and was simulated by the forming simulation program.
[0280]
The method of virtually dividing the material into a plurality of shell elements does not differ between the calculation of back stress and the determination of the material characteristic values, but for the sake of convenience of explanation, the number of each shell element is determined in the determination of the material characteristic values. Is represented by j, while it is represented by k in the calculation of the backstress α.
[0281]
Thereafter, in S303, the back stress α of the current shell element k is calculated based on the read initial stress-strain relationship. Further, data that associates the calculated back stress α with the number k of the current shell element is stored in the memory 22. Therefore, by referring to the data, the springback analysis program analyzes the springback of the material by using not only the material characteristic value but also the backstress α calculated for each shell element. Can do.
[0282]
Here, a method of calculating the back stress α will be described.
[0283]
In the press forming simulation, the material to be press formed is virtually divided into a plurality of shell elements, and a plane stress state is assumed for each shell element. Therefore, in each shell element, all stress components in the plate thickness direction are ignored as shown by the equation (101) in FIG. In this formula (101), “σ” means normal stress and “τ” means shear stress.
[0284]
In each shell element assumed in this manner, the direction of the deviation stress vector (σ′xx, σ′yy, σ′xy) and the deviation component (α′xx, α′yy, α) of the back stress α to be obtained Assuming that the direction of “xy” is the same, the equation (102) is established. In this formula (12), “A” means a constant.
[0285]
On the other hand, the equation (103) is established for the deviation stress vector σ′ij. In this equation (103), “σij” means stress (not only normal stress σ but also shear stress τ), “δij” means Kronecker delta, and “σm” means stress σij under hydrostatic pressure. The stress σm can be obtained from the equation (104).
[0286]
Further, with respect to the back stress α, the formula (105) is established. Therefore, the constant A is obtained by the equation (106) in FIG. In this equation (106), “ε” selectively represents the vertical strain ε and the shear strain γ. Specifically, when “σ” on the right side of Expression (106) represents the vertical stress σ, it represents the vertical strain ε, whereas when it represents the shear stress τ, the shear strain γ is represented.
[0287]
On the other hand, the relationship represented by the equation (107) in FIG. 42 is established by the above equation (102). Further, for the deviation component α ′ of the back stress α, the equation (108) is established in the same manner as the deviation stress vector σ ′. In this equation (108), “αij” means back stress, “δij” means Kronecker delta, and “αm” means back stress αij under hydrostatic pressure. The stress αm can be obtained from the equation (109).
[0288]
Therefore, the back stress α is calculated by using the equations (106) to (109). FIG. 43 is a graph showing the equivalent value of the back stress α to be obtained in relation to the equivalent value of the stress σ.
[0289]
After the back stress α is calculated as described above, in S304, it is determined whether or not the number k is equal to or larger than the maximum value kMAX, that is, whether or not the back stress α is calculated for all the plurality of elements of the material. Is done. If it is assumed that the number k is not greater than or equal to the maximum value kMAX this time, the determination is no, the number k is incremented by 1 in S305, and then the process proceeds to S302. Thereafter, S302 to S305 are executed for the new element.
[0290]
As a result of repeating S302 to S305 several times, if the determination in S304 is YES, one execution of the backstress calculation program is completed.
[0291]
FIG. 44 is a graph showing the effect of analyzing a springback by performing a molding simulation using the shell element model of the present embodiment. This graph shows the calculated value of the springback defined as shown in FIG. 45 and the calculation time for that. FIG. 44 shows experimental values of springback.
[0292]
FIG. 44 further shows an example using a solid element model and an example using a conventional shell element model for comparison. If the shell element model of the present embodiment is used, the calculation accuracy that is sufficiently comparable to the calculation accuracy when the solid element model is used is realized in substantially the same calculation time as when the conventional shell element model is used. .
[0293]
As is clear from the above description, in the present embodiment, step S4 in FIG. 2 constitutes an example of the “molding simulation step” in the item (1).
[0294]
Further, in the present embodiment, S501 to S514 in FIG. 13 and FIG. 14 constitute an example of the “provisional value analysis step” in the item (2) or (6), and S515 to S518 are the “final value” in the same term. It constitutes an example of “analysis process”.
[0295]
Further, in the present embodiment, S515 to S518 in FIG. 14 constitute an example of the “final value determining step” in the item (3).
[0296]
Further, in the present embodiment, S518 in FIG. 14 constitutes an example of the “step” in the item (4) or (8).
[0297]
Further, in the present embodiment, S516 in FIG. 14 constitutes an example of the “step” in the item (5) or (9).
[0298]
Further, in the present embodiment, S516 in FIG. 14 constitutes an example of the “stress approximate function specifying step” in the item (7), and S517 constitutes an example of the “in-mask stress determination step” in the term, S515 and S518 cooperate with each other to constitute an example of the “final value determining step” in the same term.
[0299]
Further, in the present embodiment, S502 and S503 in FIG. 13 and the strain approximation function specifying routine in FIG. 20 cooperate with each other to constitute an example of the “strain analysis step” in the item (12).
[0300]
Furthermore, in the present embodiment, S603 in FIG. 20 constitutes an example of the “strain approximate function specifying step” in the item (16), and S502 to S518 in FIG. 13 and FIG. It constitutes an example of “process”.
[0301]
Further, in the present embodiment, step S7 in FIG. 2 constitutes an example of the “spring back analysis step” in the item (17).
[0302]
Further, in the present embodiment, step S1 in FIG. 2 constitutes an example of the “experimental value acquisition step” in the above item (18), and step S2 constitutes an example of the “material property value acquisition step” in the same term. -ing
[0303]
As mentioned above, although one embodiment of the present invention was described in detail based on a drawing, this is an exemplification, including the mode described in the above section (Means for Solving the Problems and Effects of the Invention). It is possible to implement the present invention in the form of various modifications and improvements based on the knowledge of those skilled in the art.
[Brief description of the drawings]
FIG. 1 is a system diagram showing a springback analysis apparatus suitable for implementing a springback analysis method according to an embodiment of the present invention.
FIG. 2 is a process diagram showing the springback analysis method.
FIG. 3 is a front sectional view showing an example of a press molding machine capable of bending a plate material into a U-shape.
FIG. 4 is a graph for explaining how experimental values related to stress-strain are acquired for each pre-strain in the springback analysis method.
5 is a graph for explaining a material hardening model used for the springback analysis in step S4 in FIG. 2 and material characteristic values necessary for the springback analysis.
FIG. 6 is another graph for explaining material characteristic values necessary for the springback analysis.
FIG. 7 is a graph for explaining a material hardening model used in the springback analysis while comparing it with an isotropic hardening model and an actual stress-strain relationship.
FIG. 8 is a flowchart conceptually showing the contents of a material characteristic value acquisition program executed by a computer for acquiring material characteristic values in step S2 in FIG.
9 is a diagram showing several equations for explaining the theory and principle of obtaining material property values by the material property value obtaining program of FIG. 8. FIG.
FIG. 10 is a flowchart conceptually showing details of S106 in FIG. 8 as a material softening coefficient calculation routine.
11 is a diagram showing, in a tabular form, a state in which the relationship between pre-strain and material property value is obtained discretely by the material property value obtaining program of FIG.
12 is a flowchart conceptually showing the contents of a molding simulation program executed by a computer to perform the molding simulation in step S4 in FIG.
FIG. 13 is a part of a flowchart conceptually showing details of S403 in FIG. 12 as an element stress calculation routine.
FIG. 14 is another part of a flowchart showing the element stress calculation routine.
15 is a diagram showing several equations for explaining the theory and principle used in the element stress calculation routine of FIGS. 13 and 14. FIG.
FIG. 16 is a graph for explaining the theory and principle used in the element stress calculation routine of FIGS. 13 and 14;
FIG. 17 is a diagram showing a displacement-strain relational expression used for a conventional shell element.
18 is a diagram showing a strain distribution expressed by the displacement-strain relational expression in FIG. 17 together with bending deformation of a plate material.
FIG. 19 is a graph showing how the relationship between the cross-sectional rotation angle θy and the distance ds from the bending inner surface of the plate changes according to the bending severity Lb.
20 is a flowchart conceptually showing the contents of a strain approximation function specifying routine incorporated in a forming simulation program for executing the forming simulation in step S4 in FIG.
FIG. 21 is a graph for verifying the accuracy of the value calculated using the solid element model for the relationship between the cross-sectional rotation angle θy and the distance ds from the bent inner surface of the plate material, in comparison with the actual measurement value.
22 is a graph showing some examples of execution results of the strain approximation function specifying routine of FIG.
FIG. 23 shows three types of calculated values for the relationship between the cross-sectional rotation angle θy and the distance ds from the bent inner surface of the plate material, calculated using the solid element model and calculated using the shell element model of the above embodiment. It is a graph which shows what consists of a value and the calculated value using the shell element model of a prior art example.
FIG. 24 is a diagram showing the distribution of the shear strain γxz calculated using the shell element of the above embodiment and the distribution of the shear strain γxz calculated using the shell element of the conventional example, together with the bending deformation of the plate material.
FIG. 25 is a diagram showing a displacement-strain relational expression used in the shell element model of the embodiment.
FIG. 26 is a diagram showing the distribution in the plate thickness direction of the normal strain εx and the shear strain γxz calculated using the shell element of the embodiment, together with the bending deformation of the plate material.
FIG. 27 is a diagram showing the distribution of the normal stress σz calculated using the shell element of the embodiment and the distribution of the normal stress σz calculated using the shell element of the conventional example, together with the bending deformation of the plate material.
FIG. 28 is a diagram showing a part of several equations for explaining the theory and principle used for specifying the stress approximation function in S516 of FIG. 14;
FIG. 29 is a diagram showing another part of several equations for explaining the theory and principle used to specify the stress approximation function.
FIG. 30 is a diagram showing still another part of some equations for explaining the theory and principle used for specifying the stress approximation function.
FIG. 31 is a diagram showing still another part of some equations for explaining the theory and principle used to specify the stress approximation function.
FIG. 32 is a diagram showing still another part of some equations for explaining the theory and principle used for specifying the stress approximation function.
FIG. 33 is a diagram showing still another part of some equations for explaining the theory and principle used to specify the stress approximation function.
FIG. 34 is a diagram showing still another part of some equations for explaining the theory and principle used to specify the stress approximation function.
FIG. 35 is a diagram showing several equations for explaining the relationship between normal stress σx, σz not considering vertical stress σz and normal stress σx ′, σy ′ considering normal stress σz.
FIG. 36 is a graph for explaining the theory and principle used to move the neutral plane in S518 in FIG. 14;
FIG. 37 is a graph for explaining how the vertical stress σx differs between the case where the normal stress σz is not taken into consideration and the case where the vertical stress σz is taken into consideration, as well as the case where the cross-sectional force balance is taken into consideration.
FIG. 38 shows a case where the relationship between the bending direction stress σx and the distance ds from the bending inner surface of the plate material is calculated using the solid element model, the case where the shell element of the above embodiment is calculated, and the conventional shell element. It is a graph each shown about the case where it calculates using.
FIG. 39 is a flowchart conceptually showing the contents of a material characteristic value determination program executed by a computer for determining the material characteristic value in step S5 in FIG.
FIG. 40 is a flowchart conceptually showing a back stress calculation program executed by a computer to perform the back stress calculation in step S6 in FIG.
41 is a diagram showing a part of several formulas for explaining the theory and principle used for calculating the back stress by the back stress calculation program of FIG. 40. FIG.
FIG. 42 is a diagram showing another part of several equations for explaining the theory and principle used to calculate the backstress.
43 is a graph for explaining the execution contents of the back stress calculation program of FIG. 38;
FIG. 44 shows the calculated value and calculation time when calculating the springback of a plate using the shell element model of the above embodiment, using the solid element model and the conventional shell element model. It is a graph shown in comparison with the case.
45 is a front sectional view for explaining the definition of the springback in FIG. 44. FIG.
[Explanation of symbols]
12 Computer
16 External storage device
30 CD-ROM
32 MO
100 press molding machine
110 punches
112 dice
114 pads

Claims (14)

The plate material is pressed by bringing the punch into contact with the bent inner surface which is the other surface while supporting the bent outer surface which is one of both surfaces of the plate material to be bent and deformed by at least one of the die and the pad. A method of analyzing the deformation of the plate material during molding by computer simulation,
The shell element model modeled on the computer is used as an aggregate of a plurality of shell elements in which the plate members are arranged in the plane direction thereof, and the thickness of each shell element is considered in consideration of the thickness direction stress of each shell element. Including a molding simulation step of analyzing in-plane stress, which is stress in the plane direction , and
The molding simulation process is
For each shell element, a provisional value analysis step of analyzing a provisional value of the in-plane stress assuming a plane stress state in which the plate thickness direction stress is zero;
A forming simulation including a final value analysis step for analyzing the final value of the in-plane stress by taking account of the balance of forces in the plate thickness direction based on the analyzed provisional value for each shell element analysis method. In the final value analysis step, the final value is set such that a cross-sectional force that is a sum of the in-plane stresses on a cross section of each shell element substantially matches a cross-sectional force of a shell element adjacent to the cross section. The molding simulation analysis method according to claim 1 , further comprising a final value determining step for determining the value. In the final value determining step, on the cross-section of each shell element, the cross-sectional force that is the sum of the in-plane stresses considering the thickness direction stress is substantially the cross-sectional force that is the sum of the in-plane stresses in the plane stress state. The molding simulation analysis according to claim 2 , further comprising: determining a position of a neutral surface of each shell element so as to be maintained and determining the final value based on the determined neutral surface position. Method. The final value analysis step, and the analyzed provisional value, based on the contact pressure which the plate material during the press molding receives from the punch in the bending inner surface, we claim 1 comprising the step of analyzing the final value 3 The molding simulation analysis method according to any one of the above. The plate material is pressed by bringing the punch into contact with the bent inner surface, which is the other surface, while supporting the outer surface, which is one of both surfaces of the plate material to be bent, by at least one of the die and the pad. A method of analyzing the deformation of the plate material during molding by computer simulation,
The shell element model modeled on the computer is used as an aggregate of a plurality of shell elements in which the plate members are arranged in the plane direction thereof, and the thickness of each shell element is considered in consideration of the thickness direction stress of each shell element. Including a molding simulation step of analyzing in-plane stress, which is stress in the plane direction, and
The molding simulation process is
For each shell element, a provisional value analysis step of analyzing a provisional value of the in-plane stress assuming a plane stress state in which the plate thickness direction stress is zero;
For each shell element, based on the analyzed provisional value, identify a stress approximation function that approximately defines the distribution of the thickness direction stress in the thickness direction, and use the identified stress approximation function A final value analysis step of analyzing the final value of the in-plane stress;
Molding simulation analysis method including The final value analysis step
A stress approximation function specifying step for specifying the stress approximation function based on the analyzed provisional value;
Using the identified stress approximation function, an in-plane stress determination step for determining an in-plane stress established together with the plate thickness direction stress as an in-plane stress,
Based on the determined in-plane stress, the cross-sectional force, which is the sum of the in-plane stresses on the cross-section of each shell element, substantially matches the cross-sectional force of the shell element adjacent to the cross-section. The molding simulation analysis method according to claim 5 , further comprising: a final value determining step for determining the final value. In the final value determining step, on the cross-section of each shell element, the cross-sectional force that is the sum of the in-plane stresses considering the thickness direction stress is substantially the cross-sectional force that is the sum of the in-plane stresses in the plane stress state. The molding simulation analysis according to claim 6 , further comprising: determining a position of a neutral surface of each shell element so as to be maintained, and determining the final value based on the determined neutral surface position. Method. The final value analysis step, and the analyzed provisional value, the based on the contact pressure applied from the punch in the bending inner surface the sheet is during the press molding, 5 to claim comprising the step of identifying the stress approximation function 8. The molding simulation analysis method according to any one of 7 above. It said stress approximation function, the per shell element, forming simulation analysis method according to any one of claims 5 to 8 to formulate the distribution assuming the balance of forces in the plate thickness direction. Wherein each shell element is defined by a plurality of nodes, and the plurality of nodes is, according to any one of claims 1 free of the corresponding not essential nodes to identify shell elements geometrically 9 Molding simulation analysis method. The plate material is pressed by bringing the punch into contact with the bent inner surface, which is the other surface, while supporting the outer surface, which is one of both surfaces of the plate material to be bent, by at least one of the die and the pad. A method of analyzing the deformation of the plate material during molding by computer simulation,
The shell element model modeled on the computer is used as an aggregate of a plurality of shell elements in which the plate members are arranged in the plane direction thereof, and the thickness of each shell element is considered in consideration of the thickness direction stress of each shell element. Including a molding simulation step of analyzing in-plane stress, which is stress in the plane direction, and
The forming simulation process is performed on the computer as an aggregate of a plurality of solid elements in which equivalent plate materials equivalent to the plate materials are arranged in both the surface direction and the plate thickness direction with respect to the severity of bending by the press forming. The plate of the equivalent plate has a cross section rotation angle, which is an angle at which each point on the plane cross section perpendicular to the longitudinal axis of the equivalent plate before the bending deformation is rotated by the bending deformation. A forming simulation analysis method including a strain analysis step of analyzing the strain of the equivalent plate by using the analysis result of the cross-sectional rotation angle by analyzing the position so as to depend on the position in the thickness direction. The strain analysis step includes a strain approximation function specifying step for specifying a strain approximation function that approximately defines a distribution in the plate thickness direction of the strain based on the analysis result of the cross-sectional rotation angle for the equivalent plate material,
In-plane stress analysis in which the forming simulation step uses the shell element model for each shell element together with the specified strain approximation function and analyzes the in-plane stress in consideration of the plate thickness direction stress. The molding simulation analysis method according to claim 11 , comprising a step. A program executed by a computer in order to carry out the molding simulation analysis method according to any one of claims 1 to 12 . A recording medium on which the program according to claim 13 is recorded so as to be readable by a computer.

Images