JP4371985B2 - Stress analysis method, program, and recording medium - Google Patents

Description

  The present invention relates to a technique for analyzing bending deformation of a plate material by press molding by a computer, and in particular, by using a shell element model that is an aggregate of a plurality of shell elements that approximately represents the plate material. The present invention relates to a technique for analyzing stress in a shell element.

  There is press forming as a technique for forming a plate material by plastic working using a mold such as a punch or a die. This press molding is used to manufacture various products, and an example of such a product is an automobile. In the automobile manufacturing process, press molding is performed to manufacture various parts such as a body panel, a frame, a door panel, and a suspension arm.

  In this press molding, generally, a bending outer surface which is one of both surfaces of a plate material to be bent and deformed is supported by at least one of a die and a pad, and a punch is brought into contact with the bending inner surface which is the other surface in a pressurized state. Is done by letting During the press molding, stress due to bending, stress due to contact, and stress due to frictional force are generated in the plate material.

  As a cause of poor shape accuracy of a product manufactured by this press molding, there is a spring back of the product after press molding. This phenomenon of springback is a phenomenon that occurs in plastic working such as press molding, and when the material is unloaded after being loaded, the material elastically recovers and its shape changes. In this case, plastic deformation may be partially accompanied.

  The amount of this spring back can be improved by adjusting the press molding conditions including the shape of the mold. The press molding conditions include, for example, each shoulder R of the punch and the die and radial and width clearances between the punch and the die.

  Therefore, in order to improve the shape accuracy of the product by press molding, it is effective to design the die by accurately predicting the amount of spring back. The larger the springback amount of the plate material, the greater the springback amount of the plate material due to the reason that the plate thickness t of the plate material is large or the bending radius R of the plate material is large, and the result is fed back to the mold design. It is strongly requested to do.

  In order to accurately predict the amount of spring back, it has already been analyzed by numerical simulation using a computer that the plate material is deformed by press forming. The numerical simulation is also called molding simulation.

  In the 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.

  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.

  The conventional molding simulation analysis using shell elements cannot take into account the stress generated in the plate thickness direction during press molding. However, as in the case of thin plate molding, the severity of bending of the plate (plate thickness of the plate) If the ratio of the bending radius R of the plate to t is relatively weak, the deformation of the plate can be calculated with high accuracy. However, the present inventors have found that it is difficult to accurately predict the deformation of the thick plate by the molding simulation using the conventional shell element.

  Therefore, the present inventors consider the thickness direction stress of each shell element by using a shell element model modeled on a computer as an aggregate of a plurality of shell elements in which plate materials are arranged in the plane direction. In addition, the in-plane stress, which is the stress in the plane direction of each shell element, was analyzed, and a technique to analyze the deformation of the plate using the analysis result of the in-plane stress was proposed. This technique is disclosed in Patent Document 1.

  According to this disclosed technique, deformation of a plate material is analyzed using a shell element model, so that it is easy to shorten the calculation time for analyzing the deformation of a plate material by simulation as compared with the case of using a solid element model. It becomes. Moreover, according to this disclosed technique, unlike the case of analyzing the deformation of the plate material using the above-described conventional shell element model, the analysis is performed in consideration of the stress in the plate thickness direction. This makes it easier to approach when using a solid element model. Thus, according to this disclosed technique, it is easy to achieve both reduction in calculation time and improvement in analysis accuracy.

  Hereinafter, the contents of this disclosed technique will be described in more detail. For convenience of description, the direction in which the neutral surface extends in the plate material, that is, the bending direction is referred to as the x direction, and the plate thickness direction is referred to as the z direction. This direction is referred to as the y direction (also referred to as the cross-sectional 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.

  When using the shell element model, it is difficult to directly calculate the plate thickness direction stress σz, which is the stress in the plate thickness direction of the plate material. For example, assuming a microelement that forms a partial cylinder, If the balance of force in the plate thickness direction (radial direction) in the plate material to be bent and deformed is expressed by an equation, the plate thickness direction stress σz and the bending direction stress σx which is the stress in the plate bending direction (in the cross-sectional direction of the plate material) It is possible to induce a relationship with the stress in the cross-sectional direction σy, which is a stress. If the induced 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.

  Based on such knowledge, in the above-described disclosed technology, a provisional value of in-plane stress is analyzed for each shell element assuming a plane stress state in which the plate thickness direction stress is 0. Based on the analyzed provisional value, the final value of the in-plane stress is analyzed by considering the balance of forces in the thickness direction.

  Therefore, according to this disclosed technique, the in-plane stress calculated using the shell element model is effective even though the shell element model is difficult to directly calculate the thickness direction stress. By using this, it is possible to calculate the in-plane stress in consideration of the thickness direction stress.

  However, for each shell element, the provisional value of the in-plane stress, i.e., the in-plane stress of the in-plane stress, resulting in a value different from the in-plane stress calculated by ignoring the thickness direction stress in the plane stress state. When the final value is determined, the cross-sectional force that is the sum of the in-plane stresses on the cross section of each shell element may be different from the cross-sectional force of the shell element adjacent to the cross section. In this way, the situation in which the balance of force related to the cross-sectional force is disturbed between adjacent shell elements is contrary to the balance of forces in the in-plane direction at the nodes shared by the adjacent shell elements.

  Based on such knowledge, the present inventors believe 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. In addition, a technique to determine the final value of in-plane stress was proposed. This technique is also disclosed in Patent Document 1.

  According to this disclosed technique, the continuity of the cross-sectional force is maintained between adjacent shell elements even though the in-plane stress calculated in the plane stress state is changed in consideration of the thickness direction stress. Is possible.

  Furthermore, the inventors of the present invention substantially reduced the cross-sectional force, which is the sum of the in-plane stresses considering the thickness direction stress, to the cross-sectional force, which is the sum of the in-plane stresses in the plane stress state, on the cross-section of each shell element. A technique is also proposed in which the position of the neutral plane of each shell element is determined so as to be maintained, and the final value of the in-plane stress is determined based on the determined position of the neutral plane. This technique is also disclosed in Patent Document 1.

  According to this disclosed technique, the positions of the neutral surfaces are 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, the adjacent shell elements It becomes possible to maintain the continuity of the cross-sectional force between them.

Specifically, according to this disclosed technique, first, the provisional value of the in-plane stress is calculated without considering the thickness direction stress, and then the calculated provisional value of the in-plane stress is used. The final value of the in-plane stress is calculated in consideration of the plate thickness direction stress. The calculated value of the cross-sectional force (also referred to as “axial force”), which is the sum of the cross-sections of the plate material, is the final value of the calculated in-plane stress. Increase or decrease. Here, “the calculated value of the cross-sectional force increases or decreases in the bending direction” means that the calculated values of the cross-sectional force (axial force) are different between adjacent cross sections in the same plate material. However, according to this disclosed technique, an increase or decrease in the calculated value of the cross-sectional force is eliminated by moving an appropriate amount of the position of the neutral surface set in the plate material in order to define the elastic-plastic characteristics of the plate material. The
JP 2004-42098 A

  However, the present inventors have realized that there is room for improvement in this disclosed technique. This will be specifically described below.

  As a specific example of the disclosed technology, the present inventors have made it possible to calculate the amount of movement of the neutral plane in a short time, and the distribution of the in-plane stress in the plate material in the plate thickness direction by a prediction formula. A first method was developed to calculate the amount of movement of the neutral plane in a fixedly defined state.

  However, when this first method is adopted, the cross-sectional force is accurately corrected by the movement of the neutral surface in press forming in which the plate material is bent and deformed while applying a large tension, such as hat bending. There is a limit. For this reason, when this first method is adopted, the calculation accuracy of the in-plane stress is lowered, and as a result, the prediction accuracy of the springback amount may be lowered.

  Furthermore, as another specific example of the above-described disclosed technology, the present inventors have made it possible to calculate the amount of movement of the neutral surface as a variable parameter in order to allow the amount of movement of the neutral surface to be calculated with high accuracy. A second method for calculating the in-plane stress distribution after calculating the amount of movement of the neutral plane by convergence calculation for optimization calculation so that the continuity of the cross-sectional force between adjacent shell elements is maintained. The method was developed.

  When this second method is adopted, it is possible to accurately calculate the amount of movement of the neutral plane and correct the cross-sectional force with high accuracy, but this may increase the calculation time required. There is sex.

  In the background described above, the present invention uses a shell element model that is an aggregate of a plurality of shell elements that approximately represents a plate material in order to analyze the bending deformation of the plate material by press molding by a computer. Thus, in the technique of analyzing the stress in each shell element, it is an object to improve the accuracy of the analysis while shortening the calculation time required for the analysis.

  The following aspects are obtained by the present invention. Each aspect is divided into sections, each section is given a number, and is described in a form that cites other section numbers as necessary. This is to facilitate understanding of some of the technical features that the present invention can employ and combinations thereof, and the technical features that can be employed by the present invention and combinations thereof are limited to the following embodiments. Should not be interpreted. That is, it should be construed that it is not impeded to appropriately extract and employ the technical features described in the present specification as technical features of the present invention although they are not described in the following embodiments.

  Further, describing each section in the form of quoting the numbers of the other sections does not necessarily prevent the technical features described in each section from being separated from the technical features described in the other sections. It should not be construed as meaning, but it should be construed that the technical features described in each section can be appropriately made independent depending on the nature.

(1) Use of a shell element model modeled on the computer as an aggregate of a plurality of shell elements in which the plate material is arranged in the surface direction in order to analyze the bending deformation of the plate material by press molding by a computer. According to the stress analysis method of analyzing the in-plane stress that is the stress in the surface direction of each shell element after considering the thickness direction stress that is the stress in the thickness direction of each shell element,
For each shell element, a provisional value calculation step of calculating a provisional value of the in-plane stress assuming a plane stress state in which the sheet thickness direction stress is zero;
Based on the calculated provisional value of the in-plane stress, a stress approximation function specifying step for specifying a stress approximation function that approximately defines a distribution of the plate thickness direction stress in the plate thickness direction;
By using the identified stress approximation function, the in-plane stress that is mechanically established together with the plate thickness direction stress at each position in the plate thickness direction is defined as the in-plane stress in consideration of the plate thickness direction stress. A step of calculating the in-mask stress to be calculated;
Based on the distribution of the plate thickness direction stress in the plate thickness direction, the sum of the in-plane stress on the cross section of each shell element due to the consideration of the plate thickness direction stress when calculating the in-plane stress A cross-sectional force increment calculating step for calculating the amount of increase in cross-sectional force as a cross-sectional force increment;
A distribution amount determining step for determining a distribution amount for distributing the calculated cross-sectional force increment to each of a plurality of positions arranged in the plate thickness direction;
A final value calculation step of calculating a final value of the in-plane stress by correcting the calculated in-plane stress by the determined distribution amount for each position in the plate thickness direction. Method.

  In this method, a provisional value of in-plane stress is calculated for each shell element assuming a plane stress state in which the plate thickness direction stress is zero. Based on the calculated provisional value of the in-plane stress, a stress approximation function that approximately defines the distribution of the plate thickness direction stress in the plate thickness direction is specified. By using the specified stress approximation function, for each position in the plate thickness direction, the in-plane stress that is dynamically established along with the plate thickness direction stress is calculated as the in-plane stress considering the plate thickness direction stress. .

  In this method, based on the distribution of the thickness direction stress in the thickness direction, the sum of the in-plane stress on the cross-section of each shell element due to the consideration of the thickness direction stress when calculating the in-plane stress. The amount by which the cross-sectional force increases is calculated as the cross-sectional force increment. The cross-sectional force increment is generated in consideration of the thickness direction stress. The distribution amount of the cross-sectional force increment is determined so that the calculated cross-sectional force increment is distributed to a plurality of positions arranged in the plate thickness direction.

  In this method, the final value of the in-plane stress is further calculated by correcting the calculated in-plane stress with the determined distribution amount for each position in the thickness direction.

  In this method, the determination of the distribution amount of the cross-sectional force increment and the distribution amount were used so that the cross-sectional force increase due to the consideration of the plate thickness direction stress did not occur nonuniformly between the adjacent shell elements. It is possible to correct the in-mask stress.

  Therefore, according to this method, the in-plane stress can be analyzed so that the balance of forces related to the in-plane stress between the adjacent shell elements is substantially established before and after considering the thickness direction stress. It becomes easy.

  According to this method, the correction of the in-plane stress can be performed by simple calculation as compared with the case where the neutral surface of the plate material is moved according to the above-described disclosed technique. It is easy to shorten the total time required for analysis. Therefore, according to this method, it becomes easy to shorten the calculation time for the analysis while ensuring the analysis accuracy of the in-plane stress.

  One example of the use of this method is to analyze the deformation of the plate material during press molding (including the shape of the plate material after press molding) by simulation using a computer based on the analysis result of in-plane stress. This method can be applied to other uses. One example of the use of the analysis result of the simulation is to analyze the springback of the plate material based on the analysis result, but the analysis result can be used for other purposes.

  An example of the use of the analysis result of the springback is to optimize the conditions of press forming performed on the plate material, but the analysis result can be used for other purposes. The press molding conditions include, for example, the shape of a mold used for the press molding.

  The methods according to this section and the following sections are of course effective when the plate is a thick plate, but can be applied not only to the stress analysis of a thick plate but also to the stress analysis of a thin plate. is there.

  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.

(2) The final value calculating step corrects the calculated in-plane stress by subtracting the determined distribution amount from the calculated in-mask stress for each position in the plate thickness direction. The stress analysis method according to item (1).

  Regarding the in-plane stress at each position in the plate thickness direction, in terms of plastic mechanics, the value calculated taking into account the plate thickness direction stress is greater than the value calculated without considering the plate acting at the same position. Increases by the same amount as the thickness direction stress. On the other hand, the balance of forces in the bending direction or the cross-sectional direction of the plate material is established for the sum of in-plane stresses acting on the entire cross section, that is, the cross-sectional force, between the cross-sections adjacent to each other in the plate material.

  Therefore, in order to maintain the balance of force with respect to the cross-sectional force before and after considering the plate thickness direction stress, from the cross-sectional force that is the sum of each cross section of the in-plane stress considering the plate thickness direction stress. The cross-sectional force increment, which is the sum of the same thickness in each cross-section, may be subtracted.

  Correcting the position of the neutral surface of the plate according to the above-described disclosed technique is effective for maintaining the balance of forces in the bending direction or the cross-sectional direction of the plate, before and after considering the thickness direction stress. However, while the position of the neutral plane is fixed to the position of the center plane that is geometrically determined by the cross-sectional shape of the plate, the in-plane stress can be corrected by using the cross-sectional force increment. It is effective to maintain the balance of force in the thickness before and after considering the thickness direction stress.

  However, the increase in cross-sectional force is the increase in in-plane stress for the entire cross-section, not the increase in in-plane stress at each position in the plate thickness direction. The in-plane stress at each position in the direction cannot be accurately corrected.

  However, the distribution of the stress in the thickness direction in the thickness direction is based on the condition that the balance of force in the bending direction or the cross-sectional direction of the plate material is maintained at least with respect to the cross-sectional force, and that the stress in the thickness direction is actually in the thickness direction. It is possible to predict so that the condition of showing a distribution as close as possible to the distribution shown is established.

  Furthermore, in accordance with the predicted distribution of stress in the plate thickness direction, the distribution amount for allocating the cross-sectional force increment to each of the plurality of positions in the plate thickness direction is determined, and the determined distribution amount is subtracted from the in-mask stress, The in-mask stress is corrected with high accuracy.

  Based on such knowledge, in the method according to this section, the calculated in-mask stress is calculated by subtracting the determined distribution amount from the calculated in-mask stress for each position in the plate thickness direction. Is corrected.

  Therefore, according to this method, the in-mask stress can be corrected by a relatively simple calculation of subtracting the distribution amount of the sectional force increment from the in-mask stress.

(3) In the item (1) or (2), the distribution amount determination step determines the distribution amount so that the calculated cross-sectional force increments are distributed substantially evenly to the plurality of positions, respectively. The stress analysis method described.

  In the stress analysis method according to the item (1) or (2), the cross-sectional force is so distributed that the calculated cross-sectional force increments are not evenly distributed to a plurality of positions aligned in the plate thickness direction. It is possible to implement in a manner that determines the amount of incremental distribution. In this aspect, for example, the distribution amount can be determined so as to change linearly according to the position in the plate thickness direction, or the distribution amount can be determined so as to change nonlinearly.

  On the other hand, in the method according to this section, the distribution amount of the cross-sectional force increment is such that the calculated cross-sectional force increments are distributed substantially evenly to each other at a plurality of positions arranged in the plate thickness direction. It is determined.

  Therefore, according to this method, compared with the case where the calculated cross-sectional force increment is distributed so as not to be equal to each other at a plurality of positions arranged in the plate thickness direction, the distribution amount of the cross-sectional force increment is determined. Therefore, it becomes easy to simplify the calculation required for the purpose and shorten the time required for the calculation.

  By the way, for example, it is performed under the condition that the bending severity of the plate material is strong such that the bending index obtained by dividing the bending radius R [mm] of the plate material by the plate thickness t [mm] of the plate material is smaller than 3. Assuming press forming, the plate material reaches the plastic region in almost the entire region in the plate thickness direction by the press forming. On the other hand, when the gradient of the graph representing the relationship between the strain and the thickness direction stress is compared between the plastic region and the elastic region, the plastic region is more gradual than the elastic region. This means that the dependence of the plate thickness direction stress on the plate thickness direction position is smaller in the plastic region than in the elastic region.

  Therefore, when it is necessary to analyze the in-plane stress of a plate material obtained by press molding performed under the condition that the bending of the plate material is severe as described above, it is necessary to analyze the stress in the plate thickness direction. It is reasonable to approximate the distribution in the thickness direction with a uniform distribution. Therefore, it is appropriate to determine the distribution amount of the above-described cross-sectional force increments so as to be uniform in the plate thickness direction.

  Based on the knowledge described above, in the method according to this section, the calculated cross-sectional force increment is so distributed that the calculated cross-sectional force increments are substantially evenly distributed to each other at a plurality of positions aligned in the plate thickness direction. The distribution amount is determined.

(4) The stress analysis method according to any one of (1) to (3), wherein the distribution amount determination step determines the distribution amount so as not to depend on a position in the plate thickness direction.

  In the stress analysis method according to the item (1) or (2), the calculated cross-sectional force increment depends on the position in the plate thickness direction at each of a plurality of positions arranged in the plate thickness direction, that is, the position It is possible to implement in such a manner that the distribution amount of the sectional force increment is determined so as to be distributed differently depending on the condition. In contrast, in the method according to this section, the calculated cross-sectional force increment is distributed to a plurality of positions aligned in the plate thickness direction so as not to depend on the positions in the plate thickness direction. An incremental allocation amount is determined.

  Therefore, according to this method, compared with the case where the calculated cross-sectional force increment is distributed to a plurality of positions arranged in the plate thickness direction so as to depend on the position in the plate thickness direction, the cross-sectional force increment is increased. It becomes easy to simplify the calculation necessary for determining the allocation amount and reduce the time required for the calculation.

(5) The cross-sectional force increment calculating step calculates the cross-sectional force increment by substantially integrating the stress distribution function with respect to the position in the plate thickness direction. The stress analysis method described.

  In the method according to any one of the items (1) to (4), the stress approximation function in the item (1) is used for the distribution in the plate thickness direction of the plate thickness direction stress in order to calculate the sectional force increment. Instead, it can be implemented in a manner that defines and calculates the cross-sectional force increment based on the distribution thus defined.

  On the other hand, in the method according to this section, the stress approximation function is used not only for calculating the in-mask stress but also for calculating the sectional force increment. That is, in this method, the sectional force increment is calculated by substantially integrating the stress distribution function with respect to the position in the plate thickness direction. The calculation of the sectional force increment is performed for each shell element, for example.

  Therefore, according to this method, since it is possible to calculate the in-mask stress and the cross-sectional force increment using a function common to them, compared to the case of calculating using different functions, the calculation Design burdens such as program design.

(6) The press molding is performed under the condition that the bending index obtained by dividing the bending radius R [mm] of the plate material by the plate thickness t [mm] of the plate material is smaller than 3. (1) to (5 ) The stress analysis method according to any one of items.

  As a result of experiments, the inventors of the present invention can obtain the method according to any one of (1) to (5) by dividing the bending radius R [mm] of a plate by the thickness t [mm] of the plate. It was confirmed that even a stress generated in a plate material during press forming performed under a condition where the bending index is smaller than 3 can be analyzed with high accuracy.

  Based on the results of this experiment, the method according to this section is performed to analyze the stress generated in the plate material during press forming performed under the condition that the bending index is smaller than 3.

(7) The stress analysis method according to any one of (1) to (6), wherein the press forming includes at least one of thick plate forming, hemming, and hat bending.

  The “thick plate forming” in this section is performed, for example, for forming a plate material having a plate thickness of about 2 [mm] or more. In addition, the “hemming process” is performed, for example, in order to connect an inner panel and an outer panel of an automobile door to each other at edges, and generally, a portion protruding from the edge of the inner panel is sandwiched between the inner panels. Is done to curl. “Hat bending” is an example of deep drawing.

(8) A program executed by a computer to implement the method according to any one of (1) to (7).

  If this program is executed by a computer, the same effect as the method according to any one of the above items (1) to (7) can be realized.

  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.

  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.

(9) A recording medium on which the program according to item (8) is recorded so as to be readable by a computer.

  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 (7) can be realized.

  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 of unremovable storage and the like can be employed.

  Hereinafter, one of more specific embodiments of the present invention will be described in detail with reference to the drawings.

  FIG. 1 shows a springback analysis device (hereinafter simply referred to as “analysis device”) suitable for carrying out a springback analysis method including a stress analysis method according to an embodiment of the present invention.

  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.

  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.

  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.

  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.

  FIG. 2 is a process diagram showing this springback analysis method.

  If it demonstrates roughly, this springback analysis method is implemented in order to analyze the springback which arises in the board | plate material, after pressing the thick plate 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.

  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.

  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. It arrange | positions so that W may be pinched | interposed.

  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.

  In this springback analysis method, 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, and the acquired values are obtained. Based on a plurality of experimental values, a material characteristic value necessary for analyzing the springback of the material is 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.

  Further, in this springback analysis method, in order to analyze the stress-strain relationship remaining in the material immediately after the press molding, a molding simulation for simulating deformation of the material by press molding is performed by using the finite element method. This is done by executing a molding simulation program.

  In this springback analysis method, after the execution of the forming simulation program, the springback of the material is analyzed by executing the springback analysis program. In this analysis, the finite element method is used,

(1) a material characteristic value acquired by executing the material characteristic value acquisition program;
(2) Residual stress and residual strain calculated by execution of the forming simulation program, a neutral surface position for each element of the plate, and a plate thickness for each element of the plate,
(3) Back stress calculated by executing the back stress calculation program

Is used.

  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.

  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.

  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.

  Next, in step S2 in FIG. 2, the computer 12 executes the material characteristic value acquisition program, whereby a 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.

  Thereafter, in step S3, data necessary for executing the molding simulation, that is, data for numerical analysis is created by the operator. The data for numerical analysis includes data for virtually dividing the plate material to be press-molded into a plurality of shell elements, and the load applied to the plate material during press molding (the load applied to the plate material from the punch 110 and the pad 114). Including a boundary condition regarding a plurality of shell elements.

  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.

  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.

  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.

  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.

  Thereafter, in step S5, the material characteristic value determination program is executed by the computer 12, thereby using the plurality of material characteristic values acquired in step S2 to determine the plasticity associated with the molding of each shell element of the plate material. 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.

  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.

  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 forming simulation program, neutral surface position r ₀ and sheet thickness t ₀ ;
(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.

  An example of this springback analysis program is manufactured by the Japan Research Institute and sold under the name “JOH-NIKE3D”.

  Thus, one execution of the springback analysis method shown in FIG. 2 is completed.

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, how the plastic strain epsilon ^P and the stress σ is changed it is indicated by a parallelogram.

Here, if various symbols in FIG. 5 are described, “Y ₀ ” indicates an initial yield stress that is a yield stress when the material starts the first plastic deformation. “Y ₁ ” indicates a re-yield stress that is a stress when the material re-yields. “H” indicates a plastic hardening coefficient. The plastic hardening coefficient H is expressed by tan φ when the gradient of a straight line graph (third straight line 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 φ. “Β” indicates a material softening coefficient that defines the degree to which the re-yield stress Y ₁ decreases with respect to the initial yield stress Y ₀ . “Ε ^P ₀ ” indicates an initial strain that is a plastic strain ε ^P at the start of unloading.

  The springback analysis program can select an appropriate material hardening model from the isotropic hardening model, the kinematic hardening model, and the composite hardening model as the material hardening model that approximates the stress-strain relationship.

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 Also, 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 Further, when 0 <β <1, the 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

  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.

On the other hand, the stress - region Bauschinger effect appears among the graph showing strain relationship is plastic from the unloading start point P _S and the elastic range up again yield point P _Y, to reverse load end point from the re yield point P _Y It is an area.

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 acquired geometric characteristics, the initial yield stress Y ₀ , the re-yield stress Y ₁ , the plastic hardening coefficient H, and the material directly on the specimen, indirectly on the material. A softening factor β is obtained.

Specifically, re-yield stress Y ₁ is obtained as a stress σ at the break point of the polygonal line. By using the fact that the third straight line portion 54 is parallel to the fourth straight line portion 56, the plastic hardening coefficient H is obtained by obtaining the gradient of the fourth straight line portion 56 at an angle φ and calculating its tan φ. To be acquired. The material softening coefficient β is an initial yield stress Y ₀ , a re-yield stress Y ₁ , a plastic hardening coefficient H, and a pre-strain ε ^P _g .

Y ₁ = −Y ₀ + H · (1-2β) · ε ^P _g

It is obtained by substituting into the following equation and solving for β. The initial yield stress Y ₀ is obtained based on the experimental values at the positive load.

This springback analysis program selects material elastic modulus E (Young's modulus E and shear elastic modulus G) in addition to initial yield stress Y ₀ , re-yield stress Y ₁ , plastic hardening coefficient H and material softening coefficient β. It is designed to analyze the springback by using That is, in this springback analysis program, these initial yield stress Y ₀ , re-yield stress Y ₁ , plastic hardening coefficient H, material softening coefficient β, and elastic coefficient E constitute “material characteristic values”. .

However, what it is taken on the horizontal axis of the graph in FIG. 5, a rather total strain epsilon is the sum of the elastic strain and plastic strain, a plastic strain epsilon ^P. 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.

  Therefore, when the gradient of the first straight line portion 50 in FIG. 6 is expressed by the 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.

  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.

The stress-strain relationship changes depending on the initial strain ε ^P ₀ of the material (corresponding to the pre-strain ε ^P _g of the test piece). The geometrical feature of the parallelogram that approximately represents the stress-strain relationship changes with the initial strain ε ^P ₀ , and thus the material property value also changes with the initial strain ε ^P ₀ . Therefore, in this embodiment, as described above, the tension-compression test is repeated each time the prestrain ε ^P _g of the test piece is changed.

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 is defined so that the shape changes according to the prestrain ε ^P _g . Note that the parallelogram representing the linear kinematic hardening model used in the present embodiment is omitted in FIG.

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 shape depends on the expected initial strain ε ^P ₀ for each shell element of the material. Different material hardening models are used. However, even if the initial strain ε ^P ₀ expected for each shell element is different, the material hardening model is maintained expressed as a parallelogram.

  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 and the region that undergoes plastic deformation coincide with each other with sufficiently high accuracy.

  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.

  As shown in FIG. 7, the broken line is acquired so as to coincide as much as possible with a graph representing an actual stress-strain relationship during unloading and reverse loading. 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.

  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.

  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.

  Here, the contents of the material property value acquisition program will be described in detail.

  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 an experimental value is taken into the memory 22 of the computer 12 from the MO 32.

Next, in S102, the number i to be given to the pre-strain ε ^P _g in order is set to 1. Thereafter, in S103, data representing the stress-strain relationship related to the current pre-strain ε ^P _g (i) is fetched from the memory 22.

Subsequently, in S104, a graph representing the stress-strain relationship related to the current prestrain ε ^P _g (i) is assumed based on the captured data, 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 slope 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 φ, so that the plastic hardening coefficient corresponding to the current pre-strain ε ^P _g (i) is obtained. H (i) (= tanφ) is calculated.

Subsequently, in S106, the material softening coefficient β (i) corresponding to the current pre-strain ε ^P _g (i) is calculated by using the above-described equation represented by Equation (1) in FIG. 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.

When the material softening coefficient β is calculated using the above formula (1), the first condition shown by the formula (2) in FIG. 9 and the condition shown by the formula (3) are satisfied for all the prestrains ε ^P _g. The second condition is adopted as the constraint condition. In addition, “Y ₁ ” and “H” in the above formula (1) are obtained by experimental values.

As for the re-yield stress Y ₁ , it has been confirmed by various experiments that it shows a monotonic increase with respect to the pre-strain ε ^P _g , and when the pre-strain ε ^P _g = 0,

Y ₁ = −Y ₀

And for all prestrains ε ^P _g ,

Y ₁ > −Y ₀

The relationship is established.

  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).

  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.

In the present embodiment, it is determined whether or not the condition expressed by the above formula (7) is satisfied. If the condition is not satisfied, the current pre-strain ε ^P _g (i) is corrected so as to be satisfied. . For example, this correction may be performed by changing the current prestrain ε ^P _g (i) to

1.01 × (−Y ₁ / H)

Or by selecting the smallest of the experimental values that satisfy the condition shown in the above equation (7).

If the condition shown by the above equation (7) is satisfied, the present prestrain ε ^P _g (i) corrected as necessary, and the initial yield stress Y ₀ and re-yield stress Y ₁ obtained from experimental values. Then, the material softening coefficient β (i) of this time is calculated by substituting the plastic hardening coefficient H (i) into the above formula (1).

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. When the first condition is satisfied, the calculated material softening coefficient β (i) is set as a 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 also corrected accordingly.

The material softening coefficient β (i) and the initial yield stress Y ₀ are corrected, for example, when the original material softening coefficient β (i) is 0 or less, the material softening coefficient β (i) is set to a value larger than 0. For example, while changing to 0.1, the initial yield stress Y ₀ when the material softening coefficient β (i) is 0.1 is calculated using the above equation (1). In that case the calculated initial yield stress Y ₀ is satisfies a condition represented by Formula (3), the latest value beta (i) and the initial yield stress Y ₀ of the material softening coefficient beta is adopted as it is as a final value, respectively .

On the other hand, when the initial yield stress Y ₀ under the changed material softening coefficient β (i) does not satisfy the condition expressed by the above formula (3), the material softening coefficient β (i) is, for example, 0. In the same way, a new initial yield stress Y ₀ is calculated. Such a process is repeated within a range in which the material softening coefficient β does not exceed 1, so that the material softening coefficient β (i) becomes a re-yield stress Y ₁ based on experimental values and a plastic hardening coefficient H (i). Is corrected so as to satisfy all of the above formulas (1) to (3).

  As described above, the correction method when the original material softening coefficient β (i) is 0 or less has been described.

  FIG. 10 is a flowchart showing details of S106 as a material softening coefficient calculation routine.

First, in S251, by using the second measured point 62 the extracted, re-yield stress Y ₁ is determined. Next, in S252, the determined re-yield stress Y ₁ , the calculated plastic hardening coefficient H (i), and the current pre-strain ε ^P _g (i) are expressed by equation (7) in FIG. It is determined whether or not the conditions shown are satisfied.

This time, assuming that the conditions shown by the formula (7) is satisfied, the determination is YES in S252, in S253, the re-yield stress Y ₁ determined above, the plastic hardening coefficient H (i the calculated ) and by substituting the initial yield stress Y ₀ obtained by the experimental value to the equation (1), the material softens coefficient of this beta (i) is calculated. Thereafter, the process proceeds to S254.

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 _g (i) is as described above. Thus, it is corrected so as to satisfy the condition shown in the equation (7).

Thereafter, in S253, the current material softening coefficient β (i) is calculated by using the equation (1) under the corrected pre-strain ε ^P _g (i). Subsequently, the process proceeds to S254.

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 expressed by the equation (2) is satisfied, the determination is YES, and the calculated material softening coefficient β (i) and the initial yield stress Y ₀ based on the experimental value remain as they are. Adopted as final value. This completes one execution of this routine.

On the other hand, if it is assumed that the condition expressed by the expression (2) is not satisfied this time, the determination in S254 is NO, and the current material softening coefficient β (i) is corrected as described above in S256. Is done. Thereafter, in S257, under the material softened coefficients modified beta (i), by using the formula (1), the initial yield stress Y ₀ is calculated.

Subsequently, in S258, whether the initial yield stress Y _0, thus calculated, is greater than 0 is judged. If it is assumed that the value is larger than 0 this time, the determination is YES, and 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.

Also, this time, assuming that the initial yield stress _{Y 0,} thus calculated is less than or equal to 0, next the determination of S258 is NO, the flow 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 under the corrected material softening coefficient β (i). Calculated.

Subsequently, in S258, whether or not the calculated initial yield stress Y0 is greater than ₀ is determined in the same manner as described above. S256 to result execution of S258 is repeated several times, if a YES determination in S258, finally modified material softened coefficient β (i) and the initial yield stress Y ₀ is employed as the final value, respectively The This completes one execution of this routine.

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 γ. As tan γ, Young E (i) corresponding to the current pre-strain ε ^P _g (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.

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 executed in relation to the new pre-strain ε ^P _g (i).

S103 to result execution of S109 is repeated several times, if a YES determination in S108, in S110, the change region of the pre-strain epsilon ^P _g is between the minimum and maximum values set, constant Divided into multiple sections at intervals.

At this time, when the material characteristic value corresponding to the pre-strain epsilon ^P _g corresponding to each segment does not exist, a plurality of pre-strain epsilon ^P _g proximate to prestrain epsilon ^P _g material characteristic value is not present A necessary material property value is generated by interpolating the material property value. As a result, the relationship between the pre-strain ε ^P _g and the material property value is discretely acquired for a plurality of representative values of the pre-strain ε ^P _g as shown in a table form in FIG.

  This completes one execution of the material property value acquisition program.

  Next, the contents of the molding simulation program will be described in detail.

  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. The contents of the molding simulation program will be described below. Prior to that, the premise theory and principle will be described.

  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.

  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. .

  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.

  In the present embodiment, the material behavior accompanying one operation of the punch 110 during press molding is divided into a plurality of calculation steps each having a minute time, and a series of molding is performed by the completion of the plurality of calculation steps. The simulation ends. Furthermore, in this embodiment, a dynamic explicit method is employed as the incremental calculation method.

  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.

  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.

  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 θx, θy, and θz of each component.

  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.

The strain increment is separated into an isotropic strain increment ε _m and a deviation strain increment e _ij . Its Isotropic strain increment epsilon _m is calculated using the strain increment epsilon _ij for each component, the strain increment epsilon _ij is calculated using an approximate function strain described later.

From the isotropic strain increment ε _m , a hydrostatic pressure increment 3Kε _m (K: bulk modulus) is obtained. The hydrostatic pressure component σ _m is updated by adding the hydrostatic pressure component σ _ij ⁰ up to the previous calculation step to the hydrostatic pressure increment 3Kε _m . In other words, the updated hydrostatic pressure component σ _m ^{t + Δt} is calculated, and this is shown in FIG. 15 by equation (11).

By subtracting the isotropic strain increment ε _m from the strain increment ε _ij for each component, the deviation strain increment e _ij for each component is obtained.

A deviation stress component 2Ge _ij (G: shear elastic modulus) is obtained from the deviation strain increment e _ij . Therefore, the deviation stress component 2Ge _ij is obtained in the elastic region in the plane stress state. The trial stress (trial value of the deviation stress S) S _ij ^Trial is obtained by adding the deviation stress component S _ij ⁰ up to the previous calculation step to the deviation stress component 2Ge _ij . That is, the updated trial stress S _ij ^Trial is calculated, and this is shown by Equation (12).

By substituting this trial stress S _ij ^Trial into equation (13), a trial solution σ ^{Trial of} equivalent stress is obtained.

When the trial solution σ ^Trial of this equivalent stress is in the elastic range and the material of each shell element has not yet yielded, the above-described trial stress S _ij ^Trial is left as it is, and the updated stress component (the current stress component) ) S _ij ^Final .

On the other hand, as shown in FIG. 16A representing the stress space, when the trial solution σ ^Trial of the equivalent stress reaches the plastic region, it is determined that the material of each shell element has yielded. In this case, as shown in FIG. 15B, the trial stress S _ij ^Trial is pulled back onto the yield surface by the radial return algorithm (radius pullback method) shown in FIG.

In the radial return algorithm, the equivalent stress σ representing the radius of the yield surface and the provisional value of the final solution S _ij ^Final are calculated from the equivalent plastic strain increment ε ^P. Furthermore, compared them both calculated values, the provisional value of the final solution _S ^{ij Final} when matching equivalent stress σ is determined as the final value of the final solution _S ^{ij Final.}

In addition, when the above-described von Mises yield condition is adopted, for example, the increment Δε ^P is calculated by using the equation (16) in FIG. 15, and then the equation (14) is used. It is possible to calculate the final solution S _ij ^Final .

When the material of each shell element yields, the final value of the pulled back final solution S _ij ^Final is obtained as an updated stress component.

Then, as shown by the equation (15) in FIG. 15, the stress σ _ij ^{t + Δt} updated for each component by the updated hydrostatic pressure component σ _m ^{t + Δt} and the updated deviation stress component S _ij ^Final. Is calculated.

  One calculation step of the stress update described above is sequentially executed for all shell elements constituting the plate material.

  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.

  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.

  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.

  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.

  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.

  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. The bending severity Lb is a bending index R / t that is a value obtained by dividing the bending radius R [mm] by the plate thickness t [mm], and the smaller the value, the stronger the bending. It is expressed in

  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.

  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.

  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.

  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). The component displacements u, v, w) and the rotation angle θ (rotation angles θx, θ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. .

  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.

  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.

  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.

  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.

  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.

  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.

  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.

  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

  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.

  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). .

  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.

  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.

  As shown in FIG. 27, in the conventional shell element, as described above, 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.

  The vertical stress σz in the plate thickness direction is a contribution σzb (hereinafter referred to as “vertical stress σzb”) due to bending deformation of the thick plate, and a contribution σzc (hereinafter, “ It can be considered to be the sum of the vertical stress σzc ”. The normal stress σzc is caused by shear deformation generated in the thick plate by the contact pressure p.

  The vertical stress σzb 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”.

  In the present embodiment, first, the vertical stress σx is tentatively analyzed ignoring the vertical stress σz, and based on the provisional value and using the above equation (41), the vertical stress σzb 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.

  The specification of the stress approximation function for the normal stress σzb can be specifically performed as follows.

  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.

Equation (42) in FIG. 28 is an example of the stress approximation function. 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. . When the z-direction position of the neutral plane is represented by r ₀ and the plate thickness is represented by t ₀ , “a” is represented by Expression (43), and “b” is represented by Expression (44).

  Therefore, the normal stress σzb represented by the above formula (42) is 0 on both surfaces of the thick plate, which indicates that the contact pressure 110 received by the punch 110 from the punch 110 is not considered. Yes.

  If the vertical stress σzb is expressed by a stress approximation function f (z) expressed by the equation (42), the equation (41) is transformed into the equation (45).

  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.

  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.

  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.

  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 ′.

  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.

  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.

  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.

  If the vertical stress σzb 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.

  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.

  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.

  FIG. 33 shows that the vertical stress σz in the thickness direction corresponds to the sum of the vertical stress σzb due to bending deformation of the thick plate and the vertical stress σzc due to the contact pressure p using the weighting function in the equation (69). It is represented. In the same figure, a stress approximation function that approximately represents the distribution of the normal stress σzc in the plate thickness direction is expressed by equation (70). When this stress approximation function is employed, the vertical stress σzc is equal to the contact pressure p at the bending inner surface (z = a) of the thick plate, while the bending outer surface (z = b), as represented by the equation (71). ), The normal stress σzc becomes zero.

In FIG. 33, an example of the weighting function represented by Expression (69) is represented by Expression (72). In this example, the weighting factor for vertical stress σzb (corresponding to “k ₁ ” in equation (69)) and the weighting factor for vertical stress σzc (corresponding to “k ₂ ” in equation (69)). Are set as fixed values. Furthermore, in this example, the weighting factors and are set equal to each other.

In FIG. 33, another example of the weighting function is represented by Expression (73). “Α” in the equation (73) is represented by the equation (74). In this example, the weighting factor for vertical stress σzb (corresponding to “k ₁ ” in equation (69)) and the weighting factor for vertical stress σzc (corresponding to “k ₂ ” in equation (69)). Are set as variable values defined using a quadratic function.

In FIG. 33, still another example of the weighting function is represented by Expression (75). “X” in the equation (75) is represented by the equation (76). In this example, the weighting factor for vertical stress σzb (corresponding to “k ₁ ” in equation (69)) and the weighting factor for vertical stress σzc (corresponding to “k ₂ ” in equation (69)). Are set as variable values defined using a saturated function (for example, Brillouin function).

  In the present embodiment, the vertical stress σx (= the bending direction stress that is the first in-plane stress) and the vertical stress σy (= the cross-sectional stress that is the second in-plane stress) acting on the plate material constitute the shell element. The analysis is performed in consideration of the normal stress σz (= stress in the plate thickness direction) despite being used. Here, the vertical stress σz is expressed as F (σzb, σzc)) by the equation (69) in FIG.

  Specifically, in the present embodiment, first, a plane stress state in which the plate thickness direction stress σz is 0 is assumed, and the bending direction stress σx and the cross-sectional direction stress σy are obtained. Next, the distribution of the thickness direction stress σz in the thickness direction is defined from the obtained bending direction stress σx and the sectional direction stress σy, or the bending direction stress σx and the sectional direction stress σy obtained in the previous calculation step. In order to do so, the stress approximation function is identified. Thereafter, by using the identified stress distribution function, the plate thickness direction stress σz is obtained, and on the assumption of the obtained plate thickness direction stress σz, that is, considering the plate thickness direction stress σz, bending is performed. The directional stress σx ′ and the cross-sectional direction stress σy ′ are obtained using the aforementioned radial return algorithm.

  According to the Prandtl-Reuss flow law, the equation (81) in FIG. 34 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 an x-direction normal stress (bending direction stress) and a y-direction normal stress (cross-sectional direction stress) without considering the normal stress σz, respectively, “Σx ′” and “σy ′” in 83) mean an x-direction normal stress and a y-direction normal stress when the normal stress σz is considered.

  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.

  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.

  In FIG. 35, during the press forming, each shell element constituting the plate material is actually superimposed with the bending of each shell element and the elongation due to the axial force F parallel to the neutral surface of each shell element. It is conceptually expressed to work. FIG. 35 also conceptually shows that such material mechanical behavior is theoretically separated into the bending and elongation of each shell element. Here, since the axial force F exists for each vertical stress σx, σy and is the sum of the vertical stresses σx, σy on the cross-section of the shell element, it is also referred to as a cross-sectional force.

In the present embodiment, the vertical stresses σx ′ and σy ′ are corrected for each position in the thickness direction so that the cross-sectional force F considering the vertical stress σz is maintained at the cross-sectional force in the plane stress state. The In other words, in the present embodiment, the normal stresses σx ′ and σy ′ are corrected without changing the neutral surface position r ₀ of the plate material, so that the final values of the vertical stresses σx and σy become the normal stress σx ′. Induced as', σy ''.

  As described above, the in-plane stresses σx and σy at each position in the plate thickness direction are calculated without considering the values calculated in consideration of the plate thickness direction stress σz in terms of plastic mechanics. From this value, the plate thickness direction stress σz acting on each same position increases by the same amount. On the other hand, the balance of forces in the bending direction and the cross-sectional direction of the plate material is established for the sum of in-plane stresses acting on the entire cross section, that is, the cross-sectional force F, between the cross-sections adjacent to each other in the plate material.

  Therefore, in order to maintain the balance of force with respect to the sectional force F before and after considering the plate thickness direction stress σz, each of the in-plane stresses σx ′ and σy ′ considering the plate thickness direction stress σz What is necessary is just to subtract the cross-sectional force increment Fz which is the sum total in the same cross section of sheet thickness direction stress (sigma) z from the cross-sectional force F which is the sum total in a cross section.

  The sectional force increment Fz is represented by the formula (91) in FIG. In this equation (91), “σz” means the stress approximation function described above.

  In the present embodiment, the plate material is approximately expressed by a plurality of shell elements arranged in the plane direction thereof, and each shell element is a plane connecting only the center of the plate thickness. However, the contact treatment is performed on the inner and outer surfaces in consideration of the plate thickness. As shown in FIG. 37, the axial forces and moments of the elements are calculated at in-plane integration points Ip1 to Ip4 (the number of integration points varies depending on the elements used, but is typically four). The node is determined by using the shape function. Each integration point in the plane has an integration point (indicated by “x” in FIG. 37) in several layers in the thickness direction. A plurality of discrete representative points representing continuous regions inside each shell element are set as integration points, and these set integration points can be used as evaluation points to be noted inside each shell element. It is.

  In this case, by calculating the normal stress σz for each of a plurality of integration points arranged in the thickness direction of each shell element among the integration points, and summing the calculated values of the normal stress σz for these integration points. , The above-described cross-sectional force increment Fz is calculated.

  Correcting the in-plane stresses σx ′ and σy ′ using the cross-sectional force increment Fz is effective in maintaining the balance of the force in the bending direction and the cross-sectional direction of the plate material before and after considering the plate thickness direction stress σz. is there. However, the sectional force increment Fz is an increment of the in-plane stress σx, σy for the entire section, and is not an increment of the in-plane stress σx, σy at each position in the plate thickness direction. However, it is impossible to accurately correct the in-mask stresses σx ′ and σy ′ for each position in the thickness direction.

  However, the distribution of the plate thickness direction stress σz in the plate thickness direction is based on the condition that the balance of the forces in the bending direction and the cross section direction of the plate material is maintained at least with respect to the cross section force F, and It is possible to predict so that the condition of showing a distribution as close as possible to the distribution shown in the thickness direction is satisfied.

  Furthermore, according to the predicted distribution of the plate thickness direction stress σz, a distribution amount σzm for allocating the sectional force increment Fz to each of a plurality of positions in the plate thickness direction is determined, and the determined distribution amount σzm is used as the in-mask stress σx. By subtracting from “, σy”, the in-mask stresses σx ′, σy ′ are corrected with high accuracy.

  By the way, for example, the bending index R [mm] of the plate material is divided by the plate thickness t [mm] of the plate material, and the bending index obtained by dividing the bending radius R [mm] of the plate material is smaller than 3. Assuming press forming to be performed, the plate material reaches the plastic region in almost the entire region in the plate thickness direction by the press forming. On the other hand, when the gradients of the graphs representing the relationship between the strain and the plate thickness direction stress σz are compared with each other for the plastic region and the elastic region, the plastic region is more gradual than the elastic region. This means that the dependence of the plate thickness direction stress σz on the plate thickness direction position is smaller in the plastic region than in the elastic region.

  Therefore, when it is necessary to analyze the in-plane stresses σx and σy of a plate material obtained by press molding performed under the condition that the bending severity Lb of the plate material is strong as described above, It is appropriate to approximate the distribution of the thickness direction stress σz in the thickness direction with a uniform distribution. Therefore, it is appropriate to determine the distribution amount σzm of the aforementioned sectional force increment Fz to be uniform in the thickness direction.

Based on such knowledge, in this embodiment, as represented by the formula (92) in FIG. 36, the calculated cross-sectional force increment Fz is divided by the plate thickness t ₀ , thereby obtaining an average in the plate thickness direction. The average value is determined as the distribution amount σzm of the cross-sectional force increment Fz. Specifically, in the example shown in FIG. 37, the cross-sectional force increment Fz is divided by the number of integration points arranged in the thickness direction, whereby the cross-sectional force increment Fz is applied to each of the integration points. They will be allocated equally to each other.

  Further, as represented by the equation (93) in FIG. 36, the calculated distribution amount σz is subtracted from the bending direction stress (in-plane stress) σx ′ taking into account the plate thickness direction stress σz, whereby the bending direction A final value σx ″ of the stress σx is derived. Similarly, as represented by the equation (94), the calculated distribution amount σz is subtracted from the cross-sectional direction stress (in-plane stress) σy ′ considering the plate thickness direction stress σz, thereby obtaining the cross-sectional direction stress. A final value σy ″ of σy is derived.

  In FIG. 38, the bending direction stress σx is taken as an example, and a plurality of types of calculated values for the in-plane stress are represented by a plurality of graphs, respectively. These calculated values are

(1) A bending direction stress σx calculated using a shell element model and assuming that the thickness direction stress σz is 0;
(2) Bending direction stress σx ′ calculated using the shell element model and considering the thickness direction stress σz;
(3) The final value σx ″ (= σx′−σzm) of the bending direction stress σx calculated by subtracting the distribution amount σzm from the bending direction stress σx ′;

  Is included. In FIG. 38, the value of the plate thickness direction stress σz is also represented by a graph.

  38, the final value σx ″ is a value obtained by subtracting the distribution amount σzm from the bending direction stress σx ′, but the graph representing the final value σx ″ is a bending direction stress σx ′. In this example, the distribution amount σzm is a compressive stress, and its sign is negative, which corresponds to a graph obtained by translating the graph representing σ in the positive direction of the stress by an amount equal to the absolute value of the distribution amount σzm. Because there is.

  Therefore, according to the present embodiment, the sectional force, which is the sum of the vertical stress σx, does not change depending on whether or not the vertical stress σz is considered. 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.

  The contents of the molding simulation program shown in FIG. 12 have been schematically described above, but will be specifically described below.

  In this molding simulation program, first, in S401, data necessary for molding analysis is read from the MO 32 and the analysis is initialized.

  Next, in S402, boundary conditions regarding the load acting on the plate material are applied. In the present embodiment, the boundary condition regarding the load acting on the plate material from the pad 114 is applied.

  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.

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, a strain increment ε _ij is calculated for each component based on the strain approximation function. Here, the strain increment ε _ij includes both normal strain and shear strain.

Thereafter, in S503, based on the calculated strain increment epsilon _ij, isotropic strain increment epsilon _m is calculated. Subsequently, in S504, based on the incremental epsilon _m strain isotropic that is the calculation, hydrostatic pressure increment 3Keiipushiron _m is calculated. Thereafter, in S505, the updated hydrostatic pressure σ _m ^{t + ΔT} is calculated as described above based on the stress σ _ij ⁰ for each component and the hydrostatic pressure increment 3Kε _m in the previous calculation step.

Subsequently, in S506, the deviation strain increment e _ij is calculated for each component based on the calculated strain increment ε _ij and isotropic strain increment ε _m for each component. Subsequently, in S507, in plane stress, and, in the elastic region, deviatoric stress increment _{2GE ij} is calculated.

Thereafter, in S508, the updated trial stress S _ij ^{t + ΔT} is calculated as described above based on the deviation stress S _ij ⁰ for each component in the previous calculation step and the deviation stress increment 2Ge _ij . Subsequently, in S509, the trial solution σ ^{Trial of} equivalent stress is calculated as described above based on the calculated value of the trial stress S _ij ^{t + ΔT} .

  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.

In S511, by using the equation (16) in FIG. 15, increment [Delta] [epsilon] ^P of equivalent plastic strain increment epsilon ^P is calculated. Thereafter, in S512, the equivalent stress σ corresponding to the calculated increment Δε ^P is set as the corresponding equivalent stress σ, and the final solution S _ij ^Final corresponding to the corresponding equivalent stress σ and the trial stress S _ij ^Trial is Calculation is made using equation (14) in FIG. In the present embodiment, the final solution S _ij ^Final is an example of the “provisional value” in the item (1).

  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.

For example, when each round of execution of S511, equivalent plastic strain increment epsilon ^P is increased at a constant small increment [Delta] [epsilon] ^P. In contrast, during each round of execution of S512, the equivalent stress σ is calculated as the corresponding equivalent stress σ corresponding to the provisional value of the equivalent plastic strain increment epsilon ^P (calculated by the execution of S510). Further, a provisional value of the final solution S _ij ^Final corresponding to the corresponding equivalent stress σ and the trial stress S _ij ^Trial is calculated.

In S512, it is further determined whether or not the corresponding equivalent stress σ and the provisional value of the final solution S _ij ^Final sufficiently match each other. If they match well in the provisional value of the equivalent plastic strain increment epsilon ^P provisional value and the final solution S _ij ^Final is the final value, respectively. The final value of the final solution S _ij ^Final means an updated deviation stress component. On the other hand, if the corresponding equivalent stress σ and the provisional value of the final solution S _ij ^Final are not sufficiently coincident with each other, the process returns to S511, and the provisional value of the equivalent plastic strain increment ε ^P increases with a constant minute increment Δε ^P. Be made.

In any case, if the final value of the equivalent plastic strain increment ε ^P and the updated deviation stress component S _ij ^Final are obtained by the radial return algorithm, then, in S513, the updated deviation stress component is updated. Based on S _ij ^Final and the updated hydrostatic pressure component σ _m ^{t + Δt} , the updated stress σ _ij ^{t + Δt} is calculated for each component by using the equation (15) in FIG.

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, based on the trial stress S _ij ^Trial updated up to the previous calculation step and the updated hydrostatic pressure component σ _m ^{t + Δt} , by using the equation (17) in FIG. The updated stress σ _ij ^{t + Δt} is calculated.

If the plate yields on at least one of its both surfaces, after execution of S513, in S515, the final solution S _ij ^Final (corresponding to the in-plane stress σx in the plane stress state) or the stress in the previous calculation step. For example, by using the formula (41) to the formula (55) shown in FIG. 28, the formula (70) shown in FIG. 33, and the formula (73) or (75) shown in FIG. 33, for example. As described above, the stress approximation function that approximately defines the thickness direction stress σz (= F (σzb, σzc)) is specified. Specifically, the specified stress approximation function includes a stress approximation function that defines a distribution in the plate thickness direction of the plate thickness direction stress σzb accompanying bending deformation, and a plate thickness direction of the plate thickness direction stress σzc due to the contact pressure p. Stress approximation function that defines the distribution at.

  Thereafter, in S516, by using the stress approximation function specified as described above, the bending direction stress σx ′ and the cross-sectional direction stress σy ′ satisfying the current thickness direction stress σz are obtained based on the radial return algorithm. Calculated. The bending direction stress σx ′ and the cross-sectional direction stress σy ′ are calculated for each position in the plate thickness direction, that is, for each integration point described above. In the present embodiment, the calculated bending direction stress σx ′ and cross-sectional direction stress σy ′ are examples of “temporary in-plane stress” in the item (1).

  Subsequently, in S517, the cross-section associated with considering the thickness direction stress σz as described above by using the above-described stress approximation function and the calculation formula represented by Expression (91) in FIG. A force increment Fz is calculated. This cross-sectional force increment Fz is calculated for each shell element. If the shell element has only one in-plane integration point, this cross-sectional force increment Fz is calculated one for each shell element.

  Thereafter, in S518, the distribution amount σzm is calculated as described above by using the calculated sectional force increment Fz and the calculation formula represented by the formula (92) in FIG. This distribution amount σzm is also calculated for each shell element.

  Subsequently, in S519, by using the calculated distribution amount σzm and the calculation formula represented by the formula (93) in FIG. 36, for each position in the plate thickness direction (for example, each For each integration point), the bending direction stress σx ′ is corrected to obtain the final value σx ″ of the bending direction stress. In S519, by using the calculated distribution amount σzm and the calculation formula represented by the formula (94) in FIG. 36, for each position in the thickness direction as described above (for example, For each integration point), the cross-sectional stress σy ′ is corrected to obtain the final cross-sectional stress σy ″.

  In other words, in the present embodiment, the acquired final values σx ″ and σy ″ constitute an example of the “final value” in the item (1). In addition, the final values σx ″ and σy ″ 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.

Thereafter, in S520, the updated plastic strain ε ^P _ij ^{t + Δt} is calculated by using 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. Is done. “Δλ” in the equation (19) is obtained from the equations (18) and (19). The plastic strain ε ^P _ij ^{t + Δt} when it is determined in S409 in FIG. 12 that a series of analysis has been completed means the residual strain. The residual strain is stored in the memory 22 in association with each shell element.

Subsequently, in S521, based on the updated plastic strain ε ^P _ij ^{t + Δt} , the plate thickness t ₀ of the current attention shell element 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.

  In any case, in S522, it is determined whether or not the execution of S501 to S521 has been performed for all the shell elements of the plate material. If it is assumed that the execution of S501 to S521 is not performed for all shell elements, the determination is NO, and in S523, the next shell element is set as a new focused shell element, and then the process returns to S501. On the other hand, if the execution of S501 to S521 is performed for all shell elements, the determination in S522 is YES, and in S524, the finite element method is executed for the shell element model that approximates the plate material.

  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.

  Thereafter, in S525, 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 S524 are executed using the new overall stiffness matrix. On the other hand, when the solution of the stiffness equation has converged, the determination in S525 is YES, and one execution of this element stress calculation routine is completed.

  Thereafter, in S404 of FIG. 12, the contact force acting on the plate material from the punch 110 is calculated by considering the deformation of the plate material and the shape of the punch 110 under the boundary condition.

  Subsequently, in S405, the acceleration of each node of each shell element is updated by considering the boundary condition regarding the machining speed of the punch 110. Thereafter, in S406, the acceleration of each node of each shell element is updated by considering the analysis condition regarding the contact between the punch 110 and the plate material. Subsequently, in S407, the speed of each node is updated by performing time integration on these node accelerations, and further, the displacement of each node is updated by performing time integration on the calculated speed. Is done.

  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 a phenomenon, the punch 110 enters a small distance into the die 112. Corresponding to the action to be performed.

  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.

  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.

  Next, the contents of the material property value determination program will be described in detail.

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 ε ^P ₀ (j) corresponding to the current shell element j is read from the memory 22.

Thereafter, in S203, it is determined which of the plurality of pre-strain categories the initial strain ε ^P ₀ (j) of this time belongs to, and further, the material characteristic value corresponding to the pre-strain category determined to belong to And the data that associates the current shell element number j with each other are 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.

Subsequently, in S204, necessary correction is performed on the current initial strain ε ^P ₀ (j). Specifically, in consideration of the fact that the pre-strain ε ^P _g strictly corresponding to the material property value selected in S203 does not necessarily exactly match the initial strain ε ^P ₀ (j) of this time, the agreement If not, the current initial strain ε ^P ₀ (j) is corrected to match. Further, data that associates the number j of the current shell element with the corrected initial strain ε ^P ₀ (j) 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 corrected for each shell element.

  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.

  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.

  Next, the contents of the back stress calculation program will be described in detail.

  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.

  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 α.

  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.

  Here, a method of calculating the back stress α will be described.

  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.

  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.

  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).

  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.

  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).

  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 σ.

  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.

  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.

  In order to verify the effect of the springback analysis method according to the present embodiment, the inventors experimentally performed a hat bending process on a certain plate material as the workpiece W using a press molding machine shown in FIG. Through this experiment, experimental values described later were obtained.

  As shown in FIG. 44, the press molding machine corresponds to the press molding machine 100 shown in FIG. 3 with a blank holder 126 added. The blank holder 126 is a position facing the upper surface of the die 112 in order to prevent wrinkles from occurring during press molding at the edge of the plate material that is supported by the upper surface of the die 112. Is arranged. The blank holder 126 holds the edge portion of the plate material with the die 112 during press molding.

  44 shows the punch 110 in addition to the die 112, the workpiece W and the blank holder 126. In FIG. 44, the punch 110 is shown by a solid line in a state at the initial position, and is shown by a broken line in a state at a full stroke position. In FIG. 44, the die shoulder R that is the round radius of the shoulder of the die 112 is indicated by “Rd”, and the punch shoulder R that is the round radius of the shoulder of the punch 110 is indicated by “Rp”. ing.

  Using the press molding machine shown in FIG. 44, the hat bending process was performed under the following conditions.

(1) Conditions related to the plate material • Plate thickness t ₀ : 2.3 [mm]
・ Bending radius R: 3 [mm]

(2) Conditions related to press molding machines • Die shoulder Rd: 3 [mm]
-Punch shoulder Rp: 5 [mm]
・ Shoulder clearance: 0.23 [mm]

  The inventors selected the opening amount and the curvature of the vertical wall portion of the plate material as factors for evaluating the springback generated in the plate material after the hat bending process.

  As shown in FIG. 45, the opening amount B is defined as the amount of warping of the vertical wall portion of the plate material in the width direction (X direction) from the ideal shape (the inner shape of the die 112) in which no springback occurs. Has been.

  Specifically, in order to calculate the opening amount B, the point of interest PX, which is 15 [mm] apart from the bottom surface in the axial direction (Z direction) of the plate material on the bending outer surface of the plate material, is similarly 80 [ mm] separated points of interest PC are set. Further, an attention point PY that is 65 [mm] apart from the attention point PX is set along a straight line (indicated by a broken line in FIG. 45) connecting the attention point PX and PC.

  A distance d at which the attention point PY is separated from the press center line (indicated by a one-dot chain line in FIG. 45) in the width direction (X direction) is measured. The opening amount B is calculated by subtracting the width direction distance A (for example, 52.53 [mm]) of the inner surface of the die 112 from the press center line from the measured distance d.

  On the other hand, a curvature means the curvature of the curve formed by the vertical wall part of the board | plate material, when the bending outer surface of a board | plate material is cut | disconnected by the plane containing the centerline of the board | plate material, as shown in FIG. Is defined as

  Specifically, in order to calculate this curvature, 20 [mm], 50 [mm] and 50 [mm] are spaced apart from the bottom surface in the axial direction (Z direction) on the outer surface of the vertical wall portion of the plate material. Only three attention points PA, PB and PC points are set. A single perfect circle passing through these points of interest PA, PB and PC is assumed, and the curvature is calculated as the reciprocal of the radius ρ of the perfect circle.

  The inventors of the present invention analyzed the springback, which is expected to occur in the plate material when the hat bending process is performed on the plate material equivalent to the above-described experimental product under the same conditions, by three types of molding simulations. In FIG. 46, the three types of analysis results obtained by these three types of molding simulations are represented by a graph from the viewpoint of the amount of opening together with the experimental values. Further, in FIG. 47, three types of analysis results by these three types of molding simulations are represented by a graph from the viewpoint of curvature (1 / curvature radius ρ) together with experimental values.

  These three types of forming simulations are simulations according to the present embodiment, and the force balance in the bending direction and the cross-sectional direction due to considering the plate thickness direction stress σz without moving the neutral plane of the plate material. Includes things to rectify. 46 and 47, the result of this type of molding simulation is shown with the label “Axial force correction”.

  These three types of forming simulations further include correcting the unbalance of forces in the bending direction and the cross-sectional direction associated with considering the plate thickness direction stress σz by moving the neutral surface of the plate material. In FIGS. 46 and 47, the result of this type of molding simulation is shown with a label “middle plane movement”.

  These three types of forming simulations further include those that do not correct the force imbalance in the bending direction and the cross-sectional direction that accompanies the plate thickness direction stress σz. 46 and 47, the result of this type of molding simulation is shown with a label “no axial force correction”.

  As shown in FIG. 46, according to the molding simulation according to the present embodiment, the experimental value of the opening amount is reproduced with sufficiently high accuracy, and the reproduction accuracy is similar to that of the molding simulation with neutral plane movement. is there.

  As shown in FIG. 47, according to the molding simulation according to the present embodiment, the experimental value of the curvature is reproduced with sufficiently high accuracy, and the reproduction accuracy is higher than that of the molding simulation accompanied with the neutral plane movement.

  As is clear from the above description, in the present embodiment, S512 in FIG. 14 constitutes an example of the “provisional value calculation step” in the item (1), and S515 is a “stress approximation function specifying step” in the same term. S516 constitutes an example of the “in-mask stress calculating step” in the same term, S517 constitutes an example of the “cross-sectional force increment calculating step” in the same term, and S518 represents the “distribution amount” in the same term. An example of “decision process” is configured, and S519 constitutes an example of “final value determination process” in the same section.

  Further, in the present embodiment, S519 in FIG. 14 constitutes an example of the “final value calculation step” in the item (2), and S518 represents an example of the “distribution amount determination step” in the item (3) and the ( An example of the “distributed amount determining step” in the item 4) is configured, and S517 configures an example of the “sectional force increment calculating step” in the item (5).

  Further, in the present embodiment, the molding simulation program shown in FIG. 12 constitutes an example of the “program” according to the above item (8), and the memory 22 storing the molding simulation program, the CD-ROM 30 or Each of the MOs 32 constitutes an example of a “recording medium” according to the item (9).

  As described above, one of the embodiments of the present invention has been described in detail with reference to the drawings. It is possible to implement the present invention in other forms in which various modifications and improvements are made based on the above.

It is a systematic diagram which shows the springback analysis apparatus suitable for implementing the springback analysis method including the stress analysis method according to one Embodiment of this invention. It is process drawing which shows the said springback analysis method. It is front sectional drawing which shows an example of the press molding machine which can bend | fold a board | plate material in U shape. It is a graph for demonstrating a mode that the experimental value of stress-strain relation is acquired for every pre-strain in the said springback analysis method. It is a graph for demonstrating the material characteristic value required for the material hardening model used for the springback analysis of step S4 in FIG. 2, and its springback analysis. It is another graph for demonstrating the material characteristic value required for the said springback analysis. It is a graph for demonstrating the material hardening model used for the said springback analysis, contrasting with an isotropic hardening model and an actual stress-strain relationship. FIG. 3 is a flowchart conceptually showing the contents of a material property value acquisition program executed by a computer for performing material property value acquisition in step S <b> 2 in FIG. 2. It is a figure which shows some formulas for demonstrating the theory and principle with which a material characteristic value is acquired by the material characteristic value acquisition program of FIG. 9 is a flowchart conceptually showing details of S106 in FIG. 8 as a material softening coefficient calculation routine. It is a figure which shows a mode that the relationship between pre-strain and a material characteristic value is discretely acquired by the material characteristic value acquisition program of FIG. FIG. 3 is a flowchart conceptually showing the contents of a molding simulation program executed by a computer to perform a molding simulation in step S4 in FIG. 2. 13 is a part of a flowchart conceptually showing details of S403 in FIG. 12 as an element stress calculation routine. It is another part of the flowchart showing the said element stress calculation routine. It is a figure which shows some formulas for demonstrating the theory and principle utilized for the element stress calculation routine shown to FIG. 13 and FIG. It is a graph for demonstrating the theory and principle utilized for the element stress calculation routine shown in FIG. 13 and FIG. It is a figure which shows the displacement-strain relational expression utilized for the shell element of a prior art example. It is a figure which shows the strain distribution represented by the displacement-strain relational expression of FIG. 17 with the bending deformation of a board | plate material. It is a graph which shows a mode that the relationship between cross-section rotation angle (theta) y and the distance ds from the bending inner surface of a board | plate material changes according to the bending | striction Lb. 3 is a flowchart conceptually showing the contents of a strain approximation function specifying routine incorporated in a forming simulation program for executing a forming simulation in step S4 in FIG. It 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, against the actual measurement value. It is a graph which shows some examples of the execution result of the distortion approximation function specific routine of FIG. Three types of calculated values regarding the relationship between the shear strain γxz and the distance ds from the bent inner surface of the plate, calculated using the solid element model, calculated using the shell element model of the above embodiment, and the conventional example It is a graph which shows what consists of a calculated value using a shell element model. It is a figure which shows the distribution of the shear strain (gamma) xz calculated using the shell element of the said embodiment, and the distribution of the shear strain (gamma) xz calculated using the shell element of a prior art example with the bending deformation of a board | plate material. It is a figure which shows the displacement-strain relational expression utilized for the shell element model of the said embodiment. It is a figure which shows the distribution in the plate | board thickness direction of the normal strain (epsilon) x and shear strain (gamma) xz calculated using the shell element of the said embodiment with the bending deformation of a board | plate material. It is a figure which shows distribution of normal stress (sigma) z calculated using the shell element of the said embodiment, and distribution of normal stress (sigma) z calculated using the shell element of a prior art example with the bending deformation of a board | plate material. It is a figure which shows a part of some formula for demonstrating the theory and principle utilized in order to specify a stress approximation function in S516 of FIG. It is a figure which shows another part of some formulas for demonstrating the theory and principle utilized in order to specify the said stress approximation function. It is a figure which shows another part of some formulas for demonstrating the theory and principle utilized in order to specify the said stress approximation function. It is a figure which shows another part of some formulas for demonstrating the theory and principle utilized in order to specify the said stress approximation function. It is a figure which shows another part of some formulas for demonstrating the theory and principle utilized in order to specify the said stress approximation function. It is a figure which shows another part of some formulas for demonstrating the theory and principle utilized in order to specify the said stress approximation function. It is a figure which shows some formulas for demonstrating the relationship between normal stress (sigma) x, (sigma) z which does not consider normal stress (sigma) z, and normal stress (sigma) x ', (sigma) y' which considered normal stress (sigma) z. In S517 to S519 in FIG. 14, the axial force correction is performed to correct the unbalance of the forces in the bending direction and the cross-sectional direction due to considering the plate thickness direction stress σz without moving the neutral plane of the plate material. It is a figure for demonstrating notionally the principle. It is a figure which shows some formulas for demonstrating the principle of the said axial force correction | amendment. It is a perspective view which shows a shell element in order to demonstrate the said axial force correction | amendment. In order to explain the axial force correction, an approximate value of the plate thickness direction stress σz is shown, and the bending direction stress σx is calculated when the plate thickness direction stress σx is not taken into account, and when calculated, It is a graph which shows the calculated value at the time of performing the said axial force correction | amendment. FIG. 3 is a flowchart conceptually showing the contents of a material property value determination program executed by a computer for performing material property value determination in step S5 in FIG. 2. FIG. 3 is a flowchart conceptually showing a back stress calculation program executed by a computer for performing the back stress calculation in step S <b> 6 in FIG. 2. It is a figure which shows a part of some formula for demonstrating the theory and principle utilized in order to calculate a back stress by the back stress calculation program of FIG. FIG. 5 is a diagram showing another part of some equations for explaining the theory and principle used to calculate the back stress. It is a graph for demonstrating the execution content of the back stress calculation program of FIG. It is side surface sectional drawing which shows the principal part of the press molding machine used in order to implement a hat bending process in order to demonstrate the effect of the said embodiment. It is side surface sectional drawing for demonstrating the opening amount and curvature of the vertical wall part of the board | plate material which are the factors for evaluating the springback which generate | occur | produces in a board | plate material by the said hat bending process. It is a graph for demonstrating the analysis precision of the springback by the said embodiment with three comparative examples regarding the opening amount in FIG. It is a graph for demonstrating the analysis precision of the springback by the said embodiment with three comparative examples regarding the curvature in FIG.

Explanation of symbols

12 Computer 16 External storage device 30 CD-ROM
32 MO
100 Press molding machine 110 Punch 112 Die 114 Pad 126 Blank holder

Claims (9)

In order to analyze the bending deformation of the plate material due to press molding by a computer, by using a shell element model modeled on the computer as an aggregate of a plurality of shell elements in which the plate material is arranged in the plane direction thereof, A stress analysis method for analyzing in-plane stress, which is stress in the surface direction of each shell element, in consideration of plate thickness direction stress, which is stress in the thickness direction of each shell element,
For each shell element, a provisional value calculation step of calculating a provisional value of the in-plane stress assuming a plane stress state in which the plate thickness direction stress is zero;
Based on the calculated provisional value of the in-plane stress, a stress approximation function specifying step for specifying a stress approximation function that approximately defines a distribution of the plate thickness direction stress in the plate thickness direction;
By using the identified stress approximation function, the in-plane stress that is mechanically established together with the plate thickness direction stress at each position in the plate thickness direction is defined as the in-plane stress in consideration of the plate thickness direction stress. A step of calculating the in-mask stress to be calculated;
Based on the distribution of the plate thickness direction stress in the plate thickness direction, the sum of the in-plane stress on the cross section of each shell element due to the consideration of the plate thickness direction stress when calculating the in-plane stress A cross-sectional force increment calculation step of calculating the amount of increase in cross-sectional force as a cross-sectional force increment;
A distribution amount determining step for determining a distribution amount for distributing the calculated cross-sectional force increment to each of a plurality of positions arranged in the plate thickness direction;
A final value calculating step of calculating a final value of the in-plane stress by correcting the calculated in-plane stress with the determined distribution amount for each position in the plate thickness direction. Method.   The final value calculating step corrects the calculated in-plane stress by subtracting the determined distribution amount from the calculated in-surface stress for each position in the plate thickness direction. 2. The stress analysis method according to 1.   3. The stress analysis method according to claim 1, wherein in the distribution amount determination step, the distribution amount is determined such that the calculated cross-sectional force increments are distributed substantially evenly to the plurality of positions, respectively.   The stress analysis method according to any one of claims 1 to 3, wherein the distribution amount determination step determines the distribution amount so as not to depend on a position in the plate thickness direction.   The stress analysis method according to claim 1, wherein the cross-sectional force increment calculating step calculates the cross-sectional force increment by substantially integrating the stress distribution function with respect to a position in the plate thickness direction.   6. The press molding is performed according to any one of claims 1 to 5, wherein a bending index obtained by dividing a bending radius R [mm] of the plate material by a plate thickness t [mm] of the plate material is smaller than 3. The stress analysis method described.   The stress analysis method according to claim 1, wherein the press forming includes at least one of thick plate forming, hemming, and hat bending.   A program executed by a computer to perform the method according to claim 1.   The recording medium which recorded the program of Claim 8 so that computer reading was possible.

Images