CN102608917B - Method for simulating high-temperature phase change heat storage value with hole under microgravity condition - Google Patents

Abstract

The invention provides a method for simulating a high-temperature phase change heat storage value with a hole under the microgravity condition. Under the condition of simultaneously considering the phase change, the hole and the radiation, the value calculation is carried out on a complex phase change heat transfer problem, so that technical personnel can acquire the variation of the melting rate of a solid/liquid phase change material in a high temperature heat storage container on site and the distribution of a flow field, a temperature field and a liquid phase in the phase changing process by utilizing a computer, and therefore, the important reference basis is provided for optimizing the design of the heat storage container, inhibiting the hole, improving the heat storage efficiency and reducing the phenomena of hot spots and hot loosening. According to the method provided by the invention, by a calculating program on the basis of an enthalpy method and a finite control volume method, the complex phase change problem of melting/solidification in the high temperature heat storage container with the hole under the microgravity condition is simply, conveniently and efficiently solved. The method has important practical value.

Description

With cuniculate high-temperature phase change heat accumulation method for numerical simulation under a kind of microgravity

Technical field

The present invention relates to a kind of method in Fluid Mechanics Computation field, be specifically one kind under microgravity with cuniculate high-temperature phase change heat accumulation method for numerical simulation, i.e., by calculating enthalpy method energy hole equation group.

Background technology

For a long time, people have ignored solid-liquid variable density in the research process to Phase-change Problems and the influence that hole is brought produced in container.Under normal circumstances, density of the material in solid-state and liquid is unequal.For phase-change material (PCM, Phase Change Material) in container for the certain PCM of quality, the change of PCM density correspondingly causes PCM Volume Changes in phase transition process, causes the residual volume beyond the volume that PCM is occupied in PCM containers to increase or reduce.The change of residual volume causes the change of thermal resistance, and then influences the progress of phase transition process.For LiF, cubical contraction during solidification is about 23%, so high cubical contraction, and the influence to phase transition process is very big, it is impossible to simply ignored, the change of solid-liquid density is must take into consideration in phase transition process.

PCM (such as LiF, CaF₂Deng) solidification when volume contraction consequence be to form hole in PCM containers.The presence in hole increases thermal resistance, and produces influence to the direction that phase transition process is carried out.The presence in hole is also possible to make PCM containers produce " hot spot " and " heat is released " phenomenon, and both phenomenons can all cause the destruction of PCM containers.

Kerslake (Kerslake T W, Ibrahim M B.Analysis of thermal energy storage material with change-of phase volumetric effects.Trans.ASME, Journal of Solar Energy Engineering, 1993,115 (2)：22~One Dimension Analysis result 31) shows, because the hole that PCM solidification shrinkages are produced makes PCM container wall peak temperatures rise 200K.But, the two side of PCM containers will be played a significant role in diabatic process, and substantial amounts of heat can be passed to inside PCM containers via side wall, so as to weaken influence of the hole to diabatic process.To quantify hole influence, Kerslake (Kerslake T W, Ibrahim M B.Two-dimensional model of a space station freedom thermal energy storage canister.Trans.ASME, Journal of Solar Energy Engineering, 1992,114 (5)：114~121) by taking the PCM containers in SSF heat dumps as an example, two-dimentional (r, z) calculation procedure has been write, numerical solution has been carried out to microgravity state periodical phase-change heat transfer process.Heat transfer process includes the radiation heat transfer between steam heat conduction and each surface in shape and position of the hole model given hole in PCM containers, hole.Kerslake assumes that hole is located at outer wall, and cavity volume keeps invariable.But, in actual moving process, with PCM fusing and solidification alternately, the volume in hole is also reduced or increased therewith in container, and is not maintained as definite value.

Remove Jiang great Peng (the design portion master thesis of " CALCULATION OF THERMAL of space solar power set heat absorption-thermal storage device phase-change material cavity volume " National Space Industry Corporation the 5th 01 in the country, 1991.7) once the solid-liquid phase change heat-accumulating process containing hole has been carried out beyond Two-Dimensional Thermal Analysis, do not carry out deeper into research, there is no correlation document can be for reference.And the heat analysis that Jiang great Peng is done also simply can be seen that in result of calculation and consider hole, but it is not specified how calculate and handle the Volume Changes in hole in phase transition process.It is domestic also to there is blank with the research in terms of void nucleation and the solid-liquid phase-change process of development under microgravity.

The content of the invention

The purpose of the present invention is to solve the shortcomings of the prior art, there is provided under a kind of microgravity condition with cuniculate high-temperature phase change heat accumulation method for numerical simulation, phase transformation is considered at the same time, under conditions of hole and radiation, numerical computations are carried out to complicated phase-change heat transfer problem, so that technical staff just can obtain solid-liquid phase-change material (PCM in live high-temperature heat accumulation container using computer, Phase Change Material) melting rate change and phase transformation when flow field, temperature field, liquid phase is distributed, so as to be designed for optimization heat storage container, suppress hole, improve heat storage efficiency and reduce " hot spot " and " heat is released " phenomenon and important reference frame is provided.

To realize above-mentioned target, the present invention is integrated using Finite Volume Method for Air to control volume, and phase transition process is handled using enthalpy method.

The present invention is specially for the technical scheme taken of its technical problem of solution：

With solid-liquid phase transformation method for numerical simulation in cuniculate high-temperature heat accumulation container under a kind of microgravity condition, it is characterised in that the method for numerical simulation is comprised the following specific steps that：

(1) geometrical model according to high-temperature heat accumulation container divides liquid and solid zoning, and mesh generation is carried out to above-mentioned zoning using triangle unstructured grid；

(2) defining operation condition：R is to by Action of Gravity Field, and reference pressure is an atmospheric pressure；Define the natural convection model assumed based on BOUSSINESQ；

(3) unstable state fusing/SOLIDIFICATION MODEL is defined, discrete-ordinates radiation model is defined；

(4) material properties of the phase-change material under gas, solid, liquid three-phase are defined, the material properties specifically include density, specific heat capacity, thermal conductivity, viscosity, heat of fusion, liquidus curve and solidus temperature, absorptivity and emissivity；

(5) following partial differential energy hole equation group is set up：

$\frac{\∂\left(\mathrm{\ρe}\right)}{\∂t}+\frac{\∂\left(\mathrm{\ρue}\right)}{\∂z}+\frac{1}{r}\frac{\∂\left(\mathrm{r\ρve}\right)}{\∂r}=k{\▿}^{2}T---\left(1\right)$

In formula：E is specific enthalpy, J/kg；ρ is density, kg/m³；T is time, s；T is temperature, K；U, v are respectively z directions, the velocity component in r directions；The influence of phase-change material phase transformation is taken into account in the form of specific enthalpy e in formula, thus, the expression formula is applied to whole domain；The form that specific enthalpy e and temperature T relational expression can be expressed as：

$T=\left\{\begin{array}{cc}{T}_m+e/c,& e\&le;0\\ T_m,& 0<e<\mathrm{\&Delta;}{H}_m\\ T_m+(e-\mathrm{\&Delta;}{H}_m)/c,e& \&GreaterEqual;\mathrm{\&Delta;}{H}_m\end{array}\right.---\left(2\right)$

In formula：T_mFor phase transition temperature, K；ΔH_mFor the latent heat of phase change of material unit mass,

J/kg；

(6) boundary condition is defined：It is hot-fluid border that the inwall of high-temperature heat accumulation container, which is defined, for Equations of The Second Kind heat transfer border, and outer wall and side wall are adiabatic boundary, Gu other walls are stream/coupling boundary；Wherein, the hot-fluid border be periodicity hot-fluid border, setting sunshine period and the shade phase time, property write cycle hot-fluid SQL, wherein：

A) sunshine period heat flow density q_sunSpecific functional relation is as follows：

q_sun=15424-24.88 (T_wall- 1020), W/m²；

B) shade phase heat flow density q_shadowSpecific functional relation is as follows：

q_shadow=-12314-24.88 (T_wall- 1020), W/m²；

Wherein, T_wallFor the inner wall temperature of high-temperature heat accumulation container；

(7) primary condition is defined：Whole zoning is initialized, setting initial temperature and initial velocity；

(8) setting monitoring phase-change material temperature and liquid phase fraction distribution, setting monitoring phase-change material melting rate change；

(9) discretization is carried out to the partial differential energy hole equation in step (5), and the boundary condition and primary condition that are defined using step (6), (7) are closed and solved；

(10) whole zoning is initialized, setting time step-length and iterations, iterative calculation is repeated to the Algebraic Equation set in zoning, untill the iteration precision set by satisfaction, completed under microgravity with solid-liquid phase transformation numerical simulation in cuniculate high-temperature heat accumulation container；

(11) result of calculation is post-processed, draws out cloud atlas and correlation curve.

Further, the present invention carries out discrete region using interior nodes method.

Further, the present invention exports discrete equation using control volume integral method.

Further, the present invention carries out discretization solution using complete explicit form.

Further, in step (6), sunshine period is set as 54min, and the shade phase is 36min.

Further, in step (9), cell node P (i, j, k) place of any non-wall, the discrete form of above-mentioned partial differential energy hole equation is

$E_P^{t+\mathrm{\&Delta;t}}-E_P^t=[{\left(\mathrm{kr}\right)}_n(T_N^t-T_P^t)-{\left(\mathrm{kr}\right)}_s(T_P^t-T_S^t)]\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T}-k_b\frac{T_P^t-T_B^t}{r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B})\mathrm{\&Delta;r\&Delta;z\&Delta;t}+(k_e\frac{T_E^t-T_P^t}{\mathrm{\&Delta;}{z}_E}-k_w\frac{T_P^t-T_W^t}{\mathrm{\&Delta;}{z}_W})r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;t};$

Receiving solar radiation hot-fluid q for inwall node, outside it, (θ, t), inwall node discrete are turned to

$E_P^{t+\mathrm{\&Delta;t}}-E_P^t={\left(\mathrm{kr}\right)}_s(T_S^t-T_P^t)\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}+q(\mathrm{\&theta;},t)r_i\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T}-k_b\frac{T_P^t-T_B^t}{r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B})\mathrm{\&Delta;r\&Delta;z\&Delta;t}+(k_e\frac{T_E^t-T_P^t}{\mathrm{\&Delta;}{z}_E}-k_w\frac{T_P^t-T_W^t}{\mathrm{\&Delta;}{z}_W})r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;t},$

Wherein：

r_iRepresent the external diameter of PCM container inner walls；

E, W, N, S, T, B represent node P six adjacent nodes respectively：(i+1, j, k), (i-1, j, k), (i, j+1, k), (i, j-1, k), (i, j, k+1) and (i, j, k-1)；

(kr)_n、(kr)_s、k_e、k_w、k_t、k_bFor：

${\left(\mathrm{kr}\right)}_n=\frac{1}{\frac{\mathrm{\&Delta;}{r}_N^{-}}{k_N(r_N-\frac{\mathrm{\&Delta;}{r}_N^{-}}{2})}+\frac{\mathrm{\&Delta;}{r}_P^{+}}{k_P(r_P+\frac{\mathrm{\&Delta;}{r}_P^{+}}{2})}},$ ${\left(\mathrm{kr}\right)}_s=\frac{1}{\frac{\mathrm{\&Delta;}{r}_S^{+}}{k_S(r_S+\frac{\mathrm{\&Delta;}{r}_S^{+}}{2})}+\frac{\mathrm{\&Delta;}{r}_P^{-}}{k_P(r_P-\frac{\mathrm{\&Delta;}{r}_P^{-}}{2})}},$

$k_e=\frac{\mathrm{\&Delta;}{z}_E}{\frac{{\left(\mathrm{\&Delta;}{z}_E\right)}^{+}}{k_E}+\frac{{\left(\mathrm{\&Delta;}{z}_E\right)}^{-}}{k_P}},$ $k_w=\frac{\mathrm{\&Delta;}{z}_W}{\frac{{\left(\mathrm{\&Delta;}{z}_W\right)}^{+}}{k_P}+\frac{{\left(\mathrm{\&Delta;}{z}_W\right)}^{-}}{k_W}}$

$k_t=\frac{r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T}{\frac{{\left(r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T\right)}^{+}}{k_T}+\frac{{\left(r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T\right)}^{-}}{k_P}},$ $k_b=\frac{r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B}{\frac{{\left(r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B\right)}^{+}}{k_P}+\frac{{\left(r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B\right)}^{-}}{k_B}}.$

Further, in step (10), the Limit of Stability of time step is：

$\mathrm{\&Delta;t}\&le;\frac{(\mathrm{\&rho;}{c}_P/k)r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;z}}{(\frac{1}{\frac{\mathrm{\&Delta;r}}{r_N-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_P+\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}}+\frac{1}{\frac{\mathrm{\&Delta;r}}{r_P-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_S+\frac{\mathrm{\&Delta;r}}{4}}})2\mathrm{\&Delta;\&theta;\&Delta;z}+\frac{2r_P\mathrm{\&Delta;r\&Delta;\&theta;}}{\mathrm{\&Delta;z}}+\frac{2\mathrm{\&Delta;r\&Delta;z}}{r_P\mathrm{\&Delta;\&theta;}}}.$

Further, the present invention is used as the PCM undergone phase transition quality using the difference Δ m for calculating the solid-state PCM of former and later two moment t and t+ Δ t of each time step quality in PCM containers.

Further, the gross mass of PCM in current time container will be checked in each time step, if there is deviation, drift correction is carried out to Δ m, so as to ensure the conservation of mass；If Δ m is not 0, illustrate that PCM there occurs phase transformation.

Further, if revised Δ m ＞ 0, solid-state PCM quality increase is illustrated, PCM is solidified；If revised Δ m ＜ 0, illustrate solid-state PCM Mass lost, PCM is melted；If revised Δ m=0, illustrate that solid-state PCM quality does not change, PCM is in explicit neither endothermic nor exothermic state.

Further, if Δ m is not 0, illustrate that PCM there occurs phase transformation, the volume change Δ V of PCM in each time step is obtained by Δ m, the PCM volume change Δ V obtained according to calculating are adjusted to the cavity volume of PCM containers：If Δ V ＞ 0, illustrate PCM volume increase, be reduced by the volume shared by hole；If Δ V ＜ 0, illustrate PCM volume-diminished, correspondingly increase cavity volume.

Further, the calculation of radiation heat transferring in PCM holes be based on it is assumed hereinafter that：(1) all cavity surfaces are diffusing reflection grey body surface；(2) radiation of PCM Surface absorptions all wavelengths；(3) hole inner vapor is not involved in radiation heat transfer；Radiation heat transfer between the interface of hole is concentrated mainly on radially, and the radiation heat transfer of axial and circumferential compares much weaker, to simplify calculating, negligible axial and circumferential radiation heat exchange.

Further, hole effective heat transfer coefficient k_veffDetermine as the following formula

$k_{\mathrm{veff}}=k_{\mathrm{PCM}}+k_{\mathrm{reff}}=k_{\mathrm{PCM}}+\ln \left(\frac{r_w}{r_v}\right)\&CenterDot;\frac{r_v\mathrm{\&sigma;}(T_w\left(t\right)+T_v\left(t\right))(T_w^2\left(t\right)+T_v^2\left(t\right))}{\frac{1}{{\mathrm{\&epsiv;}}_{\mathrm{PCM}}}+\frac{r_v}{r_w}(\frac{1}{{\mathrm{\&epsiv;}}_w}-1)},$

K in formula_reffFor the equivalent equivalent heat conductivity of radiation heat transfer between the interface of hole, σ is stefan boltzmann's constant (σ=5.67 × 10^-8W·m^-2·K^-4), T_WFor outside wall temperature, T_VIt is PCM in the temperature of hole interface, ε_PCMFor PCM emissivity, ε_WFor chamber wall emissivity, v represents hole interface.

The inventive method passes through the calculation procedure based on enthalpy method and finite volume algorithm so that under microgravity with the fusing in cuniculate high-temperature heat accumulation container/solidify this complicated Phase-change Problems obtained it is easy, efficient solve, with important practical value.

The discretization method of partial differential governing equation in the present invention and the method for solving of algebraic equation are introduced below.

To carry out numerical computations, first have to carry out discrete region, i.e., domain is divided into the subregion of many non-overlapping copies, representational point is the unknown variable value on node to represent the unknown variable of consecutive variations in region in subregion.According to the difference of node position in subregion, discrete region method can be divided into two classes：Exterior node method and interior nodes method.Because the situation that interior nodes method handles physical property change is more convenient, therefore present invention preferably employs interior nodes method.

After the discretization of domain, it is possible to by control differential equation on each node it is discrete, set up corresponding discrete equation for each node.The discrete method of control differential equation has Taylor series expansion method, polynomial fitting method, control volume integral method, balancing method and calculus of variations etc..Control volume integral method is also known as Finite Volume Method for Air, is widely used method in numerical computations of conducting heat.The clear physics conception of this method derivation, derivation result has clear and definite physical significance, and the conservation property of discrete equation can be guaranteed.Preferably, the present invention exports discrete equation using control volume integral method.

For the selection of difference scheme, Thibault has ever done detailed research.By solving three dimentional heat conduction problem, the result of calculation for comparing three kinds of explicit forms, four kinds of ADI forms and two kinds of implied formats and the error accurately solved, easy programming, required calculating time and the requirement to calculating amount of storage, most preferably two kinds of ADI forms are pointed out, complete explicit form comes the 3rd.The major defect of complete explicit form is limited by stability, and time step can not obtain excessive.For the research of the present invention, because the density of PCM in phase transition process changes, there is hole in container, with the progress of phase transition process, cavity volume increase or diminution, this just determines that time step can not obtain too big, to avoid result of calculation distortion.When time step is smaller, the advantage of implied format is also just not present, and due to needing iteration, calculates the time longer than display format on the contrary.Therefore, discretization solution is carried out present invention preferably employs complete explicit form.

By PCM containers by r × θ × z (i.e. cylindrical coordinates mode, r represents radially, θ represent circumferential, z represent it is axial) grid division unit.Any cell node P (i, j, k) control volume is as shown in Figure 1.Axial cross section where node P is as shown in Figure 2.

Energy hole equation is integrated on node P (i, j, k) control volume (referring to Fig. 1)：

$\&Integral;\&Integral;\&Integral;\&Integral;\frac{\&PartialD;\left(\mathrm{\&rho;e}\right)}{\&PartialD;t}\mathrm{dVdt}=\&Integral;\&Integral;\&Integral;\&Integral;\frac{1}{r}\frac{\&PartialD;}{\&PartialD;r}\left(\mathrm{kr}\frac{\&PartialD;T}{\&PartialD;r}\right)\mathrm{dVdt}\&Integral;\&Integral;\&Integral;\&Integral;+\frac{1}{r}\frac{\&PartialD;}{\&PartialD;\mathrm{\&theta;}}\left(\frac{k}{r}\frac{\&PartialD;T}{\&PartialD;\mathrm{\&theta;}}\right)\mathrm{dVdt}+$

$\&Integral;\&Integral;\&Integral;\&Integral;\frac{\&PartialD;}{\&PartialD;z}\left(k\frac{\&PartialD;T}{\&PartialD;z}\right)\mathrm{dVdt}---\left(3\right)$

Expression formula (3) equal sign left side transient terms are distributed using stairstepping over time and space, i.e.,

${\&Integral;}_r^{r+\mathrm{\&Delta;r}}{\&Integral;}_{\mathrm{\&theta;}}^{\mathrm{\&theta;}+\mathrm{\&Delta;\&theta;}}{\&Integral;}_z^{z+\mathrm{\&Delta;z}}{\&Integral;}_t^{t+\mathrm{\&Delta;t}}\frac{\&PartialD;\left(\mathrm{\&rho;e}\right)}{\&PartialD;t}\mathrm{rdrd\&theta;dzdt}=\mathrm{\&rho;}(e^{t+\mathrm{\&Delta;t}}-e^t)r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;z}=E_P^{t+\mathrm{\&Delta;t}}-E_P^t---\left(4\right)$

E in formula_PRepresent the enthalpy of PCM in node P (i, j, k) places control volume, m_PThe quality of PCM in volume is controlled where node P (i, j, k)；

Diffusion term is distributed using stairstepping in time on the right of expression formula (3) equal sign, is spatially distributed, then had using piecewise linearity

${\&Integral;}_r^{r+\mathrm{\&Delta;r}}{\&Integral;}_{\mathrm{\&theta;}}^{\mathrm{\&theta;}+\mathrm{\&Delta;\&theta;}}{\&Integral;}_z^{z+\mathrm{\&Delta;z}}{\&Integral;}_t^{t+\mathrm{\&Delta;t}}\frac{1}{r}\frac{\&PartialD;}{\&PartialD;r}\left(\mathrm{kr}\frac{\&PartialD;T}{\&PartialD;r}\right)\mathrm{rdrd\&theta;dzdt}=[{\left(\mathrm{kr}\right)}_n(T_N^t-T_P^t)-{\left(\mathrm{kr}\right)}_s(T_P^t-T_S^t)]\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}---\left(5\right)$

${\&Integral;}_r^{r+\mathrm{\&Delta;r}}{\&Integral;}_{\mathrm{\&theta;}}^{\mathrm{\&theta;}+\mathrm{\&Delta;\&theta;}}{\&Integral;}_z^{z+\mathrm{\&Delta;z}}{\&Integral;}_t^{t+\mathrm{\&Delta;t}}\frac{1}{r}\frac{\&PartialD;}{\&PartialD;\mathrm{\&theta;}}\left(\frac{k}{r}\frac{\&PartialD;T}{\&PartialD;\mathrm{\&theta;}}\right)\mathrm{rdrd\&theta;dzdt}=(k_t\frac{T_T^t-T_P^t}{r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T}-k_b\frac{T_P^t-T_B^t}{r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B})\mathrm{\&Delta;r\&Delta;z\&Delta;t}---\left(6\right)$

${\&Integral;}_r^{r+\mathrm{\&Delta;r}}{\&Integral;}_{\mathrm{\&theta;}}^{\mathrm{\&theta;}+\mathrm{\&Delta;\&theta;}}{\&Integral;}_z^{z+\mathrm{\&Delta;z}}{\&Integral;}_t^{t+\mathrm{\&Delta;t}}\frac{\&PartialD;}{\&PartialD;z}\left(k\frac{\&PartialD;T}{\&PartialD;z}\right)\mathrm{rdrd\&theta;dzdt}=(k_e\frac{T_E^t-T_P^t}{\mathrm{\&Delta;}{z}_E}-k_w\frac{T_P^t-T_W^t}{\mathrm{\&Delta;}{z}_W})r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;t}---\left(7\right)$

In formula, E, W, N, S, T, B represent node P six adjacent nodes respectively：(i+1, j, k), (i-1, j, k), (i, j+1, k), (i, j-1, k), (i, j, k+1) and (i, j, k-1).(kr)_n、(kr)_s、k_e、k_w、k_t、k_bMeaning is as follows.

${\left(\mathrm{kr}\right)}_n=\frac{1}{\frac{\mathrm{\&Delta;}{r}_N^{-}}{k_N(r_N-\frac{\mathrm{\&Delta;}{r}_N^{-}}{2})}+\frac{\mathrm{\&Delta;}{r}_P^{+}}{k_P(r_P+\frac{\mathrm{\&Delta;}{r}_P^{+}}{2})}},$ ${\left(\mathrm{kr}\right)}_s=\frac{1}{\frac{\mathrm{\&Delta;}{r}_S^{+}}{k_S(r_S+\frac{\mathrm{\&Delta;}{r}_S^{+}}{2})}+\frac{\mathrm{\&Delta;}{r}_P^{-}}{k_P(r_P-\frac{\mathrm{\&Delta;}{r}_P^{-}}{2})}},$

$k_e=\frac{\mathrm{\&Delta;}{z}_E}{\frac{{\left(\mathrm{\&Delta;}{z}_E\right)}^{+}}{k_E}+\frac{{\left(\mathrm{\&Delta;}{z}_E\right)}^{-}}{k_P}},$ $k_w=\frac{\mathrm{\&Delta;}{z}_W}{\frac{{\left(\mathrm{\&Delta;}{z}_W\right)}^{+}}{k_P}+\frac{{\left(\mathrm{\&Delta;}{z}_W\right)}^{-}}{k_W}}$

$k_t=\frac{r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T}{\frac{{\left(r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T\right)}^{+}}{k_T}+\frac{{\left(r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T\right)}^{-}}{k_P}},$ $k_b=\frac{r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B}{\frac{{\left(r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B\right)}^{+}}{k_P}+\frac{{\left(r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B\right)}^{-}}{k_B}}$

The discrete form for obtaining partial differential energy equation (1) by formula (4), (5), (6), (7) is as follows

$E_P^{t+\mathrm{\&Delta;t}}-E_P^t=[{\left(\mathrm{kr}\right)}_n(T_N^t-T_P^t)-{\left(\mathrm{kr}\right)}_s(T_P^t-T_S^t)]\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T}-k_b\frac{T_P^t-T_B^t}{r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B})\mathrm{\&Delta;r\&Delta;z\&Delta;t}+(k_e\frac{T_E^t-T_P^t}{\mathrm{\&Delta;}{z}_E}-k_w\frac{T_P^t-T_W^t}{\mathrm{\&Delta;}{z}_W})r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;t}---\left(8\right)$

Receiving solar radiation hot-fluid q for inwall node, outside it, (θ, t), node discrete are turned to

$E_P^{t+\mathrm{\&Delta;t}}-E_P^t={\left(\mathrm{kr}\right)}_s(T_S^t-T_P^t)\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}+q(\mathrm{\&theta;},t)r_i\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T\mathrm{\&Delta;}{\mathrm{\&theta;}}_T}-k_b\frac{T_P^t-T_B^t}{r_B\mathrm{\&Delta;}{\mathrm{\&theta;}}_B})\mathrm{\&Delta;r\&Delta;z\&Delta;t}+(k_e\frac{T_E^t-T_P^t}{\mathrm{\&Delta;}{z}_E}-k_w\frac{T_P^t-T_W^t}{\mathrm{\&Delta;}{z}_W})r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;t}---\left(9\right)$

R in formula_iRepresent the external diameter of PCM container inner walls.

When the heat exchange amount with adjacent cells is calculated if the adjacent cells of certain unit are hole unit, in formula (8), (9), it should be calculated using hole-recombination heat transfer coefficient by the heat exchange of the non-cavitated unit adjacent with hole unit.

${\left(\mathrm{kr}\right)}_n=\frac{2k}{\frac{\mathrm{\&Delta;r}}{r_N-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_P+\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}},$ ${\left(\mathrm{kr}\right)}_s=\frac{2k}{\frac{\mathrm{\&Delta;r}}{r_P-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_S+\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}},$

Arrange

$T_P^{t+\mathrm{\&Delta;t}}=T_P^t+\frac{1}{\mathrm{\&rho;}{c}_P{r}_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;z}}(\frac{T_N^t-T_P^t}{\frac{\mathrm{\&Delta;r}}{r_N-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_P+\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}}+\frac{T_S^t-T_P^t}{\frac{\mathrm{\&Delta;r}}{r_P-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_S+\frac{\mathrm{\&Delta;r}}{4}}})2\mathrm{k\&Delta;\&theta;\&Delta;z\&Delta;t}+$

$\frac{1}{\mathrm{\&rho;}{c}_P{r}_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;z}}[(\frac{T_T^t-T_P^t}{r_P\mathrm{\&Delta;\&theta;}}-\frac{T_P^t-T_B^t}{r_P\mathrm{\&Delta;\&theta;}})\mathrm{k\&Delta;r\&Delta;z\&Delta;t}+\frac{k{r}_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;t}}{\mathrm{\&Delta;z}}(T_E^t+T_W^t-2T_P^t)]$

$=\frac{\mathrm{k\&Delta;t}}{\mathrm{\&rho;}{c}_P{r}_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;z}}[(\frac{T_N^t}{\frac{\mathrm{\&Delta;r}}{r_N-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_P+\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}}+\frac{T_S^t}{\frac{\mathrm{\&Delta;r}}{r_P-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_S+\frac{\mathrm{\&Delta;r}}{4}}})2\mathrm{\&Delta;\&theta;\&Delta;z}+\frac{\mathrm{\&Delta;r\&Delta;z}}{r_P\mathrm{\&Delta;\&theta;}}(T_T^t+T_B^t)$

$+\frac{r_P\mathrm{\&Delta;r\&Delta;\&theta;}}{\mathrm{\&Delta;z}}(T_E^t+T_W^t)]+\{1-\frac{\mathrm{k\&Delta;t}}{\mathrm{\&rho;}{c}_P{r}_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;z}}[(\frac{1}{\frac{\mathrm{\&Delta;r}}{r_N-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_P+\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}}+\frac{1}{\frac{\mathrm{\&Delta;r}}{r_P+\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_S+\frac{\mathrm{\&Delta;r}}{4}}})2\mathrm{\&Delta;\&theta;\&Delta;z}$

$+\frac{2\mathrm{\&Delta;r\&Delta;z}}{r_P\mathrm{\&Delta;\&theta;}}+\frac{2r_P\mathrm{\&Delta;r\&Delta;\&theta;}}{\mathrm{\&Delta;z}}\left]\right\}{T}_p^t---\left(10\right)$

Ensureing the condition of stability is

Coefficient be more than or equal to 0, i.e.,

$\mathrm{\&Delta;t}\&le;\frac{(\mathrm{\&rho;}{c}_P/k)r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;z}}{(\frac{1}{\frac{\mathrm{\&Delta;r}}{r_N-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_P+\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}}+\frac{1}{\frac{\mathrm{\&Delta;r}}{r_P-\frac{\mathrm{\&Delta;<figure-callout\; id="4"\; label="r"\; filenames="HDA0000136673010000012.png"\; state="\{\{state\}\}">r</figure-callout>}}{4}}+\frac{\mathrm{\&Delta;r}}{r_S+\frac{\mathrm{\&Delta;r}}{4}}})2\mathrm{\&Delta;\&theta;\&Delta;z}+\frac{2r_P\mathrm{\&Delta;r\&Delta;\&theta;}}{\mathrm{\&Delta;z}}+\frac{2\mathrm{\&Delta;r\&Delta;z}}{r_P\mathrm{\&Delta;\&theta;}}}---\left(11\right)$

This formula is the Limit of Stability of step-length.

The calculating that cavity volume changes is introduced with processing below.

The present invention is the phase-change heat transfer process changed in accurate simulation container along with cavity volume, introduces a variable --- cavity volume fraction f_v, represent the relative size of cavity volume in unit.f_vFor 0 when, without hole in representative unit；f_vFor 1 when, is full of by hole in representative unit；f_vWhen between 0 and 1, representative unit is with the presence of hole.By introducing cavity volume fraction f_vThe present invention proposes the algorithm for calculating cavity volume change and hole adjustment, the heat transfer process mathematical modeling of the high-temperature phase change heat accumulation container with void nucleation and development is established, the CALCULATION OF THERMAL of the phase-change heat transfer process changed under microgravity along with cavity volume is completed.

Due to there is hole in PCM containers, this is required after to governing equation discretization,, it is necessary to being judged in the unit and adjacent cells with the presence or absence of hole when calculating the heat exchange amount of each unit and adjacent cells, to be computed correctly PCM quality in the unit and the heat exchange amount with adjacent cells.

In the presence of having hole, PCM quality should be in PCM cell

m_{I, j, k}=ρ_{I, j, k}V_{I, j, k}(1-f_{Vi, j, k})

I=1,2 ... ..., II；J=1,2 ... ..., JJ；K=1,2 ... ..., KK

                                              (12)

I, j, k represent (i, j, k) individual unit in formula, and II, JJ, KK represent the unit number that each coordinate direction of PCM containers is divided, f_vFor cavity volume fraction, concrete meaning is as follows：f_v=0, no hole；0 ＜ f_v＜ 1, element memory is in partial holes；f_v=1, all holes in unit.

In addition, with the generation of PCM phase transformations, the PCM undergone phase transition density changes therewith, it is necessary to consider the change of cavity volume in some units caused by PCM Volume Changes, i.e. f_vChange therewith.PCM volume change can be obtained by the PCM undergone phase transition Mass Calculation in container.In each time step, the PCM undergone phase transition quality is the difference of former and later two moment PCM containers internal solid or liquid PCM quality.Due to PCM volumetric expansions or contraction in phase transition process, hole in liquid PCM meeting cartridges, or extracted out from certain unit, so that the unit is changed into even complete empty containing partial holes, it is inconsistent before and after the liquid PCM so calculated in each time step quality, and solid-state PCM quality is constant.To avoid the result of calculation of mistake, the present invention is using former and later two moment for calculating each time step in PCM containers：The difference Δ m of the solid-state PCM of t and t+ Δ ts quality as the PCM undergone phase transition quality, i.e.,

$\mathrm{\&Delta;m}=\underset{i,j,k}{\mathrm{\&Sigma;}}[(1-f_{\mathrm{li},j,k}^{t+\mathrm{\&Delta;t}})\&CenterDot;(1-f_{\mathrm{vi},j,k}^{t+\mathrm{\&Delta;t}})\&CenterDot;{\mathrm{\&rho;}}_{i,j,k}^{t+\mathrm{\&Delta;t}}-(1-f_{\mathrm{li},j,k}^t)\&CenterDot;(1-f_{\mathrm{vi},j,k}^t)\&CenterDot;{\mathrm{\&rho;}}_{i,j,k}^t]\&CenterDot;V_{i,j,k}---\left(13\right)$

F in formula_lFor PCM liquid phase volume fractions.

To ensure that in the conservation of mass, each time step the gross mass of PCM in current time container will be checked, if there is deviation, drift correction is carried out to Δ m.Judge revised Δ m value, if Δ m ＞ 0, illustrate solid-state PCM quality increase, PCM is solidified；If Δ m ＜ 0, illustrate solid-state PCM Mass lost, PCM is melted；If Δ m=0, illustrate that solid-state PCM quality does not change, PCM is in explicit neither endothermic nor exothermic state.If Δ m is not 0, illustrate that PCM there occurs phase transformation, by the difference Δ m of former and later two moment solid-state PCM quality, it is possible to obtain the volume change Δ V of PCM in each time step, i.e.,

Δ V=Δs m (1/ ρ_s-1/ρ_l)          (14)

In formula：ρ_sRepresent solid-state PCM density, ρ_lRepresent liquid PCM density.The PCM volume change Δ V obtained according to calculating, are adjusted to the cavity volume of PCM containers.If Δ V ＞ 0, illustrate PCM volume increase, be reduced by the volume shared by hole；If Δ V ＜ 0, illustrate PCM volume-diminished, correspondingly increase cavity volume.Because assuming that hole is located at outer wall, therefore reduce the volume shared by hole, carried out by order from inside to outside, and increase cavity volume and carry out in the opposite order.

The principle that should follow of adjustment cavity volume is：If f_v＜ 1, while mf_lThe non-full sky of ＞ 0, i.e. unit, with the presence of liquid PCM, then can therefrom cut down PCM volume, that is, increase cavity volume；If f_v＞ 0, i.e. unit underfill can then fill liquid PCM into unit, that is, reduce cavity volume.

The calculating and processing to hole effective heat transfer coefficient are introduced below.

The consideration that volume contraction and heat dump are designed when being solidified due to PCM, always with the presence of a certain proportion of hole in PCM containers.For LiF, hole accounts for PCM volumes of a container percentage and is between 9%~21%, is changed with PCM fusing and solidification.PCM and the chamber wall Volume Changes as caused by thermal expansion are ignored.

The distribution in hole there is no the theory of maturation to quote under microgravity.The space flight test (Namkoong David, Jacqmin D and Szaniszlo.Effect of microgravity on material undergoing melting and freezing-the TES experiment.AIAA 95-0614) of PCM containers has been carried out in STS-62 aerial missions.PCM containers after experiment have carried out tomography research, and experiment is as shown in Figure 3 with PCM containers and photograph result.As can be seen from the figure, solid-state LiF after solidification concentrates on the one end of PCM containers close to radiative heat rejection device rear portion (in Fig. 3 at position 9, the temperature of PCM containers in this place is minimum), and the other end also have accumulated a part of LiF due to the effect of wetability.The hole that container is internally formed is centrally located at the slightly higher region of temperature in container, and tends to towards the maximum direction of heating hot-fluid (in Fig. 3 directly over container).For the PCM containers that the present invention is studied, because heat is taken away from inwall by cycle fluid, therefore when in the shade phase, temperature is relatively low at inwall, and PCM solidifies at inwall first, and hole is finally radially formed at outer wall gradually to external expansion.Therefore can be assumed that hole is located at container outer wall, form a column annular space.

It is assumed that being full of PCM steams in hole, steam pressure is very low, and LiF steam pressures only have 0.933Pa during 1121K, it is believed that LiF steams are substantially filled with hole.Because steam pressure is very low, hole inner vapor quality can be ignored.Therefore the heat exchange for passing through hole includes the radiation heat transfer between hole heat conduction and hole interface.Axial-temperature gradient very little in hole, therefore hole inner vapor Temperature Distribution is by the determination of radial direction steady state heat transfer equation

$\frac{1}{r}\frac{\&PartialD;}{\&PartialD;r}\left(r\frac{\&PartialD;T}{\&PartialD;r}\right)=0---\left(15\right)$

The solution of equation (15) is as follows

T (r)=Alnr+B (16)

In formula

$A=\frac{T\left(r_0\right)-T\left(r_v\right)}{\ln ({\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000136673010000021.png"\; state="\{\{state\}\}">r</figure-callout>}}_0/r_v)},$ $B=\frac{T\left(r_0\right)-\ln [T\left(r_0\right)-T\left(r_v\right)]}{\ln ({\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000136673010000021.png"\; state="\{\{state\}\}">r</figure-callout>}}_0/r_v)}$

Subscript 0 represents the inner surface of PCM container outer walls, and v represents the cavity surface with PCM intersections.

Calculation of radiation heat transferring in hole be based on it is assumed hereinafter that：(1) all cavity surfaces are diffusing reflection grey body surface；(2) radiation of PCM Surface absorptions all wavelengths；(3) hole inner vapor is not involved in radiation heat transfer.In the case of two and three dimensions, transformation interface and hole interface be generally in the shape of distortion, it is irregular, can also be blocked mutually between interface, the radiation heat transfer ascent calculated between the interface of hole is extremely difficult.Although more accurate result can be obtained to calculate radiation heat transfer using methods such as Monte Carlo method, Huo Teer (Hottel) field methods and discrete transfer methods, but these methods are required for substantial amounts of calculating time and computer amount of storage in itself, the calculating time for causing whole phase transition process is greatly prolonged, therefore is uneconomic.In addition, the radiation heat transfer between the interface of hole is concentrated mainly on radially, the radiation heat transfer of axial and circumferential compares much weaker, to simplify calculating, negligible axial and circumferential radiation heat exchange.

To quantify the relative size of radiation heat transfer between hole inner vapor heat conduction and hole interface, the present invention is compared to hole interface radiation heat transfer under PCM container radial direction one-dimensional cases with LiF steam heat conduction.

The thermal conductivity factor k of LiF steams in hole_LiFIt is calculated as follows according to kinetic theory of gases

$k_v=1.457\&times;{10}^{-3}\sqrt{T}---\left(17\right)$

Radiation heat transfer between the interface of hole two is combined in LiF steam heat conduction, that is, considers the heat conduction through hole and radiation heat transfer effect, hole effective heat transfer coefficient k_veffIt can be calculated as follows：

$k_{\mathrm{veff}}=k_{\mathrm{LiF}}+k_{\mathrm{reff}}=k_{\mathrm{LiF}}+\ln \left(\frac{r_w}{r_v}\right)\&CenterDot;\frac{r_v\mathrm{\&sigma;}(T_w\left(t\right)+T_v\left(t\right))(T_w^2\left(t\right)+T_v^2\left(t\right))}{\frac{1}{{\mathrm{\&epsiv;}}_{\mathrm{PCM}}}+\frac{r_v}{r_w}(\frac{1}{{\mathrm{\&epsiv;}}_w}-1)}---\left(18\right)$

K in formula_reffFor the equivalent equivalent heat conductivity of radiation heat transfer between the interface of hole, σ is stefan boltzmann's constant (σ=5.67 × 10^-8W·m^-2·K^-4), T_WFor outside wall temperature, T_VIt is PCM in the temperature of hole interface, ε_PCMFor PCM emissivity, ε_WFor chamber wall emissivity, v represents hole interface.

The hole effective heat transfer coefficient k that Fig. 4 obtains for calculating_veffWith PCM steam thermal conductivity factors k_LiFThe ratio between k_veff/k_LiFChange within the whole orbital period, it can be seen that k_veff/k_LiFChanged between 2.29~4.57, can therefrom obtain k_veff/k_LiFArithmetic mean of instantaneous value in whole cycle is 3.407, i.e. 2.407 times of hole radiative heat exchange amount average out to hole heat conduction amount.In view of the conductive force of PCM container side walls, pressed in three-dimensional computations than more conservative estimation, it is 2 times of LiF steam heat conduction amounts, i.e. hole-recombination heat transfer coefficient k to take hole radiative heat exchange amount_eff=3.0*k_LiF。

Fig. 5 is shown under microgravity of the invention with cuniculate high-temperature phase-change heat exchange calculating program frame chart.

Brief description of the drawings

Fig. 1 is any cell node P (i, j, k) control volume schematic diagram.

Fig. 2 is the axial sectional diagrammatical view illustration where node P (i, j, k).

Fig. 3 is experiment PCM containers and photograph result.

The hole effective heat transfer coefficient k that Fig. 4 obtains for calculating_veffWith PCM steam thermal conductivity factors k_LiFThe ratio between k_veff/k_LiFChange within the whole orbital period.

Fig. 5 is shown under microgravity with cuniculate high-temperature phase-change heat exchange calculating program frame chart.

Fig. 6 is the geometric shape schematic diagram of PCM containers.

Fig. 7 is the axial sectional diagrammatical view illustration of PCM containers.

Fig. 8 calculates result and U.S. NASA scheme high temperature phase-transition heat-storage unit thermal cycle comparison of computational results after stabilization for the present invention.

Embodiment

Technical scheme for a better understanding of the present invention, is further described below in conjunction with accompanying drawing and instantiation, and embodiment is the support to the technology of the present invention feature, rather than limits.

The space station solar dynamic power system system heat pipe receiver high-temperature phase change heat accumulation unit that NASA of physical model system of the present invention (NASA, National Aeronautics and Space Administration) is studied.PCM container profiles are as shown in Figure 6, the axial cross section of container is as shown in Figure 7, PCM heat storage container outer annular diameters 89mm, annular diameters 51mm, inwall, outer wall, sidewall thickness 1.5mm, wall surface material use Hayness188, axial length is 25mm, phase-change material is LiF, and the orbital period is sunshine period 54min, shade phase 36min.

Because the experimental physics model is simple, data are detailed and representative, the present invention carries out numerical computations to the NASA schemes, and result of calculation is compared with NASA result of calculations, to verify its accuracy, with very strong convincingness.Fig. 5 is shown under microgravity with cuniculate high-temperature phase-change heat exchange calculating program frame chart.It mainly realizes that step is as follows：

(1) zoning, non-individual body, border and grid are divided；

(2) setting operation condition：Reference pressure is an atmospheric pressure；

(3) laminar flow, unstable state fusing/SOLIDIFICATION MODEL, discrete-ordinates radiation model are defined；

(4) " material " is selected to be used for setting the medium in flow field in " definition " one column：The liquid defined in flow field is phase-change material LiF, and wall surface material is Hayness188, and it is thin LiF steam to define in hole；

(5) set up the partial differential governing equation group based on energy equation (referring to expression formula (1))；

(6) boundary condition is defined：The inwall for defining high-temperature heat accumulation container is hot-fluid border, and outer wall and side wall are adiabatic boundary, Gu other walls are stream/coupling boundary；Wherein, the hot-fluid border is periodicity hot-fluid border, property write cycle hot-fluid SQL, sunshine period 54min, shade phase 36min；Sunshine period heat flow density is q_sun=15424-24.88 (T_wall- 1020), W/m²；Shade phase heat flow density is q_shadow=-12314-24.88 (T_wall- 1020), W/m²；

(7) initialize, definition initial temperature is 950K, and initial velocity is 0m/s；

(8) internal face temperature change, phase-change material area temperature and liquid phase fraction change are monitored；

(9) the partial differential governing equation group closed to boundary condition contained above carries out discrete (referring to expression formula (8), (9)), and to its iterative；

(10) whole zoning is initialized, iteration time step-length is set to 60s, and iteration precision reaches 10e-6, just it is considered that the result calculated is close to actual effect；

(11) drawn by above-mentioned solution procedure after the solution in each calculate node, it is necessary to be indicated by modes such as line value figure, polar plot, isoline figure, motion pattern, cloud atlas to result of calculation.

Fig. 8 show the present invention and calculates result and U.S. NASA scheme high temperature phase-transition heat-storage unit thermal cycle comparison of computational results after stabilization.By with can be seen that after NASA scheme comparison of computational results：Range of temperature and trend comparison obtained by heat pipe wall, heat storage container outside wall surface range of temperature and trend and NASA schemes are calculated in a cycle obtained by numerical computations of the present invention is close, the temperature value and NASA schemes calculated value of corresponding position each point at most difference are no more than 16K, and this demonstrates the reasonability and accuracy with cuniculate high-temperature phase change heat accumulation unit computational methods, model and program under the microgravity condition that the present invention is set up.

Claims (10)

1. with solid-liquid phase transformation method for numerical simulation in cuniculate high-temperature heat accumulation container under a kind of microgravity condition, it is characterised in that the method for numerical simulation is comprised the following specific steps that：

（1）Geometrical model according to high-temperature heat accumulation container divides liquid and solid zoning, and mesh generation is carried out to above-mentioned zoning using triangle unstructured grid；

（2）Defining operation condition：R is to by Action of Gravity Field, and reference pressure is an atmospheric pressure；Define the natural convection model assumed based on BOUSSINESQ；

（3）Unstable state fusing/SOLIDIFICATION MODEL is defined, discrete-ordinates radiation model is defined；

（4）Material properties of the phase-change material under gas, solid, liquid three-phase are defined, the material properties are specifically made up of density, specific heat capacity, thermal conductivity, viscosity, heat of fusion, liquidus curve and solidus temperature, absorptivity and emissivity；

（5）Set up following partial differential energy hole equation group:

$\frac{\&PartialD;\left(\mathrm{\&rho;e}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;ue}\right)}{\&PartialD;z}+\frac{1}{r}\frac{\&PartialD;\left(\mathrm{r\&rho;ve}\right)}{\&PartialD;r}={k\&dtri;}^{2}T---\left(1\right)$

In formula：E is specific enthalpy, J/kg；ρ is density, kg/m³；T is time, s；T is temperature, K；U, v are respectively z directions, the velocity component in r directions；The influence of phase-change material phase transformation is taken into account in the form of specific enthalpy e in formula, thus, expression formula（1）Suitable for whole domain;The form that specific enthalpy e and temperature T relational expression can be expressed as:

$T=\left\{\begin{array}{cc}{T}_m+e/c,& e\&le;0\\ T_m,& 0<e<{\mathrm{\&Delta;H}}_m\\ T_m+(e-{\mathrm{\&Delta;H}}_m)/c,& e\&GreaterEqual;{\mathrm{\&Delta;H}}_m\end{array}---\left(2\right)\right.$

In formula：T_mFor phase transition temperature, K；ΔH_mFor the latent heat of phase change of material unit mass,

J/kg;

（6）Define boundary condition：It is hot-fluid border that the inwall of high-temperature heat accumulation container, which is defined, for Equations of The Second Kind heat transfer border, and outer wall and side wall are adiabatic boundary, Gu other walls are stream/coupling boundary；Wherein, the hot-fluid border be periodicity hot-fluid border, set sunshine period as 54min, the shade phase be 36min, property write cycle hot-fluid SQL, wherein：

a）Sunshine period heat flow density q_sunSpecific functional relation is as follows：

q_sun=15424-24.88(T_wall-1020),W/m²；

b）Shade phase heat flow density q_shadowSpecific functional relation is as follows：

q_shadow=-12314-24.88(T_wall- 1020), W/m²；

Wherein, T_wallFor the inner wall temperature of high-temperature heat accumulation container；

（7）Define primary condition：Whole zoning is initialized, setting initial temperature and initial velocity；

（8）Setting monitoring phase-change material temperature and liquid phase fraction distribution, setting monitoring phase-change material melting rate change；

（9）To step（5）In partial differential energy hole equation carry out discretization, and utilize step（6）、（7）The boundary condition and primary condition of definition are closed and solved；

（10）Whole zoning is initialized, setting time step-length and iterations, iterative calculation is repeated to the Algebraic Equation set in zoning, untill the iteration precision set by satisfaction, completed under microgravity with solid-liquid phase transformation numerical simulation in cuniculate high-temperature heat accumulation container；

（11）Result of calculation is post-processed, cloud atlas and correlation curve is drawn out.

2. method for numerical simulation according to claim 1, it is characterised in that above-mentioned steps（9）Middle use interior nodes method carries out discrete region, discrete to the progress of partial differential energy hole equation using control volume integral method export, and discretization solution is carried out using complete explicit form.

3. method for numerical simulation according to claim 1, it is characterised in that step（9）In, the cell node P of any non-wall（I, j, k）Locate, the discrete form of above-mentioned partial differential energy hole equation is

$(k_t\frac{T_T^t-T_P^t}{r_T{\mathrm{\&Delta;\&theta;}}_T}-k_b\frac{T_P^t-T_B^t}{r_B{\mathrm{\&Delta;\&theta;}}_B})\mathrm{\&Delta;r\&Delta;z\&Delta;t}+(k_e\frac{T_E^t-T_P^t}{{\mathrm{\&Delta;z}}_E}-k_w\frac{T_P^t-T_W^t}{{\mathrm{\&Delta;z}}_W})r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;t};$

Receiving solar radiation hot-fluid q for inwall node, outside it, (θ, t), inwall node discrete are turned to

$E_P^{t+\mathrm{\&Delta;t}}-E_P^t={\left(\mathrm{kr}\right)}_s(T_S^t-T_P^t)\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}+q(\mathrm{\&theta;},t)r_i\mathrm{\&Delta;\&theta;\&Delta;z\&Delta;t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T{\mathrm{\&Delta;\&theta;}}_T}-k_b\frac{T_P^t-T_B^t}{r_B{\mathrm{\&Delta;\&theta;}}_B})\mathrm{\&Delta;r\&Delta;z\&Delta;t}+(k_e\frac{T_E^t-T_P^t}{{\mathrm{\&Delta;z}}_E}-k_w\frac{T_P^t-T_W^t}{{\mathrm{\&Delta;z}}_W})r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;t},$

Wherein：

I=1,2 ... ..., II；J=1,2 ... ..., JJ；K=1,2 ... ..., KK；I, j, k represent（I, j, k）Individual unit, II, JJ, KK represent the unit number that each coordinate direction of PCM containers is divided；

r_iRepresent phase-change material（PCM）The external diameter of container inner wall；

E, W, N, S, T, B represent node P six adjacent nodes respectively：（I+1, j, k）、（I-1, j, k）、（I, j+1, k）、（I, j-1, k）、（I, j, k+1）With（I, j, k-1）；

(kr)_n、(kr)_s、k_e、k_w、k_t、k_bFor：

${\left(\mathrm{kr}\right)}_n=\frac{1}{\frac{\mathrm{\&Delta;r}\stackrel{-}{N}}{k_N(r_N-\frac{\mathrm{\&Delta;r}\stackrel{-}{N}}{2})}+\frac{\mathrm{\&Delta;r}\stackrel{+}{P}}{k_P(r_P+\frac{\mathrm{\&Delta;r}\stackrel{+}{P}}{2})}},$ ${\left(\mathrm{kr}\right)}_s=\frac{1}{\frac{\mathrm{\&Delta;r}\stackrel{+}{S}}{k_S(r_S+\frac{\mathrm{\&Delta;r}\stackrel{+}{S}}{2})}+\frac{\mathrm{\&Delta;r}\stackrel{-}{P}}{k_P(r_P-\frac{\mathrm{\&Delta;r}\stackrel{-}{P}}{2})}},$

$k_e=\frac{{\mathrm{\&Delta;z}}_E}{\frac{{\left({\mathrm{\&Delta;z}}_E\right)}^{+}}{k_E}+\frac{{\left({\mathrm{\&Delta;z}}_E\right)}^{-}}{k_P}},$ $k_w=\frac{{\mathrm{\&Delta;z}}_W}{\frac{{\left({\mathrm{\&Delta;z}}_W\right)}^{+}}{k_P}+\frac{{\left({\mathrm{\&Delta;z}}_W\right)}^{-}}{k_W}}$

$k_t=\frac{r_T{\mathrm{\&Delta;\&theta;}}_T}{\frac{{\left(r_T{\mathrm{\&Delta;\&theta;}}_T\right)}^{+}}{k_T}+\frac{{\left(r_T{\mathrm{\&Delta;\&theta;}}_T\right)}^{-}}{k_P}},$ $k_b=\frac{r_B{\mathrm{\&Delta;\&theta;}}_B}{\frac{{\left(r_B{\mathrm{\&Delta;\&theta;}}_B\right)}^{+}}{k_P}+\frac{{\left(r_B{\mathrm{\&Delta;\&theta;}}_B\right)}^{-}}{k_B}}.$

4. method for numerical simulation according to claim 3, it is characterised in that step（10）In, the Limit of Stability of time step is：

$\mathrm{\&Delta;t}\&le;\frac{({\mathrm{\&rho;c}}_P/k)r_P\mathrm{\&Delta;r\&Delta;\&theta;\&Delta;z}}{(\frac{1}{\frac{\mathrm{\&Delta;r}}{r_N-\frac{\mathrm{\&Delta;r}}{4}}+\frac{\mathrm{\&Delta;r}}{r_P+\frac{\mathrm{\&Delta;r}}{4}}}+\frac{1}{\frac{\mathrm{\&Delta;r}}{r_P-\frac{\mathrm{\&Delta;r}}{4}}+\frac{\mathrm{\&Delta;r}}{r_S+\frac{\mathrm{\&Delta;r}}{4}}})2\mathrm{\&Delta;\&theta;\&Delta;z}+\frac{{2r}_P\mathrm{\&Delta;r\&Delta;\&theta;}}{\mathrm{\&Delta;z}}+\frac{2\mathrm{\&Delta;r\&Delta;z}}{r_P\mathrm{\&Delta;\&theta;}}}.$

5. method for numerical simulation according to claim 1, it is characterised in that using calculating phase-change material（PCM）The solid phase change material of former and later two moment t and t+ Δ t of each time step in container（PCM）The difference Δ m of quality be used as the phase-change material undergone phase transition（PCM）Quality, Δ m expression formula is：

$\mathrm{\&Delta;m}=\underset{i,j,k}{\mathrm{\&Sigma;}}[(1-f_{\mathrm{li},j,k}^{t+\mathrm{\&Delta;t}})\&CenterDot;(1-f_{\mathrm{vi},j,k}^{t+\mathrm{\&Delta;t}})\&CenterDot;{\mathrm{\&rho;}}_{i,j,k}^{t+\mathrm{\&Delta;t}}-(1-f_{\mathrm{li},j,k}^t)\&CenterDot;(1-f_{\mathrm{vi},j,k}^t)\&CenterDot;{\mathrm{\&rho;}}_{i,j,k}^t]\&CenterDot;V_{i,j,k}$

F in formula_lFor phase-change material（PCM）Liquid phase volume fraction, f_vFor cavity volume fraction, i=1,2 ... ..., II；J=1,2 ... ..., JJ；K=1,2 ... ..., KK；I, j, k represent（I, j, k）Individual unit；f_v=0, unit is interior without hole；0<f_v<1, element memory is in partial holes；f_v=1, all holes in unit.

6. method for numerical simulation according to claim 5, it is characterised in that phase-change material in current time container will be checked in each time step（PCM）Gross mass, if there is deviation, to Δ m carry out drift correction, so as to ensure the conservation of mass；If Δ m is not 0, illustrate phase-change material（PCM）It there occurs phase transformation.

7. the method for numerical simulation according to claim 5 or 6, it is characterised in that if revised Δ m>0, illustrate solid-state PCM quality increase, phase-change material（PCM）Solidified；If revised Δ m<0, illustrate solid phase change material（PCM）Mass lost, phase-change material（PCM）Melted；If revised Δ m=0, illustrates solid phase change material（PCM）Quality do not change, phase-change material（PCM）It is in explicit neither endothermic nor exothermic state.

8. method for numerical simulation according to claim 7, it is characterised in that phase-change material in each time step is obtained by Δ m（PCM）Volume change Δ V, Δ V expression formula is：

ΔV=Δm·(1/ρ_s-1/ρ_l)

In formula：ρ_sRepresent solid phase change material（PCM）Density, ρ_lRepresent liquid state phase change material（PCM）Density；

The phase-change material obtained according to calculating（PCM）Volume change Δ V, to phase-change material（PCM）The cavity volume of container is adjusted：If Δ V>0, it is reduced by the volume shared by hole；If Δ V<0, correspondingly increase cavity volume.

9. method for numerical simulation according to claim 1, it is characterised in that phase-change material（PCM）Calculation of radiation heat transferring in hole be based on it is assumed hereinafter that：（1）All cavity surfaces are diffusing reflection grey body surface；（2）Phase-change material（PCM）The radiation of Surface absorption all wavelengths；（3）Hole inner vapor is not involved in radiation heat transfer；Radiation heat transfer between the interface of hole is concentrated mainly on radially, and the radiation heat transfer of axial and circumferential compares much weaker, to simplify calculating, negligible axial and circumferential radiation heat exchange.

10. method for numerical simulation according to claim 9, it is characterised in that hole effective heat transfer coefficient k_veffDetermine as the following formula：

$k_{\mathrm{veff}}=k_{\mathrm{PCM}}+k_{\mathrm{reff}}=k_{\mathrm{PCM}}+\ln \left(\frac{r_w}{r_v}\right)\&CenterDot;\frac{r_v\mathrm{\&sigma;}(T_w\left(t\right)+T_v\left(t\right))(T_w^2\left(t\right)+T_v^2\left(t\right))}{\frac{1}{{\mathrm{\&epsiv;}}_{\mathrm{PCM}}}+\frac{r_v}{r_w}(\frac{1}{{\mathrm{\&epsiv;}}_w}-1)},$

In formula：k_reffFor the equivalent equivalent heat conductivity of radiation heat transfer between the interface of hole, k_PCMFor phase-change material（PCM）Thermal conductivity factor, σ be Stefan-Boltzmann constant（σ=5.67×10^-8W·m^-2·K^-4）, r_wRepresent phase-change material（PCM）The external diameter of container outer wall, r_vRepresent phase-change material（PCM）The external diameter of hole interface, T_WFor outside wall temperature, T_VFor phase-change material（PCM）Temperature in hole interface, ε_PCMFor phase-change material（PCM）Emissivity, ε_WFor chamber wall emissivity, v represents hole interface.

Images