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

Белгород 1982

МИНИСТЕРСТВО ЧЕРНОЙ МЕТАЛЛУРГИ СССР Управление горного производства

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

УТВЕРКДлХ): Директор института

И.Ф.Оксанич

27 оентября 1982 г.

МЕТОДИЧЕСКИЕ УКАЗАНИЯ ПО ПРИМЕНЕНИЮ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ ПЛАНОВО-ПРОСТРАНСТВЕННЫХ ЗАДАЧ ЖЛЬТРАДОИ В ТРЕХСЛОЙНОЙ СРВДВ КА ЭЦВМ

Белгород 1982

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

2.    Модуль печати исходной информации состоит из программы FEMID02. . Для его работы необходимы следупдие ресурсы: устройство печати (АЦПУ), устройство прямого доступа (МД), устройство связи о Оператором (ПМ), IOGK байтов оперативной памяти.

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

На АЦПУ распечатывается следующая информация: шифр решаемой задачи, полное время расчета, шаг расчета по времени, таг печати результатов, ширина ленты матрицы системы, количество узлов разбивки, количество элементов разбивки, требуемая точность счета.

.Кроме того, для каждого элемента разбивки печатается: номер элемента, его площадь, мощность, коэффициент фильтрации, коэффициент влагоемкости и начальный градиент фильтрации слабопроницаемого слоя, коэффициенты фильтрации верхнего и нижнего слоев, коэффициенты упругой водоотдачи верхнего и нижнего слоев, отметки кровли и подошвы верхнего и нижнего слоев, номера узлов-, координаты узлов (X иУ ), начальные условия в узлах, тип граничного условия в узлах, граничные условия в узлах.

3.    Модуль решенья задачи плансво-пространственной фильтрации в трехслойиой среде состоит из головной программ Ft И СЛ0*2. и подпрограмм FE.MMS02 ,FEMSA02 ,FEMRQ02» PEMPR02. Весь процесс решения распределен между программами модуля следующим об -разом:FEMCS02 - общее управление решением; FEMMS02 - формиро -вание матрицы системы, учет граничных условий и приведение мат -рииы к верхнему треугольному виду; FEMGA02 - решение системы методом Гаусса; FEMRQ02 - вычисление расходов; FtMDR62 - печать результатов счета.

Для работы модуля необходимы следующие ресурсы: устройство печати (АЦПУ), устройство прямого доступа (МД), устройство связи с Оператором (ПМ), 286К байтов оперативной памяти.

10

Алгоритм работы модуля следующий:

1)    устанавливаются в начальное положение счетчики шага расчета и шага печати результатов;

2)    формируется вектор начальных значений напоров и нулевой вектор промежуточных значений напоров;

3)    формируется матрица системы;

4)    учитываются граничные условия, и матрица приводится к верхнему треугольному виду;

5)    производится решение системы методом Гауоса;

6)    полученные значения напоров сравниваются со значениями вектора промежуточных значений напоров и если отклонение больше заданного, то вектору промежуточных значений напоров присваива -ется величина, равная среднему арифметическому между старым и полученным решением системы,и все повторяется с пункта 3;

7)    производится анализ - нужна ли печать результатов и если да, то вычисляются значения узлов расходов и печатаются результаты счета;

8)    вектору начальных значений напоров присваиваются вычисленные значения напоров, формируется нулевой вектор промежуточны х значений напоров;

9)    производится анализ на конец расчетного интервала времени и, если он не кончился, то повторяется о пункта 3;

10)    конец счета.

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

читаются результаты; заданное значениэ точности счета; номер узла и значения напора и расхода в нем.

Все программные модули должны быть оттранслированы, отредак -тированы и записаны в системную библиотеку загрузочных модул е й SYSILINKLIB или в личную библиотеку пользователя на магнитном диске. Задание на вызов модулей и выполнение составляется на языке управления заданиями ( ЯУЗ) операционной системы ОС ЕС в соответствии о требованиями,предъявлявшми ЯУЗ,

3. П0ДГ0Т0ЖА ИСХОДНОЙ ИНФОРМАЦИИ

При решении задач планово-пространственной фильтрации в трех-олойной среде на ЭЕМ очень важное значение имеет этап подготовки исходной информации. От того, насколько качественно выполнен этот этап, в огромной мере зависит правильность полученных poll

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

Веоь процесс подготовки задачи к решению на ЭШ можно под -разделить на следующие этапы: гидрогеологическая постановка задачи; конечноэлементная схематизация; запиоь походных данных в виде, предназначенном для перфорации; составление пакета перфокарт с данными для счета. Каждый из этих этапов следует вштол -нить со всей возможной тщательностью, так как от этого завиоит успех всей работы в целом. Ниже будет рассмотрен каждый из этапов подготовки данных в отдельности.

3.1. 1Ъдрогеологическая постановка задачи

В большинстве случаев моделированием задач фильтрации на ЭЕМ занимается инженер-гидрогеолог, решапций конкретную практичео -кую задачу. Обычно в его распоряжении имеются данные о гидрогеологическом состоянии, полученные в результате натурных наблю -дений за интересующей облаотыо фильтрации. К ним относятся сведения о действует их и проектируемых водозаборах, хвостохранвли-щах и отстойниках, данные в геологической структуре рассматри -ваемой области, г^идрогеологических параметрах и т.д. Задачей инженера-гидрогеолога на этом этапе является систематизация всей имеющейся информации с целью схематизации области фильтрации.

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

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

коэффициентом фильтрации меньшим, чем у залегающего на нем водоносного горизонта, в IC00 и более раз [4] .

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

Установление всех границ (внешних и внутренних) филътрацион -ной модели и граничных условий на них - заключительный этап схематизации. Расчетная фильтрационная модель должна быть по возможности более простой по сравнению с природной схемой области фильтрации, но она должна функционально соответствовать натурным условиям в сочетании с искусственными сооружениями. Необходимо определить также промежуток времени, для которого будет производиться расчет и выбрать шаг расчета задачи по времени и интервалы, для которых следует выдать результаты на печать.

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

а)    для водоносных горизонтов: отметки кровли и подошвы; коэффициенты фильтрации, упругой и гравитационной водоотдачи; нача -льные и граничные условия; сведения о режиме уровней подземных вод;

б)    для раз делящего слоя: коэффициенты фильтрации и упругой водоотдачи; начальный градиент фильтрации.

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

3*2. Конечноэлементная схематизация области исследований

Следующим очень важным этапом подготовки данных является этап конечноэлементной схематизации области исследований. Шчинаетоя он с разбивки области на конечные элементы (прямые треугольные призмы). Рекомендуется производить ее на гидрогеологической карте исследуемой области.

При разбивке области на конечные элементы следует руководот-воватьол следующими общими правилами: число узловых точек разби-

13

вки по верхнему слою не должно быть более 500, число треугольных элементов не должно быть более 800, построены треугольные элементы не должны иметь тупьа углов, так как это приводит к большим погрешностям счета; на участках области, представляодих паиболь -ший интерес, необходимо разбивку области сделать более мелкой, в то время как на участках, не представляющих особого интереса,элементы разбивки можно выбрать достаточно крупными; желательно,чтобы имеющиеся скважины попадали в узловые точки разбивки. В резу -льтате разбивки на карту исследуемой области,будет нанесена ?тка конечных элементов. Все полученные треугольные элементы разбивки необходимо пронумеровать в произвольном порядке по возрастанию , начиная о I, без пропусков.

Следующим шагом конечноэлементной схематизации является при -своение номеров узловым точкам (вершинам треугольников) разбивки. Для рассматриваемой задачи фильтрации в трехслойной среде при нумерации узловых точек необходимо руководствоваться следующим: номера присваиваются узлам верхнего слоя, причем нечетными цифрами, начиная с I, без пропуоков (номера узлам нижнего слоя присваиваются автоматически, например,узлу 1-го верхнего олоя будет соответствовать узел 2-го нижнего слоя и т.д.); при нумерации узлов необходимо стремиться к тому, чтобы максимальная разность между номерами соседних узлов была как можно меньше (так как от нее зависит объем требуемой оперативной памяти и время счета на ЭШ); предельно допустимое значение этой разности определяется по фор -муле    _

М^р =-^--11

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

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

В результате конечноэлементной схематизации должно быть получено следующее:    таблица номеров треугольных элементов и соответст

вующих им номеров вершин треугольников (только верхнего слоя); таблицы X и V координат узлов разбивки, таблица признаков грани-

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

3.3. Подготовка исходных данных к перфорации

В результате гидрогеологической постановки задачи и конечноэлементной схематизации получаются таблицы исходных данных, которы е надо теперь записать в виде, приго.чном для перфорации. При этом следует руководствоваться следующими рекомендациями: для записи исходных данных использовать специальные бланки или бумагу в кле -точку; данные писать четко, без помарок, по одному символу в кле -точке; числа записывать в десятичной системе счисления; целая часть отделяется от дробной точкой; знак числа ставится перед первой цифрой ("+" можно опускать); данные отделяются друг от друга запятой (например, 1,2.5,-3.65,0.24,...); общее количество сичио -лов, записываемых в одной строке (включая знаки чисел, точки, от -деляющие целые части от дробных, и запятые, разделяющие лонные) iте должно быть больше 80; переносить дачные на новую строку можн о только после запятой, разделяющей данные; каждую новую таблицу ма-чинать записывать с новой строки; в конце каждой таблицы дачных необходимо (обязательно!) ставить символ "точка с запятой".

При задании исходной информации часто встречаются случаи, когда несколько подряд следующих данных имеют одно и то же значение. При этом допускается следующая запись: в-круглых скобках записывается целое число, соответствующее количеству подряд идущих одинаковых данных, а затем значение повторяющегося данного (например, запись (3)-2.5 эквивалента записи - ?.5,-2.5,-2.5).

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

После того, как все ошибки перфорации выявлены и исиравлеиы.оо-

15

УДК (622.001.57:5I)+681.3.06

Настоящие методические указания предназначены для численного моделирования на ЭВМ методом конечных элементов (МКЭ) нестационарных планово-пространственных задач фильтрации в трехслойной среде. Программы составлены на алгоритмическом языке ПЛ/I и предназначены для выполнения под управлением операционной системы X ЕС на ЭШ о оперативной памятью не менее 5I2K. Приведены указания по подготовке исходной информации и работе с прог -рамшми на ЭВМ ЕС-1032. Указания рассчитаны для использования работниками НИИ, производственных и проектных организаций. Програмш составлены В. В.Васильевым. В отладке и апробации программ принимали учаотие Т.Е.Беняш, Ю.И.Волков, В.М.Орел.

Работа составлена В.В.Васильевым и утверждена НТС 22 июля 1982 г. в качестве методических указаний.

_ Всесоюзный нвучно-исслвдомгельский и пропггнп-кг>чс_Тругтор_Ск_ИЙ институт по осушению месторождений полезных ископаемых, сгонкапьнмм горным работам, рудничной геология я морк шейдпрсчому л<._;»_у (ВИОГЕМ), 10*12-

ВВЕДЕНИЕ

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

Наибольшее распространение в настоящее время получили два численных метода решения краевых задач механики сплошных сред: метод конечных разностей (МКР) ч метод конечных элементов (МКЭ), По оценке ряда авторов [5,7,8] МКЭ в олучаях, когда процесс описывается оложными системами дифференциальных уравнений, является гораздо эффективнее МКР. Особенно преимущество МКЭ сказывается, когда надо учитывать сложную конфигурацию области решения, что ха -рактерно для реальных задач осушения месторождений полезных ископаемых.

Раиработанное програлмное обеспечение поззоляет моделировать на ЭДЕМ задачи фильтрации в трехслойном пласте. Практика его ис -пользования в институте ВИОГЕМ показала хорошую согласованность получаемых результатов о данными натурных наблюдений.

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

Рассматривается задача о движении подземных вод в слоистых водоносных горизонтах с чередующимися хорошо и слаболронидаемыми слоями. При этом очятается, что выполняется гилотева Н.К. Гирине — кого и А.Н.Митяева, согласно которой в слабопроницаемых олсях горизонтальные составляющие о корост и фильтрации принимаются нас то — <:ько малыми, что ими можно пренебречь и считать, что фильтрация в нех происходит по вертикали. Кроме того, согласно этой же гипотезе движение в хорошо проницаемых слоях принимается горизонталь — кым. Исходя из этих предположений фильтрационный поток такого вида можно определить как планово-пространственный •

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

Рио.1. Схема планово-пространственной в трехслойной среде

4

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

где h (l,y,t) - напор (уровень) грунтовых вод;№ - интенсивность площадного питания (инфильтрации); t - время;^. - коэффициент гравитационной водоотдачи для безнапорного движения или коэффициент упругой водоотдачи для напорного движения;Т = Km для напор -ной фильтрации иТ =k (h    )    - для безнапорной фильтрации;К -

коэффициент фильтрации; m - мощность пласта; hnoq- отметка подошва водоносного слоя; L₀ - начальный градиент в прослойке; 0,1,2 -индексы, определяющие олабопроницаемый слой, верхний и нижний сдои соответственно,

В общем случае решение системы уравнений (I) должно удовлетворять условиям:

1)    начальное условие hf = Н°(2,Ц);    *    H° fotj),

где Й?и Н$ - известные функции координат;

2)    условие 1-го рода h,(i,y,t)l = Nj(x,y,t>; ho(x,y,t)|^ Н₂ (х,y,t);

3)    условие 2-гс рода

t) =-Т_г|^    Ц;

где    -    расход потока на участке границы Г& для верхнего и

нижнего слоя соответственно; 1_Х1Ц - направлящие кооинусы шещ-ней нормали границы

Разбивая область фильтрации на конечные элементы - прямые треугольные призмы (ом.рио.1) и рассматривая один такой элемент, получим следующую конечноэлементную интерпретацию системы (i)*) • ^{1 2}

2. ОПИСАНИЕ ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ

Программное обеспечение решения задач планово-пространственной фильтрации в трехелойней среде предназначено для эксплуатации под управлением операционной системы ОС ЕС на ЭВМ, имеющих объем оперативной памяти не менее 5I2K. Вое программы написаны на алгоритмическом языке ПЛ/I и объединены по своему функциональному назначению в три неэависишх модуля: модуль ввода, контроля и записи на магнитный диок походных данных; модуль печати исходных данных и модуль счета планово-проотронственной задачи фильтрации в трехслойной среде. Это сделано для того,чтобы повысить надежность работы воего программного комплекса и предотвратить проведение расчетов о заведомо неверными исходными данными. Пользователь полу -чает возможность исправить ошибки в данных до того, как они будут использоваться в модуле снега, а поскольку ооновное машинное время расходуется именно при работе модуля счета, тем самым сокрнщ-ется общее время решения задачи на ЭШ и, следовательно, стоимость решения.

Таким образом, предусматривавтоя следующий порядок решения на ЭВМ:

1)    обработка исходной информации модулем ввода и запись ее на магнитный диск. На этом этапе проиоходит контроль правильности перфорации и частичный контроль доотовернооти вводимых данных. Обнаруженные ошибки распечатываются на алфавитно-цифровом печатающем устройстве (АЩУ), что облегчает их поиок и исправление в массиве исходных данных;

2)    раопечагка введенных данных на АЦПУ модулем печати исходных данных. Эго предусмотрено о целью предоставления пользователю возможности визуального контроля исходных данных и выявления ошибок, которые не были обнаружены на этапе ввода данных;

3)    решение задачи планово-пространотвенной фильтрации в трех -олойной среде, осуществляемое модулем счета.

Ниже приводится краткое описание программных модулей. Схема решения задачи планово-проотранственной фильтрации в трехолойной среде приведена на рис.2.

I. Модуль ввода, контроля и записи на магнитный диск исходных данных оостоит из головной программы FEMID02 и двух вспомогала -льных подпрограмм INFX и INFL. Для работы этого модуля необхо-

димн следующие ресурсы: устройство ввода с перфокарт (УБК); устройство печати (АЦПУ), устройство прямого доступа (МД); устройство связи с Оператором (ПМ), I30K байтов оперативной памяти.

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

Работа модуля начинается о ввода исходной информации с перфокарт. Параллельно о вводом, как было указано выше, осуществляется частичный контроль вводимой информации. На этом этапе контролируется правильность перфорации цифровой информации (т.е.чтобы не было символов, отличных от    ",^п,0+9),а    также    ее    пол

нота (например, чтобы количество введенных координат узлов в точности равнялось чио^ узлов точек разбивки и т.д.).

В случае обнаружения ошибки на АЦПУ печатается сообщение, неверные данные не принимаются в обработку, а работа модуля на этом заканчивается. При отсутствии ошибок ввода, для каждого треугольного элемента вычисляется его площадь и проверяется ши -рина ленты матрицы элемента (разница соседних номеров узлов элемента). Тем самым контролируется правильность нумерации узлов разбивки. Если встречается элемент, площадь которого получается равной нулю (это может быть вызвано, например, неверным заданием координат узлов треугольника) или максимальная разность между номерами вершин (включая и вершины треугольника нижнего слоя) окажется больше заданной ширины ленты матрицы системы, то на АЦПУ печатается сообщение об обнаруженной ошибке и работа модуля прекращается. Формирование набора данных на магнитном диске производится только при отсутствии ошибок в исходных данных.

Такал организация процесоа ввода позволяет на начальном этгае решения задачи на ЭВМ обнаружить и устранить ошибки в исходных

9

1

*) Разработка алгоритма и програш по реализации метода конечных элементов для решения планово-пространственных задач движения подземных вод и теплопроводности. Отчет/ВЙ01ЕМ.руководитель работы В.В.Васильев.Шифр работы 1.5-2 (79)-С-&-7$И>, Белгород, 1962, 130 о.

2