ASHRAE IJHVAC 17-6-2011 HVAC&R RESEARCH An International Journal of Heating Ventilating Air-Conditioning and Refrigerating Research.pdf

1、HVAC Ei(Engineering Information, Inc.) Compendex and EngineeringIndex; ISI (Institute for Scientific Information) Web Scienceand Research Alert; BSRIA (Building Services Research ACS (American Chemical Society) Chem-ical Abstracts Service and Scientific and Technical Informa-tion Network; CSA: Guide

2、 to Discovery CSA Materials Re-search Database with METADEX, CSA Engineering ResearchDatabase, and CSA High Technology Research Database withAerospace;IIR(InternationalInstituteofRefrigeration)Bulletinof the IIR and Fridoc; and Thomson Gale. Current contents arein ISI Engineering, Computing Online I

3、SSN: 1938-5587Institutional Subscribers: $270, 150, 216.Personal Subscribers: $175, 97, 140.Production and Advertising Office: Taylor Shonderetal.2000).Morerecently,thepresentauthors(Spitleretal.2009)have been involved in a comparison of energy simulations of GHX. In both cases, notable differencesw

4、ere reported which prompted some changes in some methods and indicate the usefulness of suchexercises. However, there is a need to standardize these tests to cover the full spectrum of conditionsmuch like what has been done for building energy simulation programs with BESTEST.Finally, the guest edit

5、ors would like to thank the reviewers; without the help of about 50 reviewers itwould not have been possible to produce this high-quality special issue. Also we would like to speciallythank the Editor-in-Chief, Reinhard Radermacher, and the Journal Editorial Assistant, Mary Baugher, fortheir continu

6、ed support through this whole process.ReferencesShonder, J. A., V. Baxter, J. Thornton, and P. Hughes. 1999. A new comparison of vertical ground heat exchanger design methodsfor residential applications.” ASHRAE Transactions 102(2): 117988.Shonder, J. A., V. D. Baxter, P. J. Hughes, and J. W. Thornt

7、on. 2000. A comparison of vertical ground heat exchanger design softwarefor commercial applications.” ASHRAE Transactions 106(1): 83142.Downloaded by T accepted July 7, 2011Johan Claesson is Professor. Goran Hellstrom is Associated Professor.or thermally enhanced grouts with higher thermalconductivi

8、ty. In Northern Europe, the boreholes areusually filled with groundwater to a few meters be-low the ground surface.The heat transfer between the heat carrier fluidand the surrounding ground depends on the ar-rangement of the flow channels, the convective heattransferinsidethechannels,andthethermalpr

9、oper-ties of the materials involved. The thermal processin the ground surrounding the borehole is usuallytreated as transient heat conduction unless substan-tial groundwater flow has to be accounted for. Theprocess inside the borehole is described as a bore-holethermalresistancebetweentheheatcarrier

10、fluidand the borehole wall. Axial heat conduction is895HVAC Claesson and Bennet 1987). The methodis also applied to district heating pipes (Wallenten1991). The algorithms have recently been revisedand further developed (Claesson forthcoming). Thecondensed version presented here is based on thisnew v

11、ersion. A main change is the discussion aboutaverage radial temperature. The condition of con-stant temperature at a circle outside the boreholeis replaced by a condition on the mean tempera-ture around the wall of the borehole. This simplifiesthe formulas, and the problem of choosing a suit-able ra

12、dius is avoided. Another improvement is thatthe analysis starts with prescribed heat fluxes fromthe pipes and not prescribed fluid temperatures. Thediscussion of thermal networks is now separatedfrom the basic solutions for prescribed heat fluxesmore clearly.Theconsideredsteady-statesolutionsforacro

13、ss-section perpendicular to the borehole axis do notconsider rapid fluid temperature variation or axialeffects. But the analysis and, in particular, the ther-mal resistances represent a key performance factorfor the BHE. The resistances determine the tem-perature differences between fluids and groun

14、d justoutside the borehole, which are required to sustainany desired heat injection or extraction.The multipole method provides an exact algo-rithm to compute the steady-state borehole thermalresistance. The method has been implemented inseveral GSHP design programs, such as Earth En-ergyDesigner(EE

15、D;Blombergetal.2008)andGH-LEPRO (International Ground Source Heat PumpAssociation IGSHPA 2007).Heat flow from pipes in aborehole to groundFigure 1 shows the most common case of twopipes,aU-pipeinaborehole(N = 2).Inthegeneralcase, there are N pipes with midpoints at (xn,yn)inaborehole(N = 1,2, .).The

16、thermalconductivityin the ground outside the borehole is (W/(m,K),and it is bfor the grout in the borehole outside thepipes. The borehole radius is rb. The outer radiusof pipe n is rpn, and Rpnis the thermal resistancefrom pipe fluid to the grout adjacent to the pipe (permeter pipe, K/(W/m).Theheatf

17、lowfrompipenisdenotedbyqn(W/m).ThecorrespondingfluidtemperatureisTfn.Thegen-eral objective is to present a very precise methodto calculate the relations between heat flows andfluidtemperatures.Therelationswillberepresentedgraphically as a thermal network.Thermal problemSteady-state heat conduction f

18、or a cross-sectionperpendicular to the borehole axis and its surround-ing ground with N heat-carrier pipes is considered.The primary input data areb,rb, N, J;n = 1,.,N: xn, yn,rpn, Rpn,qn. (1)Here, J is the number of multipoles at each pipethat is used in the solution. For J = 0, only linesources at

19、 the pipes are used. The primary outputis the fluid temperatures Tfnin the pipes for any setof prescribed heat flows qn. A linear equation sys-tembetweentheseheatflowsandfluidtemperaturesDownloaded by T2rb () dTavdrvextendsinglevextendsinglevextendsinglevextendsingler=rb= qb=Nsummationdisplayn=1qn.

20、(16)The (average) heat flow qbthrough the boreholewall must be equal to the sum of the heat flows qn.As the outer boundary condition in the ground,the average temperature at the borehole wall cannow be prescribed:Tb,av=12integraldisplayTrad(rb,)d. (17)General temperature fieldThetotaltemperaturefiel

21、dwillbeobtainedbythesuperposition of heat sources and so-called multi-poles at the centers of the pipes. All these compo-nents satisfy the heat equation.Heat sources and multipolesAfirststepistointroduceheatsourcesandmulti-poles inCartesian coordinates, in polar coordinates,and in complex form. The

22、following relations areused for complex logarithms and powers of z:ln(x + i y) = ln(z) = ln(|z|) + i arg(z),Re(ln(z) = ln(|z|), Im(ln(z) = arg(z)ln(r ei) = ln(r)+i ;1z=zz z=x i yx2+ y2,1zj=eijrj=cos(j ) i sin(j )rj. (18)A heat source with the strength qnat (0,0) in theborehole with the thermal condu

23、ctivity bof thefilling material gives the temperatureT(x, y) =qn2blnparenleftBiggr0radicalbigx2+ y2parenrightBigg=qn2blnparenleftBigr0rparenrightBig=qn2bRebracketleftbigglnparenleftbiggr0zparenrightbiggbracketrightbigg.(19)Here,r0isanyradiussothattheargumentofthelog-arithm is dimensionless. Any cons

24、tant temperaturemay be added. In particular, the right-hand complexform will be used. The corresponding solution for aheat source at the center of pipe n isT(x, y) =qn2b Rebracketleftbigglnparenleftbiggr0z znparenrightbiggbracketrightbigg. (20)General angular-dependent solutions for two-dimensional,

25、 steady-state heat flow in the infiniteregionoutsideapipearealsoneeded.FromCarslawDownloaded by Tln(r), 1. (21)With these solutions (and ln(r) and a constant), anyheat flow problem outside a circle r = r0with giventemperature T(r0, ) may be solved using a Fourierexpansion of T(r0, ). A complex combi

26、nation isused in the following way:T(z) =cos(j ) i sin(j )rj=eijrj=1parenleftbigr eiparenrightbigj=1zj. (22)Here, both the real and the imaginary part of thecomplex functionT(z) are solutions to the heatequation.Anycombinationofthetwosolutionsmaybe used and the center moved to zn:T(x, y) = Rebracket

27、leftbiggPn,j(z zn)jbracketrightbigg,Pn,j= Re(Pn,j) + i Im(Pn,j),n = 1,. N, j = 1, 2, (23)These solutions are multipoles of order j.ThecomplexnumberPn,jisthestrengthofthemultipole.Heat sources and multipoles for thecomposite regionFrom Equations 20 and 23, one has the complexheat source and multipole

28、 solutions from pipe n ofthe following types:T(x, y) = Rebracketleftbigglnparenleftbiggrbz znparenrightbiggbracketrightbigg,parenleftbiggrpnz znparenrightbiggj,j = 1, 2, (24)All of these functions (real and imaginary partsfor multipoles) are solutions in the infinite regionoutside the pipe center z

29、= znfor constant ther-mal conductivity. However, the considered problemconcerns the composite region (Equation 4). Theheat sources and multipoles must be extended tocorresponding solutions for the composite region.These solutions are studied in Claesson (2011). Thederivations are not reported here s

30、ince they requiretoo much space.The basic solutions for a heat source and for amultipole of order j at z = znbecomeW0(z, zn)=lnparenleftbiggrbz znparenrightbigg+ lnparenleftbiggr2br2b z znparenrightbigg|z| rb(1 + )lnparenleftbiggrbz znparenrightbigg+ 1 + 1 lnparenleftbiggrbzparenrightbigg|z| rb, =b

31、b+ .z = x i y, (25)Wj(z, zn)=parenleftbiggrpnz znparenrightbiggj+ parenleftbiggrpn zr2b z znparenrightbiggj|z| rb(1 + )parenleftbiggrpnz znparenrightbiggj|z| rb,j = 1, 2, (26)Here, is a dimensionless parameter for the twothermal conductivities. It is zero for equal thermalconductivites, and it lies

32、between 1 and +1. Thereal part of W0, and both real and imaginary partsof Wjfor any j, satisfythe heat equation of Equation6 and the boundary conditions in Equation 9 at theborehole wall. The condition in Equation 17 withzeroaveragetemperatureattheboreholewallisalsosatisfied. The following set of so

33、lutions is obtained:U0(x, y) = ReW0(z, zn),URe,j(x, y) = ReWj(z, zn),UIm,j(x, y) = ImWj(z, zn). (27)The second logarithm in Equation 25, top, is aheat source at point z = r2b/znoutside the boreholewith the relative strength . This mirror point liesonthesameradiallinefromthecoordinateoriginaszn, and

34、the product of the distances from the originis |zn| r2b/|zn| = r2b. The second and third lines ofEquation 25 are the solution outside the borehole.The relative strength of the original line source is1 + . Another line source with a suitable relativeDownloaded by T only the state vectors T(t) will be

35、different in the simulations.Horizontal domain decomposition for distinctsoil regions is also used. These sub-zones are con-centrically structured so as to employ different timesteps, which generally leads to faster calculation.The simulation time step is adopted for the firstsub-zone,representedast

36、heinnermostcircleinFig-ure 2b, while larger time steps can be used for othersub-zones. Although the work reported here is forsingle boreholes, this horizontal sub-structuring isalso applicable to bore fields, as seen in Figure 2c.A fine unstructured mesh is used in the hor-izontal direction, as regu

37、lar polar and Cartesianmeshes cannot precisely describe the circular ge-ometries found in U-tubes and boreholes. The De-launay triangulation has been adopted in this work.An example of the mesh, generated using Gmsh(www.geuz.org/gmsh), is shown in Figure 3. Themesh generator delivers mesh informatio

38、n (i.e., amesh file, .msh) so that the grid information can beintegrated in the code with an additional identifica-tion procedure.Different numerical methods are adopted in thevertical and horizontal directions. In the horizon-tal direction, the finite-element method (FEM) hasbeen applied using an u

39、nstructured mesh, and thefinite-volume method (FVM) is adopted for the ver-ticaldirection.Theseselectionshavebeenmotivatedby numerical considerations: the FEM associated tounstructured mesh guarantees reciprocal relationsof heat transfer between nodes, and extrusion in thevertical direction generate

40、s orthogonal faces, alsoassuring reciprocity even if the FVM is used. TheFVM was selected simply because of its ease ofuse and implementation. Reciprocity in the equa-tions leads to the symmetry of state matrices of themodels in each sub-domain. Then, the symmetryof the matrices ensures real-type ei

41、genvalues andeigenvectors in the state matrices.State model representationThe governing equation for the conductive heattransfer can be expressed asctT div(T) = 0, (2)where is the density, c the specific heat, and thethermal conductivity of a conductive medium.Whether itisthe FEM orFVMbeing used,Equ

42、a-tion2can beconverted intoalgebraicform,thestatemodel, which expresses the energy conservation ateach node of the mesh:CT(t) = AT(t) + BprimeU(t), (3)where T(t), the approximated state vector (dimen-sion: n, 1), is the assembly of the time-dependenttemperature nodes of the system; that is, T(t) con

43、-tains all the node temperatures, the locations ofwhicharedefinedbythespatialdiscretizationT(t)=T1(t), T2(t), T3(t), ., Tn(t), where 1, 2, 3,., n indicate the node numbers in the mesh gen-eration. Matrix C, the square capacitance matrixFigure 3. Example of horizontal meshes generated by Gmsh: (a) fr

44、ontier between sub-zones n 1andn and (b) magnification of thehalf-circle of the first sub-zone (color figure available online).Downloaded by TMenezo et al. 2002; Gao et al. 2008) over the last 30years. One of the first formulations was proposedby Marshall (1966) in the field of automatic con-trol. T

45、he following paragraphs describe the processof one of the reduction techniques; this process hasbeen carried out for each sub-zone obtained by thehorizontal and vertical decomposition.The main idea of size reduction is to use onlya few modes considered “dominant” among X(t)inFigure 4. Responses of m

46、odes in Marshalls method: (a) slowand (b) fast.reconstructing the observed valuesY and truncatingthe other non-dominant modes. Equation 5 can berearranged to represent the dominant Xm(t) and thenon-dominant Xn-m(t) modes, as shown in Equation6:braceleftbiggXm(t)Xnm(t)bracerightbigg=bracketleftbiggWm

47、00 WnmbracketrightbiggbraceleftbiggXm(t)Xnm(t)bracerightbigg+bracketleftbiggBmBnmbracketrightbiggU(t)Y(t) =bracketleftbigOmega1mOmega1nmbracketrightbigbraceleftbiggXm(t)Xnm(t)bracerightbigg+ DU(t).(6)In the “Marshall” method (Marshall 1966), non-dominant modes are regarded as the fastest modesof the

48、 system. The essence of the reduction tech-nique can be represented schematically in Figure4. Since the eigenvalues of Wmare associated withlarge time constants (slow mode) of the system, atypical response of any element of Xm(t)toastepchange in U(t) may be represented schematically ascurve a in Fig

49、ure 4, while the response of any ele-ment of Xn-m, which is associated with the smallertime constants (fast mode), may be represented ascurve b in Figure 4. Since this response reaches itssteady state after a short time, it may be approxi-mated by an instantaneous step change as shown bycurve c in Figure 4.The“dominantmodes”canalsobeselectedusi