CN105526729A - Simulation analysis method of light-weight miniature coaxial pulse tube refrigerator based on CFD technology - Google Patents

Abstract

The invention discloses a simulation analysis method of a light-weight miniature coaxial pulse tube refrigerator based on a CFD technology. The simulation analysis method comprises the following steps: (1) establishing the geometric model of the light-weight miniature coaxial pulse tube refrigerator; (2) dividing and setting control areas; (3) setting boundary conditions; (4) dividing grids; (5) establishing unsteady governing equations; (6) conducting the unsteady numerical solution; and (7) analyzing the simulation result, and predicting the refrigeration performance of the light-weight miniature coaxial pulse tube refrigerator. The two-dimensional computer simulation of the light-weight miniature coaxial pulse tube refrigerator is realized through commercial computational fluid software Fluent, the refrigeration performance of the light-weight miniature coaxial pulse tube refrigerator is predicted within a shorter time, and the development cost is reduced greatly; and the operating mechanism and the fluid characteristic in a pulse tube are known deeply by analyzing the distribution of the temperature field, the pressure field and the speed filed in the light-weight miniature coaxial pulse tube refrigerator, and the simulation analysis method has important significance in the design optimization of the light-weight miniature coaxial pulse tube refrigerator.

Description

Based on the analog analysing method of the light weight mini-coax pulse tube refrigerating machine of CFD technology

Technical field

The present invention relates to refrigeration & cryogenic engineering field, particularly a kind of analog analysing method of the light weight mini-coax pulse tube refrigerating machine based on CFD technology.

Background technology

Pulse tube refrigerating machine is a significant innovation to philip refrigerator.It has fully phased out the displacer of cold junction, and pressure wave needed for realizing freezing by the running of hot junction phase modulating mechanism and the phase difference of mass flow.On the one hand, fully phasing out of cold junction moving component, makes pulse tube refrigerating machine can realize the low vibration of cold junction, low interference and without wearing and tearing; On the other hand, through a series of important improvement, at some typical warm areas, its actual efficiency has also reached the peak of philip refrigerator.These remarkable advantages make pulse tube refrigerating machine become a big hot topic of regenerating type low-temperature machinery cold research over nearly 30 years, all obtain a wide range of applications in Aero-Space, low-temperature electronics, superconduction industry and cryosurgery industry etc.

In recent years, ground tactics and space application all propose more and more harsher requirement to the volume of Cryo Refrigerator, weight and cooling rate, thus facilitate developing rapidly of microminiaturized Refrigeration Technique.Conventional JT refrigeration machine and miniature Stirling refrigerating machine are the representatives of the small-sized Cryo Refrigerator realizing microminiaturization and fast-refrigerating, but the JT efficiency of refrigerator of dividing wall type is low, working life is extremely short, and miniature Stirling refrigerating machine is in vibration, wearing and tearing, the aspect such as reliability and working life is also more and more difficult to meet application demand, thus correlative study has focused on the microminiaturization of pulse tube refrigerating machine again, attempt to increase substantially operating frequency and the blowing pressure, the energy loss brought because working medium volume reduces is made up in the mode increasing energy density, thus realize needed for volume little, the object of lightweight and fast-refrigerating.

Pulse tube refrigerating machine is mainly divided into linear pattern, coaxial type and U-shaped three kinds of arrangements, and wherein the structure of coaxial type is the compactest, is the optimal selection of minisize pulse tube refrigerator.But the loss of this arrangement is large, and the complexity that influences each other between parts, actual design difficulty is very large.CFD is the abbreviation of Fluid Mechanics Computation, by commercial CFD software Fluent in practice, can carry out calculating and the forecast analysis in flow field in the short period of time, significantly use manpower and material resources sparingly.Intuitively can show a series of flow field characteristic such as temperature field, pressure field simultaneously, more deep understanding be produced to the operating mechanism of control domain inner fluid, important directive function is played to optimal design and experiment test.The domestic research about light weight mini-coax pulse tube refrigerating machine at present is still in the starting stage, its outstanding behaviours is that R&D costs are high, design cycle is long, and lacks for the operating mechanism of mini refrigerating machine inside and fluid behaviour and understand systematicly, and causes actual model-performance undesirable.

Summary of the invention

In view of the deficiencies in the prior art, the present invention proposes a kind of analog analysing method of the light weight mini-coax pulse tube refrigerating machine based on CFD technology.

The object of the invention is to, business Fluid Computation software Fluent is adopted to set up the two-dimensional axial symmetric computation model of light weight mini-coax pulse tube refrigerating machine, carry out Non-Steady Numerical Simulation, the refrigeration performance of prediction refrigeration machine, analyze the temperature field of refrigerating device inner, pressure field and velocity field distribution, the operating mechanism of profound understanding pulse tube inside and fluid behaviour, for the design of practicality and high efficiency type provides theoretical foundation.

This analog analysing method comprises the following steps:

Step one: the geometrical model setting up light weight mini-coax pulse tube refrigerating machine.This model is two-dimensional axial symmetric model, comprises grade aftercooler 1, regenerator 2, cool end heat exchanger 3, pulse tube 4, hot end heat exchanger 5, first paragraph inertia tube 6, second segment inertia tube 7 and air reservoir 8.

Step 2: divide and setup control region.The internal control region at level aftercooler 1, regenerator 2, cool end heat exchanger 3, hot end heat exchanger 5 place is set as porous media region, and other internal control regions are set as gas zones.

Step 3: boundary condition is set.The left boundary arranging grade aftercooler 1 is pressure entrance border, is write by UDF file:

p _in＝p ₀+Δpsin2πft(1)

In formula: p _infor the inlet pressure of level aftercooler 1, p ₀for the blowing pressure, Δ p is dynamic pressure amplitude, and f is running frequency, and t is the time.

The outer boundary arranging grade aftercooler 1, hot end heat exchanger 5, first paragraph inertia tube 6, second segment inertia tube 7 and air reservoir 8 is isothermal border; Regenerator 2 is set, the outer boundary of pulse tube 4 is adiabatic boundary; The outer boundary of cool end heat exchanger 3 has two kinds of methods to set up, can measure the zero load cryogenic temperature of refrigeration machine when being set to adiabatic boundary, can measure the refrigerating capacity of refrigeration machine under this cold junction temperature when being set to isothermal border.

Step 4: grid division.All structured grid is adopted to the control domain of this model, near border, carries out Local grid refinement simultaneously.

Step 5: set up non steady control equation.With the Working medium gas helium of refrigerating device inner for object, respectively governing equation is set up to gas zones and porous media region:

First governing equation is set up to gas zones:

Continuity equation:

$\frac{\∂{\mathrm{\ρ}}_f}{\∂t}+\▿\·\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\right)=0---\left(2\right)$

In formula: ρ _fwith be respectively density and the speed of gas.

The equation of momentum:

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\right)+\▿\·\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\stackrel{\→}{v}\right)=-\▿p+\▿\·\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\right)---\left(3\right)$

In formula: p is static pressure, for pressure tensor.

Energy equation:

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_f{E}_f\right)+\▿\·\[\stackrel{\→}{v}\left({\mathrm{\ρ}}_f{E}_f+p\right)\]=\▿\·\[k_f\▿T+\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\·\stackrel{\→}{v}\right)\]---\left(4\right)$

In formula: k _ffor gas conduction rate, T is temperature, E _ffor gas total energy.

Afterwards governing equation is set up to porous media region:

Continuity equation:

$\frac{\∂{\mathrm{\ϵ\ρ}}_f}{\∂t}+\▿\·\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\right)=0---\left(5\right)$

The equation of momentum:

$\frac{\∂}{\∂t}\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\right)+\▿\·\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\stackrel{\→}{v}\right)=-\mathrm{\ϵ}\▿p+\▿\·\left(\mathrm{\ϵ}\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\right)+S_i---\left(6\right)$

In formula: ε is the porosity of porous media.Porous media region is filled by silk screen, and therefore porosity ε is by meshcount m and string diameter d _wdetermine:

$\mathrm{\ϵ}=1-\frac{{\mathrm{\πmd}}_w}{4\×0.0254}---\left(7\right)$

The hydraulic diameter d of silk screen _hcan be calculated by following formula:

$d_h=\frac{\mathrm{\ϵ}}{1-\mathrm{\ϵ}}{d}_w---\left(8\right)$

S _ifor the source item in the equation of momentum, lose item two parts by viscosity loss item and inertia and form, for isotropic medium, S _ican be expressed as:

$S_i=-\left(\frac{\mathrm{\μ}}{\mathrm{\α}}\stackrel{\→}{v}+\frac{C_2}{2}{\mathrm{\ρ}}_f\left|\stackrel{\→}{v}\right|\stackrel{\→}{v}\right)---\left(9\right)$

In formula: μ is the dynamic viscosity of gas, α is permeability, C ₂for inertial resistance coefficient.α and C ₂can calculate by the following method:

The axial pressure drop in porous media region can be expressed as:

$\frac{dp}{dx}=-\frac{f_{osc}}{2d_h}{\mathrm{\ρ}}_f{v}^2---\left(10\right)$

In formula: f _oscfor the friction factor under concussion stream mode, rule of thumb formula:

$f_{osc}=\frac{129}{\mathrm{Re}}+2.91{\mathrm{Re}}^{-0.103}---\left(11\right)$

In formula: Re is the Reynolds number of porous media region inner fluid, can be expressed as:

$\mathrm{Re}=\frac{{\mathrm{\ρ}}_f\left|\stackrel{\→}{v}\right|d_h}{\mathrm{\μ}}---\left(12\right)$

Convolution (9), formula (10) and formula (11), can determine permeability α and inertial resistance coefficient C ₂:

$\mathrm{\α}=\frac{d_h^2}{64.5}---\left(13\right)$

$C_2=\frac{2.91}{d_h}{\mathrm{Re}}^{-0.103}---\left(14\right)$

Energy equation:

$\frac{\∂}{\∂t}\[{\mathrm{\ϵ\ρ}}_f{E}_f+\left(1-\mathrm{\ϵ}\right){\mathrm{\ρ}}_s{E}_s\]+\▿\·\[\stackrel{\→}{v}\left({\mathrm{\ρ}}_f{E}_f+p\right)\]=\▿\·\[k_{eff}\▿T+\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\·\stackrel{\→}{v}\right)\]---\left(15\right)$

In formula: ρ _sand E _sbe respectively density and the total energy of solid, k _efffor the effective thermal conductivity of porous media.

Due to the impact of thermal contact resistance, the effective thermal conductivity k of porous media _effbe far smaller than its average conduction, can revise as follows:

$k_{eff}=k_f\mathrm{\ϵ}+k_s(1-\mathrm{\ϵ}){\left(\frac{k_s}{k_f}\right)}^{-0.835}\[\frac{3\left(\frac{k_s}{k_f}-\mathrm{\ϵ}\right)+\left(2+\frac{k_s}{k_f}\right)\mathrm{\ϵ}}{3\left(1-\mathrm{\ϵ}\right)+\left(2+\frac{k_s}{k_f}\right)\mathrm{\ϵ}}\]---\left(16\right)$

Step 6: carry out Unsteady Numerical and solve.According to the light weight mini-coax pulse tube refrigerating machine CFD model that above-mentioned steps is set up, given solver controling parameters, discrete scheme and residual error convergence, initialize flow field, start to carry out Unsteady Numerical and solve.When the external boundary of cool end heat exchanger 3 is set to adiabatic boundary condition, monitoring objective is the face mean temperature on this border, when the changing value of this temperature within continuous 100 cycles is all less than 1%, thinks to calculate and stablizes.When the external boundary of cool end heat exchanger 3 is set to boundary condition, monitoring objective is total rate of heat flow on this border, when the changing value of this rate of heat flow within continuous 100 cycles is all less than 1%, thinks to calculate and stablizes.

Step 7: analysis mode result.After calculating is stable, by the analysis to analog result, the refrigeration performance of refrigeration machine can be predicted.When the external boundary of cool end heat exchanger 3 is set to adiabatic boundary condition, under stable state, the face mean temperature on this border is the zero load cryogenic temperature that refrigeration machine can reach; When the external boundary of cool end heat exchanger 3 is set to boundary condition, under stable state, the time equal rate of heat flow on this border is the refrigerating capacity that refrigeration machine can obtain under setting cold junction temperature.In addition, by the distribution situation of observation refrigerating device inner temperature field, pressure field and velocity field, the operating mechanism of refrigeration machine and the flow regime of fluid can be analyzed further.

The invention has the advantages that:

(1) by this analog analysing method, the refrigeration performance of light weight mini-coax pulse tube refrigerating machine can be predicted in the short period of time;

(2) by this analog analysing method, a series of analog results such as the temperature field of light weight mini-coax pulse tube refrigerating machine inside, pressure field and velocity field can intuitively be shown, the operating mechanism of profound understanding pulse tube inside and fluid behaviour.

(3) the present invention significantly reduces the R&D costs of light weight mini-coax pulse tube refrigerating machine, shortens the R&D cycle.

Above-mentioned advantage can realize the fast prediction of light weight mini-coax pulse tube refrigerating machine refrigeration performance, the operating mechanism of profound understanding pulse tube inside and fluid behaviour, and the design optimization for light weight mini-coax pulse tube refrigerating machine is significant.

Accompanying drawing explanation

Fig. 1 is the geometrical model schematic diagram of light weight mini-coax pulse tube refrigerating machine.

Fig. 2 is the refrigeration machine temperature lowering curve of simulation and forecast.

Fig. 3 is the program flow diagram of analog analysing method.

Wherein: 1 is level aftercooler; 2 is regenerator; 3 is cool end heat exchanger; 4 is pulse tube; 5 is hot end heat exchanger; 6 is first paragraph inertia tube; 7 is second segment inertia tube; 8 is air reservoir.

Detailed description of the invention

Below in conjunction with drawings and Examples, the specific embodiment of the present invention is described in further detail; but be not limited thereto; everyly technical solution of the present invention is modified or equivalent to replace; and do not depart from the spirit and scope of technical solution of the present invention, all should be encompassed in protection scope of the present invention.

In view of the deficiencies in the prior art, the present invention proposes a kind of analog analysing method of the light weight mini-coax pulse tube refrigerating machine based on CFD technology.

The object of the invention is to, business Fluid Computation software Fluent is adopted to set up the two-dimensional axial symmetric computation model of light weight mini-coax pulse tube refrigerating machine, carry out Non-Steady Numerical Simulation, the refrigeration performance of prediction refrigeration machine, analyze the temperature field of refrigerating device inner, pressure field and velocity field distribution, the operating mechanism of profound understanding pulse tube inside and fluid behaviour, for the design of practicality and high efficiency type provides theoretical foundation.

In the present embodiment, the light weight mini-coax pulse tube refrigerating machine of sunykatuib analysis object to be an operating frequency be 128Hz.

This analog analysing method comprises the following steps:

Step one: the geometrical model setting up light weight mini-coax pulse tube refrigerating machine.Fig. 1 is the geometrical model schematic diagram of light weight mini-coax pulse tube refrigerating machine.This model is two-dimensional axial symmetric model, comprises grade aftercooler 1, regenerator 2, cool end heat exchanger 3, pulse tube 4, hot end heat exchanger 5, first paragraph inertia tube 6, second segment inertia tube 7 and air reservoir 8.

Step 2: divide and setup control region.The internal control region at level aftercooler 1, regenerator 2, cool end heat exchanger 3, hot end heat exchanger 5 place is set as porous media region, the shade overlay area namely shown in Fig. 1.Other internal control regions are set as gas zones.

Step 3: boundary condition is set.The left boundary arranging grade aftercooler 1 is pressure entrance border, is write by UDF file:

pin＝p0+Δpsin2πft(1)

In formula: p _infor the inlet pressure of level aftercooler 1, p ₀for the blowing pressure, Δ p is dynamic pressure amplitude, and f is running frequency.In the present embodiment, the blowing pressure 3.6MPa, dynamic pressure amplitude 0.377MPa, running frequency is 128Hz.

The outer boundary arranging grade aftercooler 1, hot end heat exchanger 5, first paragraph inertia tube 6, second segment inertia tube 7 and air reservoir 8 is isothermal border, and temperature is 300K; Regenerator 2 is set, the outer boundary of pulse tube 4 is adiabatic boundary; The outer boundary of cool end heat exchanger 3 has two kinds of methods to set up, can measure the zero load cryogenic temperature of refrigeration machine when being set to adiabatic boundary, can measure the refrigerating capacity of refrigeration machine under this cold junction temperature when being set to isothermal border.In the present embodiment, the outer boundary of cool end heat exchanger 3 is set to adiabatic boundary condition.

Step 4: grid division.All structured grid is adopted to the control domain of this model, near border, carries out Local grid refinement simultaneously.

Step 5: set up non steady control equation.With the Working medium gas helium of refrigerating device inner for object, respectively governing equation is set up to gas zones and porous media region:

First governing equation is set up to gas zones:

Continuity equation:

$\frac{\∂\mathrm{\ρ}f}{\∂t}+\▿\·\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\right)=0---\left(2\right)$

In formula: ρ _fwith be respectively density and the speed of gas.

The equation of momentum:

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\right)+\▿\·\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\stackrel{\→}{v}\right)=-\▿p+\▿\·\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\right)---\left(3\right)$

In formula: p is static pressure, for pressure tensor.

Energy equation:

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_f{E}_f\right)+\▿\·\[\stackrel{\→}{v}\left({\mathrm{\ρ}}_f{E}_f+p\right)\]=\▿\·\[k_f\▿T+\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\·\stackrel{\→}{v}\right)\]---\left(4\right)$

In formula: k _ffor gas conduction rate, T is temperature, E _ffor gas total energy.

Afterwards governing equation is set up to porous media region:

Continuity equation:

$\frac{\∂{\mathrm{\ϵ\ρ}}_f}{\∂t}+\▿\·\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\right)=0---\left(5\right)$

The equation of momentum:

$\frac{\∂}{\∂t}\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\right)+\▿\·\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\stackrel{\→}{v}\right)=-\mathrm{\ϵ}\▿p+\▿\·\left(\mathrm{\ϵ}\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\right)+S_i---\left(6\right)$

In formula: ε is the porosity of porous media.Porous media region is filled by silk screen, and therefore porosity ε is by meshcount m and string diameter d _wdetermine:

$\mathrm{\ϵ}=1-\frac{{\mathrm{\πmd}}_w}{4\×0.0254}---\left(7\right)$

The hydraulic diameter d of silk screen _hcan be calculated by following formula:

$d_h=\frac{\mathrm{\ϵ}}{1-\mathrm{\ϵ}}{d}_w---\left(8\right)$

S _ifor the source item in the equation of momentum, lose item two parts by viscosity loss item and inertia and form, for isotropic medium, S _ican be expressed as:

$S_i=-\left(\frac{\mathrm{\μ}}{\mathrm{\α}}\stackrel{\→}{v}+\frac{C_2}{2}{\mathrm{\ρ}}_f\left|\stackrel{\→}{v}\right|\stackrel{\→}{v}\right)---\left(9\right)$

In formula: μ is the dynamic viscosity of gas, α is permeability, C ₂for inertial resistance coefficient.α and C ₂can calculate by the following method:

The axial pressure drop in porous media region can be expressed as:

$\frac{dp}{dx}=-\frac{f_{osc}}{2d_h}{\mathrm{\ρ}}_f{v}^2---\left(10\right)$

In formula: f _oscfor the friction factor under concussion stream mode, rule of thumb formula:

$f_{osc}=\frac{129}{\mathrm{Re}}+2.91{\mathrm{Re}}^{-0.103}---\left(11\right)$

In formula: Re is the Reynolds number of porous media region inner fluid, can be expressed as:

$\mathrm{Re}=\frac{{\mathrm{\ρ}}_f\left|\stackrel{\→}{v}\right|d_h}{\mathrm{\μ}}---\left(12\right)$

Convolution (9), formula (10) and formula (11), can determine permeability α and inertial resistance coefficient C ₂:

$\mathrm{\α}=\frac{d_h^2}{64.5}---\left(13\right)$

$C_2=\frac{2.91}{d_h}{\mathrm{Re}}^{-0.103}---\left(14\right)$

Energy equation:

$\frac{\∂}{\∂t}\[{\mathrm{\ϵ\ρ}}_f{E}_f+\left(1-\mathrm{\ϵ}\right){\mathrm{\ρ}}_s{E}_s\]+\▿\·\[\stackrel{\→}{v}\left({\mathrm{\ρ}}_f{E}_f+p\right)\]=\▿\·\[k_{eff}\▿T+\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\·\stackrel{\→}{v}\right)\]---\left(15\right)$

In formula: ρ _sand E _sbe respectively density and the total energy of solid, k _efffor the effective thermal conductivity of porous media.

Due to the impact of thermal contact resistance, the effective thermal conductivity k of porous media _effbe far smaller than its average conduction, can revise as follows:

$k_{eff}=k_f\mathrm{\ϵ}+k_s(1-\mathrm{\ϵ}){\left(\frac{k_s}{k_f}\right)}^{-0.835}\[\frac{3\left(\frac{k_s}{k_f}-\mathrm{\ϵ}\right)+\left(2+\frac{k_s}{k_f}\right)\mathrm{\ϵ}}{3\left(1-\mathrm{\ϵ}\right)+\left(2+\frac{k_s}{k_f}\right)\mathrm{\ϵ}}\]---\left(16\right)$

In the present embodiment, the porous media area filling 635 order stainless steel cloth of regenerator 2 inside, the porous media area filling 100 order copper mesh of level aftercooler 1, cool end heat exchanger 3 and hot end heat exchanger 5 inside.Table 1 summarizes by each parameter of the above-mentioned silk screen calculated.

Table 1 porous media intra-zone silk screen parameter

Region	Material	m	d w(μm)	ε	d h(μm)	α(m 2)	C 2(m -1)
								Regenerator 2	SS304	635	20	0.607	31.4	1.5×10 -11	81460
Level aftercooler 1	Copper	100	100	0.691	222.6	7.7×10 -10	12173
								Cool end heat exchanger 3	Copper	100	100	0.691	222.6	7.7×10 -10	12173
Hot end heat exchanger 5	Copper	100	100	0.691	222.6	7.7×10 -10	12173

Step 6: carry out Unsteady Numerical and solve.According to the light weight mini-coax pulse tube refrigerating machine CFD model that above-mentioned steps is set up, given solver controling parameters, discrete scheme and residual error convergence, initialize flow field, start to carry out Unsteady Numerical and solve.When the external boundary of cool end heat exchanger 3 is set to adiabatic boundary condition, monitoring objective is the face mean temperature on this border, when the changing value of this temperature within continuous 100 cycles is all less than 1%, thinks to calculate and stablizes.When the external boundary of cool end heat exchanger 3 is set to boundary condition, monitoring objective is total rate of heat flow on this border, when the changing value of this rate of heat flow within continuous 100 cycles is all less than 1%, thinks to calculate and stablizes.

The present embodiment selects laminar model to calculate, and the PISO algorithm of solver selection pressure and speed coupling, discrete scheme adopts Second-order Up-wind form, and residual error convergence energy term is 10 ^-6, other parameters are 10 ^-3, initial temperature is 300K, time step 5 × 10 ^-5s, each walks maximum iteration 50 times.Monitoring objective is interface, the outside mean temperature of cool end heat exchanger 3, when the changing value of this temperature within continuous 100 cycles is all less than 1%, thinks to calculate and stablizes.

Step 7: analysis mode result.After calculating is stable, by the analysis to analog result, can predict the refrigeration performance of refrigeration machine, when the external boundary of cool end heat exchanger 3 is set to adiabatic boundary condition, under stable state, the face mean temperature on this border is the zero load cryogenic temperature that refrigeration machine can reach; When the external boundary of cool end heat exchanger 3 is set to boundary condition, under stable state, the time equal rate of heat flow on this border is the refrigerating capacity that refrigeration machine can obtain under setting cold junction temperature.In the present embodiment, the external boundary of cool end heat exchanger 3 is set to adiabatic boundary condition, therefore can predict the zero load cryogenic temperature of refrigeration machine.In addition, by the distribution situation of observation refrigerating device inner temperature field, pressure field and velocity field, the operating mechanism of refrigeration machine and the flow regime of fluid can be analyzed further.

Fig. 2 is the refrigeration machine temperature lowering curve of simulation and forecast, i.e. interface, the outside mean temperature monitoring curve of cool end heat exchanger 3.When computing time is after about 100s, calculate and reach stable state, this temperature stabilization is at about 64.5K, and namely the zero load temperature of refrigeration machine is 64.5K.

Claims (1)

1., based on the analog analysing method of the light weight mini-coax pulse tube refrigerating machine of CFD technology, it is characterized in that:

Business Fluid Computation software Fluent is adopted to set up the two-dimensional axial symmetric computation model of light weight mini-coax pulse tube refrigerating machine, carry out Non-Steady Numerical Simulation, the refrigeration performance of prediction refrigeration machine, analyze the temperature field of refrigerating device inner, pressure field and velocity field distribution, the operating mechanism of profound understanding pulse tube inside and fluid behaviour; Specifically comprise the following steps:

Step one: the geometrical model setting up light weight mini-coax pulse tube refrigerating machine, this model is two-dimensional axial symmetric model, comprises a grade aftercooler (1), regenerator (2), cool end heat exchanger (3), pulse tube (4), hot end heat exchanger (5), first paragraph inertia tube (6), second segment inertia tube (7) and air reservoir (8);

Step 2: divide and setup control region, the internal control region at level aftercooler (1), regenerator (2), cool end heat exchanger (3), hot end heat exchanger (5) place is set as porous media region, and other internal control regions are set as gas zones;

Step 3: boundary condition is set, the left boundary arranging a grade aftercooler (1) is pressure entrance border, is write by UDF file:

p _in＝p ₀+Δpsin2πft(1)

In formula: p _infor the inlet pressure of level aftercooler (1), p ₀for the blowing pressure, Δ p is dynamic pressure amplitude, and f is running frequency, and t is the time;

The outer boundary arranging a grade aftercooler (1), hot end heat exchanger (5), first paragraph inertia tube (6), second segment inertia tube (7) and air reservoir (8) is isothermal border; Regenerator (2) is set, the outer boundary of pulse tube (4) is adiabatic boundary; The outer boundary of cool end heat exchanger (3) has two kinds of methods to set up, can measure the zero load cryogenic temperature of refrigeration machine when being set to adiabatic boundary, can measure the refrigerating capacity of refrigeration machine under this cold junction temperature when being set to isothermal border;

Step 4: grid division, all adopts structured grid to the control domain of this model, carries out Local grid refinement near border simultaneously;

Step 5: set up non steady control equation, with the Working medium gas helium of refrigerating device inner for object, respectively governing equation is set up to gas zones and porous media region:

First governing equation is set up to gas zones:

Continuity equation:

$\frac{\∂{\mathrm{\ρ}}_f}{\∂t}+\▿\·\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\right)=0---\left(2\right)$

In formula: ρ _fwith be respectively density and the speed of gas.

The equation of momentum:

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\right)+\▿\·\left({\mathrm{\ρ}}_f\stackrel{\→}{v}\stackrel{\→}{v}\right)=-\▿p+\▿\·\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\right)---\left(3\right)$

In formula: p is static pressure, for pressure tensor;

Energy equation:

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_f{E}_f\right)+\▿\·\[\stackrel{\→}{v}\left({\mathrm{\ρ}}_f{E}_f+p\right)\]=\▿\·\[k_f\▿T+\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\·\stackrel{\→}{v}\right)\]---\left(4\right)$

In formula: k _ffor gas conduction rate, T is temperature, E _ffor gas total energy;

Afterwards governing equation is set up to porous media region:

Continuity equation:

$\frac{\∂{\mathrm{\ϵ\ρ}}_f}{\∂t}+\▿\·\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\right)=0---\left(5\right)$

The equation of momentum:

$\frac{\∂}{\∂t}\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\right)+\▿\·\left({\mathrm{\ϵ\ρ}}_f\stackrel{\→}{v}\stackrel{\→}{v}\right)=-\mathrm{\ϵ}\▿p+\▿\·\left(\mathrm{\ϵ}\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\right)+S_i---\left(6\right)$

In formula: ε is the porosity of porous media, and porous media region is filled by silk screen, and therefore porosity ε is by meshcount m and string diameter d _wdetermine:

$\mathrm{\ϵ}=1-\frac{{\mathrm{\πmd}}_w}{4\×0.0254}---\left(7\right)$

The hydraulic diameter d of silk screen _hcan be calculated by following formula:

$d_h=\frac{\mathrm{\ϵ}}{1-\mathrm{\ϵ}}{d}_w---\left(8\right)$

S _ifor the source item in the equation of momentum, lose item two parts by viscosity loss item and inertia and form, for isotropic medium, S _ican be expressed as:

$S_i=-\left(\frac{\mathrm{\μ}}{\mathrm{\α}}\stackrel{\→}{v}+\frac{C_2}{2}{\mathrm{\ρ}}_f\left|\stackrel{\→}{v}\right|\stackrel{\→}{v}\right)---\left(9\right)$

In formula: μ is the dynamic viscosity of gas, α is permeability, C ₂for inertial resistance coefficient, α and C ₂can calculate by the following method:

The axial pressure drop in porous media region can be expressed as:

$\frac{dp}{dx}=-\frac{f_{osc}}{2d_h}{\mathrm{\ρ}}_f{v}^2---\left(10\right)$

In formula: f _oscfor the friction factor under concussion stream mode, rule of thumb formula:

$f_{osc}=\frac{129}{\mathrm{Re}}+2.91{\mathrm{Re}}^{-0.103}---\left(11\right)$

In formula: Re is the Reynolds number of porous media region inner fluid, can be expressed as:

$\mathrm{Re}=\frac{{\mathrm{\ρ}}_f\left|\stackrel{\→}{v}\right|d_h}{\mathrm{\μ}}---\left(12\right)$

Convolution (9), formula (10) and formula (11), can determine permeability α and inertial resistance coefficient C ₂:

$\mathrm{\α}=\frac{d_h^2}{64.5}---\left(13\right)$

$C_2=\frac{2.91}{d_h}{\mathrm{Re}}^{-0.103}---\left(14\right)$

Energy equation:

$\frac{\∂}{\∂t}\[{\mathrm{\ϵ\ρ}}_f{E}_f+\left(1-\mathrm{\ϵ}\right){\mathrm{\ρ}}_s{E}_s\]+\▿\·\[\stackrel{\→}{v}\left({\mathrm{\ρ}}_f{E}_f+p\right)\]=\▿\·\[k_{eff}\▿T+\left(\stackrel{\‾}{\stackrel{\‾}{\mathrm{\τ}}}\·\stackrel{\→}{v}\right)\]---\left(15\right)$

In formula: ρ _sand E _sbe respectively density and the total energy of solid, k _efffor the effective thermal conductivity of porous media;

Due to the impact of thermal contact resistance, the effective thermal conductivity k of porous media _effbe far smaller than its average conduction, can revise as follows:

$k_{eff}=k_f\mathrm{\ϵ}+k_s\left(1-\mathrm{\ϵ}\right){\left(\frac{k_s}{k_f}\right)}^{-0.835}\[\frac{3\left(\frac{k_s}{k_f}-\mathrm{\ϵ}\right)+\left(2+\frac{k_s}{k_f}\right)\mathrm{\ϵ}}{3\left(1-\mathrm{\ϵ}\right)+\left(2+\frac{k_s}{k_f}\right)\mathrm{\ϵ}}\];---\left(16\right)$

Step 6: carry out Unsteady Numerical and solve, according to the light weight mini-coax pulse tube refrigerating machine CFD model that above-mentioned steps is set up, given solver controling parameters, discrete scheme and residual error convergence, initialize flow field, start to carry out Unsteady Numerical and solve.When the external boundary of cool end heat exchanger (3) is set to adiabatic boundary condition, monitoring objective is the face mean temperature on this border, when the changing value of this temperature within continuous 100 cycles is all less than 1%, thinks to calculate and stablizes.When the external boundary of cool end heat exchanger (3) is set to boundary condition, monitoring objective is total rate of heat flow on this border, when the changing value of this rate of heat flow within continuous 100 cycles is all less than 1%, thinks to calculate and stablizes;

Step 7: analysis mode result, after calculating is stable, by the analysis to analog result, can predict the refrigeration performance of refrigeration machine.When the external boundary of cool end heat exchanger (3) is set to adiabatic boundary condition, under stable state, the face mean temperature on this border is the zero load cryogenic temperature that refrigeration machine can reach; When the external boundary of cool end heat exchanger (3) is set to boundary condition, under stable state, the time equal rate of heat flow on this border is the refrigerating capacity that refrigeration machine can obtain under setting cold junction temperature, in addition, by the distribution situation of observation refrigerating device inner temperature field, pressure field and velocity field, the operating mechanism of refrigeration machine and the flow regime of fluid can be analyzed further.