етодические указания по расчету на ЭЦВМ температурных попей при замораживании горных пород

Белгород 1980

МИНИСТЕРСТВО ЧЕРНОЙ ШШЗФГВВ ССОР управление горного правэводева

Всесоюзный иаучно^исояедованешшнй в проевшо-ношз^рувБорсвдШ вногигу® по ооушенив ыесэорэадевйй полоза® нокопа8вж,впеаЕаа&Ш4а гориш рабогаи,рудничной гаоаахта i маргвдйдарадезд логу В 0 О I* S й

МЕТОДИЧЕСКИЕ ШШШ

ПО РАСЧЕТУ НА ЭЦВМ ЭДООДЕИДО ПОЛЕЙ ПРИ ЗАМОРАЖИВАНИИ РОРШ ИООД

Белгород I960

ТМАХ ” иаг по времени для выдачи на печать результа -sob расчета (тмдх * К*тЫ» где К - целое чиоло шагов по времени, к * 1| 2, ... а. )| 45

Ткон ~ полное время расчета, 45 МЫ - чиоло узлов}

ММ - чиоло элементов}

тт - гемпература фазового перехода, °0{

ТФт Qop- теплота фазового перехода, ккал/м³;

Нв “ помер варианта»

Ортер подготовки информационной карты.

ТЫ * Ю Ч| ТМАХ* 20 4} ТКОН* 500 ч; ЫЫ * 230} ММ* 403; ТТ* 0°С} ТФП* 25000 ккал/м³; НВ= I.

	+	2	i	0
*	+	2	2	0
+	+	3	5	0	0
*	4*	3	23	0	0
+	*	3	403	0	0
*	*	0	0
*		5	25	0	0
4*	+	I	I

2.    Массивы Xx(NN)* YV (NN) - масоивы прямоугольных коор » д'шат х и у рассматриваемой области. Количество данных в каждом маосиве равно количеству узлов МЫ. Массивы данных записываются на бланке для передачи на ВЦ в нормализованной форме.

3.    Маосив рр ( МЫ) - значение граничных уоловий,формнруетоя по

узлам сетки, причен если в узле задана температура, то рр (I)    =

Т (I ), в противном олучае - рр (I ) « о (I), где (I ) - тепло-вой поток.

А. Массив ун (NH ) узлам сетки).

5.    Маооив ТТГ (мм) по числу элементов

6.    Маооив тмг(мм) ется по элементам.

7.    Массив ЛТ(ММ) витоя по элементам.

8.    Массив ЛЗ (ММ) витоя по элементам*

9.    Массив мр (ЫМ) - признаки граничных условий» вводится фор ватным вводов но формату Ц . При подготовке исходных данных запя-оавается на бланках ФОРТРАН. При ЫР(1, ) =* 0 з узел с номером ь задается тепловой поток величиной РР(ъ)« о при МР ( I ) «= I г ;,У85 узел задаетоя значение температура рр( I).

10.    Массив Им( 3» ММ ) - номера узлов, окружающих данный треугольный элемент. Массив записывается па бланках. ФОРТРАН по фор -мату 39.

Йля расчета по программе ХОЛОД I комплектуется колода перфокарт с директивами.0СЧИ220» которые вызывают программу» записан -ную на магнитной ленте. После вызова программы вводятся походные данные. Для выдача узловых значений ffeanepasypa на перфокарты необходимо на Н.8Л 8ЦВМ ШЗНЙЗ набрать произвольный код, отличный от нуля. Выданная колода перфокарт ставится в читающее устройство и нажатием клавши "Пуск” производится фиктивный ввод о контрольным суммированием. При несовпадении контрольной суммы нажат нем ^иПуок^й вывод повторяется. Массивы 2, - 8. аапвоываются яа бланках в нормализованной форме. Яйле приводятся распечатка программы холод I.

program ХОЛОД!

DIMENSION ХХ*УУП7Н(500)

333    FORMAT(1 ОХ/ШИРИНА ПОЛОШ РАША%ХЮ)

шшаь чтзу,зтау*у1Ш>Ф

РР*2ИУР \ У Q* гшга; КМ- 2HNM; НР= 2HNP; № 2НЛТ; ЛЗ =*2НЛЗ ;ТТГ=3 НТТГ;

CALL L0A»§O«2.S?,bPP»l*5OQ,yQ_tl_?5OO,NM_}l_?8OO,NP,i,5OO^a-,l_t 800,ЛЗ i

«beOO^nr^^OOtWf^tiOOjSMii^S&OO^HRMO^y)

CALL LOADGO Oft*УУ,»,РР ^NIS^N,TMAX,TKOH_tM,Ш_}TT_fM\*

ТМГ,ТФП»ЛТ

* »ЛЗ ,ra,NB,4T3y ₅ЗЙЗУ *1Н»ВОДТ)

T=TN

CALL LOADGO( NM _# Ш, NN, KW , ЧТЭУ * SHPOLOiA)

№TMAX;T«TN

KHm

GOTO 10078 10086 CONTINUE

CALL LOADGO(XX,yy_>¥Q,VH_iSM_>TN,KW,NN*ММ*ТТ|ТТГ,ТЫГ₈ТФН»ЛТ^ ЛЗ,ЧТЗ

*У _эЭПЗУ _&6HMATRIX)

DO 10150 Ь=1,Ш

ШЧТЗУ(№,1) 310150,10149 >10150

10149    И^ЧТЗУ(РР,I);call зпзу(чтзу<уа,1)+ндо,1)

10150    CONTINUE

CALL LOADGO(SM,NN_#YQ_tKW _tPP,NP,4T3y,ЗПЗУ,6НМОДИФ13 CALL LOADGO(YQ^SMjKWjNNДН,ЧТЭУ |ЗЛЗУ* 5ЙГАУСС)

кн*кк^1;т-тт

XP(I-TMAX)10086,10086,10087

10087    T=T-TN

10078 CALL LOADGO(KK,T_>YH,NN,yni3^,6HnE4001)

IF(KK)1009 9» 10086₁Ю099 10059 TMAX=TMAX+T П T=T+TN

IF(TMAX-TKOH)10086 j10086 * 10088

10088    STOP

END

subroutine вводТ(хх,уу,вр,рр,тн,нм*тн,тиАх,ткон,ня,им,тт

ттг,тмг,т

*ФП,ЛТ,ЛЭ,ТП,НВ,ЧТЗУ,ЗПЗУ)

DIMENSION ХХ,УУДН(1),»Ф(60) ,NR(t2)

10006 FORMAT С 5бХ» "ВАРИАНТ* ,33C,I4///)

10003 FORMAT!'ШАГ СЧЕТА ПО ВРЕМЕНИ' ,F8.1/*MAF П01АТИ ПО ВРЕМР"’ F0.1/'H

*СШЯНОЕ РАСЧЕТ. ВРНЯ* ,F8.1/'tlJ4CAO УЗЛОВ’,15/'ЧИСЛО ТТ*' ГОЛЕИИКОВ*

^w j15/"S ЗШПЕРАТУРА ТАЯНИЯ ГРУНТА,ГРАД С' ,F7.2/‘ТЕПЛОТА

вого я®

* ВХОДА* ,F10.l/)

770    FORMAT Ш)

771    FORMATilf)

772    FORMAT Ш0)

CALL ВВОД!(flf/jTMAX ДКОН,Ш,ИМ,ТТ,ТФП,ИВ)

PRINT 10003    ,ТИОН,НИ,ТГ    ,Т#П

CALL ШСЩШМН)

CALL ВВОДМ(УУ,ЮО CALL

DO 2221*1,N#

222    CALL ЗПЗУ<УИШ,ГМ)

CALL ВВОДИ(YH.NN)

1ПЕ=1

39    1Р(1ПЕ-5>30,63,63

38    IF(MM-60)40,41,41

40    K«»HM;00 TO 42

41    M*60

42    Ы

13    CALL ИВСДМ(**,К#);О0 TOtSO,61,611,62),1ПЕ

60    DO 6 <7=1,КФ

6    CALL ЭПЭУ(ФФ(<7),ТТГ,1+<7-1)

T    IF(t+60-MM)11,11,12

11    IF(I+120-MM)15,15,16

15    I*I+60;GO TO 13

!6    1=1+6 0(K#4U=!*i

00 ТЯ 13

13

34	12	1ПЕ=1ПВ+1 ;GQ ТО 39
35	61	DO 53 J=l,K<j>
36	53	CALL ЗПЗУ(*Ф(0),ТМГ,1+0-1)
37		GO ТО 7
38	6U	DO 533 «М,К®
39	533	CALL ЗПЗУ(**№),ЛТ,1+<1-1)
40		QO ТО 7
41	62	DO 54J=l,KO
42	54	сам. 3n3y(«(j),jja,i+j-i)
43		СЗО ТО 7
44	63	1=1
45	47	READ 770,(jNH(J),J=1,72)
46		PRINT 772,(NR<Ll),Ll=i,Tg5
47		DO 44 J=l,72
48	44	CALL 3tI3y(NR(J),NP8I+J“1)
49		IF<I+72-NN) 45,45,46
50	45	1=1+72}QO TO 47
51	46	in
52	5 \	READ 771,(NR(J),J=l,a)
53		DO 48 J=l,8
54	48	CALL 3n3y(RR(J),NM,I+J-l)
55		Ш1+8-ММ) 49,49,50
5-6	49	I=I+8jGO TO 51
57	50	PRINT 10006,NB
58		RETURN
59		END
1		SUBROUTINE М0ДИФ1(SM,NN, YQ,KW,PP,NP,ЧТЗУ,ЗПЗУ)
г		DIMENSION PPU)
3		DO 1 1=1,NN
4		IF(ЧТЗУ(NP,I))1,1,2
5	г	L=(I-1)*KW+1
ё		CALL ЗПЭУ(ЧТЗУ(SM,L)*10**10,SM,L)
7		CALL ЭПЗУСЧТЗУ(ЗМ,Ь)*РР(1),У0,1)
8	1	CONTINUE
9		RETURN; END

SUBROUTINE MATWX{XX,yy^Q,W»SM,TN,IW,NN,NM,»W,TT,T'Tr_! T'<' Т«!,ЛХ,

»ла,чтэу,8П8У)

DIMENSION Х,У,Н,МВ(Э),А,АА,В(3,Э),Ю,РР,ХХ,УУ,¥Н(1)

10002 FCmTUX,'ПЛОЩАДЬ ЭЛЕМЕНТА*,2Х»14»‘РАВНА НУЛЮ*///) KW=KW*NN 10086 Ш 11125 I“1,NN шгз саы. эпэу(о.,те,1)

DO 10027 1*1,К?

10027    САЫ. ЭПЗУ(0.,8М,1)

DO 10148 1*1,ММ

10085 Ш1*ЧТЭУ(КЦ,1)

NH(1)*NM1/I00000d NRl=NMl-NH(l)f1000000 NR(2)“NRl/1000 HR(3)*NB1-NH(2)*1000 DO 10020 J=l,3 NPN=NR(J)

H(J)*yH(NPN)

X(J)=XX(NPN)

10028    y(J)=yy(HPH) x2i*x(2)>x(i)}Xi3*x(i)-x{3hW2=x{3)‘-x(2)?yi2*y(t)-?{?] угз=7(2)-У(3)|У31“У(3)“У(1) в^(хп*У31‘-хп*У1г)*г

IF(S)10029,10030,10129

10029    S“-S

GO TO 10129

10030    PRINT 10002,1 STOP

10129 HCP«(H(l)+H(2)+H(3))/3 ШНСР-ТТ)Ш8,1978,1977

1977    0Т*ЧТЗУ<ТТГ,1)|ЛТ1»ЧТЗУ(ЛТ,1)

ФК1«ЛТ1*2/3;D»QT*S/14 4/TN

GO TO 1979

1978    Ш»ЧКГУ(ТМГ,1)1ЛЭ1*ЧТ9У(Л8,1)

#K1«ЛЭ1* 2/3 5 D“(QM+ »tl/( 1+(TT-HCP))*■*2) »S/144 /ТН

1979    CONTINUE

15

Аа_>2)*$К1*<Ш*У31+ХЭ2*Х13)/3 А С1*3 ) «Ш*( У23 *У12+Х32*Х21) /3 А{1,1)=~А(1_#2)~А(1,3) А(2,3)*ФК1МУ31*У12+Х13*Хг1)/3 A(2,1)*AU,2)

A(2_>2)*-<A(2_flHA(2_>3))

A(3,1>»A(1_s3)

AO»2)^KA<2*3>

A(3|3)^B-(A<3,1)+A<3*2>)

AA<l,l)*A(l_tl)*14*D

AA(2_t2)»A(2,2)*H*P

AA(3,3>^BA(3»3)+U*X>

АА(М)*АА(2,1}*А(1,2)+5*В

AA(l,3)^BAA(3,l)^BA<l,3)+5*I>

AA(2 ,3)^(3,2)*A(2„3)+5*B

В(1Д>*14*Р~А(1_?1)/2

B(2,2)*14*B-A(2*2)/2

B<3,3)=14*D-A(3,3)/2

B(1,2)^(2 _tl)*5*B-A(l_t2)/2

B(l_#3)*B(3_sl)^a5«D-A(l,3)/2

ВСг^^вО^^в-А^зТ/г

ВО 10042 Jn,3

NPN«NR(J)

Ш 10042 K=l,3 №>II-NR(K)

10042    CAUL ЗПЗУ(ЧТЭУ(¥Q,NPN)+B(G *К)*УН(MPM),YQ,WPH}

10043    DO 10148 J*l,3

DO 10148 К=^,3 N2-NR(K)

IF(N1~N2)10044 _#1Q044₍1004 5 i b044 L={N1~1 )*KW+N2~Nm GOTO 10046 1004 5 L*(N2-l)*KSr+Nl~N2 + l 10046 Sl=4T3y(Stt,L)

10И8 CALL 3n3y{Sl+AA(J_JK)_#SM_>L)

RETURN; END

U9 к»-т

Ьг*Ы+и 1 №b*KlMKW~t) *i+i 51    014(Ъ2) *ЧТЗУ{SM, JJ)

goto юо 200    №= 1

YH(KN)*yyQ(t)/0Ma)

5Ю 202

YlHs4T3y(YQ >Ш-1)

00 203 J*1,ХЛ? JJ*KW*(NN-l-l)+t+J

203    YlH=YlH*4T3y(SM,JJ)*YH(J-I+HN) YH0W-I)«Y1H TF(KW-LR-1)202,202,204

204    LB=LR+1 202 CONTINUE

REFORM

ть

SUBROUTINE ?OLOSmM*lW>HN,KV^qray)

DIMENSION NB(3)

К|Г=0

DO 10124 1*1,MM

MAX = 0; MINrl 000 ;NJ4l=4T3y(fOfi, I); NR( 1) =N111/100 0000

NR1=NM1-NR{ 1) *1000000

NR(2)“NB1/1000;RR(3)*NfU-NR(2)*l 00 0

DO 10021 J-1,3

NA*NR(J)

IF(NA-MINH0019,10020,10020

10019    MIN“NA

10020    IF(N4-MAX)10022,10022,10021

10021    MAX=NA

10022    IF(МАХ-МШ-KW) 10024,10023,10023

10023    KW=MAX-MIN+1 10124 CONTINUE

RETURN

END

1		ттттж ввч«нйк,т,?н»И1,у!ш**)
8		М9Ш91Ш *И<1>
3	new	ммм*(и,'гмт* на ерщр т*1 .рю.г/ часов1///)
4	10008	wjw»f( m,' начальное г&тиятюъ гшг-)
5	10003	УвНЙЯТ(1Х,5(15,4Х,РП.2,1Х))
б	10000	РОШИКIHt, 118{18= )///1Х ,5(2?Н: УЗЕЛ : ТЕМПЕРАТУРА ) 3X»im///il
1		»§(№))
8		IFOBO10079 »»Ов7в,100 7 9
9	10018	РИШТ 10008)009010108
10	10079	PRIRT Н00Й,Т
11	10108	РИЙТ 10080
12		1*1
13	10081	m кда( I) | вг*ш s+ i), из «тпкх+г) i w *лис i+э) ? К5 *лн < т * о
14		11*1+1;12»1+2;13*1+3;Т4*1+4
15		print юооз,s,И1 ,н,нг,хг„щ,хз,Н4,и ,т
16		S«I+5;XP(l-mr)l0061t1008i,10062
И	10082	сдатшл
18		call упёрфскзу)
19		гр«кзу)1,г,1
20	1	CALL ПЕРФ(1йС1),7Н(гда)>
21	г	ВЕТОК»
22		т>

Литература

1.    Методические указания so пряыеневив метода конечных впеиен-sob дня решения плановых задач фяльтрапки подземных под на ЭЦВМ. Белгород, ВИОМ, 1979, 50 с.

2.    Сегерливд Л. Приивненйе метода конечных элементов. "Мир", М., 1979, 392 С,

3.    Тихонов А,Н», Сакарский А,А. Уравнения математической физики. М., "Наука”, 1972, 736 с.

Настоящие аетедячесйке указания праанаамадвы дня чш~ данного аодеайрош» нотодок извечных ътттт (ВК8) двумзршк еаяач йвоадионаршей тепаонроводяоогг ш ЭЦШ« Приведена шсвдкцш т айшгауатацщ ирогршш’ Ш0Д 2, аапиоанной га язнге WPTPAH-2 ванд.§из,"Ш$,наук Ваеаяа «> ввын В.А. и яащигйхн.ааук сеаеаневык Н.А. Яра отвадке и апробация програюш прашмаяк участке Овадчяй й.Ф.* Иу* равад» ТЛ» ш Веня® Т*Е«

Работа утверждена секцией шуч&о-явштттт совета ввстягута ВЙОШ В яюая 1980 г» й качестве вееодочве » вас укааашй*

Всесоюзный научно-исследовательский и проектно-конструкторский ия^ пгитут по осушению месторождений полезных ископаемых* специальным гор -ным работам, рудничной геологии и маркшейдерскому делу (ВИОГЕМ)* 1980,

6028?148КВ

Введение . . . , , , . . . . * » . . , , . , , . , , .    3

1.    Постановка задачи .................

2. Харктериотиве прогрмбш . * , . » . „ , „ _е , . »    9

3. Инструкция но подготовке исходной информаци* . . .    9

МЕТОДИЧЕСКИ! УКАЗАНИЯ 80 РАСЧЕТУ НА ЭЦШ ТЩЕРАТШЫХ ПОЛЕЙ ПРИ ЗАМОРАЖИВАНИЙ ГОРНЫХ ПОРОД

Научный редактор канд.техи.наук С.Г.Аксенов

Литературный редактор Л.А.Порубай Технический редактор А.Г.Ворнцова Корректор й.А.Соянр

Подписано к печати 2? июня 1980 г.

Объем 1,1 уч.-изд.а. Тира ISO ев8, 8акаа 1 §89» Ротапринт ВВОШ, г_аБенгород, ул.Б.Хйваьнйцкого,, 86. Цена 18 коп.

ВВЕДЕНИЕ

Искусственное замораживание грунтов находит широкое применение при проходке горних выработок в слоеных гидрогеологических уоло -ваях, а также при сооружении промышленных объектов глубокого за ~ ложения. Решение инженерных задач при этом требует звания распределения температуры в промерзающих грунтах для расчета скорости промерзания и оценка толщины формируемой ледопородной стенки. Задача теплопроводности с учетом фазового перехода, характерная для промерзающих грунтов, при произвольном размещении замораживающих колонок и жешшфизической неоднородности горных пород, не имеет аналитического решения» Сложность проблемы и в то же время ее исключительная важность обусловили широкое применение методов математического моделирования тепдофизических процессов»

Среди различных методов моделирования в настоящее время полу -чает наибольшее развитие метод вычислительного эксперимента, основанный на использовании вычислительной математики и средств вычислительной техники, причем для влажных задач наиболее эффективный оказывается метод конечных элементов (МКЗ) [I, 2].

Алгоритм теплофизичеокого процесса и разработанная на его ос -ново программа для ЭЦВМ фигурируют в качестве постоянно действу -кидай математической модели, которую несложно приспособить для конкретного объекта. Как показала практика использования МКЭ в институте ВйОГЕМ [l], постоянно действующие модели нестационарного температурного поля позволяют успешно решать сложные теплофи -зичесяие задачи замораживания. Цель настоящей работы - изложить методику использования МКЭ в двумерных задачах промерзания горных пород.

I. ПОСТАНОВКА ЗАДАЧИ

йетод конечных элементов успешно применяемся при моделирований на ЭЦВМ процеооа двумерной нестационарной теплопроводности. Рас -прострпнввие тепла в промерзающих или протаивающих грунтах описываемся дифференциальным уравнением пораболяческого мила

где Т - температура; время! - координаты» Сэ- зффек -тивная теплоемкость? X ~ коэффициент теплопроводности.

В общем случае теплоемкость и коэффициент теплопроводности являются функциями координат и температуры» следовательно» уравнение (I) - нелинейное. Для решения уравнения задают краевые уело -вия;

1.    Начальное распределение температуры внутри    рассматриваемой

плоской области

Т(*,у,0)« То (х,1|)    ,    ₍₂)

где ">о(Х,у) - известная функция координат,

2.    Граничное условие первого рода, состоящее в задании темпе -ратуры на некоторой части границы области в любой момент времени

^т(М,‘^с)^в‘ ^т"(*»ул) ,    (3)

где Т_п~ температура на границе •

3.    Если на границе обяаоти Sg задан поток тепла q (условие П рода) или конвективный теплообмен (условие Ш рода), то граничное условие записывается в следующем виде [2}.:

A(^rad    T.it)+q+*(T-T_c^^eO ,    W

где с( — козффициент теплообмена; Т - температура на границе (неизвестная)! 1^- температура окружающей среды или теплоноси -теля в колонках! К - внешняя нормаль к поверхности.

Поток тепла q считается положительным» если тепло теряется телом, причем q, и конвективная потеря тепла о( (T-Тс) не задаются на одной участке границы одновременно. Если задано условие И ро -да» то q = О, если задано q, то на этом участке <*. = о. Если одновременно q - 0 и о( о о, то условие (А) сводится к соотноие -

НИИ

И

Эт

которое выдавав? уоловив теплоизоляция.

Объединение S₄« в Sj образует полную границу рассматриваемой области. Ори решении задачи методой конечных элементов изучаемая область разбивается на конечное число МИ треугольных елеыентов произвольной формы. Все треугольники области нумеруются от I до ММ« а вое вершины ах (узлы) - от I до МЫ (рисунок). Для про *■ невольного злемен.а с номерами вершин 1₍|,к основное уравнение метода конечных элементов имеет вид [I]

М * {т} ⁺ [rn]« |c{^THf} I

причем элементы матриц определяются выразенияии

%i~ J^A (Ix^ ⁺ if |f-)dxd^ +

Д    $2:

тц = jCaMiMj dxd^ ;

A

ji * J c^rtids - J<* Tc Mids ₅

$a    5*

где Ы_п=2д(ап ⁺ 6_г,х + с_ггу)₎ п = ц^,к . Функции Ь\_п называются функциями формы [2]

a_L = x^_K-x_K^j    a_K=Xi.y_rxjy_Li

*,

с 1 = Х_К-Х_а-; Cj=X_L~X_K- C_K=Xj-li,;

Д= ^ [ (XL'XjKy^y^)-(Xi,~X_K)(tj_l--ур] .

Система (б) представляет собой систему трех обыкновенных дифференциальных уравнений» Рассматривая расчетный интервал времени

д

и заменяя производную от температуры по времени конечно-разностныы отношением, о помощью метода Галернина получим следующую систему линейных алгебраических уравнений:

(к ["*!    ⁺

f^fl«) ¹

Рис. Пример траангуляажи области:

I. 2, 3, ... - номера узлов;ф,@»(|), ... - камера элементов; • - замораживавдда кояовкв.

где векторы {т} в {т}₀- влечения температуры в вершинах элемента ва моменты вреиени 'С = д'О в *С * О соответственно.

Матрица теплоемкости [т] в развернутой виде имеет структуру

["Ml

Суммируя уравнения (II) по воем элементам области, подучим результирующую систему, которая имеет вид г

[*H^tU"[^bM^t}ct,,>⁺M ■ _т

Матрица [А] являетоя комбинацией матриц [ср и trr\], зависит от шага но времени Д'С. Система линейных алгебраических уравне н и й обладает свойствами, которые делают ее весьма удобной для решения, а именно, матрица системы САЗ симметрична и положительно опреде -ленная, причем имеет ленточную структуру. Последнее обстоятельство позволяет значительно сократить объем намята ЭЦВМ.

Конечно-разностная схема, используемая для временной аппроксимации уравнения (6), являетоя безусловно устойчивой. Безусловная устойчивость означает, что волн распределение температуры по времени преобразовать по Фурье в частотную область., то коэффициент усиления для каждой частотной компоненты будет затухать во времени [23. Однако при этом обычно возникают колебания числовых зна -чений искомых величин температуры. Размах колебаний зави сит от теплофизичеоких свойств грунтов, размеров используемых треугольных элементов, величины шага по времени &V и значений фурье-ком-понента температурного распределения, соответствующих началу Поскольку теплофизические свойства грунтов обычно заданы (в пря -мой задаче), переменными, которыми можно варьировать, будут только размеры элемента и шаг по вреиени . Одновременное уменьше -ние размеров элемента и Д'С существенно снижает размах колебаний, тогда как изменение только одной из этих величин при фиксирован -ном значении другой не всегда уменьшает колебания. Не рекоменду -етоя сочетание грубой разбивки области на элементы и малого шага по времени. Опыт показывает, что такая комбинация может привести к результатам, противоречащим физичеокому смыслу гадачи. Алгоритм решения системы (13) прост и заключается в следующем: при известном тепловом состоянии грунта в рассматриваемой области на момент t - 0 вычисляются матрицы [А] и [В], а затем вектор правой части.

7

После этого» зная начальное распределение температуры {T}_ctap пуней решения системы (13)» вычисляется распределение температуры на расчетный промежуток времени {т}^. Этот процесс продолжает в я от исходного состояния до «оборо заданного момента времена. Яри учете фазового перехода» следует иметь в виду то оботоятельот -во» что корда температура приближается к температура фазового перехода» эквивалентная теплоемкость Сэ претерпевает резкое импульсное изменение» подобно дельта-функции, вояедотвии второ при за -давни аффективной объемной теплоемкости широко применяется метод усреднения» согласно которому распределение теплоты фазового т « рехода Qq, задается на некотором интервале температур ГзЗ

С_э-См +(ТЬ    <»)

где С_м - объемная теплоемкость замороженного грунта? |(т)~ неко -торая функция температуры.

Для талого грунта аффективная теплоемкость принимается равной объемной теплоемкости талого грунта

Сэ ⁼ Ст

Функцию j(т) обычно подбирают так» что<& распределение теплота фазового перехода происходило в небольшом интервале температур, Однако при малых АТ появляется опасность проскоков пиковых значений эффективной теплоемкости» что при численном решении задачи га ЭЦВМ приводит к осложнениям. Этот недостаток устраняется, если теплота фазового перехода распределяется по интервалу Tjfej^Tn.» где т₃ - температура замерзания воды» а Т_п~ температура на гра -нице (на стенках замораживающих колонок) на конечный момент вре -пени. При этом функция f(т) должна удовлетворять условию

где 'f - некоторый коэффициент, зависящий от вида функции j (Т) и температуры Тп • Очевидно» можно подобрать много функций, удов -летворяющих условию (15). Как показала практическая реализация такого подхода«удобна та-

кая функция ^7>=(Д-*)! ’ (16) где f = \ + ——г . т 1 Т,-Тп (П)

При указанном методе определения эффективной теплоемкости» уравнена» (X) можно реиать на ШШМ численным методом как обычное уравнение теплопроводности о переменными коэффициентами»

8а основании описанного алгоритма ооохавлена для ЭЦВМ программа ХОЛОД X.

2» ХАРАКТЕРИСТИКА ПРОГРАММЫ

Программе ХОЛОД I предназначена для моделирования яеохациоаар -йот температурного поля в горных породах и сосхоих иа ооновн о й програмш и иеохи подпрограмм. Весь пакех написан и отлажен на алгоритмическом языке ФОРТРАН, транслятор ф-20 для ЭЦВМ БЭСМ-дм.Оригинал хранится не перфокартах (ГОСТ €198-75) а фонде алгоритмов я программ отдела осушения инсхяхуха ВИОГВМ . Транслированная про -грамма вапиоана на магнитной ленте»

Максимальное число угловых точек ЫН = 500, максимальное число элементов-треугольников - 800»

Ввод исходной информации осуществляется подпрограммой ВВОДТ.

В результате расчета по программе печатается?

1)    номер варианта?

2)    число узловых точек области!

3)    число элементов?

4)    максимальная разность номеров узлов в треугольниках (ширина полосы матрицы теплопроводности)!

5)    начальное распределение температуры?

6)    расчетное время?

?) вначение температуры на расчетный момент времени для каждого узла.

Выдача информации по пп. б) и 7) повторяется.

Предусмотрена промежуточная выдача результата расчета значений узловых температур на перфокарты. При продолжении расчета полученный маосив данных используется в качестве начального распределения температур.

3. ИНСТРУКЦИЯ ПО ПОДГОТОВКЕ исходной ИНФОРМАЦИИ

Для расчетов по программе ХОЛОД I готовится оледующая информа -

ция. _    _____

I. Информационная карта состоит из следующей информации: тН- расчетный интервал времени, ч;

9