CN102760183A - Bidirectional evolutionary structure topological optimization based disk type flywheel optimized design method - Google Patents

Abstract

The invention discloses a bidirectional evolutionary structure topological optimization based disk type flywheel optimized design method. The disk type flywheel is made of a homogeneous material. The design method comprises the following steps: performing finite element grid division on an initial design domain of the disk type flywheel and giving a target volume, a volume evolution rate and a filtering radius of a structure; performing finite element analysis on the flywheel structure and calculating a stress sensitivity number of each unit; obtaining the sensitivity number of the current iteration and the structure target volume of the next iteration; and deleting and adding units until the target volume and the condition of convergence are satisfied at the same time. Through the bidirectional evolutionary structure topological optimization based disk type flywheel optimized design method, disclosed by the invention, practical flywheel optimized design can be guided rationally, and therefore, mass of the flywheel can be reduced on the basis of a certain intensity.

Description

A kind of dish-type flywheel Optimization Design based on two-way progressive structure topological optimization

Technical field

The present invention relates to mechanical industry energy storing structure design field, relate in particular to a kind of homogeneous material dish-type flywheel structure method of topological optimization design.

Background technology

Structure optimization not only can effectively improve existing structure, but also important, valuable theoretical direction can be provided for the conceptual design that exceeds deviser's experience.Structural Topology Optimization is optimization objects with the distribution of material, in the equally distributed design space of material, searches best distribution scheme.Topological optimization is the advanced stage of structural optimization problems; Its difficulty will be far above dimensionally-optimised and Shape optimization; But along with the lifting of computer computation ability and the development of finite element technique, the structural Topology Optimization technology not only becomes possibility, and has grown up gradually.At present; Representative Topology Optimization Method mainly contains homogenising (Homogenization) method, SIMP (Solid Isotropic Material with Penalization) method, level set (level set) method, progressive structure optimization (Evolutionary Structural Optimization, ESO) method and ICM (Independent Continuous Method) method etc.Compare with other several method, the progressive structure Topology Optimization Method has that notion is simple, convergence is fast, sharpness of border, the characteristics that are easy to realize, at present, in engineering design, has obtained application more and more widely.

For homogeneous material dish-type flywheel, under certain constraint condition, only depend on Shape optimization often can not obtain the most economical and the most practical flywheel structure, can not obtain through Shape optimization like common disc formula, spoke type dish-type flywheel in the life.Can what shape where " hole " dig, dig just make result optimal for disc or spoked flywheel? This has just related to the structure topology optimization problem.Therefore,, flywheel structure is carried out topological structure optimization,, improve the flywheel performance, have actual directive significance to reduce flywheel weight to the problem and shortage that traditional flywheel design method exists.

Summary of the invention

The technical matters that the present invention will solve is: to the technical matters that prior art exists, the present invention proposes the Topology Optimization Method that a kind of convergence to the dish-type flywheel structure is faster, be easier to realize.

For addressing the above problem, the present invention adopts following technical scheme:

A kind of based on two-way progressive structure topological optimization (Bi-directional Evolutionary Structural Optimization, flywheel Optimization Design BESO) is characterized in that, step is:

(1) finite element grid is carried out in flywheel initial designs territory and divide, and give the target volume V of fixed structure ^*, volume evolution rate and filter radius;

(2) flywheel structure is carried out finite element analysis, calculate the stress sensitivity number of each unit and the sensitivity number of this iteration;

(3) confirm the structural object volume of next iteration; The target volume of supposing structural design is for give the target volume V of fixed structure ^*, then in process of topology optimization, deleting of unit added and can whether be satisfied target volume V according to structure residual volume V ^*Control;

(4) repeating step (2) ~ (3) are until target volume V ^*And condition of convergence formula is met simultaneously.

Adopt the performance judge index of stress mean square deviation in the said step (1) as two-way progressive structure topological optimization algorithm, through will be real with dummy cell by the unified way of lining up of sensitivity number size, realized that the unit deletes carrying out synchronously of adding.

Adopt the sensitivity filtering method effectively to improve the topological optimization result of BESO algorithm in the said step (2).

The sensitivity filtering method that in said step (2), adopts comprises the steps:

At first, the sensitivity number at each node place, i unit being asked on average, is the sensitivity number of this unit with this mean value:

${\mathrm{\α}}_i^e=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_i^n}{N}$

In the formula, α _i ^eBe the average sensitivity number at each node place, i unit, N is the node number of i unit;

Secondly, be the center of circle with the barycenter of i unit, with filter radius r _MinBe that radius makes a circle, establish and include K unit in this border circular areas, then after weighted filtering, the sensitivity numerical table of i unit is shown:

${\mathrm{\α}}_i=\frac{\underset{j=1}{\overset{K}{\mathrm{\Σ}}}w\left(r_{\mathrm{ij}}\right){\mathrm{\α}}_j^e}{\underset{j=1}{\overset{K}{\mathrm{\Σ}}}w\left(r_{\mathrm{ij}}\right)}$

W (r in the formula _Ij) be weight factor.

For quadrilateral units N=4, for hexahedral element N=8.

Weight factor w (r _Ij) be calculated as follows:

w(r _ij)＝r _min-r _ij,j＝1,2,…,K

R in the formula _IjDistance for barycenter barycenter of i of j unit in the border circular areas to the unit.

In order to guarantee that iterative process has better convergence, carry out once historical iteration again and ask average sensitivity number, promptly

${\mathrm{\α}}_i=\frac{{\mathrm{\α}}_i^k+{\mathrm{\α}}_i^{k-1}}{2}$

sensitivity number when being the k time iteration in i unit in the formula.

In the said step (3), when residual volume V greater than target volume V ^*The time, the BESO algorithm is main with delete cells mainly, makes it satisfy target volume V gradually to reduce V ^*, but in the optimizing process, the dummy cell of high sensitivity number still can be converted into real unit to be added in the structure; When residual volume V less than target volume V ^*The time, the BESO algorithm is main with adding device mainly, makes it satisfy target volume V gradually to increase V ^*But, still can from structure, the deleting of the real unit of the muting sensitive number of degrees.

In the evolutionary process, the volume of structure is expressed as:

V _k+1＝V _k(1±ER),(k＝1,2,3,…)

In the formula, ER is the volume evolution rate; K is an iterations.

In case structural volume V _kReached target volume V ^*, the volume of structure will be retained as V ^*Constant; The sensitivity that obtains each unit through step (3) is counted α _iAfter, ranked according to the size of sensitivity number from high to low or from low to high in all real unit and dummy cell; Suppose that deletion threshold values and interpolation threshold values are respectively

Then for satisfying

Real unit carry out deletion action, promptly become dummy cell, for satisfying

Dummy cell add operation, promptly become real unit; Threshold values

The step of confirming the face of pressing calculate:

1) for the k time iteration, if sequence in all unit of team front M ₁The volume of individual unit with equal V _K+1, if M ₁+ 1 unit sensitivity number is α ^Th, then

2) calculating is satisfied

The volume AR of dummy cell, if AR is less than the predefined maximum volume AR that adds _Max, then skip the 3rd) and the step, otherwise, change the 3rd) step, recomputate

3) all dummy cells are lined up by sensitivity number order from big to small, if front M ₀The volume of individual dummy cell with equal AR _Max, then with M ₀The sensitivity number of+1 dummy cell is made as

In order to guarantee the volume V of the k+1 time iteration _K+1, in like manner,

Can and equal (V through the delete cells volume _k-V _K+1+ AR _Max) method confirm.

Description of drawings

Fig. 1 is the schematic flow sheet of topological optimization of the present invention;

Fig. 2 is 1/6 model of dish-type flywheel initial configuration of the present invention;

Fig. 3 is the region of search synoptic diagram;

Fig. 4 a-d is the topologies change (=30%) of dish-type flywheel 1/12 structure in the iterative process of the present invention;

Fig. 5 a-b is the volume fraction and the stress mean square deviation convergence curve (=30%) of dish-type flywheel 1/12 structure;

Flywheel topological structure (1/12 model) after Fig. 6 a-c different target volume is optimized down;

Fig. 7 a-b illustrates the flywheel topological structure (=30%) after the different rotary cycle is optimized down.

Embodiment

Below will combine Figure of description and specific embodiment that the present invention is done further detailed description.

The dish-type flywheel of high speed rotating mainly bears the effect of centrifugal load.In application of practical project, on the flywheel usually along circumferentially being uniform-distribution with circular hole of the same size, to reduce the quality of flywheel.The flywheel initial configuration of the embodiment of the invention (1/6 model with the dish-type flywheel is an example) as shown in Figure 2; Wherein the inner of dish-type flywheel and outer end (zone 1 and zone 3) is non-design section; This is in order to guarantee that the dish-type flywheel can have certain moment of inertia simultaneously with the axle assembling; Zone 2 parts (except the circular hole) in the middle of only but be design section, the material parameter of dish-type flywheel is seen table 1, dimensional parameters is seen table 2.When having dug hole in advance on the flywheel, the flywheel optimal design has just belonged to the three-dimensional topology optimization problem.Fig. 2, according to plane symmetry, only needs to calculate the half the of 1/6 model for 1/6 model (swing circle is 60 °) of dish-type flywheel when optimal design, promptly 1/12 of flywheel gets final product.

The material parameter of table 1 dish-type flywheel

The basic size of table 2 dish-type flywheel

(1) finite element grid is carried out in flywheel initial designs territory and divide, and give target volume, volume evolution rate and the filter radius of fixed structure;

(2) flywheel structure is carried out finite element analysis, calculate the stress sensitivity number of each unit and the sensitivity number of this iteration;

It is limited unit grid that non-individual body is dispersed, and with causing the discontinuous of sensitivity number, this is the basic reason that topology optimization problem produces " gridiron pattern " phenomenon.In addition, the result of topological optimization often also can be different because of the thickness that unit grid is divided, promptly to the dependency problem of unit grid.For fear of producing " gridiron pattern " phenomenon, overcome dependence simultaneously to finite element grid, the present invention adopts the sensitivity filtering method, and roughly step is following:

At first, the sensitivity number at each node place, i unit being asked on average, is the sensitivity number of this unit with this mean value:

${\mathrm{\α}}_i^e=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_i^n}{N}---\left(1\right)$

In the formula (1), α _i ^eBe the average sensitivity number at each node place, i unit, N is the node number of i unit, for quadrilateral units N=4, for hexahedral element N=8.

Secondly, be the center of circle with the barycenter of i unit, with filter radius r _MinFor radius is made a circle, as shown in Figure 3, establish and include K unit in this border circular areas, then after weighted filtering, the sensitivity numerical table of i unit is shown:

${\mathrm{\α}}_i=\frac{\underset{j=1}{\overset{K}{\mathrm{\Σ}}}w\left(r_{\mathrm{ij}}\right){\mathrm{\α}}_j^e}{\underset{j=1}{\overset{K}{\mathrm{\Σ}}}w\left(r_{\mathrm{ij}}\right)}---\left(2\right)$

W (r in the formula _Ij) be weight factor, it can be calculated as follows:

w(r _ij)＝r _min-r _ij,j＝1,2,…,K (3)

R in the formula _IjDistance for barycenter barycenter of i of j unit in the border circular areas to the unit.

In order to guarantee that iterative process has better convergence, the sensitivity number can also carry out once historical iteration again to be asked on average, promptly

${\mathrm{\α}}_i=\frac{{\mathrm{\α}}_i^k+{\mathrm{\α}}_i^{k-1}}{2}---\left(4\right)$

sensitivity number when being the k time iteration in i unit in the formula.

(3) confirm the structural object volume of next iteration; The target volume of supposing structural design is V ^*, then in process of topology optimization, deleting of unit added and can whether be satisfied target volume V according to structure residual volume V ^*Control.When residual volume V greater than target volume V ^*The time, the BESO algorithm is main with delete cells mainly, makes it satisfy target volume V gradually to reduce V ^*, but in the optimizing process, the dummy cell of high sensitivity number still can be converted into real unit to be added in the structure; When residual volume V less than target volume V ^*The time, the BESO algorithm is main with adding device mainly, makes it satisfy target volume V gradually to increase V ^*But, still can from structure, the deleting of the real unit of the muting sensitive number of degrees.Therefore, in the evolutionary process, the volume of structure can be expressed as:

V _k+1＝V _k(1±ER), (k＝1,2,3,…) (5)

In the formula, ER is the volume evolution rate; K is an iterations.

In case structural volume V _kReached target volume V ^*, the volume of structure will be retained as V ^*Constant.The sensitivity that obtains each unit through step (3) is counted α _iAfter, with all real unit and dummy cell according to the size of sensitivity number from high to low (or from low to high) rank.Suppose that the deletion threshold values is respectively

with the interpolation threshold values and then carries out deletion action (becoming dummy cell) for the real unit that satisfies

, add operation (becoming real unit) for the dummy cell that satisfies .Confirming of threshold values

is also fairly simple, can calculate by following step:

1) for the k time iteration, if sequence in all unit of team front M ₁The volume of individual unit with equal V _K+1, if M ₁+ 1 unit sensitivity number is α ^Th, then

2) calculating is satisfied The volume AR of dummy cell, if AR is less than the predefined maximum volume AR that adds _Max(desirable structure collectivity long-pending 1%) then skips the 3rd) step, otherwise, change the 3rd) step, recomputate

3) all dummy cells are lined up by sensitivity number order from big to small, if front M ₀The volume of individual dummy cell with equal AR _Max, then with M ₀The sensitivity number of+1 dummy cell is made as

In order to guarantee the volume V of the k+1 time iteration _K+1, in like manner,

Can and equal (V through the delete cells volume _k-V _K+1+ AR _Max) method confirm.

What be worth explanation is, in this algorithm, all realities and dummy cell have all been participated in ordering, and therefore, the deleting to add simultaneously of unit carried out, and the method for deleting, adding separate operation of unit is compared, and this algorithm is simply much effective.

(4) repeat (2) ~ (3) step, until target volume V ^*And condition of convergence formula is met simultaneously.

In statistics, the dispersion degree of data acquisition can be come constant with its mean square deviation, and the big more declarative data of mean square deviation is overstepping the bounds of propriety looses, otherwise data are concentrated more.Obvious reason, whether the structural stress level is even, certainly judges with the mean square deviation of stress.Therefore, the performance index PI of BESO algorithm may be defined as:

$\mathrm{PI}=\sqrt{\frac{\underset{i=1}{\overset{N_k}{\mathrm{\Σ}}}{\left(\right|{\mathrm{\σ}}_i^{\mathrm{VM}}|-{\stackrel{\‾}{\mathrm{\σ}}}_k^{\mathrm{VM}})}^2}{N_k}}---\left(6\right)$

In the formula, N _kBe the k time iteratively-structured unit sum;

The mean stress of the k time each unit of iteration structure

So far, the topology optimization problem of improvement BESO algorithm may be defined as:

$\begin{array}{cc}\mathrm{Min}& \mathrm{PI}=\sqrt{\frac{\underset{i=1}{\overset{N_k}{\mathrm{\Σ}}}{\left(\right|{\mathrm{\σ}}_i^{\mathrm{VM}}|-{\stackrel{\‾}{\mathrm{\σ}}}_k^{\mathrm{VM}})}^2}{N_k}}\end{array}$

s.t.KU=f

$V^{*}-\underset{i=1}{\overset{N_k}{\mathrm{\Σ}}}{V}_i{x}_i=0,$ x _i＝0?or?1

(7)

In the formula, K, U, f are respectively integral rigidity matrix, motion vector and the outer force vector in the finite element analysis; V _iIt is the volume of i unit.As front and back n ₀The relative error error of the PI value sum of inferior iteration is less than given error value epsilon _eThe time, can be considered algorithm convergence, therefore, ε _eBe defined as:

$\mathrm{error}=\frac{|\underset{i=1}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{PI}}_{k-i+1}-\underset{i=1}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{PI}}_{k-n_0-i+1}|}{\underset{i=1}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{PI}}_{k-i+1}}\≤{\mathrm{\ϵ}}_e---\left(8\right)$

Fig. 4 a-d has provided the change procedure of dish-type flywheel 1/12 structural Topology Optimization in the iterative process, wherein, and the target volume V of flywheel remainder ^*Be controlled to be 30%, the unit is 87 * 29 * 15 hexahedral elements.Fig. 4 a shows the situation of k=0, and Fig. 4 b shows the situation of k=50, and Fig. 4 c shows the situation of k=110, and Fig. 4 d shows the situation of k=195, from Fig. 4 a-d, can obviously see the process of adding of deleting of unit.Fig. 5 a-b has provided the target volume of dish-type flywheel 1/12 structure and the convergence situation of stress mean square deviation PI in the optimizing process respectively.What deserves to be mentioned is that the part of stress mean square deviation PI thermal agitation causes owing to deletion " rod member " in the structure causes structural topology obviously to change among Fig. 5 b.

In addition, the target volume V of flywheel structure ^*Difference, optimum topological structure can be obviously different.Fig. 6 a-c has provided to the flywheel topological structure after optimizing under the different target volume of flywheel 1/12 model, and Fig. 6 a illustrates V ^*=30% topological structure, Fig. 6 b illustrate V ^*=40% topological structure, Fig. 6 c illustrate V ^*=50% topological structure can be found out along with target volume V ^*Constantly reduce, the optimum topological structure of dish-type flywheel becomes spoke type from the disc formula gradually, it is simple that structure is tending towards.In addition, swing circle is different, and the optimum topological structure of periodic symmetry structure also can be obviously different.Fig. 7 a-b has provided to volume fraction V ^*=30% different rotary is the flywheel topological structure after (1/8 and 1/18 model) optimization under the cycle; Fig. 7 a illustrates the flywheel topological structure of 1/8 model and the flywheel topological structure that Fig. 7 b illustrates 1/18 model, can find out the increase along with the revolution issue; Hole location is inwardly radial direction extension gradually; Simultaneously, the member structure in the dish-type flywheel also increases gradually, and structure is also complicated gradually.In the processing of reality, should take all factors into consideration the performance requirement and the production and processing cost thereof of structure, design, the more economical and practical topological structure of selection.

Claims (10)

1. dish-type flywheel Optimization Design based on two-way progressive structure topological optimization is characterized in that step is:

(1) finite element grid is carried out in flywheel initial designs territory and divide, and give the target volume V of fixed structure ^*, volume evolution rate and filter radius;

(2) flywheel structure is carried out finite element analysis, calculate the stress sensitivity number of each unit and the sensitivity number of this iteration;

(3) confirm the structural object volume of next iteration; The target volume of supposing structural design is for give the target volume V of fixed structure ^*, then in process of topology optimization, deleting of unit added and can whether be satisfied target volume V according to structure residual volume V ^*Control;

(4) repeating step (2) ~ (3) are until target volume V ^*And condition of convergence formula is met simultaneously.

2. the dish-type flywheel Optimization Design based on two-way progressive structure topological optimization according to claim 1; It is characterized in that: adopt the performance judge index of stress mean square deviation in the said step (1) as two-way progressive structure topological optimization algorithm; Through will be real with dummy cell by the unified way of lining up of sensitivity number size, realized that the unit deletes carrying out synchronously of adding.

3. the dish-type flywheel Optimization Design based on two-way progressive structure topological optimization according to claim 1 and 2 is characterized in that: adopt the sensitivity filtering method effectively to improve the topological optimization result of BESO algorithm in the said step (2).

4. the dish-type flywheel Optimization Design based on two-way progressive structure topological optimization according to claim 3 is characterized in that: the sensitivity filtering method that in said step (2), adopts comprises the steps:

At first, the sensitivity number at each node place, i unit being asked on average, is the sensitivity number of this unit with this mean value:

${\mathrm{\α}}_i^e=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_i^n}{N}$

In the formula, α _i ^eBe the average sensitivity number at each node place, i unit, N is the node number of i unit;

Secondly, be the center of circle with the barycenter of i unit, with filter radius r _MinBe that radius makes a circle, establish and include K unit in this border circular areas, then after weighted filtering, the sensitivity numerical table of i unit is shown:

${\mathrm{\α}}_i=\frac{\underset{j=1}{\overset{K}{\mathrm{\Σ}}}w\left(r_{\mathrm{ij}}\right){\mathrm{\α}}_j^e}{\underset{j=1}{\overset{K}{\mathrm{\Σ}}}w\left(r_{\mathrm{ij}}\right)}$

W (r in the formula _Ij) be weight factor.

5. the dish-type flywheel Optimization Design based on two-way progressive structure topological optimization according to claim 4 is characterized in that: for quadrilateral units N=4, for hexahedral element N=8.

6. the dish-type flywheel Optimization Design based on two-way progressive structure topological optimization according to claim 4 is characterized in that: weight factor w (r _Ij) be calculated as follows:

w(r _ij)＝r _min-r _ij,j＝1,2,…,K

R in the formula _IjDistance for barycenter barycenter of i of j unit in the border circular areas to the unit.

7. according to the described dish-type flywheel Optimization Design of one of claim 4-6, it is characterized in that: have better convergence in order to guarantee iterative process, carry out once historical iteration again and ask average sensitivity number, promptly based on two-way progressive structure topological optimization

${\mathrm{\α}}_i=\frac{{\mathrm{\α}}_i^k+{\mathrm{\α}}_i^{k-1}}{2}$

sensitivity number when being the k time iteration in i unit in the formula.

8. according to the described dish-type flywheel Optimization Design of one of claim 1-7, it is characterized in that based on two-way progressive structure topological optimization: in the said step (3), when residual volume V greater than target volume V ^*The time, the BESO algorithm is main with delete cells mainly, makes it satisfy target volume V gradually to reduce V ^*, but in the optimizing process, the dummy cell of high sensitivity number still can be converted into real unit to be added in the structure; When residual volume V less than target volume V ^*The time, the BESO algorithm is main with adding device mainly, makes it satisfy target volume V gradually to increase V ^*But, still can from structure, the deleting of the real unit of the muting sensitive number of degrees.

9. the dish-type flywheel Optimization Design based on two-way progressive structure topological optimization according to claim 8, it is characterized in that: in the evolutionary process, the volume of structure is expressed as:

V _k+1＝V _k(1±ER), (k＝1,2,3,…)

In the formula, ER is the volume evolution rate; K is an iterations.

10. the dish-type flywheel Optimization Design based on two-way progressive structure topological optimization according to claim 9 is characterized in that: in case structural volume Vk has reached target volume V ^*, the volume of structure will be retained as V ^*Constant; The sensitivity that obtains each unit through step (3) is counted α _iAfter, ranked according to the size of sensitivity number from high to low or from low to high in all real unit and dummy cell; Suppose that deletion threshold values and interpolation threshold values are respectively

Then for satisfying

Real unit carry out deletion action, promptly become dummy cell, for satisfying

Dummy cell add operation, promptly become real unit; Threshold values The step of confirming the face of pressing calculate:

1) for the k time iteration, if sequence in all unit of team front M ₁The volume of individual unit with equal V _K+1, if M ₁+ 1 unit sensitivity number is α ^Th, then

2) calculating is satisfied The volume AR of dummy cell, if AR is less than the predefined maximum volume AR that adds _Max, then skip the 3rd) and the step, otherwise, change the 3rd) step, recomputate

3) all dummy cells are lined up by sensitivity number order from big to small, if front M ₀The volume of individual dummy cell with equal AR _Max, then with M ₀The sensitivity number of+1 dummy cell is made as

In order to guarantee the volume V of the k+1 time iteration _K+1, in like manner,

Can and equal (V through the delete cells volume _k-V _K+1+ AR _Max) method confirm.

Images