A Phase Change Heat Sink Design Method Based on Topology Optimization

Abstract

The present invention discloses a phase change heat sink design method based on topological optimization. The phase change heat sink includes phase change material, fins, shell, and substrate. The design method includes the following steps: determine the heat sink shell size, material, and phase change material type according to the design working conditions; establish a finite element analysis model of the phase change heat sink heat transfer, set optimization goals and constraints, and establish a topological optimization model; use Python language to solve the topological optimization model to preliminarily obtain the size, structure, and layout of the fins; use a binary processing method to reconstruct the topologically optimized fins, and check the structural strength and heat transfer characteristics of the heat sink. The method is used to optimize the design of the phase change heat sink, and the maximum heat transfer performance can be achieved under certain materials. The method provides theoretical guidance for the structural design of the phase change heat sink.

Description

Technical Field

[0001] The present invention belongs to the field of semiconductor technology, and in particular relates to a phase change heat sink design method based on topology optimization. Background Art

[0002] With the advancement of semiconductor technology, the number of transistors integrated on a chip per unit area is increasing, and the heat generation power of electronic device chips is also increasing. In addition, the heat generation power of electronic components does not remain stable during operation, and usually encounters transient loads, causing a sudden change in chip power density. This will cause uneven temperature distribution on the chip surface, which in turn generates huge thermal stress, affecting the reliability of the chip and even causing permanent damage. According to statistics, more than 50% of integrated circuit damage is caused by heat dissipation. For every 10°C increase in the node temperature of an electronic chip, its reliability will be reduced by half. It is particularly important to use appropriate cooling methods to keep electronic chips running at a safe temperature.

[0003] Phase change material (PCM) has good application prospects in non-steady-state thermal management due to its characteristics of large phase change latent heat and small temperature fluctuation. However, PCM usually has low thermal conductivity and poor heat transfer efficiency, and high thermal conductivity material (HCM) fins are needed to enhance heat transfer. The optimization of fins in the study mainly focuses on the structure, shape, size, position and number of fins, which belongs to the category of shape optimization and size optimization. The optimized results are heavily dependent on the experience of researchers and lack effective theoretical guidance. Summary of the invention

[0004] The purpose of the present invention is to optimize the design of the phase change heat sink by using the topology optimization method, so as to reduce the weight of the high thermal conductivity fins as much as possible and increase the filling rate of the phase change material while meeting the heat dissipation performance. The specific scheme is as follows:

[0005] The phase change heat sink design method based on topology optimization is implemented through the following steps:

[0006] Step 1: Determine the shell size, material and phase change material type of the phase change heat sink to be optimized;

[0007] Step 2: Establish a finite element analysis model of heat transfer of phase change heat sink;

[0008] Step 3: Set the physical geometric dimensions, initial conditions and corresponding boundary conditions of the computational domain;

[0009] Step 4: Set the objective function of the topology optimization model. Global integral variables or local feature point variables can be used as the expression form of the optimization objective;

[0010] Step 5: Establish a mapping system from design variables to physical models. Use a (0, 1) continuous distribution function to represent the presence or absence of high thermal conductivity materials in the computational domain, and set a corresponding penalty factor to complete the mapping of design variables to physical models, where the independent variable in the interpolation function is used as the design variable;

[0011] Step 6: Set non-equality constraints. Constrain the volume percentage of high thermal conductivity materials, temperature variation range, etc.;

[0012] Step 7: Iteratively solve the topology optimization model to preliminarily obtain the size, structure and layout position of the fins; wherein, the design variables are filtered during the solution process and projected after the solution is completed, and the iterative solution method is used as the solution algorithm;

[0013] Step 8: Reconstruct the topology optimization model. Determine the final topology optimization configuration shape according to the evolution law of the geometric configuration during the topology optimization process, extract the key geometric points of the geometric configuration, establish a parametric model, complete the reconstruction of the topology optimization configuration, and test and verify the structural strength and heat transfer characteristics of the heat sink.

[0014] Among them, in step 2, the sigmoid function is used to smooth the function H(T) of enthalpy changing with temperature, and T is further differentiated to obtain the smooth equivalent heat capacity function c(T). The θ difference method is used to discretize the time domain of the energy conservation equation in the phase change heat storage process, and the four-node quadrilateral unit and the eight-node quadrilateral unit are used to discretize the spatial domain of the equation to construct the thermal stiffness matrix, temperature change matrix, heat source vector and latent heat vector of the heat transfer equation.

[0015] In step three, optimization objectives such as minimum entropy production, minimum fire accumulation dissipation, minimum average temperature, and minimum local high temperature can be used as objective functions.

[0016] In step five, the SIMP interpolation model of the variable density method is used to interpolate the thermal conductivity, density, and specific heat capacity of the material, and a mapping relationship from the design variables to the physical property parameters is established.

[0017] In step seven, the global moving asymptote algorithm is used to iterate the design variables, and the Helmholtz partial differential equation is used to filter the design variables obtained after each iteration to determine the influence of different filtering radii on the convergence of heat transfer and flow equations, so as to correct the optimization model; the grayscale filtering function is used to project the final topology optimization results, and the filtering function threshold value suitable for the final topological configuration is obtained through numerical calculation, making the structure clearer and highly feasible.

[0018] The Helmholtz partial differential equation has the form:

[0019]

[0020] Where r is the filter radius, ρ f is the filtered design variable, and ρ is the original design variable.

[0021] The projection function form is:

[0022]

[0023] Where: fina is the projected design variable, β is the steepness value, and θ is the threshold value.

[0024] In step eight, a binarization method is used to convert the topology optimization configuration into a model without grayscale units by adjusting the threshold size; the binarized matrix is ​​traversed and the boundary density matrix is ​​extracted, and then the matrix element positions are converted into rectangular coordinates for storage, and the topology optimization structure is converted into a parameterized model; the contour boundary points of the topology optimization result are tracked and sorted using chain codes, and all contour boundary points are traversed in a clockwise or counterclockwise order to obtain the chain code geometry of the contour boundary information, and the coordinate information of the boundary points is saved in order; the differential chain code is used instead of the discrete curvature of the contour to express the geometric information of the contour, and the contour chain code value is weighted averaged; by analyzing the differential chain code graph, the topological contour segmentation points are determined according to the contour curvature mutation points; the B-spline interpolation method is used to complete the reconstruction of the contour boundary.

[0025] The beneficial effects of the present invention are as follows: (1) the present invention takes into account the latent heat change during the phase change process, adopts the equivalent heat capacity method to establish the relationship between the phase change latent heat and the temperature, and smoothes the enthalpy function H(T) through the Sigmoid function, thereby ensuring the convergence and stability of the topological optimization of the phase change problem;

[0026] (2) The present invention establishes a binary processing method for extracting the contour of the topological optimization configuration, converts the topological optimization configuration into a parameterized model, and reconstructs the contour using a B-spline interpolation method, thereby reducing the processing cost while ensuring the heat transfer advantage of the topological optimization configuration;

[0027] (3) The layout, quantity and shape of the high thermal conductivity fins can be preliminarily determined in the design stage based on the heat dissipation intensity, size and heat dissipation boundary conditions. While ensuring the heat dissipation performance, the weight of the high thermal conductivity material can be reduced and the filling rate of the phase change material can be increased. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 Design flow chart for topological optimization of phase change heat sink configuration;

[0029] Figure 2 Schematic diagram of the phase change heat sink structure;

[0030] Figure 3 The topological rib phase change heat sink model when the heat flux is 1.5W / cm2 and the HCM volume accounts for 15%;

[0031] Figure 4 To reconstruct the rib phase change heat sink model using parametric modeling;

[0032] Figure 5 This is a trend diagram of the objective function changing with the iteration process. DETAILED DESCRIPTION

[0033] To make the purpose, technical solution and advantages of the embodiments of the present invention more clear, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0034] Example

[0035] In this embodiment, Figure 2 The phase change heat sink structure is a model, which includes phase change material 1, fins 2 and shell 3. Phase change material 1 and fins 2 are located inside shell 3. The size of shell 3 matches the size of substrate 4. Phase change material 1 is placed on the surface of substrate 4. The heat emitted by the heating element is transferred to the interior of the phase change material 1 heat sink through heat conduction and converted into latent heat of PCM. Substrate 4 is usually an electronic chip. Figure 2 The phase change heat sink structure including the base 4 has an overall model length L of 40 mm, a width W of 40 mm, and a height H of 10 mm.

[0036] The maximum temperature threshold allowed by the electronic chip is set to 80°C. Considering the influence of the inherent thermal resistance between the phase change material 1 and the chip, the phase change temperature of the phase change material 1 should be lower than 80°C. Usually, the phase change material 1 is an organic phase change material with a melting point range of 30 to 90°C. Myristic acid is selected as the phase change material 1, and its phase change temperature is between 55 and 60°C, and the phase change enthalpy is about 200 kJ / kg.

[0037] The fin 2 is made of high thermal conductivity material (HCM). In order to enhance the heat transfer of the phase change material 1, the fin 2 is selected to be 6061 aluminum alloy.

[0038] Correspondingly, the size of the shell 3 should match the size of the phase change material 1, and the material of the shell 3 is a high thermal conductivity metal material such as aluminum alloy or copper.

[0039] In order to accommodate the volume change of the phase change material 1 during the phase change process, a certain volume margin should be left in the shell 3 , and the reserved space volume is 5-10% of the total volume of the shell 3 .

[0040] The fins 2 in this embodiment are optimized by using a topology optimization method, and the phase change heat sink structure model can be obtained by wire cutting or 3D printing.

[0041] Figure 2 The phase change heat sink structure model is topologically optimized through the following steps. See the flowchart for Figure 1 :

[0042] Step 1: Establish a three-dimensional unsteady-state heat transfer model for the phase change process. The expression is as follows:

[0043]

[0044]

[0045]

[0046] Where ρ(γ) is the distribution function of material density; c(γ) is the distribution function of material specific heat capacity; k(γ) is the distribution function of material thermal conductivity; γ is the design variable, representing material distribution; H L is the latent heat of the phase change material, kJ / kg.

[0047] Step 2: Set the objective function. Take the system minimum fire dissipation as the objective function, and the expression is as follows:

[0048]

[0049] Where: k is the thermal conductivity of the material, W / (m·K); is the temperature gradient, ℃; A is the design area, m 2 .

[0050] Step 3: Use the variable density method to complete the mapping of design variables to the physical model. Use the variable density method to interpolate the node density, heat capacity and thermal conductivity coefficient. The calculation formula is as follows:

[0051] ρ(γ)=ρ PCM +(ρ HCM -ρ PCM )ψ(γ) p

[0052] c(γ)=c PCM +(c HCM -c PCM )ψ(γ) p

[0053] k(γ)=k PCM +(k HCM -k PCM )ψ(γ) p

[0054] Where: ρ(γ) is the density at the node; c(γ) is the density specific heat capacity at the node; k(γ) is the thermal conductivity coefficient at the node. The subscript PCM represents phase change material, and the subscript HCM represents high thermal conductivity material. p is the penalty factor, which is set to 3.

[0055] When ρ(γ) approaches 0, the material is of low thermal conductivity, and when ρ(γ) approaches 1, the material is of high thermal conductivity.

[0056] Step 4: Set the volume percentage of the high thermal conductivity fins. For example, γ can be 6%, 9%, 12%, 15%, or 18%.

[0057] Step 5: Use the global moving asymptote algorithm to iteratively calculate the system topology optimization model. Figure 5 The objective function changes with the iteration process, and the Helmholtz partial differential equation is used to filter the design variables obtained after each iteration to determine the influence of different filter radii on the convergence of heat transfer and flow equations, so as to correct the optimization model; the grayscale filter function is used to project the final topology optimization result, and the filter function threshold value suitable for the final topology configuration is obtained through numerical calculation, making the structure clearer and highly feasible, such as Figure 3 shown.

[0058] The Helmholtz partial differential equation has the form:

[0059]

[0060] Where r is the filter radius, ρ f is the filtered design variable, and ρ is the original design variable.

[0061] The projection function form is:

[0062]

[0063] Where: fina is the projected design variable, β is the steepness value, and θ is the threshold value.

[0064] Step 6: Use the binarization method to convert the topology optimization configuration into a model without grayscale units by adjusting the threshold value; traverse the binarized matrix and extract the boundary density matrix, then convert the matrix element position into rectangular coordinates and save it, and convert the topology optimization structure into a parameterized model; use the chain code to track and sort the contour boundary points of the topology optimization result, traverse all the contour boundary points in a clockwise or counterclockwise order to obtain the chain code geometry of the contour boundary information, and save the coordinate information of the boundary points in order; use the differential chain code instead of the contour discrete curvature to express the geometric information of the contour, and perform weighted averaging on the contour chain code value; by analyzing the differential chain code graph, determine the topological contour segmentation points according to the contour curvature mutation point; use the B-spline interpolation method to complete the reconstruction of the contour boundary, such as Figure 4 shown.

[0065] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them.

Claims

1. A phase change heat sink design method based on topology optimization, the phase change heat sink comprising a phase change material, a fin, and a shell, characterized in that: The method comprises the following steps: 1) Determine the heat sink shell size, material, and phase change material type according to the design conditions; 2) Establish a finite element analysis model for heat transfer of non-steady-state phase change heat sink, set optimization objectives and constraints, and establish a topology optimization model; 3) Solve the topology optimization model to preliminarily obtain the size, structure and layout of the fins; 4) Establish a parametric model of the topologically optimized fin configuration, reconstruct the topologically optimized configuration, and obtain the phase change heat sink fin structure; The expression of the finite element analysis model established in step 2) is: Where ρ(γ) is the distribution function of material density; c(γ) is the distribution function of material specific heat capacity; k(γ) is the distribution function of material thermal conductivity; γ is the design variable, representing material distribution; H L is the latent heat of the phase change material, kJ / kg.

2. A phase change heat sink design method based on topology optimization according to claim 1, characterized in that: In step 2), the objective function adopts at least one of minimum entropy production, minimum fire dissipation, minimum average temperature, and minimum local high temperature.

3. The phase change heat sink design method based on topology optimization according to claim 1 is characterized in that: In step 2), the minimum fire dissipation of the system is taken as the objective function, and the expression is as follows: Where: k is the thermal conductivity of the material, W / (m·K); is the temperature gradient, ℃; A is the design area, m 2 .

4. The phase change heat sink design method based on topology optimization according to claim 1 is characterized in that: In step 3), the variable density method is used to interpolate the node density, heat capacity and thermal conductivity. The calculation formula is as follows: p(γ)=p PCM +(r HCM -r PCM )ψ(c) p c(γ)=c PCM +(c HCM -c PCM )ψ(c) p k(γ)=k PCM +(k HCM -k PCM )ψ(c) p Where: ρ(γ) is the density at the node; c(γ) is the density specific heat capacity at the node; k(γ) is the thermal conductivity coefficient at the node, the subscript PCM represents phase change material, the subscript HCM represents high thermal conductivity material, and p is the penalty factor, which is taken as 3; When ρ(γ) approaches 0, the node is made of low thermal conductivity material, and when ρ(γ) approaches 1, the node is made of high thermal conductivity material.

5. The phase change heat sink design method based on topology optimization according to claim 1 is characterized in that: After step 3) and before step 4), non-equality constraints on high thermal conductivity materials are set.

6. The phase change heat sink design method based on topology optimization according to claim 1, characterized in that: In step 4), the method for establishing a parameterized model of the topologically optimized fin configuration is: The global moving asymptote algorithm is used to iteratively calculate the system topology optimization model, and the Helmholtz partial differential equation is used to filter the design variables obtained after each iteration to determine the influence of different filtering radii on the convergence of heat transfer and flow equations, so as to correct the optimization model; the grayscale filtering function is used to project the final topology optimization results, and the threshold value of the filtering function suitable for the final topological configuration is obtained through numerical calculation; The Helmholtz partial differential equation has the form: ▽·(-r▽ρ f )+▽ρ f =ρ Where r is the filter radius, ρ f is the filtered design variable, ρ is the original design variable; The projection function form is: Where: fina is the projected design variable, β is the steepness value, and θ is the threshold value.

7. The phase change heat sink design method based on topology optimization according to claim 1, characterized in that: In step 4), the method for reconstructing the topology optimization configuration is: The binarization method is used to convert the topology optimization configuration into a model without grayscale units by adjusting the threshold value; the binarized matrix is ​​traversed and the boundary density matrix is ​​extracted, and then the matrix element positions are converted into rectangular coordinates for storage, and the topology optimization structure is converted into a parameterized model; the contour boundary points of the topology optimization results are tracked and sorted using chain codes, and all contour boundary points are traversed in a clockwise or counterclockwise order to obtain the chain code geometry of the contour boundary information, and the coordinate information of the boundary points is saved in order; The geometric information of the contour is expressed by using differential chain code instead of the discrete curvature of the contour, and the contour chain code value is weighted averaged. The topological contour segmentation points are determined according to the contour curvature mutation points by analyzing the differential chain code graph. The B-spline interpolation method is used to reconstruct the contour boundary.

Images