CN111832203A - Graphical method for generating heat dissipation topology by zero-deficiency grid curved surface - Google Patents
Graphical method for generating heat dissipation topology by zero-deficiency grid curved surface Download PDFInfo
- Publication number
- CN111832203A CN111832203A CN202010632026.3A CN202010632026A CN111832203A CN 111832203 A CN111832203 A CN 111832203A CN 202010632026 A CN202010632026 A CN 202010632026A CN 111832203 A CN111832203 A CN 111832203A
- Authority
- CN
- China
- Prior art keywords
- grid
- formula
- heat
- fluid
- mesh
- 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.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- 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
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- 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/08—Thermal analysis or thermal optimisation
Abstract
A graphical method for generating a heat dissipation topology by a zero-deficiency grid curved surface utilizes a double-layer grid to design a heat dissipation structure; the upper layer uses a zero-deficiency deformable triangular mesh to describe a topological geometric model, and a heat dissipation topological structure can be generated only through one component, so that the initial layout dependency is eliminated, a clear geometric boundary is obtained, the design variable is reduced, and the calculation speed is improved; the lower layer uses the isogeometric units to perform heat flow coupling steady-state heat conduction calculation and sensitivity analysis, so that the accurate description of the analysis model on the geometric model is realized, errors caused by the fact that the traditional finite element approaches the geometric model by adopting a piecewise function are avoided, the solving precision of the fluid heat dissipation problem is improved, redundant mesh division links in the classical finite element are avoided in the isogeometric analysis process, and the calculation cost is saved; the invention provides a more reliable idea for the topological optimization of the heat dissipation structure.
Description
Technical Field
The invention relates to the technical field of heat dissipation performance optimization design of cold plates, in particular to a graphical method for generating a heat dissipation topology by a zero-deficiency grid curved surface.
Technical Field
High power and high heat flux density are important problems which continuously exist in the development process of modern electronics, a plurality of electronic devices, such as electronic chips, batteries, CPUs and the like, are influenced by heat loads, most of heat is harmful, and the heat must be taken away from a working area in time to ensure the normal work and the service life of the electronic devices and tools; the heat dissipation problem of electronic equipment has been well developed in the past decades, and mainstream cold plate design methods include fluid micro-channel radiators, special-shaped channel designs, multilayer cold plate heat exchangers and the like; however, the traditional cold plate design method depends on the experience of engineers, and the heat dissipation target of complex equipment is difficult to be efficiently and pertinently completed.
Different from the traditional design method, the topological optimization design method is based on calculation, materials are reasonably distributed in a design domain to achieve optimal distribution, the optimal heat dissipation effect is obtained, and a level set method (LST) based on pixel-based variable density process (SIMP), evolutionary algorithm (ESO) and RAMP and high-dimensional space description is produced; in the variable density methods such as SIMP, a design domain is dispersed into pixel points with the same size, black and white pixels and interpolation thereof are used for representing the structure, the number of defects in the design domain can be freely changed, the concept is simple, the use threshold is low, and the variable density method is widely used in a plurality of fields, however, the pixel points of the variable density method are extremely difficult to converge, and serious gray scale and checkerboard phenomena exist, so that the design result lacks clear geometric boundaries, and the geometric shape of a mechanical structure cannot be accurately described; in addition, the design result of the variable density method lacks clear design parameters, the structure characteristic dimension is difficult to control, and great difficulty is caused to the processing and manufacturing of a mechanical structure; the level set method utilizes high-dimensional space description, realizes the variable topology function of a bottom layer structure through the movement and deformation of an upper layer grid, can clearly and accurately describe the geometric boundary of the structure by using a projection method, but has stronger initial layout dependency, the variable topology capability mainly depends on the movement of the upper layer grid and the number of components, and the free change of the deficiency in a design domain is difficult to realize; therefore, the current topology optimization method cannot meet the requirement of heat dissipation topology optimization design in many cases.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention aims to provide a graphical method for generating a heat dissipation topology by a zero-deficiency mesh curved surface.
In order to achieve the aim, the invention adopts the technical scheme that:
a graphical method for generating a heat dissipation topology from a zero-deficiency mesh surface comprises the following steps:
1) defining a design working condition:
a two-dimensional square plane structure is used as a design domain, the size parameter of the design domain is l multiplied by l, the heat load uniformly distributed in the design domain is Q, the boundaries of the design domain are all heat insulation boundaries, a fluid inlet and a fluid outlet are respectively positioned in the centers of the upper boundary and the lower boundary of the design domain, and the width of the fluid inlet and the fluid outlet is l/10;
2) constructing an initial zero-deficiency mesh curved surface:
taking n non-coincident points as vertexes of a geometric model in the zero-deficiency grid to construct a square grid, wherein the coordinate values of the vertexes are (x)i,yi) I-1, 2, …, n, the parameterized geometry further divides the square mesh into triangular meshes, forms an initial component layout described by the mesh surface, sets the thresholds for mesh refinement and splitting, L respectivelySAnd DSSetting the upper limit of the volume fraction beta of the flow passage0Maximum iteration step Loop, variable tolerance tv；
3) Initializing an isogeometric analysis grid:
3.1) constructing an isogeometric analysis grid by using a non-uniform rational B spline (NURBS) curve:
taking m control points to construct a P-th-order B-spline curve, wherein the control vector of the P-th-order B-spline curve is 0Non-decreasing sequence xi between-1 ═ { xi ═ xi1，ξ2，···，ξm+p+1B-spline basis function recursion formula:
using two sets of B-spline basis functions Ni,p(ξ)、Nj,p(xi) and a weight coefficient ωi,jObtaining a two-dimensional NURBS curved surface as a bottom layer isogeometric analysis grid and a basic function expression thereof
3.2) constructing a unit rigidity matrix of the isogeometric units:
when geometric analysis needs to realize mapping from a physical space and a parameter space to a mother space, the unit rigidity for solving the heat flow coupling problem is as follows:
in the formula (I), the compound is shown in the specification,
expressing the i-th and J-th Gaussian integration points in the u and v directions of the e-th geometric unit, J (u, v) is a two-dimensional Jacobian matrix, R(i,l)And R(k,j)Are respectively the u direction
Segment and v direction
The derivative of the segment NURBS curve equation;
3.3) applying loads and constraints:
on the computing level, according to given boundary conditionsApplying Noeman boundary condition, uniformly distributing heat load Q, applying Dirichlet boundary condition on geometric grid control points of fluid inlet and outlet, respectively constant temp. TinConstant pressure Pout(ii) a Given fluid inlet velocity Uin；
4) Calculating the volume of the flow channel:
according to the half-side data structure of the upper zero-deficiency triangular mesh, traversing all triangular units to obtain the geometric boundary of the mesh curved surface, and setting the mesh vertex on the boundary as a control vertex, namely an optimization variable, wherein the calculation method of the triangular mesh volume V is as follows:
wherein ng is the number of triangle units, lngIs the side length of the triangular unit, and t is the thickness of the triangular unit;
plthe calculation method of (2) is as follows:
5) solving the physical field:
5.1) calculating the projection value of the control point:
the topological structure described by the upper triangular mesh is mapped to the lower equal-geometry analysis layer by a specified projection rule, and the projection value of the equal-geometry unit control point is calculated by the following formula:
in the formula, NAcg，NBcgAnd NCcgThe specific calculation method is as follows:
wherein ng is the number of units in the triangular mesh curved surface, e is an euler number, and b is a maximum constant, where b is 100000;
5.2) material property interpolation model:
the heat conductivity coefficient interpolation model of the control points of the equal geometric units is as follows:
in the formula, keIs the thermal conductivity, k, of cell esIs the thermal conductivity, k, of a solid materialwIs the thermal conductivity of the fluid material itself, (x)j e,yj e) J is 1,2,3,4, which is the coordinate of the jth node of element e;
solid, fluid permeability interpolation model:
in the formula, κePermeability of unit e, κsPermeability of solid material, kwPermeability of the fluid material itself;
solid, fluid density interpolation model:
in the formula, ρeIs the density of the cell e, psIs the density of a solid material, pwThe density of the fluid material itself;
solid, fluid specific heat capacity interpolation model:
in the formula, cpeIs the specific heat capacity of unit e, cpsIs the specific heat capacity of a solid material, cpwIs the specific heat capacity of the fluid material itself;
5.3) constructing an integral rigidity matrix and solving a physical field:
the governing equation of the incompressible fluid in the heat flow coupling problem is as follows:
KpP＝fp(12)
(Kt+C(P))T＝ft(13)
in the formula, KpAs a whole penetration matrix, KtIs an integral heat-conducting matrix, P is a pressure field vector, T is a temperature field vector, fpIs a pressure load vector, ftIs the heat flow load vector, C (P) is the overall convection matrix, which is a function of the pressure field P;
integral penetration matrix KpA heat conducting matrix KtAnd the convection matrix C (P) is:
in the formula, NeFor the number of finite element elements contained after the discretization of the region omega,
is a heat conducting matrix of the cell e,
a penetration matrix of cells e, ceA convection matrix for cell e;
6) determining an optimization model:
setting the average temperature of the structure as an objective function by taking the optimal heat dispersion performance of the heat flow structure as an optimization objective; setting the energy dissipation of the heat dissipation structure not to exceed the allowable energy dissipation and the material consumption not to exceed the allowable material consumption, and taking the pressure drop between the fluid inlet and the fluid outlet and the volume ratio of the flow passage as constraint functions;
the mathematical model is optimized as follows:
wherein d ═ d (d)1,d2,d3,…,dng)TFor the coordinate vectors of all control vertices of the design variables, ng is the total number of control vertices involved in the optimization problem, I is the objective function of the optimization problem, Δ P*For a defined pressure drop constraint value, V0To design the total volume of the domain, beta0For designing the upper limit value of the volume fraction, U, of the volume, V, of the fluid cooling channels in the fielddIs the set of all values in d;
7) and (3) sensitivity analysis:
taking any design variable a as a representative, and carrying out sensitivity analysis on a target function and a constraint function of the heat dissipation structure;
7.1) objective function sensitivity:
the sensitivity calculation formula of the objective function is as follows:
in the formula:
7.2) constraint function sensitivity:
the sensitivity of the pressure drop constraint is calculated as follows:
in the formula:
the sensitivity calculation for volume constraints is as follows:
in the formula, ngp is the number of Gaussian points in the e-th equal geometric unit, and ncp is the number of nodes in the e-th equal geometric unit;
8) iterative optimization:
8.1) updating design variables by a moving asymptote optimization algorithm:
substituting the objective function, the constraint function value and the sensitivity value of the objective function and the constraint function value which are obtained by the isogeometric analysis grid calculation into a moving asymptote optimization algorithm (MMA), and iteratively updating the coordinate value of the control vertex;
8.2) updating the curved surface shape of the upper layer grid:
8.2.1) grid operation decision:
carrying out grid operation judgment on the coordinate value after the control vertex is updated; the side length of the triangular mesh satisfies Lmax>LsCarrying out grid subdivision under the condition; the distance between the unshared vertexes of the adjacent triangular meshes satisfies Dij>DsCarrying out grid splitting under the condition; performing self-intersection judgment on the triangular mesh according to a Greiner-Hormann algorithm, and performing subdivision on the triangular mesh by adopting a Constrained DelaunayTriangulation algorithm; obtaining the shape of the upper grid boundary;
8.2.2) update non-control point coordinates:
after the updated coordinates of the control vertexes are obtained, the updated coordinates of other non-control vertexes are calculated by using a stiffness-preserving algorithm, the movement of the non-control vertexes is obtained according to a minimum energy function, and a calculation formula is as follows:
8.3) obtaining an optimal solution by iteration:
circularly iterating the processes from the step 4) to the step 8.2.2) until the iteration step number k is larger than the maximum iteration step number Loop or the difference value of the design variables before and after updating is smaller than the variable tolerance tvSo far, the optimal runner heat dissipation structure meeting the material consumption and energy dissipation constraints is obtained;
9) adaptive processing: and rounding the heat flow structure layout according to the production process requirements so as to obtain the final flow channel distribution.
The invention has the beneficial effects that:
the invention uses zero-deficiency grid curved surface, only configures one initial assembly, realizes continuous topology through grid movement and refinement and control point splitting operation, overcomes the serious dependence of a topology optimization design result on the initial layout of the assembly, and fully releases the deformability of the assembly; the triangular mesh can explicitly describe the topological structure, so that the topological structure has clear and definite geometric boundaries, the defects of checkerboard, fuzzy boundaries and the like in the design result of implicit topological optimization methods such as a variable density method and the like are overcome, the characteristic size of the fluid heat dissipation structure can be effectively controlled, and great convenience is created for processing and manufacturing; according to the method, the bottom layer grid adopts the isogeometric analysis grid, the non-uniform rational B spline curve is adopted to describe the analysis model, the geometric model described by the upper layer grid can be accurately reproduced, errors caused by adopting a piecewise function to approach the geometric model in a classical finite element are avoided, meanwhile, the isogeometric analysis grid uses the division of spline model parameter domains and the mapping of a parent space and a parameter model to a physical model, the complex grid division process in the traditional finite element method is avoided, the time cost is greatly reduced, the isogeometric analysis avoids the frequent interaction process of the geometric model and the analysis model, and the self-adaptive grid refinement is realized; the method can be further expanded, design targets or constraints such as structural thermal resistance, structural corrosion resistance, heat dissipation weakness and the like are increased, the method is assisted to adapt to the actual engineering requirements, the economic benefit is increased, and the method has great significance for the optimal design of the heat flow heat dissipation structure.
Drawings
FIG. 1 is a flow chart of an embodiment of the present invention.
FIG. 2 is a schematic diagram of boundary conditions according to an embodiment of the present invention.
Fig. 3 is a schematic diagram of a projection of a dual-layer grid in an embodiment of the invention.
Fig. 4 is an optimization result of the distribution of the heat dissipation structure channels in the embodiment of the present invention.
Detailed Description
The invention is further illustrated by the following figures and examples.
As shown in fig. 1, a graphical method for generating a heat dissipation topology from a zero-genus mesh surface includes the following steps:
1) defining a design working condition:
a two-dimensional square plane structure is taken as a design domain, the size parameter of the design domain is 0.01m multiplied by 0.01m, the heat load is uniformly distributed in the design domain as Q, the boundaries of the design domain are all heat insulation boundaries, a fluid inlet and a fluid outlet are respectively positioned in the centers of the upper boundary and the lower boundary of the design domain, and the widths of the fluid inlet and the fluid outlet are both 0.001m, as shown in figure 2;
2) constructing an initial zero-deficiency mesh curved surface:
taking n non-coincident points as vertexes of a geometric model in the zero-deficiency grid to construct a square grid, wherein the coordinate values of the vertexes are (x)i,yi) I-1, 2, …, n, the parameterized geometry further divides the square mesh into triangular meshes, forms an initial component layout described by the mesh surface, sets the thresholds for mesh refinement and splitting, L respectivelySAnd DSSetting the upper limit of the volume fraction beta of the flow passage050%, maximum iteration step Loop 1000, variable tolerance tv＝0.01；
3) Initializing an isogeometric analysis grid:
3.1) constructing an isogeometric analysis grid by using a non-uniform rational B spline (NURBS) curve:
m control points are taken to construct a p-th-order B spline curve, and a non-decreasing sequence xi ═ xi { xi } between control vectors of the curve and 0-11，ξ2，···，ξm+p+1B-spline basis function recursion formula:
using two sets of B-spline basis functions Ni,p(ξ)、Nj,p(xi) and a weight coefficient ωi,jObtaining a two-dimensional NURBS curved surface as a bottom layer isogeometric analysis grid and a basic function expression thereof
3.2) constructing a unit rigidity matrix of the isogeometric units:
when geometric analysis needs to realize mapping from a physical space and a parameter space to a mother space, the unit rigidity for solving the heat flow coupling problem is as follows:
in the formula (I), the compound is shown in the specification,
expressing the i-th and J-th Gaussian integration points in the u and v directions of the e-th geometric unit and the like, (J (u, v) is a two-dimensional Jacobian matrix, R(i,l)And R(k,j)Are respectively the u direction
Segment and v direction
The derivative of the segment NURBS curve equation;
3.3) applying loads and constraints:
according to given boundary conditions, applying a Noelman boundary condition on a calculation layer, uniformly distributing heat loads Q, and applying Dirichlet boundary conditions on geometric grid control points of a fluid inlet and a fluid outlet, wherein the geometric grid control points are respectively at constant temperature TinConstant pressure Pout(ii) a Given fluid inlet velocity Uin；
4) Calculating the volume of the flow channel:
according to the half-side data structure of the upper zero-deficiency triangular mesh, traversing all triangular units to obtain the geometric boundary of the mesh curved surface, and setting the mesh vertex on the boundary as a control vertex, namely an optimization variable, wherein the calculation method of the triangular mesh volume V is as follows:
wherein ng is the number of triangle units, lngIs the side length of the triangular unit, and t is the thickness of the triangular unit;
plthe calculation method of (2) is as follows:
5) solving the physical field:
5.1) calculating the projection value of the control point:
as shown in fig. 3, the topological structure described by the upper triangular mesh is mapped to the lower isogeometric analysis layer according to a predetermined projection rule, and the projection value of the isogeometric unit control point is calculated, where the calculation formula is:
in the formula, NAcg，NBcgAnd NCcgThe specific calculation method is as follows:
wherein ng is the number of units in the triangular mesh curved surface, e is an euler number, and b is a maximum constant, where b is 100000;
5.2) material property interpolation model:
the heat conductivity coefficient interpolation model of the control points of the equal geometric units is as follows:
in the formula, keIs the thermal conductivity, k, of cell esIs the thermal conductivity, k, of a solid materialwIs the thermal conductivity of the fluid material itself, (x)j e,yj e) J is 1,2,3,4, which is the jth element of unit eCoordinates of the nodes;
solid, fluid permeability interpolation model:
in the formula, κePermeability of unit e, κsPermeability of solid material, kwPermeability of the fluid material itself;
solid, fluid density interpolation model:
in the formula, ρeIs the density of the cell e, psIs the density of a solid material, pwThe density of the fluid material itself;
solid, fluid specific heat capacity interpolation model:
in the formula, cpeIs the specific heat capacity of unit e, cpsIs the specific heat capacity of a solid material, cpwIs the specific heat capacity of the fluid material itself;
5.3) constructing an integral rigidity matrix and solving a physical field:
the governing equation of the incompressible fluid in the heat flow coupling problem is as follows:
KpP＝fp(12)
(Kt+C(P))T＝ft(13)
in the formula, KpAs a whole penetration matrix, KtIs an integral heat-conducting matrix, P is a pressure field vector, T is a temperature field vector, fpIs a pressure load vector, ftIs the heat flow load vector, C (P) is the overall convection matrix, which is a function of the pressure field P;
integral penetration matrix KpA heat conducting matrix KtAnd the convection matrix C (P) is:
in the formula, NeFor the number of finite element elements contained after the discretization of the region omega,
is a heat conducting matrix of the cell e,
a penetration matrix of cells e, ceA convection matrix for cell e;
6) determining an optimization model:
setting the average temperature of the structure as an objective function by taking the optimal heat dispersion performance of the heat flow structure as an optimization objective; setting the energy dissipation of the heat dissipation structure not to exceed the allowable energy dissipation and the material consumption not to exceed the allowable material consumption, and taking the pressure drop between the fluid inlet and the fluid outlet and the volume ratio of the flow passage as constraint functions;
the mathematical model is optimized as follows:
wherein d ═ d (d)1,d2,d3,…,dng)TFor the coordinate vectors of all control vertices of the design variables, ng is the total number of control vertices involved in the optimization problem, I is the objective function of the optimization problem, Δ P*For a defined pressure drop constraint value, V0To design the total volume of the domain, beta0For designing the upper limit value of the volume fraction, U, of the volume, V, of the fluid cooling channels in the fielddIs the set of all values in d;
7) and (3) sensitivity analysis:
taking any design variable a as a representative, and carrying out sensitivity analysis on a target function and a constraint function of the heat dissipation structure;
7.1) objective function sensitivity:
the sensitivity calculation formula of the objective function is as follows:
in the formula:
7.2) constraint function sensitivity:
the sensitivity of the pressure drop constraint is calculated as follows:
in the formula:
the sensitivity calculation for volume constraints is as follows:
in the formula, ngp is the number of Gaussian points in the e-th equal geometric unit, and ncp is the number of nodes in the e-th equal geometric unit;
8) iterative optimization:
8.1) updating design variables by a moving asymptote optimization algorithm:
substituting the objective function, the constraint function value and the sensitivity value of the objective function and the constraint function value which are obtained by the isogeometric analysis grid calculation into a moving asymptote optimization algorithm (MMA), and iteratively updating the coordinate value of the control vertex;
8.2) updating the curved surface shape of the upper layer grid:
8.2.1) grid operation decision:
carrying out grid operation judgment on the coordinate value after the control vertex is updated; the side length of the triangular unit satisfies Lmax>LsCarrying out grid subdivision under the condition; the distance between the unshared vertexes of the adjacent triangle units satisfies Dij>DsCarrying out grid splitting under the condition; performing self-intersection judgment on the unit according to a Greiner-Hormann algorithm, and performing subdivision on the grid by adopting a constrainedDelaunayTriangulation algorithm; obtaining the shape of the upper grid boundary;
8.2.2) update non-control point coordinates:
after the updated coordinates of the control vertexes are obtained, the updated coordinates of other non-control vertexes are calculated by using a stiffness-preserving algorithm, the movement of the non-control vertexes is obtained according to a minimum energy function, and a calculation formula is as follows:
8.3) obtaining an optimal solution by iteration:
circularly iterating the processes from the step 4) to the step 8.2.2) until the iteration step number k is larger than the maximum iteration step number Loop or the difference value of the design variables before and after updating is smaller than the variable tolerance tvSo far, the optimal runner heat dissipation structure meeting the material consumption and energy dissipation constraints can be obtained;
9) adaptive processing: the heat flow structure layout is rounded according to the production process requirements to obtain the final flow channel distribution, as shown in fig. 4.
Claims (1)
1. A graphical method for generating a heat dissipation topology from a zero-deficiency mesh surface is characterized by comprising the following steps:
1) defining a design working condition:
a two-dimensional square plane structure is used as a design domain, the size parameter of the design domain is l multiplied by l, the heat load uniformly distributed in the design domain is Q, the boundaries of the design domain are all heat insulation boundaries, a fluid inlet and a fluid outlet are respectively positioned in the centers of the upper boundary and the lower boundary of the design domain, and the width of the fluid inlet and the fluid outlet is l/10;
2) constructing an initial zero-deficiency mesh curved surface:
taking n non-coincident points as vertexes of a geometric model in the zero-deficiency grid to construct a square grid, wherein the coordinate values of the vertexes are (x)i,yi) I-1, 2, …, n, the parameterized geometry further divides the square mesh into triangular meshes, forms an initial component layout described by the mesh surface, sets the thresholds for mesh refinement and splitting, L respectivelySAnd DSSetting the upper limit of the volume fraction beta of the flow passage0Maximum iteration step Loop, variable tolerance tv；
3) Initializing an isogeometric analysis grid:
3.1) constructing an isogeometric analysis grid by using a non-uniform rational B spline (NURBS) curve:
m control points are taken to construct a p-th-order B spline curve, and a non-decreasing sequence xi ═ xi { xi } between control vectors of the curve and 0-11，ξ2，···，ξm+p+1B-spline basis function recursion formula:
using two sets of B-spline basis functions Ni,p(ξ)、Nj,p(xi) and a weight coefficient ωi,jObtaining a two-dimensional NURBS curved surface as a bottom layer isogeometric analysis grid and a basic function expression thereof
3.2) constructing a unit rigidity matrix of the isogeometric units:
when geometric analysis needs to realize mapping from a physical space and a parameter space to a mother space, the unit rigidity for solving the heat flow coupling problem is as follows:
in the formula (I), the compound is shown in the specification,
expressing the i-th and J-th Gaussian integration points in the u and v directions of the e-th geometric unit, J (u, v) is a two-dimensional Jacobian matrix, R(i,l)And R(k,j)Are respectively the u direction
Segment and v direction
The derivative of the segment NURBS curve equation;
3.3) applying loads and constraints:
according to given boundary conditions, applying a Noelman boundary condition on a calculation layer, uniformly distributing heat loads Q, and applying Dirichlet boundary conditions on geometric grid control points of a fluid inlet and a fluid outlet, wherein the geometric grid control points are respectively at constant temperature TinConstant pressure Pout(ii) a Given fluid inlet velocity Uin；
4) Calculating the volume of the flow channel:
according to the half-side data structure of the upper zero-deficiency triangular mesh, traversing all triangular units to obtain the geometric boundary of the mesh curved surface, and setting the mesh vertex on the boundary as a control vertex, namely an optimization variable, wherein the calculation method of the triangular mesh volume V is as follows:
wherein ng is the number of triangle units, lngIs the side length of the triangular unit, and t is the thickness of the triangular unit;
plthe calculation method of (2) is as follows:
5) solving the physical field:
5.1) calculating the projection value of the control point:
the topological structure described by the upper triangular mesh is mapped to the lower equal-geometry analysis layer by a specified projection rule, and the projection value of the equal-geometry unit control point is calculated by the following formula:
in the formula, NAcg，NBcgAnd NCcgThe specific calculation method is as follows:
wherein ng is the number of units in the triangular mesh curved surface, e is an euler number, and b is a maximum constant, where b is 100000;
5.2) material property interpolation model:
the heat conductivity coefficient interpolation model of the control points of the equal geometric units is as follows:
in the formula, keIs the thermal conductivity, k, of cell esIs the thermal conductivity, k, of a solid materialwIs the thermal conductivity of the fluid material itself,
is the coordinate of the jth node of element e;
solid, fluid permeability interpolation model:
in the formula, κePermeability of unit e, κsPermeability of solid material, kwBeing fluid material itselfPermeability of (d);
solid, fluid density interpolation model:
in the formula, ρeIs the density of the cell e, psIs the density of a solid material, pwThe density of the fluid material itself;
solid, fluid specific heat capacity interpolation model:
in the formula, cpeIs the specific heat capacity of unit e, cpsIs the specific heat capacity of a solid material, cpwIs the specific heat capacity of the fluid material itself;
5.3) constructing an integral rigidity matrix and solving a physical field:
the governing equation of the incompressible fluid in the heat flow coupling problem is as follows:
KpP＝fp(12)
(Kt+C(P))T＝ft(13)
in the formula, KpAs a whole penetration matrix, KtIs an integral heat-conducting matrix, P is a pressure field vector, T is a temperature field vector, fpIs a pressure load vector, ftIs the heat flow load vector, C (P) is the overall convection matrix, which is a function of the pressure field P;
integral penetration matrix KpA heat conducting matrix KtAnd the convection matrix C (P) is:
in the formula, NeFor the number of finite element elements contained after the discretization of the region omega,
is a heat conducting matrix of the cell e,
a penetration matrix of cells e, ceA convection matrix for cell e;
6) determining an optimization model:
setting the average temperature of the structure as an objective function by taking the optimal heat dispersion performance of the heat flow structure as an optimization objective; setting the energy dissipation of the heat dissipation structure not to exceed the allowable energy dissipation and the material consumption not to exceed the allowable material consumption, and taking the pressure drop between the fluid inlet and the fluid outlet and the volume ratio of the flow passage as constraint functions;
the mathematical model is optimized as follows:
wherein d ═ d (d)1,d2,d3,…,dng)TFor the coordinate vectors of all control vertices of the design variables, ng is the total number of control vertices involved in the optimization problem, I is the objective function of the optimization problem, Δ P*For a defined pressure drop constraint value, V0To design the total volume of the domain, beta0For designing the upper limit value of the volume fraction, U, of the volume, V, of the fluid cooling channels in the fielddIs the set of all values in d;
7) and (3) sensitivity analysis:
taking any design variable a as a representative, and carrying out sensitivity analysis on a target function and a constraint function of the heat dissipation structure;
7.1) objective function sensitivity:
the sensitivity calculation formula of the objective function is as follows:
in the formula:
7.2) constraint function sensitivity:
the sensitivity of the pressure drop constraint is calculated as follows:
in the formula:
the sensitivity calculation for volume constraints is as follows:
in the formula, ngp is the number of Gaussian points in the e-th equal geometric unit, and ncp is the number of nodes in the e-th equal geometric unit;
8) iterative optimization:
8.1) updating design variables by a moving asymptote optimization algorithm:
substituting the objective function, the constraint function value and the sensitivity value of the objective function and the constraint function value which are obtained by the isogeometric analysis grid calculation into a moving asymptote optimization algorithm (MMA), and iteratively updating the coordinate value of the control vertex;
8.2) updating the curved surface shape of the upper layer grid:
8.2.1) grid operation decision:
carrying out grid operation judgment on the coordinate value after the control vertex is updated; the side length of the triangular mesh satisfies Lmax>LsCarrying out grid subdivision under the condition; the distance between the unshared vertexes of the adjacent triangular meshes satisfies Dij>DsUnder the condition of performing grid divisionCracking; performing self-intersection judgment on the triangular mesh according to a Greiner-Hormann algorithm, and performing subdivision on the triangular mesh by adopting a Constrained DelaunayTriangulation algorithm; obtaining the shape of the upper grid boundary;
8.2.2) update non-control point coordinates:
after the updated coordinates of the control vertexes are obtained, the updated coordinates of other non-control vertexes are calculated by using a stiffness-preserving algorithm, the movement of the non-control vertexes is obtained according to a minimum energy function, and a calculation formula is as follows:
8.3) obtaining an optimal solution by iteration:
circularly iterating the processes from the step 4) to the step 8.2.2) until the iteration step number k is larger than the maximum iteration step number Loop or the difference value of the design variables before and after updating is smaller than the variable tolerance tvSo far, the optimal runner heat dissipation structure meeting the material consumption and energy dissipation constraints is obtained;
9) adaptive processing: and rounding the heat flow structure layout according to the production process requirements so as to obtain the final flow channel distribution.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010632026.3A CN111832203B (en) | 2020-07-02 | 2020-07-02 | Graphical method for generating heat dissipation topology by zero-deficiency grid curved surface |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010632026.3A CN111832203B (en) | 2020-07-02 | 2020-07-02 | Graphical method for generating heat dissipation topology by zero-deficiency grid curved surface |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111832203A true CN111832203A (en) | 2020-10-27 |
| CN111832203B CN111832203B (en) | 2022-12-09 |
Family
ID=72900525
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010632026.3A Active CN111832203B (en) | 2020-07-02 | 2020-07-02 | Graphical method for generating heat dissipation topology by zero-deficiency grid curved surface |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111832203B (en) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112308975A (en) * | 2020-11-06 | 2021-02-02 | 中国石油天然气集团有限公司 | RVM binary model analysis method and system for three-dimensional lightweight engine |
| CN112446163A (en) * | 2020-11-24 | 2021-03-05 | 西安交通大学 | Energy finite element topological optimization method based on parameterized level set |
| CN113343462A (en) * | 2021-06-07 | 2021-09-03 | 西安交通大学 | High-order isogeometric multi-oil-cavity dynamic and static pressure sliding bearing oil film characteristic simulation method |
| CN113434921A (en) * | 2021-07-05 | 2021-09-24 | 西安交通大学 | Structure equal-geometry topological optimization method considering mesoscale effect |
| CN113836651A (en) * | 2021-08-31 | 2021-12-24 | 厦门大学 | Turbine blade cascade flow channel topology design method based on fluid topology optimization |
| CN113836831A (en) * | 2021-08-12 | 2021-12-24 | 上海理工大学 | Fluid cooling channel shape optimization method based on isogeometric analysis |
| CN114117877A (en) * | 2021-11-26 | 2022-03-01 | 西安交通大学 | Topological optimization method based on isogeometric particle description |
| CN114491769A (en) * | 2022-02-17 | 2022-05-13 | 河海大学 | Free-form surface structure integrated form creation method based on isogeometric analysis method |
| CN114662375A (en) * | 2022-03-31 | 2022-06-24 | 西安交通大学 | Generation type design method for special-shaped fuel structure of fast neutron reactor core |
| CN115169274A (en) * | 2022-06-20 | 2022-10-11 | 浙江大学 | Method and device for generating geometric adaptive numerical simulation grid of electronic device assembly |
| CN116228993A (en) * | 2023-05-08 | 2023-06-06 | 中国空气动力研究与发展中心计算空气动力研究所 | Grid edge construction method |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106845021A (en) * | 2017-02-28 | 2017-06-13 | 湘潭大学 | Anisotropic material heat structure Topology Optimization Method based on mesh free RKPM |
| CN109670200A (en) * | 2018-11-13 | 2019-04-23 | 华中科技大学 | A kind of equal geometry density of material field structure Topology Optimization Method |
| CN109800507A (en) * | 2019-01-22 | 2019-05-24 | 西安电子科技大学 | A kind of pair of secondary Structural shape optimization of heat dissipation cold plate topology boundary |
-
2020
- 2020-07-02 CN CN202010632026.3A patent/CN111832203B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106845021A (en) * | 2017-02-28 | 2017-06-13 | 湘潭大学 | Anisotropic material heat structure Topology Optimization Method based on mesh free RKPM |
| CN109670200A (en) * | 2018-11-13 | 2019-04-23 | 华中科技大学 | A kind of equal geometry density of material field structure Topology Optimization Method |
| CN109800507A (en) * | 2019-01-22 | 2019-05-24 | 西安电子科技大学 | A kind of pair of secondary Structural shape optimization of heat dissipation cold plate topology boundary |
Non-Patent Citations (2)
| Title |
|---|
| 傅晓锦等: "基于等几何裁剪分析的拓扑与形状集成优化", 《振动与冲击》 * |
| 杜义贤等: "利用无网格法进行几何非线性热固耦合柔性机构拓扑优化设计", 《固体力学学报》 * |
Cited By (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112308975A (en) * | 2020-11-06 | 2021-02-02 | 中国石油天然气集团有限公司 | RVM binary model analysis method and system for three-dimensional lightweight engine |
| CN112446163A (en) * | 2020-11-24 | 2021-03-05 | 西安交通大学 | Energy finite element topological optimization method based on parameterized level set |
| CN112446163B (en) * | 2020-11-24 | 2022-12-09 | 西安交通大学 | Energy finite element topological optimization method based on parameterized level set |
| CN113343462A (en) * | 2021-06-07 | 2021-09-03 | 西安交通大学 | High-order isogeometric multi-oil-cavity dynamic and static pressure sliding bearing oil film characteristic simulation method |
| CN113434921A (en) * | 2021-07-05 | 2021-09-24 | 西安交通大学 | Structure equal-geometry topological optimization method considering mesoscale effect |
| CN113836831A (en) * | 2021-08-12 | 2021-12-24 | 上海理工大学 | Fluid cooling channel shape optimization method based on isogeometric analysis |
| CN113836651B (en) * | 2021-08-31 | 2023-09-19 | 厦门大学 | Turbine cascade runner topology design method based on fluid topology optimization |
| CN113836651A (en) * | 2021-08-31 | 2021-12-24 | 厦门大学 | Turbine blade cascade flow channel topology design method based on fluid topology optimization |
| CN114117877A (en) * | 2021-11-26 | 2022-03-01 | 西安交通大学 | Topological optimization method based on isogeometric particle description |
| CN114117877B (en) * | 2021-11-26 | 2024-04-09 | 西安交通大学 | Topological optimization method based on isogeometric particle description |
| CN114491769A (en) * | 2022-02-17 | 2022-05-13 | 河海大学 | Free-form surface structure integrated form creation method based on isogeometric analysis method |
| CN114662375A (en) * | 2022-03-31 | 2022-06-24 | 西安交通大学 | Generation type design method for special-shaped fuel structure of fast neutron reactor core |
| CN114662375B (en) * | 2022-03-31 | 2023-11-28 | 西安交通大学 | Method for designing generation type special-shaped fuel structure of fast neutron reactor core |
| CN115169274B (en) * | 2022-06-20 | 2023-05-02 | 浙江大学 | Method and device for generating geometric self-adaptive numerical simulation grid of electronic device assembly |
| CN115169274A (en) * | 2022-06-20 | 2022-10-11 | 浙江大学 | Method and device for generating geometric adaptive numerical simulation grid of electronic device assembly |
| CN116228993A (en) * | 2023-05-08 | 2023-06-06 | 中国空气动力研究与发展中心计算空气动力研究所 | Grid edge construction method |
| CN116228993B (en) * | 2023-05-08 | 2023-08-25 | 中国空气动力研究与发展中心计算空气动力研究所 | Grid edge construction method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111832203B (en) | 2022-12-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111832203B (en) | Graphical method for generating heat dissipation topology by zero-deficiency grid curved surface | |
| CN111709171B (en) | Isogeometric solving and heat dissipation topology generation method for heat flow strong coupling problem | |
| De Boer et al. | Mesh deformation based on radial basis function interpolation | |
| CN109800507B (en) | Secondary shape optimization design method for topological boundary of heat dissipation cold plate | |
| CN111709096B (en) | Design method of special-shaped fin structure for strengthening natural convection heat transfer | |
| CN114741753B (en) | Thin-wall reinforcement structure optimization method and device, electronic equipment and storage medium | |
| Zhang et al. | Kriging-based shape optimization framework for blended-wing-body underwater glider with NURBS-based parametrization | |
| CN114077802B (en) | Particle modeling method using shape function interpolation to replace kernel function approximation | |
| CN104809334A (en) | Calculating method of large-volume concrete cooling temperature field | |
| Zhang et al. | A multiple level set method for modeling grain boundary evolution of polycrystalline materials | |
| CN110276131B (en) | Wing body fusion underwater glider appearance optimization method based on polynomial response surface model | |
| Lu et al. | Three-dimensional temperature field inversion calculation based on an artificial intelligence algorithm | |
| Nambu et al. | Adjoint-based Shape Optimization of High-lift Airfoil using the NSU2D Unstructured Mesh Solver | |
| CN111128316B (en) | Thermal performance analysis method for splicing material with straight cracks or heterogeneous materials | |
| CN114117877A (en) | Topological optimization method based on isogeometric particle description | |
| CN114548526A (en) | Satellite component layout temperature field prediction method based on physical prior neural network | |
| CN113268910A (en) | Gravity-driven natural convection special-shaped heat sink structure topology optimization method | |
| Chen et al. | An improved LU-SGS scheme with faster convergence for unstructured grids of arbitrary topology | |
| CN113420392B (en) | Conjugate heat transfer radiator design method based on flow channel track optimization | |
| Sloboda et al. | Numerical approach in aeroelasticity | |
| Lassaline | A Navier-Stokes equation solver using agglomerated multigrid featuring directional coarsening and line-implicit smoothing. | |
| CN110941882A (en) | Thermal performance analysis method of composite material with curved interface | |
| CN114077800B (en) | Gridless particle method capable of realizing non-traversal solving | |
| Soleimani et al. | Effect of Micro-Fin Geometry on Liquid Heat Transfer Rate and Pressure Drop | |
| Wang et al. | Aerodynamic Optimization Design of Airfoil Shape Using Entropy Generation as an Objective |
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 |