CN116187133B - Dimension separation spring comparison method for spinning mobile grid method - Google Patents
Dimension separation spring comparison method for spinning mobile grid method Download PDFInfo
- Publication number
- CN116187133B CN116187133B CN202310098738.5A CN202310098738A CN116187133B CN 116187133 B CN116187133 B CN 116187133B CN 202310098738 A CN202310098738 A CN 202310098738A CN 116187133 B CN116187133 B CN 116187133B
- Authority
- CN
- China
- Prior art keywords
- grid
- spring
- nodes
- spinning
- encryption
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 85
- 238000009987 spinning Methods 0.000 title claims abstract description 38
- 238000000926 separation method Methods 0.000 title claims abstract description 16
- 230000007704 transition Effects 0.000 claims abstract description 35
- 238000004364 calculation method Methods 0.000 claims abstract description 20
- 239000011159 matrix material Substances 0.000 claims description 9
- 238000012935 Averaging Methods 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 3
- 238000004088 simulation Methods 0.000 abstract description 15
- 230000001133 acceleration Effects 0.000 abstract description 7
- 230000015556 catabolic process Effects 0.000 abstract description 2
- 238000006731 degradation reaction Methods 0.000 abstract description 2
- 238000012821 model calculation Methods 0.000 abstract description 2
- 238000012545 processing Methods 0.000 abstract description 2
- 229910000838 Al alloy Inorganic materials 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 235000015842 Hesperis Nutrition 0.000 description 1
- 235000012633 Iberis amara Nutrition 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 238000004904 shortening Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Algebra (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Feedback Control In General (AREA)
- Springs (AREA)
Abstract
The invention relates to a dimension separation spring comparison method for a spinning mobile grid method, belonging to the field of spinning forming processing numerical simulation; firstly, designing a size field of a workpiece grid according to the position of a rotary wheel, then generating a variable-density full-quadrilateral grid on a cylindrical surface according to the size field by a dimension separation spring comparison method, projecting the grid onto the workpiece to generate inner and outer surface grids, finally obtaining internal nodes through linear interpolation, and connecting the nodes to obtain the final variable-density grid. According to the method, the dimension separation method and the cylindrical surface degradation grid are adopted to greatly reduce the calculation amount of grid node position iteration, and the linear spring comparison method is adopted to effectively reduce the nonlinearity degree of an iteration equation, so that the smooth transition local encryption full hexahedral grid can be generated rapidly. In the variable density grid, the model calculation scale can be greatly reduced by more than 70% due to the adoption of a point encryption form, so that simulation acceleration is up to 4 times.
Description
Technical Field
The invention belongs to the field of spinning forming processing numerical simulation, relates to a dimension separation spring comparison method for a spinning movable grid method, and particularly relates to a dimension separation spring comparison method for rapidly changing the grid density of a workpiece in a barrel-type workpiece flow spinning simulation.
Background
With the rapid development of high-end equipment in the fields of aviation, aerospace and the like, particularly the development of new-generation missiles, rockets and the like, higher requirements are put forward on the long service life, light weight and high reliability of key components of the high-performance long-barrel-shaped aluminum alloy large-diameter thin-wall long-barrel-shaped components integrally formed by utilizing a flow spinning technology, and the high-performance requirements can be met. However, when the finite element method (Finite Element Method, FEM) is applied to the flow spinning forming process for process research, the local point loading characteristic and the extreme size collocation of the forming area and the workpiece of the process limit the unit size of the model, and meanwhile, the extremely small time increment step and the frequent grid repartition exacerbate the simulation time consumption, so that the process only needs tens of seconds and the simulation time of weeks to months. Therefore, the research proposes that the efficient flow spinning simulation technology has important significance for accelerating process design and shortening the research and development period of high-end aerospace equipment.
The incremental forming characteristic of spinning can be fully utilized by adopting grid density control, and the whole scale of the model can be effectively reduced by only keeping the density of the grid in a plastic deformation area contacted with the spinning roller, so that the simulation time consumption is greatly reduced. However, the adoption of the grid density control method in the spinning process also has the following two problems: (1) The structured full hexahedral mesh adopted by the workpiece to ensure the precision is difficult to control the density of high quality; (2) three-dimensional grid mesh reconstruction is relatively time consuming. Therefore, there is an urgent need for a fast high quality full hexahedral mesh density control algorithm that can rapidly generate a variable density mesh of a workpiece according to the position of a spinning wheel.
Disclosure of Invention
The technical problems to be solved are as follows:
in order to avoid the defects of the prior art, the invention provides a dimension separation spring comparison method for a spinning moving grid method, which aims to construct a smooth transition full hexahedral grid model only locally encrypted in a plastic deformation area of a flow spinning workpiece, so that the calculation scale of the model is greatly reduced, and the accurate and efficient simulation of the model is ensured. Meanwhile, the self-adaptive high-fidelity dynamic reconstruction according to the original model geometric shape in the forming process of the grid model and the field variable mapping based on the indirect data transmission scheme meet the calculation continuity of the model sequence.
The technical scheme of the invention is as follows: a dimension separation spring comparison method for spinning a mobile grid method comprises the following specific steps:
s1: determining a grid size field;
s11: determining the positions of an encryption area and a transition area;
s12: determining a cell size field;
s2: generating a smooth transition partial encryption full-quadrilateral grid on the cylindrical surface according to the cell size field;
s21: arranging control nodes with single dimension;
s22: iteratively determining the position distribution of the control nodes according to a spring comparison method;
s23: generating a smooth transition partial encryption quadrilateral mesh by the control node;
s3: constructing a variable-density full hexahedral calculation grid of the workpiece according to the cylindrical grid;
s31: extracting axial lines of the inner surface and the outer surface of the original grid of the workpiece;
s32: intercepting an axial line by utilizing an axial control node of the cylindrical grid to obtain a series of contour nodes;
s33: reconstructing circumferential lines of the inner surface and the outer surface by adopting a cubic spline curve of a non-kinking boundary condition according to the contour nodes;
s34: intercepting the reconstructed circumferential line by using a circumferential control node of the cylindrical grid, and creating nodes on the inner surface and the outer surface of the grid;
s35: based on the newly generated internal and external surface nodes, constructing internal nodes through linear interpolation, and creating hexahedral units according to topological relations among the structured grid nodes.
The invention further adopts the technical scheme that: the position of the encryption area in S11 is determined:
the position of the spinning wheel is expressed in the cylindrical coordinate system as (ρ, θ, z), assuming that the plastic deformation of the spinning wheel is localized at [ θ - θ ] per moment l ，θ+θ r ]×[z-z b ，z+z t ]Then the envelope [ theta-omega delta T-theta ] of the rotor influence region is within delta T time l ，θ+θ r ]×[z-z b ，z+vΔT+z t ]Represents the range of the encryption zone, wherein omega is the rotation speed of the core mold, v is the feeding speed of the rotary wheel and theta l And theta r The radian value z of the contact position of the rotary wheel from the left and right boundaries of the influence area t And z b Is the dimension of the contact position of the rotor from the upper and lower boundaries of the affected area.
The invention further adopts the technical scheme that: wherein θ l 、θ r 、z t And z b Determining the size of the rotor influence area, theta l And theta r The value is within 15-30 DEG, z b At least take on a value of 10mm, z t The value should be able to include material stacking and non-sticking, at least 10mm under small pressure.
The invention further adopts the technical scheme that: the position of the transition zone in S11 is determined:
the encryption zone RZ extends circumferentially and axially outwardly along the boundary by θ trans And z trans Is defined as a transition zone, and assuming that there are n cells in the transition zone, the size of the cells in the encryption zone is s0, and the magnification of the cells is m, the size of the transition zone should be [ ns ] 0 ，nms 0 ]Within the range of (1) is taken as
The invention further adopts the technical scheme that: the method for determining the cell size field in S12, where the size of the cells in the encryption area is S 0 Less than or equal to (alpha+beta) r/3, wherein alpha, beta and r are attack angle, exit angle and fillet radius of the spinning wheel respectively; the cell size in the coarse grid area is ms 0 The method comprises the steps of carrying out a first treatment on the surface of the The cell size in the transition zone is from s 0 By ms 0 Is a linear transition of (2); the distance between the current position and the encryption area is defined in a single dimension as follows:
wherein delta 1 And delta 2 For the boundary of the encryption area in the corresponding direction, l is the transition area range in the corresponding direction, and theta is respectively in the circumferential direction and the axial direction trans And z trans The method comprises the steps of carrying out a first treatment on the surface of the The cell size field for the circumferential two-way single dimension can be defined as:
the invention further adopts the technical scheme that: the method for arranging the control nodes with single dimension in the S21 is based on the grid seed number N of the two directions of the grid circumference axis T And N A Uniformly distributed equidistant control nodes along two directions are arranged; the number of seed points is determined according to the formula (3), as follows:
wherein R, T and L are respectively the inner radius, thickness and length of the cylinder blank;
the invention further adopts the technical scheme that: the method for iteratively determining the position distribution of the control node according to the spring comparison method in S22 is that the virtual spring force applied to the control node i is f assuming that the control nodes are connected by virtual linear springs i Assembled force vector F (X) = { F 1 ，f 2 ，…，f N } T The method comprises the steps of carrying out a first treatment on the surface of the When F (X) =0, the grid density distribution reaches the standard, and a normal differential equation is constructed according to the grid density distributionAnd recursively calculating X by a forward Euler differential format (k+1) ＝X (k) +ΔτF (n) Until convergence; wherein Deltaτ is the virtual time increment step of the balance equation, and 10 is taken -4 A constant value of s;
the invention further adopts the technical scheme that: the method for determining the position distribution of the control node comprises the following specific steps:
s221: initializing; let k=0, the initial control node position is X (0) ；
S222: calculating the current length and the original length of a virtual spring connected with a control node; the current length of the ith spring is calculated as (4) 1 Obtaining by the difference between the coordinates of adjacent control nodes; the original length of the ith spring is required to meet the size distribution of cells in the grid, and the cell size values at the two end points of the spring are required to be determined according to the value (4) 2 Averaging;
s223: calculating virtual spring force of the control node; the stress f of the ith spring is calculated according to the formula (5) by adopting a linear spring with better convergence i s(k) Wherein K is the stiffness coefficient of the virtual spring, 10 4 The method comprises the steps of carrying out a first treatment on the surface of the The stress of the ith control node is provided by the spring connected with the ith control node, which isWherein->
S224: applying fixed boundary conditions to cause
S225: updating control node position X (k+1) ＝X (k) +ΔτF (n) ；
S226: if I X (k+1) -X (k) || 2 Iterative convergence is performed when TOL is less than TOL; otherwise k=k+1 and S222.
The invention further adopts the technical scheme that: the specific steps of generating the smooth transition local encryption quadrilateral mesh by the control node in the step S23 are as follows:
s231: generating control points of variable density in the x-direction and the y-direction by the unit size field and forming longitudinal lines and transverse lines;
s232: calculating the intersection point coordinates of the longitudinal line and the transverse line; the position (I, j) of the I-th node on the structured grid is calculated by formula (6), where N is the number of circumferential elements; the coordinates of the node are (x i ，y j )；
S233: constructing quadrilateral units on the cylindrical surface; the position (I, j) of the ith cell on the structured grid is calculated by formula (7), the node of which is taken in the clockwise direction as { (I, j), (i+1, j+1), (I, j+1) };
the invention further adopts the technical scheme that: in S33, a cubic spline curve formed by contour nodes is written as r i ＝R i (θ i )＝a i +b i (θ-θ i )+c i (θ-θ i ) 2 +d i (θ-θ i ) 3 ，i＝0，1，…，n-1；
Let h i ＝r i+1 -r i The formula is defined by the cubic spline condition:
conversion intoBinding to non-kinking boundary conditions->And->Writing the segmented spline curve into a matrix form:
Ax＝b (9)
x＝{m 0 ，m 1 ，…，m n } T
the coefficient matrix in the formula (9) is a tri-diagonal matrix, and the coefficients of the spline curve are rapidly solved by a catch-up method.
Advantageous effects
The invention has the beneficial effects that: according to the method, the dimension separation method and the cylindrical surface degradation grid are adopted to greatly reduce the calculation amount of grid node position iteration, and the linear spring comparison method is adopted to effectively reduce the nonlinearity degree of an iteration equation, so that the smooth transition local encryption full hexahedral grid can be generated rapidly. In the variable density grid, the model calculation scale can be greatly reduced by more than 70% due to the adoption of a point encryption form, so that simulation acceleration is up to 4 times. Meanwhile, the model ensures the accuracy of plastic deformation simulation of a contact area below the rotary wheel through local encryption, and the geometric detail characteristics of the model are reserved through spline curves, so that the geometric shape deviation of the method provided by the invention and a conventional fine grid computing model in the forming process is not more than 0.01mm.
In addition, the invention is not limited to the specific geometric shapes of the FEM solver and the deformation body, so that the FEM solver and the deformation body can be expanded to the rapid numerical simulation of any FEM software platform and any local point loading incremental forming process, and have certain universality.
Experiments prove that the calculation time of each model on the Intel Core i9-10900K CPU is shown in FIG. 6, so that the calculation acceleration of the method can be 2-6 times by controlling related parameters, and the multi-Core CPU can further improve the calculation efficiency of the model on the basis of the acceleration of the method.
Drawings
FIG. 1 is a flow chart of a spin-on moving mesh method calculation using a dimension-separation spring comparison method;
FIG. 2 is a schematic diagram of a variable density cylindrical mesh generated based on a dimension separation spring comparison method;
FIG. 3 is a schematic diagram of three-dimensional mesh generation of a workpiece based on a cylindrical mesh;
FIG. 4 is a single spin, rotary spinning model generated in accordance with the method of the present invention;
FIG. 5 is a comparison of geometric accuracy of a single spin, rotary spinning model according to the method of the present invention with a conventional fine grid;
fig. 6 is a comparison of the calculation time of a single spin-turn spinning model of a conventional fine grid according to the method of the present invention.
Detailed Description
The embodiments described below by referring to the drawings are illustrative and intended to explain the present invention and should not be construed as limiting the invention.
The dimension separation spring comparison method for the spinning mobile grid method aims at constructing a smooth transition full hexahedral grid model which is only locally encrypted in a plastic deformation area of a flow spinning workpiece, so that the calculation scale of the model is greatly reduced, and accurate and efficient simulation of the model is ensured. Meanwhile, the self-adaptive high-fidelity dynamic reconstruction according to the original model geometric shape in the forming process of the grid model and the field variable mapping based on the indirect data transmission scheme meet the calculation continuity of the model sequence.
Firstly, designing a size field of a workpiece grid according to the position of a rotary wheel, then generating a variable-density full-quadrilateral grid on a cylindrical surface according to the size field by a dimension separation spring comparison method, projecting the grid onto the workpiece to generate inner and outer surface grids, finally obtaining internal nodes through linear interpolation, and connecting the nodes to obtain the final variable-density grid.
The FEM simulation overall flow adopting the method is shown in fig. 1, and the method comprises the following specific steps:
s1: and (5) determining a grid size field.
S11: the positions of the encryption Zone (RZ) and the transition Zone are determined. The position of the spinning wheel can be expressed in cylindrical coordinates as (ρ, θ, z), assuming that the plastic deformation of the spinning wheel is localized to [ θ - θ ] at each moment l ，θ+θ r ]×[z-z b ，z+z t ]The range of the encryption zone may be represented by the envelope of the spin wheel influence zone [ theta-omega delta T-thetal, theta + thetar x z-zb, z + v delta T + zt ] in deltat time. Wherein omega is the rotation speed of the core mold, v is the feeding speed of the rotary wheel and theta l And theta r The radian value z of the contact position of the rotary wheel from the left and right boundaries of the influence area t And z b Is the dimension of the contact position of the rotor from the upper and lower boundaries of the affected area. θ l 、θ r 、z t And z b Determining the size of the rotor influence area, theta l And theta r The value is 15-30 degrees, z b At least take on a value of 10mm, z t The value should be able to include material stacking and non-sticking, at least 10mm under small pressure. RZ extends circumferentially and axially, respectively, outwardly along the boundary θ trans And z trans Is defined as a transition zone, provided that there are n cells in the transition zone, the size of the cells in the encryption zone is s 0 The magnification of the cell is m, the size of the transition zone should be [ ns 0 ，nms 0 ]Within the range of (1) can be taken as
S12: the cell size field is determined. The size of the cells in the encryption zone is s 0 And less than or equal to (alpha+beta) r/3, wherein alpha, beta and r are attack angle, exit angle and fillet radius of the spinning wheel respectively. The cell size in the coarse grid area is ms 0 . The cell size in the transition zone is from s 0 By ms 0 Is a linear transition of (c). Defining the distance of the current position from the encryption area in a single dimension as
Wherein delta 1 And delta 2 For the boundary of the encryption area in the corresponding direction, l is the transition area range in the corresponding direction, and theta is respectively in the circumferential direction and the axial direction trans And z trans . The unit size field of the circumferential two-way single dimension can be defined as
S2: generating a smooth transition partial encryption full quadrilateral mesh on the cylinder according to the cell size field.
S21: control nodes of a single dimension are arranged. According to the grid seed number N of the two directions of the grid circumference axis T And N A And uniformly distributed equidistant control nodes along two directions are arranged. The number of seed points can be determined according to a formula (3), wherein R, T and L are respectively the inner radius, thickness and length of the cylinder blank.
S22: and iteratively determining the position distribution of the control nodes according to a spring comparison method. Assuming that the control nodes are connected by virtual linear springs, the virtual spring force applied to the control node i is f i Assembling the available force vectorF(X)＝{f 1 ，f 2 ，…，f N } T . When F (X) =0, the grid density distribution reaches the standard, and a normal differential equation can be constructed according to the grid density distributionAnd recursively calculating X by a forward Euler differential format （k+1） ＝X （k） +ΔτF （n） Until convergence. Where Δτ is the virtual time increment step of the equilibrium equation, which may be 10 -4 s, a constant value of s. The method comprises the following specific steps:
s221: initializing. Let k=0, the initial control node position is X (0) 。
S222: and calculating the current length and the original length of the virtual spring connected with the control node. The current length of the ith spring is calculated as (4) 1 Obtaining by the difference between the coordinates of adjacent control nodes; the original length of the ith spring is required to meet the size distribution of cells in the grid, and the cell size values at the two end points of the spring are required to be determined according to the value (4) 2 And (5) averaging.
S223: a virtual spring force of the control node is calculated. The stress f of the ith spring can be calculated according to the formula (5) by adopting the linear spring with better convergence i s(k) Where K is the stiffness coefficient of the virtual spring, preferably 10 4 . The stress of the ith control node is provided by the spring connected with the ith control node, which isWherein->
S224: applying fixed boundary conditions to cause
S225: updating control node position X (k+1) ＝X (k) +ΔτF (n) 。
S226: if I X (k+1) -X (k) || 2 Iterative convergence is performed when TOL is less than TOL; otherwise k=k+1 and S222.
S23: a smooth transition partially encrypted quadrilateral mesh (see fig. 2) is generated by the control node, comprising the following steps:
s231: variable density control points are generated by the cell size field in the x-direction and the y-direction and form longitudinal and transverse lines.
S232: and calculating the intersection point coordinates of the longitudinal line and the transverse line. The position (I, j) of the I-th node on the structured grid can be calculated by equation (6), where N is the number of circumferential elements. The coordinates of the node are (x i ，y j )。
S233: quadrilateral cells on the cylinder are constructed. The position (I, j) of the ith cell on the structured grid can be calculated by equation (7) with its node taken in the clockwise direction as { (I, j), (i+1, j+1), (I, j+1) }.
S3: and constructing a variable-density full hexahedral calculation grid of the workpiece according to the cylindrical grid. The calculation grid can be directly generated by radial scanning of the cylindrical grid due to the regular shape of the workpiece at the initial moment. The calculation grid generation at the rest time is required to be according to the steps shown in fig. 3 so as to ensure the consistency of the model before and after grid transformation.
S31: and extracting axial lines of the inner surface and the outer surface of the original grid of the workpiece.
S32: and intercepting an axial line by utilizing the axial control nodes of the cylindrical grid to obtain a series of contour nodes.
S33: and reconstructing the circumferential lines of the inner surface and the outer surface by adopting a cubic spline curve of a non-kinking boundary condition according to the contour nodes. A cubic spline curve consisting of contour nodes can be written as r i ＝R i (θ i )＝a i +b i (θ-θ i )+c i (θ - θi2+di θ - θi3, i=0, 1, …, n-1. Let hi=ri+1-ri, the formula can be defined by a cubic spline condition
Conversion intoIn combination with non-kinking boundary condition R '' 0 (θ 1 )＝R″′ 1 (θ 2 ) With R'. n-2 (θ n-1 )＝R″′ n-1 (θ n-1 ) The segmented spline curve can be written in a matrix form
Ax＝b#(9)
x＝{m 0 ，m 1， …，m n } T
The coefficient matrix in the formula (9) is a tri-diagonal matrix, and the coefficients of the spline curve can be rapidly solved by a catch-up method.
S34: and intercepting the reconstructed circumferential line by using the circumferential control nodes of the cylindrical grid, and creating the nodes of the inner surface and the outer surface of the grid.
S35: based on the newly generated internal and external surface nodes, constructing internal nodes through linear interpolation, and creating hexahedral units according to topological relations among the structured grid nodes.
Examples:
the example is a 2219 aluminum alloy cylindrical part single-spin rotary spinning fast calculation model based on ABAQUS6.16 under windows10 operating system, the geometric and technological parameters of the model are listed in table 1, and the parameters of dimension separation springs are listed in table 2 compared with those of a moving grid method.
Table 1 geometry and process parameters of single-turn wheel model
Table 2 dimension split spring compares to mobile grid parameters
The full hexahedral variable density grid model generated according to the above parameters is shown in fig. 4, and a full-encrypted fine grid (number of cells: 350 x 3 x 50) is constructed for convenience of comparison of simulation accuracy and efficiency of the model. By projecting the surface nodes of the locally encrypted mesh at each instant onto the fine mesh surface, it can be seen that the geometric deviation of the two models is very low under the conditions shown in tables 1 and 2 (fig. 5). By adjusting z of the rotor influence area t And z b Parallel computing of 2 cores (AR 3-2) and 4 cores (AR 3-4) is enabled for AR3 for 15mm (AR 1), 10mm (AR 2), and 5mm (AR 3) to compare acceleration performance compared to the serial computing model AR 3-1. The calculation time of each model on the Intel Core i9-10900K CPU is shown in FIG. 6, so that the calculation acceleration of the method can be 2-6 times by controlling the related parameters, and the multi-Core CPU can further improve the calculation efficiency of the model on the basis of the acceleration of the method.
Although embodiments of the present invention have been shown and described above, it will be understood that the above embodiments are illustrative and not to be construed as limiting the invention, and that variations, modifications, alternatives, and variations may be made in the above embodiments by those skilled in the art without departing from the spirit and principles of the invention.
Claims (6)
1. A dimension separation spring comparison method for spinning mobile grid method is characterized by comprising the following specific steps:
s1: determining a grid size field;
s11: determining the positions of an encryption area and a transition area;
and determining the position of the encryption area:
the position of the spinning wheel is expressed in the cylindrical coordinate system as (ρ, θ, z), assuming that the plastic deformation of the spinning wheel is localized at [ θ - θ ] per moment l ,θ+θ r ]×[z-z b ,z+z t ]Then the envelope [ theta-omega delta T-theta ] of the rotor influence region is within delta T time l ,θ+θ r ]×[z-z b ,z+vΔT+z t ]Represents the range of the encryption zone, wherein omega is the rotation speed of the core mold, v is the feeding speed of the rotary wheel and theta l And theta r The radian value z of the contact position of the rotary wheel from the left and right boundaries of the influence area t And z b The dimension of the contact position of the rotating wheel from the upper boundary and the lower boundary of the influence area;
and determining the position of the transition zone:
the encryption zone RZ extends circumferentially and axially outwardly along the boundary by θ trans And z trans Is defined as a transition zone, provided that there are n cells in the transition zone, the size of the cells in the encryption zone is s 0 The magnification of the cell is m, the size of the transition zone should be [ ns 0 ,nms 0 ]Within the range of (1) is taken as
S12: determining a cell size field;
s2: generating a smooth transition partial encryption full-quadrilateral grid on the cylindrical surface according to the cell size field;
s21: arranging control nodes with single dimension;
s22: iteratively determining the position distribution of the control nodes according to a spring comparison method;
the method for iteratively determining the position distribution of the control nodes according to the spring comparison method is that the virtual spring force applied to the control node i is f under the assumption that the control nodes are connected by virtual linear springs i Assembled force vector F (X) = { F 1 ,f 2 ,…,f N } T The method comprises the steps of carrying out a first treatment on the surface of the When F (X) =0, the grid density distribution reaches the standard, and a normal differential equation is constructed according to the grid density distributionAnd recursively calculating X by a forward Euler differential format (k+1) ＝X (k) +ΔτF (n) Until convergence; wherein Deltaτ is the virtual time increment step of the balance equation, and 10 is taken -4 A constant value of s;
the method for determining the position distribution of the control node comprises the following specific steps:
s221: initializing; let k=0, the initial control node position is X (0) ；
S222: calculating the current length and the original length of a virtual spring connected with a control node; the current length of the ith spring is calculated as (4) 1 Obtaining by the difference between the coordinates of adjacent control nodes; the original length of the ith spring is required to meet the size distribution of cells in the grid, and the cell size values at the two end points of the spring are required to be determined according to the value (4) 2 Averaging;
s223: calculating virtual spring force of the control node; the stress of the ith spring is calculated according to the formula (5) by adopting a linear spring with better convergenceWherein K is the stiffness coefficient of the virtual spring, 10 4 The method comprises the steps of carrying out a first treatment on the surface of the The stress of the ith control node is provided by the spring connected with the ith control node, which is +.>Wherein->
S224: applying fixed boundary conditions to cause
S225: updating control node position X (k+1) ＝X (k) +ΔτF (n) ；
S226: if I X (k+1) -X (k) || 2 <TOL is iteration converged; otherwise k=k+1 and S222;
s23: generating a smooth transition partial encryption quadrilateral mesh by the control node;
s3: constructing a variable-density full hexahedral calculation grid of the workpiece according to the cylindrical grid;
s31: extracting axial lines of the inner surface and the outer surface of the original grid of the workpiece;
s32: intercepting an axial line by utilizing an axial control node of the cylindrical grid to obtain a series of contour nodes;
s33: reconstructing circumferential lines of the inner surface and the outer surface by adopting a cubic spline curve of a non-kinking boundary condition according to the contour nodes;
s34: intercepting the reconstructed circumferential line by using a circumferential control node of the cylindrical grid, and creating nodes on the inner surface and the outer surface of the grid;
s35: based on the newly generated internal and external surface nodes, constructing internal nodes through linear interpolation, and creating hexahedral units according to topological relations among the structured grid nodes.
2. A dimension splitting spring analogy method for spinning moving mesh method according to claim 1, characterized in that: wherein θ l 、θ r 、z t And z b Determining the size of the rotor influence area, theta l And theta r The value is within 15-30 DEG, z b At least take on a value of 10mm, z t The value should be able to include material stacking and non-sticking, at least 10mm under small pressure.
3. A dimension splitting spring analogy method for spinning moving mesh method according to claim 1, characterized in that: the method for determining the cell size field in S12, where the size of the cells in the encryption area is S 0 Less than or equal to (alpha+beta) r/3, wherein alpha, beta and r are attack angle, exit angle and fillet radius of the spinning wheel respectively; the cell size in the coarse grid area is ms 0 The method comprises the steps of carrying out a first treatment on the surface of the The cell size in the transition zone is from s 0 By ms 0 Is a linear transition of (2); the distance between the current position and the encryption area is defined in a single dimension as follows:
wherein delta 1 And delta 2 For the boundary of the encryption area in the corresponding direction, l is the transition area range in the corresponding direction, and theta is respectively in the circumferential direction and the axial direction trans And z trans The method comprises the steps of carrying out a first treatment on the surface of the The cell size field for the circumferential two-way single dimension can be defined as:
4. a dimension splitting spring analogy method for spinning moving mesh method according to claim 3, characterized in that: the method for arranging the control nodes with single dimension in the S21 is based on the grid seed number N of the two directions of the grid circumference axis T And N A Uniformly distributed equidistant control nodes along two directions are arranged; the number of seed points is determined according to the formula (3), as follows:
wherein R, T and L are the inner radius, thickness and length of the cylinder blank respectively.
5. A dimension splitting spring analogy method for spinning moving mesh method according to claim 4, characterized in that: the specific steps of generating the smooth transition local encryption quadrilateral mesh by the control node in the step S23 are as follows:
s231: generating control points of variable density in the x-direction and the y-direction by the unit size field and forming longitudinal lines and transverse lines;
s232: calculating the intersection point coordinates of the longitudinal line and the transverse line; the position (I, j) of the I-th node on the structured grid is calculated by formula (6), where N is the number of circumferential elements; the coordinates of the node are (x i ,y j )；
S233: constructing quadrilateral units on the cylindrical surface; the position (I, j) of the ith cell on the structured grid is calculated by formula (7), the node of which is taken in the clockwise direction as { (I, j), (i+1, j+1), (I, j+1) };
6. a dimension splitting spring analogy method for spinning moving mesh method according to claim 5, characterized in that: in S33, a cubic spline curve formed by contour nodes is written as r i ＝R i (θ i )＝a i +b i (θ-θ i )+c i (θ-θ i ) 2 +d i (θ-θ i ) 3 ，i＝0,1,…,n-1；
Let h i ＝r i+1 -r i The formula is defined by the cubic spline condition:
conversion intoIn combination with non-kinking boundary condition R '' 0 (θ 1 )＝R″′ 1 (θ 2 ) With R'. n-2 (θ n-1 )＝R″′ n-1 (θ n-1 ) Writing the segmented spline curve into a matrix form:
Ax＝b (9)
x＝{m 0 ,m 1 ,…,m n } T
the coefficient matrix in the formula (9) is a tri-diagonal matrix, and the coefficients of the spline curve are rapidly solved by a catch-up method.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310098738.5A CN116187133B (en) | 2023-02-10 | 2023-02-10 | Dimension separation spring comparison method for spinning mobile grid method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310098738.5A CN116187133B (en) | 2023-02-10 | 2023-02-10 | Dimension separation spring comparison method for spinning mobile grid method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116187133A CN116187133A (en) | 2023-05-30 |
| CN116187133B true CN116187133B (en) | 2023-10-27 |
Family
ID=86451836
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310098738.5A Active CN116187133B (en) | 2023-02-10 | 2023-02-10 | Dimension separation spring comparison method for spinning mobile grid method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116187133B (en) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010067170A (en) * | 2008-09-12 | 2010-03-25 | Hitachi Ltd | Analyzer, and mesh generation program |
| CN103226634A (en) * | 2013-04-19 | 2013-07-31 | 华南理工大学 | Computing method for unsteady flow field of rotary jet pump based on three-dimensional dynamic mesh |
| CN103699727A (en) * | 2013-12-17 | 2014-04-02 | 华中科技大学 | Power spinning spinnability analysis value simulating method |
| CN105183996A (en) * | 2015-09-14 | 2015-12-23 | 西北工业大学 | Surface element correction and grid beforehand self-adaption calculation method |
| CN107895057A (en) * | 2017-05-26 | 2018-04-10 | 宝沃汽车（中国）有限公司 | The numerical analysis method and system of rotary press modelling |
| CN108197418A (en) * | 2018-03-14 | 2018-06-22 | 上海理工大学 | A kind of hexahedron FEA Meshing Method for simulating thread fitting |
| CN112613209A (en) * | 2020-12-17 | 2021-04-06 | 上海交通大学 | Full hexahedron unit encryption method based on finite element model |
| JP2021140285A (en) * | 2020-03-03 | 2021-09-16 | 大成建設株式会社 | Mesh model generation apparatus and mesh model generation method |
| CN114969860A (en) * | 2022-06-15 | 2022-08-30 | 苏州流场信息技术有限公司 | Automatic hexahedron non-structural grid generation method |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| FR2870621B1 (en) * | 2004-05-21 | 2006-10-27 | Inst Francais Du Petrole | METHOD FOR GENERATING A THREE-DIMENSIONALLY THREADED HYBRID MESH OF A HETEROGENEOUS FORMATION CROSSED BY ONE OR MORE GEOMETRIC DISCONTINUITIES FOR THE PURPOSE OF MAKING SIMULATIONS |
-
2023
- 2023-02-10 CN CN202310098738.5A patent/CN116187133B/en active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010067170A (en) * | 2008-09-12 | 2010-03-25 | Hitachi Ltd | Analyzer, and mesh generation program |
| CN103226634A (en) * | 2013-04-19 | 2013-07-31 | 华南理工大学 | Computing method for unsteady flow field of rotary jet pump based on three-dimensional dynamic mesh |
| CN103699727A (en) * | 2013-12-17 | 2014-04-02 | 华中科技大学 | Power spinning spinnability analysis value simulating method |
| CN105183996A (en) * | 2015-09-14 | 2015-12-23 | 西北工业大学 | Surface element correction and grid beforehand self-adaption calculation method |
| CN107895057A (en) * | 2017-05-26 | 2018-04-10 | 宝沃汽车（中国）有限公司 | The numerical analysis method and system of rotary press modelling |
| CN108197418A (en) * | 2018-03-14 | 2018-06-22 | 上海理工大学 | A kind of hexahedron FEA Meshing Method for simulating thread fitting |
| JP2021140285A (en) * | 2020-03-03 | 2021-09-16 | 大成建設株式会社 | Mesh model generation apparatus and mesh model generation method |
| CN112613209A (en) * | 2020-12-17 | 2021-04-06 | 上海交通大学 | Full hexahedron unit encryption method based on finite element model |
| CN114969860A (en) * | 2022-06-15 | 2022-08-30 | 苏州流场信息技术有限公司 | Automatic hexahedron non-structural grid generation method |
Non-Patent Citations (6)
| Title |
|---|
| Developing an Extended IFC Data Schema and Mesh Generation Framework for Finite Element Modeling;Zhao Xu et al.;《Advances in Civil Engineering》;1-19 * |
| Improvement of rib-grid structure of thin-walled tube with helical grid-stiffened ribs based on the multi-mode filling behaviors in flow forming;W. Lyu et al.;《Journal of Materials Processing Technology》;1-12 * |
| 三维有限元六面体网格几何自适应再生成方法;王忠雷等;《计算力学学报》(第01期);72-77 * |
| 基于拓扑分割的泵叶轮六面体网格划分方法;许磊等;《重庆工商大学学报(自然科学版)》(第04期);39-45 * |
| 塑性成形快速数值仿真方法的研究进展;詹梅等;《机械工程学报》;第28卷(第16期);2-20 * |
| 带环向内筋筒形件旋压成形工艺试验及缺陷分析;写旭等;《航天制造技术》(第5期);62-65 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116187133A (en) | 2023-05-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110610050B (en) | Airfoil aerodynamic drag reduction method based on improved radial basis function deformation algorithm | |
| CN112765732B (en) | Aviation blade topology optimization design method based on selective laser melting process | |
| CN116244988B (en) | High-quality quadrilateral grid conformal construction method for plate spinning multi-grid method simulation | |
| CN108319224B (en) | A kind of multiaxis NC maching spiral path generation method based on diametral curve interpolation | |
| CN111177906A (en) | Method for accurately compensating discrete die profile | |
| CN104573281A (en) | Complex space curved surface thin plate forming die face designing method taking springback compensation | |
| CN114444216B (en) | Aircraft attitude control method and system under high-altitude condition based on numerical simulation | |
| CN116187133B (en) | Dimension separation spring comparison method for spinning mobile grid method | |
| CN112286886B (en) | Multi-block structure grid data compression storage format | |
| CN105260563A (en) | Finite element prestressed mode analysis method for three-dimensional entity blade and two-dimensional axisymmetric wheel disc variable-dimensionality model of impeller structure | |
| CN112966454A (en) | Wind power plant fan wake flow dynamic coupling simulation method | |
| CN108763668B (en) | Gear model region parameterization method based on subdivision technology and boundary replacement | |
| CN117272870B (en) | Dynamic formation flight numerical simulation method based on self-adaptive overlapped grid | |
| CN109702931B (en) | Method for designing mold surface of computer-aided large-scale component precise hot forming mold | |
| CN115203997A (en) | Dot matrix-entity composite structure topology optimization method based on multivariate design | |
| CN112464531B (en) | B-spline parameterization-based reinforcement modeling and optimizing method for thin-wall structure | |
| Zheng et al. | An enhanced topology optimization approach based on the combined MMC and NURBS-curve boundaries | |
| CN112800557A (en) | Topological optimization method for transition plate of speed reducer of industrial robot | |
| CN114840926B (en) | Method for generating three-dimensional finite element grid with complex tire patterns | |
| CN116756851A (en) | Parameterized grid deformation method and system based on NFFD background grid | |
| CN111881629A (en) | Pneumatic heat-structure heat conduction coupling nonlinear reduced order model method | |
| CN109948253A (en) | The GPU accelerated method of thin plate mesh free Galerkin Constructional Modal Analysis | |
| CN104820729A (en) | Temperature succession method of whole rolling process simulation analysis of steel rail | |
| CN113158377B (en) | Cambered surface indexing cam model creation and transmission performance optimization design method and CAD (computer-aided design) optimization design system | |
| CN113093998A (en) | Space distributed storage optimization method based on geographical hash |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |