A shape optimization method for fluid cooling channels based on isogeometric analysis

Abstract

The present invention proposes a method for optimizing the shape of a fluid cooling channel based on isogeometric analysis, comprising the following steps: step S1: initialization step; step S2: performing isogeometric analysis and calculating coupled temperature field; step S3: performing sensitivity analysis calculation; step S4: substituting the parameters obtained in step S2 and step S3 into the MMA algorithm for iterative optimization, and applying the MMA algorithm to complete the optimization of the shape of the fluid cooling channel based on isogeometric analysis; step S5: judging whether the average temperature change is less than 0.001°C, if the average temperature change is less than 0.001°C, the optimization is completed, if the average temperature change is greater than or equal to 0.001°C, repeating steps S2‑S4. The present application optimizes the shape of the fluid cooling channel by means of isogeometric analysis, overcoming the limitations of the traditional finite element method in structural optimization. In isogeometric analysis, the same spline information, such as control points and spline basis functions, represents geometric figures in the CAD system and can also be used for numerical calculation analysis.

Description

Technical Field

[0001] The invention relates to the field of structural analysis and optimization, and in particular to a fluid cooling channel shape optimization method based on isogeometric analysis. Background Art

[0002] At present, complete sets of electromechanical equipment are increasingly developing towards the trend of large-scale power and small size. The problem of high heat flux density per unit area of ​​electronic equipment has become more serious, making the heat dissipation problem of electronic equipment increasingly prominent. The improvement of the heat dissipation capacity of electromechanical equipment can be achieved on the one hand by preparing heat dissipation substrate materials with high thermal conductivity, and on the other hand by designing efficient heat dissipation structures. Commonly used heat dissipation structures mainly include air-cooled fins, fluid cooling channels and heat conduction channels. The fluid cooling channel uses convection heat exchange to dissipate heat, which is suitable for occasions with high heat dissipation requirements. Under the current trend, the fluid cooling channel has a good development prospect.

[0003] In the field of structural analysis and optimization research, the finite element method (FEM) is favored by many researchers, however, these methods usually have some disadvantages. In shape optimization, the design domain needs to be cumbersomely parameterized and time-consuming meshing is required. In topology optimization, the design results are usually limited to the initial design space, and additional post-processing is required to communicate with the CAD system. These shortcomings are caused by the use of different mathematical languages ​​in design or geometric modeling and numerical analysis. In traditional finite element topology optimization methods, the boundary continuity of the obtained topology optimization structure is poor. Summary of the invention

[0004] The purpose of the present invention is to propose a new fluid cooling channel analysis calculation and shape optimization method based on isogeometric analysis, which can effectively solve the problem of poor continuity of structural boundaries.

[0005] To achieve the above object, the present invention proposes a method for optimizing the shape of a fluid cooling channel based on isogeometric analysis, comprising the following steps:

[0006] Step S1: initialization step;

[0007] Step S2: perform isogeometric analysis and calculate coupled temperature field;

[0008] Step S3: perform sensitivity analysis calculation;

[0009] Step S4: Substituting the parameters obtained in step S2 and step S3 into the MMA algorithm for iterative optimization, and applying the MMA algorithm to complete the shape optimization of the fluid cooling channel based on isogeometric analysis;

[0010] Step S5: Determine whether the average temperature change is less than 0.001°C. If the average temperature change is less than 0.001°C, the optimization is completed. If the average temperature change is greater than or equal to 0.001°C, repeat steps S2-S4.

[0011] Furthermore, in step S1, the initialization step includes the volume of the design domain, calculation of the maximum material consumption rate, material properties, setting of the initial channel configuration, and meshing to generate boundary conditions.

[0012] Furthermore, in step S2, the following steps are also included:

[0013] Step D1: simplify the NS equation to form a Darcy reduction model, discretize the Darcy reduction model, and convert the formula into a matrix form;

[0014] Step D2: perform material property interpolation to determine the material properties of each unit;

[0015] Step D3: assembling the overall permeability matrix to obtain the pressure field P;

[0016] Step D4: Assemble the overall convection matrix and obtain the temperature field T by coupling the matrix equation;

[0017] Step D5: derive isogeometric analysis formulas for the Darcy reduced-order model to achieve isogeometric analysis of the fluid cooling channel.

[0018] Furthermore, in step S3, the following steps are also included:

[0019] Step L1: Calculate the sensitivity of the objective function;

[0020] Step L2: Calculate the sensitivity of the pressure reduction constraint and the volume constraint;

[0021] Furthermore, in step L1, the objective function is

[0022] Where: H is the number of control points, ones(1,H) is an H-dimensional row vector whose elements are all 1.

[0023] Compared with the prior art, the advantages of the present invention are: it overcomes the limitations of the traditional finite element method in structural optimization. In isogeometric analysis, the same spline information, such as control points and spline basis functions, which represent geometric figures in the CAD system, can also be used for numerical calculation and analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 This is a method for optimizing the shape of a fluid cooling channel based on isogeometric analysis. The general forced convection heat transfer problem is described as follows:

[0025] Figure 2 is the mapping relationship between the physical domain, the parameter domain and the parent domain in step S2 of the present invention;

[0026] Figure 3 It is the pipeline model 1 in step S3 of the present invention;

[0027] Figure 4 It is the pipeline model 2 in step S3 of the present invention;

[0028] Figure 5 The first boundary after fitting of the pipeline model in step S3 of the present invention;

[0029] Figure 6 The boundary after fitting of the second pipeline model in step S3 of the present invention;

[0030] Figure 7 It is a slice diagram after fitting of the pipeline model 1 in step S3 of the present invention;

[0031] Figure 8 It is a slice diagram after fitting of the pipeline model 1 in step S3 of the present invention;

[0032] Fig. 9 It is a schematic flow chart of a method for optimizing the shape of a fluid cooling channel based on isogeometric analysis in the present invention;

[0033] Fig.10 A comparison diagram of the pipeline model before and after optimization in step S4 of the present invention;

[0034] Fig.11 This is a diagram of the pipeline model optimization iteration process in step S4 of the present invention;

[0035] Fig.12 This is a comparison diagram before and after optimization of pipeline model 2 in step S4 of the present invention;

[0036] Fig.13 This is a diagram of the iterative process of optimizing pipeline model 2 in step S4 of the present invention. DETAILED DESCRIPTION

[0037] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be further described below.

[0038] The present invention is described in more detail by way of example in the following paragraphs with reference to the accompanying drawings. Advantages and features of the present invention will become more apparent from the following description.

[0039] like Figure 1 This is the description of the forced convection heat transfer problem. The region Ω includes the solid region Ω s and fluid area Ω w , and in the solid region Ω s There is a heat source Ω h,Firstly, the initialization steps are performed, including the volume of the design domain, the maximum material consumption rate, ,material properties, setting the initial channel configuration, meshing, and generating ,boundary conditions.

[0040] Among them, the initial channel configuration is to simulate the pipeline fitting. There are multiple methods and approaches for pipeline fitting. The selected method has no direct impact on the test results. In this scheme, the topology optimization results based on the traditional finite element method ( Figure 3 , Figure 4 ) is used to fit a fluid cooling channel model suitable for isogeometric analysis. Since the original model is complex and has distorted boundaries, it is not conducive to analysis. Therefore, its boundaries must first be smoothed. The classic Douglas-Peucker algorithm is used in the present invention.

[0041] After smoothing, we will segment the boundary. Then we apply the standard linear least squares fitting technique to fit the boundary. The fitted boundary model is as follows: Figure 5 , Figure 6 As shown. The Coons interpolation method is selected for patch fitting this time. When parameterizing the internal control points of the surface, the surface control points are required to be input sequentially. Therefore, the parameterization direction of each patch needs to be specified. Calculate the average normal vector of the boundary line, and sort the boundary lines according to the average normal vector to make them meet the Coons interpolation requirements. After fitting the pipeline model 1, there are 66 NURBS patches in total. The degree of each patch in the U and V directions is 2, and the node vectors in these two directions are {0, 0, 0, 0.5, 1, 1, 1}. After fitting the pipeline model 2, there are 65 NURBS patches in total. The degree of each patch in the U and V directions is 2, and the node vectors in these two directions are {0, 0, 0, 0.5, 1, 1, 1}. The fitted fragments and control point diagrams are shown as follows. Figure 7 , Figure 8 shown.

[0042] In this technical solution, the energy control equation followed by the forced convection heat transfer problem can be expressed in the form of the convection-diffusion equation:

[0043]

[0044] Where: ρ is the density of the fluid; c p is the specific heat capacity of the fluid working medium; k is the thermal conductivity of the fluid working medium; u is the velocity field of the fluid working medium; T is the temperature field of the area Ω, and Q is the volume heat source intensity of the area Ω.

[0045] It can be seen from the convection-diffusion equation that calculating the fluid velocity field is the prerequisite for solving the coupled temperature field. In the analysis and calculation of the invention, the working fluid is followed as an incompressible liquid, and the NS equation and continuity equation followed in its flow are shown in equations (2) and (3):

[0046]

[0047]

[0048] The present invention approximates the flow of fluid in the forced convection heat transfer problem as an inviscid flow in a porous medium, thereby applying Darcy's law to simplify the NS equations to form a Darcy reduced-order model method. The expression of Darcy's law is as follows:

[0049]

[0050] Where: u is the velocity vector of the fluid; κ is the permeability of the porous medium; μ is the dynamic viscosity coefficient of the fluid; P is the pressure field in the fluid.

[0051] Substituting Darcy's law expression into equation (3), the flow control equation can be transformed into:

[0052]

[0053] Also based on Darcy's law, the convection-diffusion equation in equation (1) can be transformed into:

[0054]

[0055] Applying the finite element method, the matrix form obtained by the above formula after discretizing the entire region Ω is:

[0056] K p P = f p (7)

[0057] (K t +C(P))T=f t (8) Among them, K p is the permeability matrix of the entire area; K t is the thermal conductivity matrix of the entire area; P is the pressure field of the object; T is the temperature field of the object; f p is the external pressure load; f t is the external heat flux load; C(P) is the convection matrix of the entire area, which is a function of P. The pressure field and temperature field should be solved in order, that is, only after the pressure field P is calculated, can the subsequent coupled temperature field T be solved.

[0058] Step 2: Derivation of isogeometric analysis formulas for Darcy reduced-order models.

[0059] In order to calculate the above matrix by isogeometric analysis method, three domains are defined, which are the physical domain Ω and its subdomain Ω e , parameter domain and its subdomains and the mother domain The relevant mapping relationships are as follows: It is a mapping from the parent domain to the parameter subdomain. S: is mapped from the parameter subdomain to the physical subdomain. Used for mapping from parent domain to physical subdomain. The mapping relationship between the three domains is as follows: Figure 2 As shown:

[0060] Mapping S: The Jacobian transformation matrix is:

[0061]

[0062] The Jacobian transformation matrix is:

[0063]

[0064] X e : The Jacobian matrix of the mapping is:

[0065] J=J 1 *J 2 (11)

[0066] The permeation matrix K of the model p , thermal conductivity matrix K t And the integral form of the convection matrix C(P) in the entire physical space is as follows:

[0067]

[0068]

[0069]

[0070] Where N e is the number of cells or subdomains. B e The matrix is ​​defined as follows:

[0071]

[0072] In order to calculate the correlation matrix, the thermal conductivity, permeability, specific heat capacity and other parameters of each unit need to be determined. First, it is necessary to determine whether the area where each control point of the unit is located belongs to solid or liquid. The original geometric analysis control point only contains the information of (x, y, z, ω). Now the material information is added to the control point, that is, (x, y, z, ω, m). When the control point is in the solid area or the solid-liquid boundary, m = 0, and when it is in the liquid area, m = 1.

[0073] The present invention adopts an isogeometric control grid of order 2, and its thermal conductivity properties can be obtained by the following interpolation model:

[0074]

[0075] Where: k e is the thermal conductivity of element e, k s is the thermal conductivity of the solid material, k w is the thermal conductivity of the fluid material, are the coordinates of the jth node of element e.

[0076] The permeability properties can be obtained by the following interpolation model:

[0077]

[0078] Where: κ e is the permeability of unit e, κ s is the permeability of the solid material, κ w is the fluid material permeability.

[0079] The density attribute can be obtained by the following interpolation model:

[0080]

[0081] Where: e is the density of unit e, ρ s is the density of the solid material, ρ w is the density of the fluid material.

[0082] The specific heat capacity property can be obtained by the following interpolation model:

[0083]

[0084] Where: c pe is the specific heat capacity of unit e, c ps is the specific heat capacity of the solid material, c pw is the specific heat capacity of the fluid material.

[0085] In convective heat transfer, the heat exchanged by convection accounts for the majority. Therefore, when calculating the temperature field, the shape function of the unit is corrected by the following simplified upwind stability term, as shown in the following formula:

[0086]

[0087] Where: u e is the average velocity of unit e, h e is the size of the unit.

[0088] The overall pressure load vector can be obtained as follows:

[0089] f p =∫ Γ R T (n·u in )dΓ (21)

[0090] The overall heat flow load vector can be obtained as follows:

[0091]

[0092] Thus, the derivation of the isogeometric analysis formula of the Darcy reduced-order model and the calculation of the coupled temperature field were completed, realizing the isogeometric analysis of the fluid cooling channel.

[0093] Next, the sensitivity analysis calculation step is performed. The present invention selects minimization of average temperature as the design goal, and the objective function is as follows:

[0094]

[0095] Where: H is the number of control points, ones(1,H) is an H-dimensional row vector whose elements are all 1.

[0096] The volume of the fluid cooling channel in the design domain is calculated using the following expression:

[0097] V D =∫ D J 1 J 2 dξdη (24)

[0098] The volume of the fluid cooling channel in the design domain needs to set an upper limit. The volume limit of the fluid cooling channel in the optimization is written in the form of a ratio, which can be expressed as:

[0099]

[0100] Where: β 0 is the volume constraint ratio; V is the total volume of the design domain.

[0101] In engineering applications of fluid flow, the pressure drop between the inlet and outlet represents the energy required to maintain flow. The pressure drop between the inlet and outlet of the fluid cooling channel is calculated as follows:

[0102]

[0103] Where: H 2 is the inlet boundary of the fluid cooling channel Γ 2 The number of control points contained in H 1 is the outlet boundary of the fluid cooling channel Γ 1In order to ensure good fluidity in the fluid cooling channel during optimization, a constraint is added to the pressure drop ΔP between the inlet and outlet of the fluid cooling channel:

[0104]

[0105] Where: ΔP * Represents a given pressure drop constraint.

[0106] At this point, the mathematical model for the shape optimization of the fluid cooling channel in this paper can be obtained:

[0107]

[0108] The design variables for shape optimization of the present invention are selected from the control point coordinate components, and the adjoint method is used to perform sensitivity analysis on the objective function. The following function is introduced according to the Lagrangian number multiplication method:

[0109]

[0110] The above formula is derived for any design variable α,

[0111]

[0112] Eliminating the terms containing dT / dα and dP / dα in the above formula, we can obtain:

[0113]

[0114] In the formula, In formula (31),

[0115]

[0116] Where: N e is the number of units. From equations (31, 32), equation (30) can be simplified to:

[0117]

[0118] in:

[0119]

[0120] In the above formula, we need to deduce and In order to calculate these two terms, two matrices are first defined. The matrix expressions are as follows:

[0121]

[0122]

[0123] Take any unit as the research object and set the basis function vector of the unit to be R = (R 1 ,R 2 ,...,R L ), then the interpolation form of the unit geometric field is:

[0124]

[0125] in

[0126]

[0127] Taking the derivative of the left side of equation (37) with respect to (ξ,η), we can obtain:

[0128]

[0129] Taking the derivative of the right side of equation (37) with respect to (ξ,η), we can obtain:

[0130]

[0131] From formula (37, 39, 40), we can get:

[0132]

[0133] Γ, M and J 1 There is the following relationship between them:

[0134]

[0135] The derivatives of both sides of equations (41, 42) with respect to the design variables are:

[0136]

[0137]

[0138] Substituting formula (44) into formula (43), we get:

[0139]

[0140]

[0141] Through formula (45), we can get After that, we can get Next we need to calculate:

[0142]

[0143] In this way, the sensitivity information of the objective function can be obtained.

[0144] Similar to the calculation of the sensitivity of the objective function to the design variables, the following function is constructed to solve the sensitivity of the pressure drop constraint to the design variables:

[0145]

[0146] Then, the sensitivity of the pressure drop constraint to any design variable α is as follows:

[0147]

[0148] in,

[0149]

[0150] Where:

[0151]

[0152] The sensitivity of the volume constraint to the design variables is:

[0153]

[0154] According to the above mathematical optimization model and sensitivity analysis, the parameters obtained in the above steps are substituted into the MMA algorithm for iterative optimization, and the MMA algorithm is applied to complete the shape optimization of the fluid cooling channel based on isogeometric analysis.

[0155] By judging the change of the average temperature, observe whether the change is less than 0.001°C. If the average temperature change is less than 0.001°C, the optimization is completed. If the average temperature change is greater than or equal to 0.001°C, repeat the above isogeometric analysis and calculation of the coupled temperature field, sensitivity analysis calculation, substitution of parameters into the MMA algorithm for iterative optimization, and application of the MMA algorithm to complete the optimization steps until the average temperature change is less than 0.001°C.

[0156] The average temperature is used as the objective function for global optimization. The relevant constraints and parameters are shown in Table 1. The pipeline model before and after optimization is shown in Fig.10 shown.

[0157] Table 1 Fluid and solid material properties

[0158]

[0159] Fig.11 The figure shows the optimization iteration diagram of the fluid cooling channel. The average temperature before optimization (before fitting) is 1.563°C, and the average temperature after optimization is 1.287°C. The optimization efficiency is 18%.

[0160] The average temperature is used as the objective function for global optimization. The relevant constraints and parameters are shown in Table 1, except that the pressure drop constraint is changed to 120 Pa. The models of pipeline model 2 before and after optimization are as follows Fig.12 shown. Fig.13 The figure shows the optimization iteration diagram of the fluid cooling channel. The average temperature before optimization (before fitting) is 2.13°C, and the average temperature after optimization is 1.825°C. The optimization efficiency is 14%.

[0161] The above is only a preferred embodiment of the present invention and does not limit the present invention in any way. Any technician in the relevant technical field, without departing from the scope of the technical solution of the present invention, makes any form of equivalent replacement or modification to the technical solution and technical content disclosed in the present invention, which does not depart from the content of the technical solution of the present invention and still falls within the protection scope of the present invention.

Claims

1. A method for optimizing the shape of a fluid cooling channel based on isogeometric analysis, characterized in that: The steps include: Step S1: initialization step; Step S2: perform isogeometric analysis and calculate coupled temperature field; Step S3: perform sensitivity analysis calculation; Step S4: Substituting the parameters obtained in step S2 and step S3 into the MMA algorithm for iterative optimization, and applying the MMA algorithm to complete the shape optimization of the fluid cooling channel based on isogeometric analysis; Step S5: Determine whether the average temperature change is less than 0.001°C. If the average temperature change is less than 0.001°C, the optimization is completed. If the average temperature change is greater than or equal to 0.001°C, repeat steps S2-S4. In step S3, the following steps are also included: Step L1: Calculate the sensitivity of the objective function; Step L2: Calculate the sensitivity of the pressure reduction constraint and the volume constraint; In step L1, the objective function is Where: T1 is the object temperature field; H is the number of control points, and ones(1,H) is an H-dimensional row vector whose elements are all 1; In the formula, K p is the permeability matrix of the entire area; K t is the thermal conductivity matrix of the overall area; P is the pressure field of the object; in: In the above formula, κ is the permeability of the porous medium; μ is the dynamic viscosity coefficient of the fluid; it is necessary to derive and In order to calculate these two terms, two matrices are first defined. The matrix expressions are as follows: Take any unit as the research object and set the basis function vector of the unit to be R = (R1, R2, ..., R L ), then the interpolation form of the unit geometric field is: Next we need to calculate: In this way, the sensitivity information of the objective function can be obtained.

Images