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

После закрепления номеров за каждой зоной краевых условий верхней границы заполняется табл. 3. Каждое значение этой таблицы имеет вид:

,

где KOL=0, если во всех элементах (i, j, K), где 1£i£NI, а j, K фиксированы, внутренние источники отсутствуют;

KOL=1, если внутренние источники присутствуют хотя бы в одном из элементов (i, j, K), где 1£i£NI, а j, K фиксированы;

NP - номер типа источника, NP=1,2,...,9. Если источник отсутствует NP=1;

MD - номер зоны краевых условий верхней границы, в которой находится элемент области исследования;

MDK - номер зоны краевых условий нижней границы для соответствующего элемента области исследования.

Таблица 3

j

K

K=1

K=2

· · ·

K =NK

j=1

KL (1, 1)

KL (1, 2)

· · ·

KL (1, NK)

j=2

KL (2, 1)

KL (2, 2)

· · ·

KL (2, NK)

·

·

·

·

·

·

·

·

·

·

·

·

j=NJ

KL (NJ, 1)

KL (NJ, 2)

· · ·

KL (NJ, NK )

Для двухмерной задачи табл. 3 состоит из одного столбца для K=1, для одномерной из одного числа KL(1,1).

3.5. Физические и теплофизические свойства Расчетная область может содержать до 10 литологических разностей, характеризующихся следующими свойствами, учитывающимися алгоритмом:

- плотность сухого грунта кг/м³;

- весовая влажность грунта, доли единицы;

- удельная теплоемкость скелета, кДж/(кг·°С), ккал/(кг·°С);

- коэффициент теплопроводности грунта в талом состоянии Вт/(м´ ´°С), ккал/(м·ч·°С);

U⁰ - температура начала фазовых переходов, °С;

W_W - количество незамерзшей воды в диапазоне отрицательных температур, доли, задаваемой кривой

,

где A, B, C - коэффициенты, определяющие количество незамерзшей воды в данной литологической разности.

Все вышеуказанные свойства должны быть определены для каждой литологической разности по лабораторным данным, справочным или нормативным документам. Они остаются неизменными для всего времени счета.

Составляется табл. 4. При этом каждой литологической разности присваивается постоянный номер от 1 до NSL, где NSL £ 10.

Таблица 4

Номер литологической разности

,

кг/м³

,

ккал/(кг·°С)

,

доли

, Вт/(м·°С),

ккал/(м·ч·°С)

, Вт/(м·°С),

ккал/(м·ч·°С)

, °С

W_W

доли

1

2

·

·

·

В графе W_W помещается информация о незамерзшей воде в данном слое. Если известна кривая, то есть уже известны коэффициенты A, B, C ее задающие, то записываются значения коэффициентов. Если известны лабораторные данные, то приводятся они. Например: и так далее не менее трех значений. Если в расчете W_W=0, то коэффициенты, задающие кривую незамерзшей воды в данном слое будут равны: .

После закрепления номеров за каждой литологической разностью заполняется табл. 5, структура которой повторяет табл. 1. Таблицей 5 за каждым элементом (i, j, K) , полученным при разбивке области исследования, закрепляется номер той литологической разности, которой он принадлежит.

Таблица 5

K=1

i

j

j =1

j =2

j =NJ

i=1

N (1, 1, 1)

N (1, 2, 1)

· · ·

N (1, NJ, 1)

i=2

N (2, 1, 1)

N (2, 2, 1)

· · ·

N (2, NJ, 1)

·

·

·

·

·

·

·

·

·

·

·

·

i =Ni

N (Ni, 1, 1)

N (1, 2, 1)

· · ·

N (Ni, NJ, 1)

K=2

i

j

j =1

j =2

j =NJ

i=1

N (1, 1, 2)

N (1, 2, 2)

· · ·

N (1, NJ, 2)

i=2

N (2, 1, 2)

N (2, 2, 2)

· · ·

N (2, NJ, 2)

·

·

·

·

·

·

·

·

·

·

·

·

i =Ni

N (Ni, 1, 2)

N (Ni, 2, 2)

· · ·

N (Ni, NJ, 2)

и т.д.

K=NK

i

j

j =1

j =2

j =NJ

i=1

N (1, 1, K)

N (1, 2, K)

· · ·

N (1, NJ, NK)

i=2

N (2, 1, K)

N (2, 2, K)

· · ·

N (2, NJ, NK)

·

·

·

·

·

·

·

·

·

·

·

·

i =Ni

N (Ni, 1, K)

N (Ni, 2, K)

· · ·

N (Ni, NJ, NK)

Каждое значение N(i, j, K) в табл. 5 может быть только целым числом от 1 до NSL, где NSL£10 - число выделенных литологических разностей.

3.6. Дополнительные данные, необходимые для проведения расчета Для проведения расчета необходимо также определить значение следующих величин, задействованных алгоритмом.

Шаг по времени. Размеры элементов и шаг по времени связаны между собой следующим соотношением:

где Dt - шаг по времени, ч;

i - номер грунтовой разности;

NSL - число выделенных грунтовых разностей;

- теплоемкость мерзлого грунта i-го слоя при , кДж/(м³·°С), ккал/(м³·°С);

- количество незамерзшей воды, доли;

- удельная теплоемкость сухого грунта, кДж/(кг·°С), ккал/(кг·°С);

- удельная теплоемкость льда, кДж/(кг·°С), ккал/(кг·°С);

n- коэффициент объемного расширения, n=1,1;

- весовая влажность грунта, доли;

- коэффициент теплопроводности мерзлого грунта, Вт/(м·°С), ккал/ /(м·ч·°С);

n - мерность области исследования, 1£n£3;

Dh - размер ребра элемента.

Шаг по времени может быть изменен в процессе счета NT раз, где NT³1.

Время измерения расчетного шага от начала расчета задается одномерным массивом размерностью .

где - время измерения шага от начала счета в часах.

Отсчет происходит от моментов начала счета, т.е. от 0.

Определяется Ткон (час) - время окончания решения задачи. Задаются Т₁, Т₂, Т₃, ..., Т_Р - время выдачи результатов расчета на печать (в часах). Время последней выдачи результата T_P обязательно меньше T_КОН на один шаг по времени. В противном случае счет может прекратиться раньше, чем будет получен последний результат.

Определяются интервалы, через которые происходит очередная выдача результатов на печать.

и т.д.

Для последующего ввода в ЭВМ формируется одновременный массив HPRI размерностью INT. Размерность INT массива зависит от того, с постоянным шагом выдаются результаты расчета или нет.

INT=2 - выдача производится с постоянным шагом;

INT=P+1 - выдача производится с переменным шагом.

Если INT=2, то

Если INT=P+1, то

и т.д.

Задание печати распределения температур

Для удобства обработки информации, получаемой в результате расчета, фиксируются необходимые для анализа температурного режима «линии резервов», по которым выдается на печать распределение температур.

Положение «разрезов», параллельных плоскости ZOX, задается координатой элементов, составляющих этот разрез.

Формируется массив целых чисел KP размерностью KW, где KW - число выдаваемых на печать «разрезов» по направлению . Если по этому направлению ни одного разреза выдавать на печать не надо, то KW=1 и KP(1)=0. В противном случае и т.д. , где

Положение «разрезов», параллельных плоскости ZOY, задается координатой j элементов, составляющих этот разрез.

Формируется массив целых чисел jP размерностью jW, где jW - число выдаваемых на печать «разрезов» по направлению . Если по этому направлению ни одного разреза выдавать на печать не надо, то jW=1 и jP(1)=0. В противном случае и т.д. где

Положение «разрезов», параллельных плоскости XOY, задается координатой i элементов, составляющих этот разрез.

Формируется массив целых чисел iP размерностью iW, где iW - число выдаваемых на печать «разрезов» по направлению .

Если по этому направлению ни одного разреза выдавать на печать не надо, то и . В противном случае и т.д. , где .

4. ПРОВЕДЕНИЕ РАСЧЕТА 4.1. Вызов и загрузка программы в ЭВМ Каждая из модификаций программы PROGNOZ состоит из нескольких модулей, объединенных главной программой MAIN. Подпрограммы-модули GRUND, WNS, HF, FIND содержатся в личной библиотеке исходных модулей (или библиотеке загрузочных модулей).

Ниже приводится пример пакета вызова программы PROGNOZ-3S из библиотеки исходных модулей. Задание требует не менее 300 килобайт памяти и использует оптимизирующий транслятор (ФОРТРАН-ОП) системы ОС ЕС.

Здесь:

имя - биб - имя личной библиотеки исходных модулей;

имя - тома - имя тома, на которой размещена библиотека.

// JOB

// EXEC FORTHCLG

// FORT. SYSIN DD DDNAME=PKAR

//                DD DSN=<ИМЯ - БИБ>(GRUND),DISP=SHR,

//                          VOL=SER=<ИМЯ-ТОМА>,UNIT=SYSDA

//                DD DSN=<ИМЯ - БИБ>(WNS),DISP=SHR,

//                          VOL=SER=<ИМЯ-ТОМА>,UNIT=SYSDA

////              DD DSN=<ИМЯ - БИБ>(HF),DISP=SHR,

//                          VOL=SER=<ИМЯ-ТОМА>,UNIT=SYSDA

//                DD DSN=<ИМЯ - БИБ>(FIND),DISP=SHR,

//                          VOL=SER=<ИМЯ-ТОМА>,UNIT=SYSDA

// PKAP DD^*

<текст программы MAIN на перфокартах>

<текст подпрограммы BLOCK DATA на перфокартах>

// GO. SYSIN DD^*

<исходные данные>

/_*

//

Текст подпрограммы для модификации PROGNOZ-3S приведен в обязательном приложении 3.

4.2. Составление программы MAIN и подпрограммы BLOCK DATA Программа MAIN в модификациях PROGNOZ-3S и PROGNOZ-2S составляется для каждого варианта расчета. Пример программы MAIN для расчета контрольного примера (приложение 1) приведен в обязательном приложении 3.

Программа MAIN содержит обращение к основной подпрограмме GRUND алгоритма. Перед обращением к подпрограмме GRUND в MAIN должны быть определены все формальные аргументы подпрограммы, приведенные в табл. 6.

Таблица 6

Идентификатор

Фактическое значение

NSL

Число выделенных литологических разностей (NSL£10)

NI

Число элементов по направлению i, (NI£100)

NJ

Число элементов по направлению j, (Nj£100)

NK

Число элементов по направлению K, (NK£100)

NT

Число изменений шага по времени (NT£10)

LI

Число интервалов времени в периоде, во время которых задаются средние значения характеристик внешней среды на верхней границе (Li£36)

M

Число зон с различными краевыми условиями на верхней границе (M£9)

M1

Число зон с различными краевыми условиями на нижней границе (M1£9)

TKON

Время окончания счета (в часах)

A1(NSL), B1(NSL)

Одномерные массивы размерностью NSL

C1(NSL)

Коэффициенты, задающие кривую незамерзшей воды по слоям

N(NI,NJ,NK)

Трехмерный массив размерностью NI´NJ´NK

Номера слоев грунта в соответствующих элементах (см. табл. 5)

KL(NJ,NK)

Двухмерный массив размерностью NJ´NK

Массив номеров зон краевых условий по элементам верхней и нижней границ (см. табл. 3)

U(NI,NJ,NK)

Трехмерный массив размерностью NI´NJ´NK

Начальная температура грунта по элементам (см. табл. 1)

IW

Число выдаваемых на печать «разрезов» по направлению i (см. п.3.6)

JW

Число выдаваемых на печать «разрезов» по направлению j (см. п.3.6)

KW

Число выдаваемых на печать «разрезов» по направлению K (см. п.3.6)

INT

Размер массива HPRI (см. п.3.6.)

HPRI(INT)

Одномерный массив размерности INT

Задает интервалы выдачи на печать температуры грунта (час) (см. п.3.6.)

H(NI,NG,NK)

Трехмерный массив размерностью NI´NJ´NK

Значение энтальпии по элементам на текущем временном слое

UN(NI,NG,NK)

Трехмерный массив размерностью NI´NJ´NK

Температура грунта по элементам на следующем временном слое

HNDV(NI,NG,NK)

Трехмерный массив размерностью NI´NJ´NK

Значение энтальпии по элементам на следующем временном слое

Nsi

Nsi =1, если единица измерения количества теплоты - кДж

Nsi ¹1, если единица измерения количества теплоты - ккал

Оператором DIMENSION определяется размерность массивов H, UN, HNDV, U, A1, B1, C1. Ввод массивов N, KL, U, HPRI осуществляется операторами описания типа (INTEGER, REAL), а также операторами присваивания.

Применение операторов описания типа для ввода указанных массивов, в подавляющем большинстве проводимых на практике расчетов, уменьшает количество информации, вводимой с перфокарт.

Обращение к подпрограммам WNS (NSL, A1, B1, C1) присутствует в том случае, если коэффициенты A1, B1, C1, определяющие кривую незамерзшей воды по грунтовым разностям, неизвестны. В противном случае эти коэффициенты вводятся операторами присваивания или операторами описания типа.

Обращение к подпрограмме GRUND присутствует обязательно и имеет вид:

CALL GRUND (NSL, NI, NJ, NK, NT, LI, M, M1, TKON, A1, B1, C1, N, KL, U, IW, JW, KW, INT, HPRI, H, UN, HNOV, NSI).

Простые переменные NSL, NI, NJ, NK, NT, LI, M, M1, TKON, IW, JW, KW, INT, NSI при обращении к подпрограмме заменяются своими фактическими значениями (см. табл. 6).

Массивы H, UN, HNOV определяются только размерностью в операторе DIMENSION.

В подпрограмме BLOCK DATA (см. приложение 3) оператором DATA вводятся значения по слоям следующих физических и теплофизических характеристик грунта соответствующими идентификаторами:

ВСК (10) - плотность сухого грунта, кг/м³;

ССК (10) - удельная теплоемкость сухого грунта, кДж/(кг·°С), ккал/(кг´ ´°С);

WC (10) - суммарная влажность грунта в долях к весу абсолютно сухого грунта (доли);

UKR (10) - температура начала фазовых переходов.

Значения характеристик берутся из табл. 4.

Текст программы MAIN и подпрограммы BLOCK DATA записывается на стандартных бланках для записи текста на FORTRANe с последующей набивкой на перфокарты.

4.3. Подготовка к вводу группы данных с перфокарт Начальные данные, вводимые с перфокарт описаны в порядке ввода в табл. 7 и 8.

По табл. 7 вводятся данные, необходимые для расчета коэффициентов кривой незамерзшей воды подпрограммой WNS, для каждой выделенной грунтовой разности. Если в программе MAIN обращения к подпрограмме WNS нет, то числовые данные для ввода готовятся только по табл. 8.

Таблица 7

Порядковый номер оператора ввода

Идентификатор

Формат

Единица измерения

Значение

1

T(3,NSL)

1256.2

°C

По три значения температуры, для которых известно количество незамерзшей воды, для каждой грунтовой разности в порядке их нумерации (см. табл. 4)

2

W(3,NSL)

1256.2

доли ед.

Соответствующие указанным выше температурам количество незамерзшей воды

Таблица 8

Порядковый номер оператора ввода

Идентификатор

Формат

Единица измерения

Значение

3

TLM(NSL)

1356.2

Вт/(м·°С),

ккал/(м·ч·°С)

Одномерный массив размерностью NSL. Теплопроводность мерзлого грунта по слоям в порядке их нумерации (см. табл. 4)

TLT(NSL)

Вт/(м·°С),

ккал/(м·ч·°С)

Одномерный массив размерностью NSL. Теплопроводность талого грунта по слоям в порядке их нумерации (см. табл. 4)

WI(N)

м

Одномерный массив размерностью NI. Размер элементов по направлению I

NJ(NJ)

м

Одномерный массив размерностью NJ. Размер элементов по направлению J

HK(NK)

м

Одномерный массив размерностью NK. Размер элементов по направлению K

SHT(NT)

ч

Одномерный массив размерностью NT. Размер шага по времени (см. п.3.6)

4

TR(NT)

6F12.2

ч

Одномерный массив размерностью NT. Граница (время) измерения шага по времени (см. п.3.6)

5

TB(LI,M)

12F6.2

°С

Двухмерный массив размерностью LI´M.

Температура внешней среды на верхней границе области исследования при периодическом ее изменении по времени и зонам краевых условий в порядке их нумерации (см. табл. 2)

RB(LI,M)

м²·°С/Вт,

м²·ч·°С/ккал

Двухмерный массив размерностью LI´M.

Термическое сопротивление поверхности на верхней границе при периодическом изменении по времени и зонам краевых условий в порядке их нумерации (см. табл. 2)

AB(LI,M)

Вт/(м²·°С),

ккал/(м²·ч·°С)

Двухмерный массив размерностью LI´M.

Коэффициент конвективного теплообмена на верхней границе при периодическом изменении по времени и зонам краевых условий в порядке их нумерации (см. табл. 2)

TBIK(M1)

°С

Одномерный массив размерностью M1.

Температура внешней среды по зонам нижней границы в порядке их нумерации

TBJ

Вт/(м²·°С),

ккал/(м²·ч·°С)

Простая переменная. Теплопоток на левой боковой границе области исследования

TBJK

Вт/(м²·°С),

ккал/(м²·ч·°С)

Простая переменная. Теплопоток на правой боковой границе области исследования

TBK

Вт/(м²·°С),

ккал/(м²·ч·°С)

Простая переменная. Теплопоток на ближней боковой границе области исследования

TBKK

Вт/(м²·°С),

ккал/(м²·ч·°С)

Простая переменная. Теплопоток на дальней боковой границе области исследования

6

MR(M)

2713

б/р

Одномерный массив размерностью M.

Род краевых условий по зонам верхней границы в порядке их нумерации (см. табл. 2)

MK(5)

б/р

Одномерный массив размерностью 5.

Род краевых условий соответственно по нижней, ближней, дальней, левой и правой границам

JP(JW)

б/р

Одномерный массив размерностью JW.

Номера элементов по направлению J, задающие «разрезы», для которых выводится на печать распределение температур (см. п. 3.6). Если печать не нужна, то JP(1)=0

KP(KW)

б/р

Одномерный массив размерностью KW.

Номера элементов по направлению K, задающие «разрезы», для которых выводится на печать распределение температур (см. п. 3.6). Если печать не нужна, то KP(1)=0

IP(IW)

б/р

Одномерный массив размерностью IW.

Номера элементов по направлению I, задающие «разрезы», для которых выводится на печать распределение температур (см. п.3.6). Если печать не нужна, то IP(1)=0

Все переменные, соответствующие одному порядковому номеру оператора ввода, вводятся по единому формату.

4.4. Выходные данные При работе программы за каждый шаг расчета по времени обновляется трехмерный массив температур UN(i,j,K), содержащий значения температур в центрах элементов исследуемого массива грунта. На заданный в исходных данных (п. 3.6) момент времени может быть выведено на печать двухмерное температурное поле, представляющее собой разрез массива U, определенный исходными параметрами (см. п. 3.6). Если в некотором элементе происходят фазовые переходы, то есть температура грунта равна температуре начала фазовых переходов (U*), то для этого элемента печатается число, показывающее, какая его часть находится в талом состоянии, по формату: 7ХХ.ХХ, где 7 - условная цифра, а ХХ.ХХ - процентное содержание талого грунта в элементе.

Кроме того, на печать выводятся номер элемента, через который проведен разрез, расстояние от его центра до начала координат и время на моменты печати. Для идентификации точек температурного поля двухмерная матрица окаймляется номерами элементов и значениями расстояния от их центров до начала координат.

Пример печати выходных данных приведен на рис. 5.

Рис. 5. Пример печати выходных данных

4.5. Подготовка данных для расчета по программе PROGNOZ-2S Исходные данные готовятся согласно пп. 4.2-4.3 настоящего раздела, но следующие переменные имеют постоянные значения:

NK=1

iW=1

jW=1

KW=1

HK(1)=1

Ввод массивов jP, KP, iP отсутствует.

4.6. Подготовка данных для расчета по программе PROGNOZ-3N Вся исходная информация вводится с перфокарт. Данные готовятся согласно табл. 9. Все переменные, соответствующие одному порядковому номеру оператора ввода, вводятся по единому формату. Максимальное число расчетных элементов 3000, то есть (Ni´Nj´NK £ 3000).

Таблица 9

Порядковый номер оператора ввода

Формат

Идентификатор

Единица измерения

Значение

1

3I5

Ni

б/р

Число элементов по направлению i

Nj

б/р

Число элементов по направлению j

NK

б/р

Число элементов по направлению K

2

16I5

NSL

б/р

Число выделенных литологических разностей (NSL£10)

NT

б/р

Число изменений шага по времени (NT£10)

Li

б/р

Число интервалов времени в периоде, во время которых задаются средние значения характеристик внешней среды на верхней границе (Li£36)

M

б/р

Число зон с различными краевыми условиями на верхней границе (М£9)

M1

б/р

Число зон с различными краевыми условиями на нижней границе (М£9)

iW

б/р

Число выдаваемых на печать «разрезов» по направлению i (см. п.3.6)

jW

б/р

Число выдаваемых на печать «разрезов» по направлению j (см. п.3.6)

KW

б/р

Число выдаваемых на печать «разрезов» по направлению К (см. п.3.6)

iNT

б/р

Размер массива (см. п.3.6)

NSt

б/р

=1, если количество теплоты измеряется в кДж

¹1, если количество теплоты измеряется в ккал

NGR

б/р

Число элементов, задающих «ступеньки» на верхней границе области исследования. Примечание. Присутствует только в алгоритме PROGNOZ-L

3

12F6.2

T(3,NSL)

°C

По три значения температуры, для которых задается количество незамерзшей воды, для каждой грунтовой разности в порядке их нумерации (см. табл. 4)