CN106021680A - Prediction method for transient working parameter of water sublimator during starting process - Google Patents

Abstract

The invention provides a prediction method for a transient working parameter of a water sublimator during a starting process. The method can predict a transient parameter of the water sublimator during the starting process, and provide the necessary technical support for grasping operating characteristics of the water sublimator and developing the water sublimator. The method mainly comprises the following steps of dividing the starting process of the water sublimator into three typical stages, namely an evaporation process within a water supply chamber, an evaporation process within a porous plate, and an alternating work process of evaporation and sublimation; according to the heat transfer characteristic within the water sublimator in the three stages, respectively constructing a lumped parameter model of variable mass temperature of each phase region of the water sublimator; determining a heat and mass transfer and rarefied gas flow differential equation with a mobile phase-change interface within a porous media; combined with the above model and equation, constructing a transient model of a water sublimation working process; and after solving the transient model of the water sublimator, the change rule of the temperature, feedwater flow, phase-change interface position and phase-change mass flow parameter during the starting process of the water sublimator can be predicted.

Description

Water sublimate device start-up course transient working parameter prediction method

Technical field

The present invention relates to a kind of Forecasting Methodology, be specifically related to water sublimate device start-up course transient working parameter prediction method, Belong to spacecraft thermal control technical field.

Background technology

Water sublimate device is a kind of by carrying out the thermal controls apparatus for spacecraft heat radiation to space drainage expendable medium, and it utilizes High vacuum environment residing for spacecraft, can make working-medium water first freeze, is more directly risen Huawei's water vapour by ice and be emitted into the outer space, Thus spacecraft used heat is taken away and realizes the hot type of spacecraft thermal control system and dissipate.Water sublimate utensil has simple in construction, volume little, heavy Measure light feature, be that there is the biggest power consumption or the higher preferable auxiliary radiating device of spacecraft of operating ambient temperature.Water liter The thermal control of China's device spacecraft thermal control system the most abroad and extravehicular activity unit (EMU) life support system (PLSS) obtains Many successes are applied, and the extravehicular space suit of China the most once successful Application Used For Eva takes water sublimate device.

Along with what China's lunar exploration activity and survey of deep space worked progressively carries out, water sublimate device thermal control technology has become China not Carry out one of thermal control means that are essential in survey of deep space spacecraft thermal control and that must break through.But research finds, although water liter China's device has obtained many successes application in engineering reality, but the existing design about water sublimate device and research are noted mostly Heavy is the stable state heat-sinking capability of water sublimate device, and less focuses on water sublimate device transient working characteristic, and especially water sublimate device opens (the whole work process of water sublimate device is all transient changing, simply the transient state mistake of start-up course in the transient response research of dynamic process Journey is more complicated).Owing to water sublimate device can only work under vacuum conditions, and in water sublimate device work process, its working medium experiences repeatedly Phase transformation and occurring in porous media, so in water sublimate device process of experimental, transient working mistake internal to water sublimate device The observation of journey is highly prone to observe the restriction with measurement means.

Therefore, in order to grasp the transient working characteristic of water sublimate device, the method for the employing numerical modeling work to water sublimate device Make process simulation analysis to be just particularly important.It is possible not only to disclose the microscopic workings process of water sublimate device, to water sublimate The service behaviour of device indicates, and can provide theories integration to the design of later water sublimate device, therefore has the heaviest The meaning wanted.

Summary of the invention

In view of this, the present invention provides a kind of water sublimate device start-up course transient working parameter prediction method, uses the party Method can predict the parameter such as temperature, transformation interface position during water sublimate device transient working, for grasping the work of water sublimate device Characteristic, carrying out water sublimate device development provides necessary technical guarantee.

Concretely comprising the following steps of the method:

Step one: water sublimate device start-up course is divided into three phases, the first stage: after feedwater, water is given at water sublimate device Evaporator section in water cavity；Second stage: water evaporator section in porous plate；Phase III: water freeze in porous plate after steaming Send out and distillation alternation section；

Step 2: the first stage of water sublimate device start-up course is carried out HEAT EXCHANGE ANALYSIS and numerical modeling and solves, it is thus achieved that The Changing Pattern of this stage running parameter；This running parameter includes: water sublimate device is to water cavity bottom-heated surface temperature, to water cavity gas District's temperature, to water cavity pool temperature, porous plate temperature, feed-water quality flow, transformation interface position；

When judging the condition of second stage that water sublimate device work process meets setting, water sublimate device proceeds to second-order Section, enters step 3；

Step 3: the second stage of water sublimate device start-up course is carried out HEAT EXCHANGE ANALYSIS and numerical modeling and solves, it is thus achieved that The Changing Pattern of this stage running parameter；This running parameter includes: water sublimate device is to water cavity bottom-heated surface temperature, to water cavity temperature Degree, porous plate pool temperature, porous plate gas district temperature, feed-water quality flow, transformation interface position；

When judging the condition of phase III that water sublimate device work process meets setting, water sublimate device proceeds to the 3rd rank Section, enters step 4；

Step 4: the phase III of water sublimate device start-up course is carried out HEAT EXCHANGE ANALYSIS and numerical modeling and solves, it is thus achieved that The Changing Pattern of this stage running parameter；This running parameter includes: water sublimate device is to water cavity bottom-heated surface temperature, to water cavity temperature Degree, porous plate pool temperature, porous plate ice formation temperature, porous plate gas district temperature, sublimated mass flow, ice-water transformation interface position Put, distil interface location；

Step 5: when the ice sheet in water sublimate device surface porous metal plate disappears, when porous plate ice formation thickness is zero, water sublimate device Work process proceeds to the evaporation process in porous plate, is i.e. back to above-mentioned steps three；

Step 6: after reaching time of end of setting or temperature conditions, the running parameter that each stage is obtained according to Time series collects, and obtains its rule over time.

The above three stage carries out transient state phase transformation work process the most successively solve, particularly as follows: (1) builds water sublimate device The variable mass temperature lumped parameter model of each phase region；(2) have in determining porous media mobile transformation interface heat and mass and The rarefied gaseous flow differential equation；(3) combine above-mentioned model and equation, constitute water sublimate work process transient model；(4) use The water sublimate work process transient model constituted is solved by numerical solution, it is thus achieved that temperature in water sublimate device start-up course, Feedwater flow, the Changing Pattern of transformation interface position.

In step 2, the first stage to water sublimate device start-up course carries out HEAT EXCHANGE ANALYSIS and numerical modeling the tool solved Body process is:

At this stage water sublimate device heating surface, to water cavity pool, to the thermally conductive relation difference of water cavity gas district and porous plate For:

Heating surface:

To water cavity pool:

To water cavity gas district:

Porous plate:

Wherein, c_i、M_i、T_iBeing respectively the thermal capacitance of region i, quality and temperature, i=k, we, wm, p, kwe, wem, wmp are respectively Represent heating surface, to water cavity pool, to water cavity gas district, porous plate, heating surface with to interface, water cavity pool, to water cavity pool with give Water cavity gas regional boundary face, to water cavity gas district and porous plate interface, Q₀Heat is added for heating surface；R_kweFor heating surface and to water cavity water Thermal resistance between district；For evaporation of water mass flow；h_eFor evaporation of water latent heat；R_wemFor to water cavity pool and to water cavity gas Thermal resistance between district；R_wmpFor to the thermal resistance between water cavity gas district and porous plate；

In formula (1)-(4), the thermal resistance of water sublimate Qi Ge district transformation interface is obtained by following formula (5):

$R_j=\frac{L_j}{k_{ef,j}{A}_j}---\left(5\right)$

Wherein K_{Ef, j}, L_j,A_jRespectively Equivalent Thermal Conductivities, thickness and sectional area, j=kwe, wem, wmp,

When in porous media being water: k_ef,i=ε k_w+(1-ε)k_m

When in porous media being ice: k_ef,i=ε k_i+(1-ε)k_m (6)

Wherein k_m、k_wAnd k_iIt is respectively porous media solid skeletal, water and the heat conductivity of ice；ε is porous media hole Rate；

The mass change amount of feedwater intracavity water is:

Wherein:For the mass flow of the intracavity water that feeds water, m_wFor the quality of the intracavity water that feeds water, t is the time；

The position relation over time of feedwater intracavity water evaporation interface is:

Wherein, δ_eT () is the position of feedwater intracavity water evaporation interface, ρ_vFor the density of the intracavity steam that feeds water, A_rFor feedwater The equivalent disengagement area in chamber；

Evaporation of water mass flowObtained by following formula:

$\frac{{\stackrel{\·}{m}}_v}{{\mathrm{nA}}_r(P_{sat}-P_o)}=(1-F(\frac{r}{\mathrm{\λ}}\left)\right)(A\stackrel{\‾}{p}+B)+F\left(\frac{r}{\mathrm{\λ}}\right)C---\left(11\right)$

In formula, n is the number in hole on porous plate, A_p For the area in the single hole of porous plate, r is porous plate average pore radius, μ_vFor the dynamic viscosity of steam, T_eFor evaporation interface temperature Degree, δ_pFor porous plate thickness, m_mFor the quality of water vapour molecule,For stream molecule average speed, P_oFor ambient pressure；P_sat Saturated vapor pressure for evaporating surface；

After completing above-mentioned numerical modeling, the running parameter in this stage is solved, particularly as follows: give fixed pattern (1)-formula (4), The initial value of differential variable in formula (7), formula (8), the thermal resistance in each region, Equivalent Thermal Conductivities can obtain according to formula (5), formula (6), Evaporative mass flow can obtain according to formula (11), then use variable step implicit expression Runge-Kutta method solve formula (1)-formula (4), The differential equation shown in formula (7), formula (8), is simulated iterative computation to the evaporation process of feedwater intracavity；

The condition of set entrance second stage is: if T_w＞ 0 and δ_e(t) ＜ δ_w, δ_wFor to water cavity height；Then recognize For feedwater underfill to water cavity, it is unsatisfactory for entering the condition of second stage；If T_w＞ 0 and δ_e(t)=δ_w, then it is assumed that feedwater Being full of to water cavity, water sublimate device enters the second stage of start-up course.

Step 3 the second stage of water sublimate device start-up course is carried out HEAT EXCHANGE ANALYSIS and numerical modeling and solve concrete Process is:

In this stage, water moves in porous plate and evaporates, and the equation of momentum is:

$F_c+F_p-F_{sat}-F_g-F_f=\frac{d\left(mu\right)}{dt}={\mathrm{\ρ}}_w\frac{d\left({\mathrm{\δ}}_e\right(t\left)A_r\mathrm{\ϵ}u\right)}{dt}---\left(13\right)$

That is:

Wherein, F_cFor capillary force, F_pFor feed pressure, F_satFor the saturated vapor pressure of evaporating surface, F_gFor liquid in porous plate The gravity of body, F_fFor the liquid pressure loss by certain thickness porous material, δ_eT () is the distance that water enters porous plate, A_r For the sectional area of porous plate, ε is the porosity of material, and u is water flowing velocity in porous plate, ρ_wFor the density of water, m is many The quality of water in orifice plate；P_cFor capillary pressure, P_pFor feedwater pressure, P_satFor the saturated vapor pressure of evaporating surface, P_fFor liquid Pressure loss by certain thickness porous material；

Evaporation interface position δ in porous plate_eT () is determined by following formula (17):

$\frac{2\mathrm{\σ}\cos \mathrm{\θ}}{r}+P_p-P_{sat}-\frac{\mathrm{\μ}{\stackrel{\·}{m}}_v^{\′}{\mathrm{\δ}}_e\left(t\right)}{{\mathrm{K\ρ}}_w{A}_r}-\frac{{\mathrm{\μ\ϵ\δ}}_e\left(t\right)}{K}\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}={\mathrm{\ρ}}_w\frac{d}{dt}\left({\mathrm{\δ}}_e\right(t)\frac{{\stackrel{\·}{m}}_v^{\′}}{{\mathrm{\ρ}}_w{A}_r\mathrm{\ϵ}}+{\mathrm{\δ}}_e(t\left)\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}\right)---\left(17\right)$

In formula (17), steam mass flowFor:

Wherein, σ is the surface tension of water, and θ is Liquid contact angle, and μ is the dynamic viscosity of water vapour, and K is oozing of porous plate Permeability coefficient,For the average pressure of steam, d_p For the average pore size of porous plate, η is the dynamic viscosity of water vapour, and k is Boltzmann constant,

Water enters after porous plate, is all filled by water to water cavity, now water sublimate device heating surface and give the temperature of water cavity can Obtained by the transient temperature equation shown in formula (19), formula (20):

$c_k{M}_k{\stackrel{\·}{T}}_k=Q_0-(T_k-T_w^{\′})/R_{kw}^{\′}---\left(19\right)$

$c_w^{\′}{M}_w^{\′}{\stackrel{\·}{T}}_w^{\′}=(T_k-T_w^{\′})/R_{kw}^{\′}-(T_w^{\′}-T_{pw})/R_{wp}^{\′}---\left(20\right)$

Be divided into by porous plate containing water section and without water section, then porous plate is containing water section and without water section temperature equation It is respectively as follows:

$c_{pw}{M}_{pw}{\stackrel{\·}{T}}_{pw}=(T_w^{\′}-T_{pw})/R_{wp}-(T_{pw}-T_{pm})/R_{pm}-{\stackrel{\·}{m}}_v{h}_e---\left(21\right)$

$c_{pm}{M}_{pm}{\stackrel{\·}{T}}_{pm}=(T_{pw}-T_{pm})/R_{pm}---\left(22\right)$

After completing above-mentioned numerical modeling, the running parameter in this stage is solved, particularly as follows: use variable step implicit expression Runge-Kutta method solves the differential equation shown in formula (17), (19)-formula (22), is simulated the evaporation process in porous plate Iterative computation；Wherein, the thermal resistance in each region, Equivalent Thermal Conductivities can obtain according to method shown in formula (5), formula (6), evaporative mass Flow can be obtained by formula (18)；

The condition of set entrance second stage is: if T_w＞ 0, T_pw＜ 0 and δ_e(t) ＜ δ_w+δ_p, then feedwater is judged Freezing in water sublimate device porous plate, water sublimate device enters the phase III of start-up course.

In step 4, the phase III to water sublimate device start-up course carries out HEAT EXCHANGE ANALYSIS and numerical modeling the tool solved Body step is:

At this stage water sublimate device heating surface, to the heat biography in water cavity, porous plate pool, porous plate ice formation and porous plate gas district Lead relation:

$c_k{M}_k{\stackrel{\·}{T}}_k=Q_0-(T_k-T_w^{\′})/R_{kw}^{\′}---\left(23\right)$

$c_w^{\′}{M}_w^{\′}{\stackrel{\·}{T}}_w^{\′}=(T_k-T_w^{\′})/R_{kw}^{\′}-(T_w^{\′}-T_{pw})/R_{wp}^{\′}---\left(24\right)$

$c_{pw}{M}_{pw}{\stackrel{\·}{T}}_{pw}=(T_w-T_{pw})/R_{wpw}-(T_{pw}-T_{pi})/R_{pwi}-{\stackrel{\·}{m}}_{wi}{h}_f---\left(25\right)$

$c_{pi}{M}_{pi}{\stackrel{\·}{T}}_{pi}=(T_{pw}-T_{pi})/R_{pwi}-(T_{pi}-T_{pm})/R_{pim}-{\stackrel{\·}{m}}_s{h}_s---\left(26\right)$

$c_{pm}{M}_{pm}{\stackrel{\·}{T}}_{pm}=(T_{pi}-T_{pm})/R_{pim}---\left(27\right)$

Wherein,For the thawing/solidifying phase variable Rate at ice-water interface,Rate of sublimation for the interface that distils；

Can be determined by formula (28):

Thawing or freezing rate at frozen water interface are then determined by formula (29):

$\left|\frac{{\mathrm{dm}}_{wi}}{dt}\right|=\frac{Q_{pw}-Q_{wi}}{h_f}=\left|{\mathrm{\ρ}}_w{\mathrm{\ϵA}}_r\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}\right|---\left(29\right)$

Wherein, Q_pwFor by the heat to water cavity conduction to aqueous porous plate, Q_wiFor the heat by ice sheet dissipation to space outerpace Amount,

After completing above-mentioned numerical modeling, the running parameter in this stage is solved, particularly as follows: use variable step implicit expression Runge-Kutta method solves the differential equation shown in formula (23)-formula (27), and the sublimation process in porous plate is simulated iteration Calculate.

Beneficial effect:

(1) present invention breaches water sublimate device Numerical Study on Transient technology, meets in water sublimate device development process Carrying out the demand of design and Analysis of Working Performance, the development for water sublimate device provides requisite technical guarantee.

(2) water sublimate device transient state modeling method of the present invention, can start and in work process water sublimate device Heat exchange and phase transition process carry out numerical simulation calculation, provide convenience for understanding water sublimate device work process in detail.

(3) method of the present invention is used to be possible not only to obtain the temperature of each several part in water sublimate device start-up course and to water quality Amount flow change, it is also possible to obtain by experiment can not measure the parameter obtained such as transformation interface change in location and evaporation/liter China's gas mass flow change, research and application for water sublimate device provide technical guarantee.

Accompanying drawing explanation

Fig. 1 water sublimate of the present invention device work process divides schematic diagram；

Fig. 2 result of calculation of the present invention and comparison of test results.

Detailed description of the invention

Below in conjunction with the accompanying drawings and embodiment, the present invention is described in further detail.

The present embodiment proposes a kind of water sublimate device start-up course transient working parameter prediction method, it is possible to solve water sublimate device In research process, its micro physical model work process and transient response can not utilize existing steady-state behaviour analysis method to carry out indicating, also The technical barrier that laboratory facilities are observed cannot be used.

First the startup work process of water sublimate device is divided into three representative continuous print working stages carry out Analyze: after wherein first stage is feedwater, water is at the evaporator section of water sublimate device feedwater intracavity, and second stage is that water is in porous Evaporator section in plate, three phases is the evaporation after water freezes in porous plate and distillation alternation section.

To water sublimate device start work process modeling analysis need respectively the above three stage is carried out HEAT EXCHANGE ANALYSIS with Mathematical description.Described HEAT EXCHANGE ANALYSIS is set up on the basis of following several hypothesis with mathematical description:

A () is due to less to water cavity height, therefore temporarily ignore the inhomogeneities of water distribution in water sublimate device；

B () enters the flash distillation effect of water sublimate device initial time due to water, water sublimate device rapid drawdown occurs to the temperature of water cavity, After reaching triple point of water, will quickly freeze；

C () is shorter due to the water sublimate device start-up course time, therefore ignore water and enter to steam in evaporation process after water cavity To the heat convection in water cavity；

D () ignores the radiation heat transfer between water sublimate device and surrounding.

Based on this, water sublimate device start-up course carries out transient state modeling with the step analyzed is:

Step one: first stage (water is at the evaporator section of water sublimate device feedwater intracavity) is carried out HEAT EXCHANGE ANALYSIS and numerical value is built Mould, and determine that its boundary condition solves, obtain the running parameter in this stage, particularly as follows:

Water is i.e. exposed under vacuum after entering water sublimate device, and water surface pressure drops suddenly to far below its initial temperature pair The saturation pressure answered, the water of feedwater intracavity is become superheat state by initial steady statue, makes water that moment rapid evaporation, water to occur And reduce rapidly to the temperature of water cavity.Owing in evaporation process, the latent heat of water release changes far above the convection current of water with feedwater intracavity Heat, therefore ignore the heat convection of feedwater intracavity temporarily, water sublimate device heating surface, to water cavity pool, to water cavity gas district and porous plate Thermally conductive relation is respectively such as formula (1)-formula (4):

Heating surface:

$c_k{M}_k{\stackrel{\·}{T}}_k=Q_0-(T_k-T_{we})/R_{kwe}---\left(1\right)$

To water cavity pool:

$c_{we}{M}_{we}{\stackrel{\·}{T}}_{we}=(T_k-T_{we})/R_{kwe}-(T_{we}-T_{wm})/R_{wem}-{\stackrel{\·}{m}}_v{h}_e---\left(2\right)$

To water cavity gas district:

$c_{wm}{M}_{wm}{\stackrel{\·}{T}}_{wm}=(T_{we}-T_{wm})/R_{wem}-(T_{wm}-T_p)/R_{wmp}---\left(3\right)$

Porous plate:

$c_p{M}_p{\stackrel{\·}{T}}_p=(T_{wm}-T_p)/R_{wmp}---\left(4\right)$

Wherein, c_i、M_i、T_iBeing respectively the thermal capacitance of region i, quality and temperature, i=k, we, wm, p, kwe, wem, wmp are respectively Represent heating surface, to water cavity pool, to water cavity gas district, porous plate, heating surface with to interface, water cavity pool, to water cavity pool with give Water cavity gas regional boundary face, to water cavity gas district and porous plate interface, Q₀Heat is added for heating surface；R_kweFor heating surface and to water cavity water Thermal resistance between district；For evaporation of water mass flow；h_eFor evaporation of water latent heat；R_wemFor to water cavity pool and to water cavity gas Thermal resistance between district；R_wmpFor to the thermal resistance between water cavity gas district and porous plate.

Wherein, the thermal resistance in water sublimate Qi Ge district changes, water sublimate Qi Ge district with the movement of transformation interface in work process The thermal resistance of transformation interface is obtained by following formula (5):

K in formula_ef,i、L_i、A_iIt is respectively Equivalent Thermal Conductivities, thickness and sectional area, i=kwe, wem, wmp.For porous Medium, its Equivalent Thermal Conductivities should consider its solid skeletal and internal working medium, then simultaneously:

Wherein k_m、k_wAnd k_iBeing respectively porous media solid skeletal, water and the heat conductivity of ice, ε is porous media hole Rate.

The mass change amount of feedwater intracavity water is:

$\frac{{\mathrm{dm}}_w\left(t\right)}{dt}={\stackrel{\·}{m}}_g-{\stackrel{\·}{m}}_v---\left(7\right)$

In formula (7)For the mass flow of the intracavity water that feeds water, m_wFor the quality of the intracavity water that feeds water, t is the time.

The position relation over time of feedwater intracavity water evaporation interface is:

$\frac{{\mathrm{d\δ}}_e\left(t\right)}{dt}=\frac{1}{{\mathrm{\ρ}}_v{A}_r}\frac{{\mathrm{dm}}_w\left(t\right)}{dt}---\left(8\right)$

Wherein, δ_eT () is the position of feedwater intracavity water evaporation interface, ρ_vFor the density of the intracavity steam that feeds water, A_rFor feedwater The equivalent disengagement area in chamber.

Gas flowing in minute yardstick and nanoscale systems, generally uses Knudsen number (K_n) carry out flow mechanism Judge.K_nDefinition be: the ratio of characteristic size Λ of fluid molecule mean free path λ and system:

$K_n=\frac{\mathrm{\λ}}{\mathrm{\Λ}}---\left(9\right)$

$\mathrm{\λ}=\frac{{\mathrm{kT}}_v}{\sqrt{2}{\mathrm{\πd}}_v^2{P}_{sat}}---\left(10\right)$

Wherein: k is Boltzman constant, T_vFor gas temperature, P_satFor saturation vapour pressure, d_vFor gas molecule diameter. Gas by the relation of the specific mass flow of capillary tube with pressure and temperature is:

$\frac{{\stackrel{\·}{m}}_v}{{\mathrm{nA}}_r(P_{sat}-P_o)}=(1-F(\frac{r}{\mathrm{\λ}}\left)\right)(A\stackrel{\‾}{p}+B)+F\left(\frac{r}{\mathrm{\λ}}\right)C---\left(11\right)$

In formula, n is the number in hole on porous plate hole, A_pFor the area in the single hole of porous plate, r is porous plate average pore size, μ_vFor the dynamic viscosity of steam, T_eFor evaporation interface temperature Degree, δ_pFor porous plate thickness, m_mFor the quality of water vapour molecule,For stream molecule average speed, P_oFor ambient pressure；

P_satFor the saturated vapor pressure of evaporating surface, T is evaporating temperature:

log₁₀(P_sat/10³)=8.42926609-1.82717843 (10³/(T+273.15))-0.071208271(10³/ (T+273.15))² (12)

After completing above-mentioned numerical modeling, the running parameter in this stage is solved, particularly as follows: give fixed pattern (1)-formula (4), The initial value of differential variable in formula (7), formula (8), the thermal resistance in the most each region of remaining parameter, Equivalent Thermal Conductivities can according to formula (5), Formula (6) obtains, and evaporative mass flow can obtain according to formula (9), (10), (11), (12), then uses variable step implicit expression Runge- Kutta method solves the differential equation shown in formula (1)-formula (4), formula (7), formula (8), is simulated the evaporation process of feedwater intracavity repeatedly Generation calculate, it is possible to obtain water sublimate device start-up course is in feedwater intracavity evaporation process in water cavity bottom-heated face temperature Degree T_k, to water cavity gas district temperature T_wm, to water cavity pool temperature T_we, porous plate temperature T_p, feed-water quality flowTransformation interface Position δ_e(t) isoparametric Changing Pattern.

During above-mentioned calculating, if T_we> 0 and δ_e(t) ＜ δ_w(δ_wFor to water cavity height), then it is assumed that feedwater underfill To water cavity, being unsatisfactory for entering the condition of second stage, feedwater goes successively to, to water cavity, continue solving, directly of differential equation group To T_we≤0；If T_we≤ 0 and δ_e(t)=δ_w, then it is assumed that feedwater is full of to water cavity, meets the condition entering second stage, enters Enter following according to step 2.

Step 2: water sublimate device proceeds to second active section (feedwater evaporator section in porous plate), and to second rank Duan Jinhang HEAT EXCHANGE ANALYSIS, numerical modeling, determine boundary condition and solve, obtaining the running parameter in this stage.Particularly as follows:

Due to water cavity height δ_wLess, water generally has little time to freeze in being full of to the short period of water cavity, therefore, gives Water in water cavity often continues evaporation after entering porous plate, is finally reached freezing temperature and freezes.In the process, feedwater Moving in porous plate and evaporate, being analyzed liquid stress relation in moving process in porous plate can be such as formula (13) The shown equation of momentum:

$F_c+F_p-F_{sat}-F_g-F_f=\frac{d\left(mu\right)}{dt}={\mathrm{\ρ}}_w\frac{d\left({\mathrm{\δ}}_e\right(t\left)A_r\mathrm{\ϵ}u\right)}{dt},---\left(13\right)$

In formula, F_cFor capillary force, F_pFor feed pressure, F_satFor the saturated vapor pressure of evaporating surface, F_gFor liquid in porous plate The gravity of body, F_fFor the liquid pressure loss by certain thickness porous material, δ_eT () is the position of water evaporation interface, A_rFor The sectional area of porous plate, ε is the porosity of material, and u is water flowing velocity in porous plate, ρ_wFor the density of water, m is porous The quality of water in plate.

Formula (13) is arranged, can obtain:

$(P_c+P_p-P_{sat}-P_f)-F_g/\left({\mathrm{\ϵA}}_r\right)={\mathrm{\ρ}}_w\frac{d\left({\mathrm{\δ}}_e\right(t\left)u\right)}{dt}---\left(14\right)$

In formula, P_cFor capillary pressure, P_pFor feedwater pressure, P_satFor the saturated vapor pressure of evaporating surface, P_fLead to for liquid Cross the pressure loss of certain thickness porous material.

It is d for average pore size_pPorous plate, pressure produced by surface tension component on liquid moving direction For:R is porous plate average pore radius,θ is Liquid contact angle, and σ is the surface tension of water.Water is many When moving in orifice plate, pressure produced by the viscosity of water and the structure of porous plate is reduced to: For water Mass flow, K is porous plate permeability coefficient, under the conditions of space microgravity, F_g=0.

At evaporation interface, can obtain according to the seriality at interface:

${\mathrm{\ρ}}_w{A}_r\mathrm{\ϵ}(u-\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt})={\stackrel{\·}{m}}_v^{\′}---\left(15\right)$

In formula,For the translational speed of evaporation interface,For steam mass flow now；Can be obtained by formula (15):

$u=\frac{{\stackrel{\·}{m}}_v^{\′}}{{\mathrm{\ρ}}_w{A}_r\mathrm{\ϵ}}+\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}---\left(16\right)$

Formula (16) substitution formula (14) can be obtained the equation about evaporation interface position in porous plate, such as formula (17):

$\frac{2\mathrm{\σ}cos\mathrm{\θ}}{r}+P_p-P_{sat}-\frac{\mathrm{\μ}{\stackrel{\·}{m}}_v^{\′}{\mathrm{\δ}}_e\left(t\right)}{{\mathrm{K\ρ}}_w{A}_r}-\frac{{\mathrm{\μ\ϵ\δ}}_e\left(t\right)}{K}\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}={\mathrm{\ρ}}_w\frac{d}{dt}\left({\mathrm{\δ}}_e\right(t)\frac{{\stackrel{\·}{m}}_v^{\′}}{{\mathrm{\ρ}}_w{A}_r\mathrm{\ϵ}}+{\mathrm{\δ}}_e(t\left)\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}\right)---\left(17\right)$

The evaporation interface position in porous plate is can determine that by formula (17), in formula (17),Determined by following formula (18):

$\frac{{\stackrel{\·}{m}}_v^{\′}}{{\mathrm{nA}}_r(P_{sat}-P_o)}=(1-F(\frac{r}{\mathrm{\λ}}\left)\right)(A^{\′}\stackrel{\‾}{p}+B^{\′})+F\left(\frac{r}{\mathrm{\λ}}\right)C^{\′}---\left(18\right)$

Wherein n is the number of porous plate hole,For the average pressure of steam,

Water enters after porous plate, owing to all being filled by water to water cavity, therefore water sublimate device heating surface and give the temperature of water cavity Can be obtained by the transient temperature equation shown in formula (19), formula (20):

$c_k{M}_k{\stackrel{\·}{T}}_k=Q_0-(T_k-T_w^{\′})/R_{kw}^{\′}---\left(19\right)$

$c_w^{\′}{M}_w^{\′}{\stackrel{\·}{T}}_w^{\′}=(T_k-T_w^{\′})/R_{kw}^{\′}-(T_w^{\′}-T_{pw})/R_{wp}^{\′},---\left(20\right)$

Be divided into by porous plate containing water section and without water section, then porous plate pool is with porous plate gas district temperature equation respectively For:

$c_{pw}{M}_{pw}{\stackrel{\·}{T}}_{pw}=(T_w^{\′}-T_{pw})/R_{wp}-(T_{pw}-T_{pm})/R_{pm}-{\stackrel{\·}{m}}_v{h}_e,---\left(21\right)$

$c_{pm}{M}_{pm}{\stackrel{\·}{T}}_{pm}=(T_{pw}-T_{pm})/R_{pm}---\left(22\right)$

Wherein, subscript pw, pm refers to gentle district, porous plate pool, R respectively_wpRefer to the interface resistance in porous plate gas district and pool, After completing above-mentioned numerical modeling, the running parameter in this stage is solved, particularly as follows: use variable step implicit expression Runge- Kutta method solves the differential equation shown in formula (17), (19)-formula (22), and the evaporation process in porous plate is simulated iteration meter Calculate.Wherein, the thermal resistance in each region, Equivalent Thermal Conductivities can obtain according to method shown in formula (5), formula (6), and evaporative mass flow can Obtained by formula (18), face temperatures, feedwater cavity temperature, porous plate pool temperature, porous plate gas district temperature, evaporative mass flow Initial value is obtained by the analog value obtained in step one respectively.Can be obtained in water sublimate device start-up course by above-mentioned calculating and be in In evaporation process in porous plate to water cavity bottom surface temperature, feedwater cavity temperature, porous plate pool temperature, porous plate gas district temperature Degree, feed-water quality flow, the isoparametric Changing Pattern in transformation interface position.

During above-mentioned calculating, if T_w＞ 0, T_pw＜ 0 and δ_e(t) ＜ δ_w+δ_p(T_pwRefer to porous plate pool temperature), then Judging that feedwater freezes in water sublimate device porous plate, water sublimate device work process proceeds to the sublimation process in porous plate, under entrance State step 3.

Step 3: water sublimate device proceeds to the 3rd active section, and (evaporation after feedwater freezes in porous plate replaces with distillation Active section), three phases is carried out HEAT EXCHANGE ANALYSIS, numerical modeling, determines boundary condition and solve, obtain the work in this stage Parameter.Particularly as follows:

In water moving process in porous plate, under evaporation, in porous plate, the temperature of water is gradually lowered.Work as porous When plate pool temperature is down under 0 DEG C, it is believed that the water in aqueous porous plate will quickly freeze.Ice sheet will be opened towards the side of vacuum Begin to distil, and will occur to freeze or melt according to the heat exchange equilibrium relation between water and ice at ice-water interface.For porous plate Interior evaporation-distillation alternation section, is reduced to the temperature as shown in Fig. 1 (b) by interstructural for water sublimate device heat transmission relation Degree lumped parameter model.Thus can set up to water sublimate device heating surface in sublimation process in porous plate, to water cavity, perforated-plate-water District, porous plate ice formation and the thermally conductive relation in porous plate gas district:

$c_k{M}_k{\stackrel{\·}{T}}_k=Q_0-(T_k-T_w^{\′})/R_{kw}^{\′}---\left(23\right)$

$c_w^{\′}{M}_w^{\′}{\stackrel{\·}{T}}_w^{\′}=(T_k-T_w^{\′})/R_{kw}^{\′}-(T_w^{\′}-T_{pw})/R_{wp}^{\′}---\left(24\right)$

$c_{pw}{M}_{pw}{\stackrel{\·}{T}}_{pw}=(T_w-T_{pw})/R_{wpw}-(T_{pw}-T_{pi})/R_{pwi}-{\stackrel{\·}{m}}_{wi}{h}_f---\left(25\right)$

$c_{pi}{M}_{pi}{\stackrel{\·}{T}}_{pi}=(T_{pw}-T_{pi})/R_{pwi}-(T_{pi}-T_{pm})/R_{pim}-{\stackrel{\·}{m}}_s{h}_s---\left(26\right)$

$c_{pm}{M}_{pm}{\stackrel{\·}{T}}_{pm}=(T_{pi}-T_{pm})/R_{pim}---\left(27\right)$

Wherein,For the thawing/solidifying phase variable Rate at ice-water interface,Rate of sublimation for the interface that distils.

Can be determined by formula (28):

${\stackrel{\·}{m}}_s={\mathrm{nA}}_r\frac{d_p^3}{3({\mathrm{\δ}}_p-{\mathrm{\δ}}_i(t)-{\mathrm{\δ}}_e(t\left)\right)}\sqrt{\frac{\mathrm{\π}m}{2{\mathrm{kT}}_{pi}}}(P_{sat}-P_o),---\left(28\right)$

Wherein, δ_iT () is porous plate ice formation thickness,

Thawing or freezing rate at frozen water interface are then determined by formula (29):

$\left|\frac{{\mathrm{dm}}_{wi}}{dt}\right|=\frac{Q_{pw}-Q_{wi}}{h_f}=\left|{\mathrm{\ρ}}_w{\mathrm{\ϵA}}_r\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}\right|,---\left(29\right)$

Wherein, Q_pwFor by the heat to water cavity conduction to aqueous porous plate, Q_wiFor the heat by ice sheet dissipation to space outerpace Amount, works as Q_pw＜ Q_wiTime, the water in porous plate continues to be frozen into ice；Work as Q_pw＞ Q_wiTime, the ice in porous plate melts at ice-water interface Turn to water.Therefore, the position δ of ice-water interface_eT () changes the most therewith.When ice sheet in distillation and melts or the effect of solidification Under when becoming 0, there is evaporation and move forward in feedwater in porous plate, starts the next distillation cycle until again freezing.

After completing above-mentioned numerical modeling, variable step implicit expression Runge-Kutta method is used to solve shown in formula (23)-formula (27) The differential equation, is simulated iterative computation to the sublimation process in porous plate.Wherein, the thermal resistance in each region, Equivalent Thermal Conductivities Can obtain according to method shown in formula (5), formula (6), ice sublimated mass flow can be obtained by formula (28), and ice-water interface location can join Close differential equation (29) to obtain, face temperatures, feedwater cavity temperature, porous plate pool temperature, porous plate ice formation, porous plate gas District's temperature, ice-water interface location initial value are obtained by the analog value obtained in step 2 respectively.Water can be obtained by above-mentioned calculating Sublimating apparatus work process is in the sublimation process in porous plate to water cavity bottom surface temperature, feedwater cavity temperature, perforated-plate-water District's temperature, porous plate ice formation temperature, porous plate gas district temperature, sublimated mass flow, ice-water transformation interface position, distillation interface The isoparametric Changing Pattern in position.

During calculating, work as T_w＞ 0, T_pw＞ 0, δ_e(t) ＜ δ_w+δ_p,δ_iT, during ()=0, the ice sheet in surface porous metal plate disappears Losing, porous plate ice formation thickness is zero, and water sublimate device work process proceeds to the evaporation process in porous plate, proceeds to above-mentioned steps two and opens Sublimation process in beginning porous plate calculates.

Step 4: after meeting time or the temperature conditions of the end of calculating and setting, each stage is obtained to water cavity bottom surface Temperature, feedwater cavity temperature, porous plate pool temperature, porous plate ice formation temperature, porous plate gas district temperature, sublimated mass flow, ice- Aqueous phase becomes the parameter value of calculation such as interface location, distillation interface location and collects according to time series, obtains its change in time Law.

In sum, these are only presently preferred embodiments of the present invention, be not intended to limit protection scope of the present invention. All within the spirit and principles in the present invention, any modification, equivalent substitution and improvement etc. made, should be included in the present invention's Within protection domain.

Claims (4)

1. water sublimate device start-up course transient working parameter prediction method, it is characterised in that

Step one: water sublimate device start-up course is divided into three phases, the first stage: after feedwater water at water sublimate device to water cavity Interior evaporator section；Second stage: water evaporator section in porous plate；Phase III: water freeze in porous plate after evaporation with Distillation alternation section；

Step 2: the first stage of water sublimate device start-up course is carried out HEAT EXCHANGE ANALYSIS and numerical modeling and solves, it is thus achieved that these rank The Changing Pattern of section running parameter；This running parameter includes: water sublimate device is to water cavity bottom-heated surface temperature, to water cavity gas district temperature Spend, to water cavity pool temperature, porous plate temperature, feed-water quality flow, transformation interface position；

When judging the condition of second stage that water sublimate device work process meets setting, water sublimate device proceeds to second stage, enters Enter step 3；

Step 3: the second stage of water sublimate device start-up course is carried out HEAT EXCHANGE ANALYSIS and numerical modeling and solves, it is thus achieved that these rank The Changing Pattern of section running parameter；This running parameter includes: water sublimate device to water cavity bottom-heated surface temperature, feedwater cavity temperature, Porous plate pool temperature, porous plate gas district temperature, feed-water quality flow, transformation interface position；

When judging the condition of phase III that water sublimate device work process meets setting, water sublimate device proceeds to the phase III, enters Enter step 4；

Step 4: the phase III of water sublimate device start-up course is carried out HEAT EXCHANGE ANALYSIS and numerical modeling and solves, it is thus achieved that these rank The Changing Pattern of section running parameter；This running parameter includes: water sublimate device to water cavity bottom-heated surface temperature, feedwater cavity temperature, Porous plate pool temperature, porous plate ice formation temperature, porous plate gas district temperature, sublimated mass flow, ice-water transformation interface position, Distillation interface location；

Step 5: when the ice sheet in water sublimate device surface porous metal plate disappears, when porous plate ice formation thickness is zero, water sublimate device works Process proceeds to the evaporation process in porous plate, is i.e. back to above-mentioned steps three；

Step 6: after reaching the time of end of setting or temperature conditions, the running parameter obtaining each stage is according to the time Sequence collects, and obtains its rule over time.

The above three stage carries out transient state phase transformation work process the most successively solve, particularly as follows: (1) builds each phase of water sublimate device The variable mass temperature lumped parameter model in district；(2) there is in determining porous media the heat and mass of mobile transformation interface and thin The gas flowing differential equation；(3) combine above-mentioned model and equation, constitute water sublimate work process transient model；(4) numerical value is used The water sublimate work process transient model constituted is solved by solution, it is thus achieved that temperature, feedwater in water sublimate device start-up course Flow, the Changing Pattern of transformation interface position.

Water sublimate device start-up course transient working parameter prediction method the most according to claim 1, it is characterised in that: step In two, first stage to water sublimate device start-up course carries out HEAT EXCHANGE ANALYSIS and numerical modeling and the detailed process that solves is:

It is respectively as follows: at this stage water sublimate device heating surface, thermally conductive relation to water cavity pool, to water cavity gas district and porous plate

Heating surface:

To water cavity pool:

To water cavity gas district:

Porous plate:

Wherein, c_i、M_i、T_iBeing respectively the thermal capacitance of region i, quality and temperature, i=k, we, wm, p, kwe, wem, wmp represent respectively Heating surface, to water cavity pool, to water cavity gas district, porous plate, heating surface with to interface, water cavity pool, to water cavity pool with to water cavity Gas regional boundary face, to water cavity gas district and porous plate interface, Q₀Heat is added for heating surface；R_kweFor heating surface with to water cavity pool it Between thermal resistance；For evaporation of water mass flow；h_eFor evaporation of water latent heat；R_wemFor to water cavity pool with to water cavity gas district it Between thermal resistance；R_wmpFor to the thermal resistance between water cavity gas district and porous plate；

In formula (1)-(4), the thermal resistance of water sublimate Qi Ge district transformation interface is obtained by following formula (5):

$R_j=\frac{L_j}{k_{ef,j}{A}_j}---\left(5\right)$

Wherein K_{Ef, j}, L_j,A_jRespectively Equivalent Thermal Conductivities, thickness and sectional area, j=kwe, wem, wmp,

When in porous media being water: k_ef,i=ε k_w+(1-ε)k_m

When in porous media being ice: k_ef,i=ε k_i+(1-ε)k_m (6)

Wherein k_m、k_wAnd k_iIt is respectively porous media solid skeletal, water and the heat conductivity of ice；ε is porosity of porous medium；Give In water cavity, the mass change amount of water is:

Wherein:For the mass flow of the intracavity water that feeds water, m_wFor the quality of the intracavity water that feeds water, t is the time；

The position relation over time of feedwater intracavity water evaporation interface is:

Wherein, δ_eT () is the position of feedwater intracavity water evaporation interface, ρ_vFor the density of the intracavity steam that feeds water, A_rFor to water cavity Equivalence disengagement area；

Evaporation of water mass flowObtained by following formula:

$\frac{{\stackrel{\·}{m}}_v}{{\mathrm{nA}}_r(P_{sat}-P_o)}=(1-F(\frac{r}{\mathrm{\λ}}\left)\right)(A\stackrel{\‾}{p}+B)+F\left(\frac{r}{\mathrm{\λ}}\right)C---\left(11\right)$

In formula, n is the number in hole on porous plate, A_pFor the area in the single hole of porous plate, r is porous plate average pore radius, μ_vFor the dynamic viscosity of steam, T_eFor evaporation interface temperature Degree, δ_pFor porous plate thickness, m_mFor the quality of water vapour molecule,For stream molecule average speed, P_oFor ambient pressure；P_sat Saturated vapor pressure for evaporating surface；

After completing above-mentioned numerical modeling, the running parameter in this stage is solved, particularly as follows: give fixed pattern (1)-formula (4), formula (7), the initial value of differential variable in formula (8), the thermal resistance in each region, Equivalent Thermal Conductivities can obtain according to formula (5), formula (6), steam Send out mass flow to obtain according to formula (11), then use variable step implicit expression Runge-Kutta method to solve formula (1)-formula (4), formula (7), the differential equation shown in formula (8), the evaporation process of feedwater intracavity is simulated iterative computation；

The condition of set entrance second stage is: if T_w＞ 0 and δ_e(t) ＜ δ_w, δ_wFor to water cavity height；Then think feedwater Underfill, to water cavity, is unsatisfactory for entering the condition of second stage；If T_w＞ 0 and δ_e(t)=δ_w, then it is assumed that feedwater be full of to Water cavity, water sublimate device enters the second stage of start-up course.

Water sublimate device start-up course transient working parameter prediction method the most according to claim 1, it is characterised in that: step The second stage of three pairs of water sublimate device start-up courses carries out HEAT EXCHANGE ANALYSIS and numerical modeling and the detailed process that solves is:

In this stage, water moves in porous plate and evaporates, and the equation of momentum is:

$F_c+F_p-F_{sat}-F_g-F_f=\frac{d\left(mu\right)}{dt}={\mathrm{\ρ}}_w\frac{d\left({\mathrm{\δ}}_e\right(t\left)A_r\mathrm{\ϵ}u\right)}{dt}---\left(13\right)$

That is:

Wherein, F_cFor capillary force, F_pFor feed pressure, F_satFor the saturated vapor pressure of evaporating surface, F_gFor liquid in porous plate Gravity, F_fFor the liquid pressure loss by certain thickness porous material, δ_eT () is the distance that water enters porous plate, A_rFor many The sectional area of orifice plate, ε is the porosity of material, and u is water flowing velocity in porous plate, ρ_wFor the density of water, m is porous plate The quality of interior water；P_cFor capillary pressure, P_pFor feedwater pressure, P_satFor the saturated vapor pressure of evaporating surface, P_fPass through for liquid The pressure loss of certain thickness porous material；

Evaporation interface position δ in porous plate_eT () is determined by following formula (17):

$\frac{2\mathrm{\σ}\cos \mathrm{\θ}}{r}+P_p-P_{sat}-\frac{\mathrm{\μ}{\stackrel{\·}{m}}_v^{\′}{\mathrm{\δ}}_e\left(t\right)}{{\mathrm{K\ρ}}_w{A}_r}-\frac{{\mathrm{\μ\ϵ\δ}}_e\left(t\right)}{K}\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}={\mathrm{\ρ}}_w\frac{d}{dt}\left({\mathrm{\δ}}_e\right(t)\frac{{\stackrel{\·}{m}}_v^{\′}}{{\mathrm{\ρ}}_w{A}_r\mathrm{\ϵ}}+{\mathrm{\δ}}_e(t\left)\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}\right)---\left(17\right)$ Formula (17) In, steam mass flowFor:

Wherein, σ is the surface tension of water, and θ is Liquid contact angle, and μ is the dynamic viscosity of water vapour, and K is the permeability of porous plate Coefficient,For the average pressure of steam, d_p For the average pore size of porous plate, η is the dynamic viscosity of water vapour, and k is Boltzmann constant,

Water enters after porous plate, is all filled by water to water cavity, now water sublimate device heating surface and give the temperature of water cavity can be by formula (19), the transient temperature equation shown in formula (20) obtains:

$c_k{M}_k{\stackrel{\·}{T}}_k=Q_0-(T_k-T_w^{\′})/R_{kw}^{\′}---\left(19\right)$

$c_w^{\′}{M}_w^{\′}{\stackrel{\·}{T}}_w^{\′}=(T_k-T_w^{\′})/R_{kw}^{\′}-(T_w^{\′}-T_{pw})/R_{wp}^{\′}---\left(20\right)$

Be divided into by porous plate containing water section and without water section, then porous plate is containing water section and without water section temperature equation difference For:

$c_{pw}{M}_{pw}{\stackrel{\·}{T}}_{pw}=(T_w^{\′}-T_{pw})/R_{wp}-(T_{pw}-T_{pm})/R_{pm}-{\stackrel{\·}{m}}_v{h}_e---\left(21\right)$

$c_{pm}{M}_{pm}{\stackrel{\·}{T}}_{pm}=(T_{pw}-T_{pm})/R_{pm}---\left(22\right)$

After completing above-mentioned numerical modeling, the running parameter in this stage is solved, particularly as follows: use variable step implicit expression Runge- Kutta method solves the differential equation shown in formula (17), (19)-formula (22), and the evaporation process in porous plate is simulated iteration meter Calculate；Wherein, the thermal resistance in each region, Equivalent Thermal Conductivities can obtain according to method shown in formula (5), formula (6), and evaporative mass flow can Obtained by formula (18)；

The condition of set entrance second stage is: if T_w＞ 0, T_pw＜ 0 and δ_e(t) ＜ δ_w+δ_p, then judge that feedwater is at water Freezing in sublimating apparatus porous plate, water sublimate device enters the phase III of start-up course.

Water sublimate device start-up course transient working parameter prediction method the most according to claim 1, it is characterised in that: step In four, phase III to water sublimate device start-up course carries out HEAT EXCHANGE ANALYSIS and numerical modeling concretely comprising the following steps of solving:

At this stage water sublimate device heating surface, close to the conduction of heat in water cavity, porous plate pool, porous plate ice formation and porous plate gas district System:

$c_k{M}_k{\stackrel{\·}{T}}_k=Q_0-(T_k-T_w^{\′})/R_{kw}^{\′}---\left(23\right)$

$c_w^{\′}{M}_w^{\′}{\stackrel{\·}{T}}_w^{\′}=(T_k-T_w^{\′})/R_{kw}^{\′}-(T_w^{\′}-T_{pw})/R_{wp}^{\′}---\left(24\right)$

$c_{pw}{M}_{pw}{\stackrel{\·}{T}}_{pw}=(T_w-T_{pw})/R_{wpw}-(T_{pw}-T_{pi})/R_{pwi}-{\stackrel{\·}{m}}_{wi}{h}_f---\left(25\right)$

$c_{pi}{M}_{pi}{\stackrel{\·}{T}}_{pi}=(T_{pw}-T_{pi})/R_{pwi}-(T_{pi}-T_{pm})/R_{pim}-{\stackrel{\·}{m}}_s{h}_s---\left(26\right)$

$c_{pm}{M}_{pm}{\stackrel{\·}{T}}_{pm}=(T_{pi}-T_{pm})/R_{pim}---\left(27\right)$

Wherein,For the thawing/solidifying phase variable Rate at ice-water interface,Rate of sublimation for the interface that distils；

Can be determined by formula (28):

Thawing or freezing rate at frozen water interface are then determined by formula (29):

$\left|\frac{{\mathrm{dm}}_{wi}}{dt}\right|=\frac{Q_{pw}-Q_{wi}}{h_f}=\left|{\mathrm{\ρ}}_w{\mathrm{\ϵA}}_r\frac{d\left({\mathrm{\δ}}_e\right(t\left)\right)}{dt}\right|---\left(29\right)$

Wherein, Q_pwFor by the heat to water cavity conduction to aqueous porous plate, Q_wiFor the heat by ice sheet dissipation to space outerpace,

After completing above-mentioned numerical modeling, the running parameter in this stage is solved, particularly as follows: use variable step implicit expression Runge- Kutta method solves the differential equation shown in formula (23)-formula (27), and the sublimation process in porous plate is simulated iterative computation.