На правах рукописи
БОБАРЫКИН Николай Дмитриевич
УДК 556.324.001.57(06)
ОПТИМАЛЬНОЕ УПРАВЛЕНИЕ РЕЖИМОМ ГРУНТОВЫХ ВОД НА ОСНОВЕ ИНВАРИАНТНОЙ НЕСТАЦИОНАРНОЙ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ПОЛЬДЕРНЫХ СИСТЕМ
Специальность:
05.13.18 - Математическое моделирование, численные методы
и комплексы программ
Автореферат диссертации на соискание ученой степени
доктора технических наук
Калининград 2007
Работа выполнена в ГОУВПО Калининградском государственном техническом университете (КГТУ)
Научный консультант:
доктор физико-математических наук, профессор, атышев К.С.
(Российский государственный университет им. И. Канта (РГУ им. И. Канта))
Официальные оппоненты:
член-корреспондент РАН,
доктор физико-математических наук, профессор Холодов А.С.
(Институт автоматизации проектирования РАН (ИАП РАН), г. Москва);
доктор физико-математических наук, профессор Дикусар В.В.
(Вычислительный центр РАН, г. Москва);
доктор физико-математических наук, профессор Альес М.Ю.
(Институт прикладной механики УрО РАН, г. Ижевск).
Ведущая организация: Институт математического
моделирования РАН (г. Москва)
Защита состоится 10 декабря 2007г. в 14 часов
на заседании диссертационного совета Д 212.065.04
в ИжГТУ по адресу: 426069, г. Ижевск, ул. Студенческая, 7.
Отзыв на автореферат, заверенный гербовой печатью,
просим выслать по указанному адресу.
С диссертацией можно ознакомиться в библиотеке института,
а с авторефератом диссертации - на официальном сайте ВАК РФ vak.ed.gov.ru.
Автореферат разослан _______ ноября 2007 г.
Ученый секретарь диссертационного совета,
доктор технических наук, профессор Б.Я. Бендерский
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Актуальность темы. Агропромышленная отрасль страны, использующая комплексную мелиорацию заболоченных и переувлажненных земель, оказывает исключительное влияние на развитие сельскохозяйственного производства в целом, без которого не могут быть решены вопросы продовольственной безопасности страны и конкурентоспособности российской сельскохозяйственной продукции на мировом рынке. Сохранение, воспроизводство и рациональное использование плодородия земель сельскохозяйственного назначения является актуальной проблемой, требующей принятия эффективных решений и действий. Улучшение плодородия сельхозугодий расценивается как важнейшая забота о национальном достоянии страны. Выполненный в 2004 г. в соответствии с Концепцией развития мелиорации в России анализ состояния агропромышленного комплекса показал, что в целом для устойчивого развития сельского хозяйства страны необходимо увеличить площадь осушаемых земель до 7 - 8 млн. га. Решение этой задачи требует создания и внедрения наукоемких технологий, совершенных мелиоративных систем, методов, форм и средств управления режимами комплексных осушительных мелиораций.
Основными средствами транспортировки воды от осушаемых массивов (ОМ) в рассматриваемых мелиоративных (польдерных) системах (ПС) служат сети проводящих открытых каналов (СПОК). Главными элементами СПОК являются: пассивные элементы - линейные участки проводящих каналов, узлы - места их соединения и активные элементы - откачивающие насосные станции, сосредоточенные в активных узлах. Для эффективного и равномерного осушения мелиорированных земель используются дренажные трубы (ДТ) различных диаметров, укладываемые на определенную глубину и расстояние между ними, перпендикулярно руслам проводящих каналам.
Наиболее широкое применение на практике проектирования и эксплуатации мелиоративных систем имеет метод водного баланса, который позволяет проводить оценки значений уровня грунтовых вод (УГВ) осушаемого массива с учетом речного стока или проводящими каналами польдерных систем.
Расчеты водного баланса по совокупности приходных и раснходных потоков влаги в их динамике и взаимосвязи весьма сложны и трудоемки. Составляющие компоненты водного баланса выступают как случайные переменные величины, и могут быть как ненпрерывными функциями, так и дискретными, и их взаимозависимость не выражается однонзначно в явном виде и является малоизученной, например, взаимосвязь между приходящими и уходящими потоками влаги на границах зоны аэрации и поверхностью грунтовых вод, что делает этот метод малопригодным для управления водным режимом осушаемых массивов в реальном масштабе времени.
Другой вид математических моделей мелиоративных систем, базирующихся на гидродинамических дифференциальных уравнениях, имеет два существенных недостатка и также не способствует построению методологии управления режимом увлажнения корнеобитаемого слоя почвы (РУКС). Во-первых, невозможность увязки в единую математическую модель фрагментарного моделирования процессов переноса воды в проводящей сети и в осушаемом массиве без учета их взаимодействия. Во-вторых, отсутствие инструментария системного подхода при моделировании сложных инженерно-технических систем, включающего семантическую модель технологического процесса, использования теории направленного графа объекта, синтез разрешающей системы дифференциальных уравнений и т.д.
Таким образом, постановка задачи математического моделирования и управления мелиоративными системами требует системного подхода и дальнейшего усовершенствования. В частности, необходимо учесть взаимодействие осушаемого массива со стоком каналов и рек, и обеспечить создание энергосберегающих технологий, обеспечивающих оптимальный режим увлажнения корнеобитаемого слоя почвы, путем управления УГВ ОМ (оптимальному положению УГВ Нор соответствует максимально возможный поток влаги Umax от поверхности грунтовых вод).
При отсутствии теоретических разработок, методик и работоспособных алгоритмов оптимального управления РУКС, основанных на модельных расчетах, а также стратегии управления УГВ, требуется создание новых концепций и постановок задачи трехмерного моделирования ПС и управления РУКС почвы осушаемых земель.
Для реализации решения данной проблемы необходимо провести анализ и обобщить опыт моделирования отдельных задач в процессе эксплуатации и конструирования мелиоративных систем, сформулировать интегральную постановку задачи математического моделирования ПС с единых позиций, где каждая прикладная задача рассматривается как частный случай общего подхода, разработать эффективные методы ее численного решения в малоисследованной области многокритериальных задач, включающих нелинейные нестационарные одномерные и двухмерные дифференциальные уравнения гидродинамики, нелинейные алгебраические соотношения и системы неравенств.
Решению данных актуальных задач и посвящена настоящая диссертационная работа.
Объектом исследования являются польдерные системы, состоящие из сетей проводящих открытых каналов с откачивающими насосными станциями, осуществляющими отвод воды с осушаемых массивов и, тем самым, управляя уровнем грунтовых вод УГВ Нор, обеспечивая эффективный режим увлажнения корнеобитаемого слоя почвы мелиорированных земель.
Предметом исследования являются методы системного анализа конструкций и технологий ПС; методы и алгоритмы численного расчета уровневого режима в СПОК, автоматически учитывающие баланс потоков воды в каналах, основанные на системе дифференциальных уравнений Сен-Венана с естественными начальными и граничными условиями; методы и алгоритмы численного решения двумерного уравнения Буссинеска, описывающего пространственно-временное распределение УГВ в осушаемых массивах; методы и алгоритмы численного расчета потоков воды в дренажных системах; методы и алгоритмы численного расчета влагообмена в зоне аэрации описываемого уравнением капиллярного потенциала влаги; стратегия и алгоритмы оптимального управления УГВ, включая целевые функции и функциональные ограничения параметров технологических процессов; выбор и реализация эффективных методов многокритериальной скалярной и векторной оптимизации; информационное обеспечение; реализованная интегральная структура математического моделирования ПС и управления РУКС.
Цель работы состоит в создании научно-обоснованной интегральной трехмерной нестационарной математической модели ПС (ИТНММ ПС) и управления РУКС почвы мелиорированных земель на основе стратегии управления УГВ, обеспечивающей энергосбережение, высокую эффективность вегетации сельскохозяйственных культур и минимизацию материальных и трудовых затрат при проектировании, строительстве и эксплуатации польдерных систем. Для реализации поставленной цели необходимо решить следующие задачи:
- анализ существующих технологий и теоретических работ по проектированию, эксплуатации и математическому моделированию польдерных систем с целью выработка научно-обоснованных решений для построения трехмерной математической нестационарной инвариантной модели ПС и управления УГВ, РУКС; разработка ее методологии, состава и структуры;
- разработка алгоритмов численного расчета потоков влаги от грунтовых вод с учетом испарения с поверхности почвы, транспирации растений и выпадающих осадков, основанных на численном решении дифференциальных уравнений капиллярного потенциала почвы и теплового баланса, обеспечивающих оптимальный режим увлажнения в корнеобитаемом слое мелиорированных почв;
- разработка вектора целевых функций, имеющего следующие компоненты: интенсивность снижения уровня грунтовых вод (ИС УГВ), время снижения уровня грунтовых вод (ВС УГВ); равномерность распределения УГВ (РР УГВ); исследование и выбор функциональных ограничений параметров технологического процесса; разработка стратегии оптимального управления УГВ: ИС УГВ max, BC УГВ min и РР УГВ min при выполнении функциональных ограничений параметров технологического процесса; выбор эффективных методов многокритериальной скалярной и векторной оптимизации;
- декомпозиция исходной постановки задачи моделирования СПОК на подзадачи с последующим синтезом системы разрешающих дифференциальных уравнений, осуществляемым на ребрах направленного графа, с заданием граничных условий в вершинах графа;
- построение алгоритма численного решения системы дифференциальных уравнений Сен-Венана, описывающих уровневый режим в СПОК, для которого автоматически выполняется закон сохранения потока воды в точках ветвления проводящих каналов;
- моделирование динамики УГВ и РУКС почвы ПС, основанного на двухмерном уравнении Буссинеска, дифференциальном уравнении переноса воды в ДТ и трехмерном уравнении для капиллярного потенциала влаги, при задании естественных граничных условий: а) взаимодействия СПОК с ОМ; б) влагообмена между УГВ и корнеобитаемым слоем почвы;
- разработка экономичных алгоритмов, основанных на консервативных и эффективных разностных схемах, численного решения нелинейных дифференциальных уравнений для дренажных систем, капиллярного потенциала влаги и уровня грунтовых вод в ОМ;
- разработка алгоритма интерпретатора полученных результатов расчета параметров ПС, преобразующего их из ЛСК в ОСК ПС;
- разработка алгоритмов, основанных на консервативных и эффективных разностных схемах повышенного порядка точности, численного решения дифференциальных уравнений в обыкновенных и частных производных параболического и гиперболического типов, необходимых для тестирования и выверки применимости более экономных разностных схем первого порядка.
Методы исследования. В работе используются методы теории оптимального управления, системного анализа конструкций и технологий польдерных систем, математической статистики, теории и практики математического моделирование сложных инженерно-технологических систем, численные методы решения нелинейных дифференциальных уравнений гидродинамики и теории фильтрации жидкостей через пористые структуры, методы решения оптимизационных многокритериальных задач, теория графов.
При моделировании и управлении РУКС почвы польдерная система представляется в виде множеств пассивных, активных элементов и точек соединений проводящих каналов, на которые накладывается структура направленного графа, соответствующая расчетной схеме ПС. Каждый элемент множеств получает свое математическое описание в рамках теории гидродинамики течений (фильтрации) жидкости.
Синтез разрешающей системы моделирующих уравнений, описывающих динамику воды в СПОК с соответствующими начальными условиями, осуществляется на ребрах, а формирование граничных условий, выражающих законы сохранения расхода воды, в вершинах графа ПС, которые являются условиями сшивания системы моделирующих уравнений вдоль ветвей обхода графа ПС. Производительность активных элементов - насосных станций, при их наличии в ветвях графа ПС, используется в качестве соответствующих граничных условий.
Моделирование режима УГВ и РУКС почвы ОМ, прилегающих к СПОК, основывается на численном решении двухмерного уравнения Буссинеска, трехмерного уравнения для капиллярного потенциала и задании естественных граничных условий взаимодействия СПОК и ОМ. При наличии дренажных систем к указанной системе дифференциальных уравнений подключается и дифференциальное уравнение переноса воды в ДТ. Таким образом, задача моделирования РУКС почвы сводится к решению трехмерной нестационарной задачи, в которой пространственно - временное распределение УГВ вычисляются в горизонтальных, а потоков влаги U- в вертикальных плоскостях.
Для численного интегрирования системы дифференциальных уравнений Сен-Венана вдоль ветвей обхода графа ПС с учетом сохранения расхода воды в точках соединения каналов, указанная система уравнений приводилась к каноническому виду. А алгоритм численного решения строился при введении новых прогоночных соотношений, учитывающих число обхода каждого проводящего канала, что автоматически обеспечивает выполнение закона сохранения потоков в точках ветвления каналов. Численные решения остальных уравнений в частных производных строились на основе консервативных разностных схем с первым порядком точности по времени и вторым - по координатам, а также проводились итерации по нелинейности и связанности дифференциальных уравнений.
Стратегия управления РУКС почвы мелиорированных земель определяется управлением потоком влаги U(Нор, h) Umax, при t min и учете технологических ограничений, в процессе которого уровень грунтовых вод ННор (критерий управления уровнем грунтовых вод УГВНор) и высота корнеобитаемого слоя h hор. Вектор целевых функций для управления УГВ формируются из трех компонент - интенсивности снижения уровня грунтовых вод, времени снижения УГВ до определенного значения и равномерности распределения УГВ.
Стратегия оптимального управления уровнем грунтовых вод, обеспечивающая энергосбережение и эффективность вегетации растений мелиорированных земель сводится к следующему: ИС УГВ max, BC УГВ min и РР УГВ min, при t min и учитываются технологические ограничения.
При эксплуатации существующих ПС управление УГВ сводится к подбору оптимальной производительности насосных станций и их режимов работы. А при проектировании новых польдерных систем, кроме управления насосными станциями, осуществляется еще и подбор оптимальных значений геометрических параметров проводящих каналов, дрен, глубины их закладки и т.д. Методы и алгоритмы скалярной и векторной многокритериальной оптимизации выбирались из библиотеки математического обеспечения пакета прикладных программ MATHCAD.
Математическое моделирование проводилось на персональном компьютере Pentium IV - 1700 c использованием среды MATHCAD, при этом программное обеспечение расчетов ПС разработано на структурно-блочном алгоритмическом языке данной системы.
Достоверность и обоснованность полученных в работе результатов и выводов подтверждается сопоставительным анализом между уже существующими и разработанными математическими моделями и методами, а также итогами практического использования интегральной трехмерной нестационарной математической модели ПС и алгоритмов управления УГВ и РУКС в ОМ при эксплуатации существующих и проектировании новых польдерных систем.
Созданная научно-обоснованная интегральная трехмерная нестационарная математическая модель ПС и управления РУКС, базируется на положениях теории гидродинамических и фильтрационных течений жидкости, системного анализа конструкций и технологий польдерных систем, математической статистики, теории и практики управления и математического моделирования, сложных инженерно-технологических систем, численных методов решения нелинейных дифференциальных уравнений, теории направленных графов, методов теории оптимального управления.
Достоверность практических и теоретических разработок ИТНММ ПС и управления РУКС почвы ПС подтверждается результатами многочисленных численных экспериментов по расчету параметров СПОК и ПС в целом, а также оптимальному управлению УГВ и РУКС в ОМ и сравнения с их натурными значениями и характерными временами управления. Все разработанные алгоритмы численного решения дифференциальных уравнений тестировались на многочисленных соответствующих тестовых задачах.
На защиту выносятся результаты проведенных численных исследований по проектированию совершенных (оптимальных) конструкций ПС и эксплуатации существующих польдерных систем в штатных и нештатных ситуациях в режиме оптимального управления УГВ, выполненных на основе созданной научно-обоснованной интегральной трехмерной нестационарной математической модели ПС и управление РУКС почвы мелиорированных земель, построенной на единой методологии моделирования ПС, учитывающей водообмен на границе между сетью проводящих каналов и осушаемых массивов, а также принципы и методы разработки этой модели и стратегии управления УГВ и РУКС, в том числе:
- обоснованы и сформулированы функциональные ограничения параметров технологического процесса увлажнения корнеобитаемого слоя на основе учета факторов, дестабилизирующих вегетацию растений, а также компоненты вектора целевых функций, характеризующего качества ПС, стратегия оптимального управления уровнем грунтовых вод и РУКС, обеспечивающая энергосбережение и оптимальную эффективность вегетации растений, методология, состав и структура интегральной трехмерной нестационарной математической модели ПС и условия практического использования математического аппарата теории гидродинамических (фильтрационных) течений жидкости и капиллярного потенциала влагопереноса, позволяющего однозначно описывать трехмерное пространственно-временное распределение воды и влаги в сети проводящих открытых каналов и осушаемых массивах;
- разработка программного и информационного обеспечения научно-обоснованной ИТНММ ПС и управления РУКС почвы мелиорированных земель, позволяющей решать конкретные задачи максимально возможного увеличения плодородия переувлажненных земель для каждого вида сельскохозяйственных культур. При этом проводится декомпозиция исходной постановки задачи моделирования ПС на подзадачи с последующим синтезом системы разрешающих дифференциальных уравнений, осуществляемым на ребрах направленного графа, с заданием граничных условий в вершинах графа;
- алгоритмы и программное обеспечение численного решения разрешающей системы дифференциальных уравнений Сен-Венана, описывающих динамику воды в СПОК, для которого автоматически выполняются законы сохранения потоков воды в точках ветвления проводящих каналов; реализация алгоритмов расчета СПОК и реальных польдерных систем, в том числе с учетом рельефа местности;
- результаты моделирования динамики УГВ и управления РУКС почвы ПС, основанные на двухмерном уравнении Буссинеска, дифференциальном уравнении переноса воды в дренажных трубах и трехмерном уравнении переноса капиллярного потенциала влаги, при задании естественных граничных условий: а) взаимодействия СПОК с ОМ; б) влагообмена между УГВ и корнеобитаемым слоем почвы;
- результаты разработки алгоритмов и программного обеспечения расчетов потоков влаги от грунтовых вод с учетом поверхностного испарения, транспирации растений и выпадающих атмосферных осадков, основанных на численном решении дифференциальных уравнений капиллярного потенциала и теплового баланса почвогрунтовой системы; решения вопросов математического моделирования сезонных потоков влаги от грунтовых вод с учетом поверхностного испарения и транспирации влаги;
- алгоритм расчета оптимальной вегетации сельскохозяйственных культур; формирование целевой функции для потока влаги U(Н, h) и функциональных ограничения параметров, дестабилизирующих вегетацию растений; построение стратегии управления РУКС;
- реализация разработанных алгоритмических средств управления РУКС и УГВ с учетом технологических ограничений в виде программного комплекса; реализация алгоритмов оптимального управления режимом грунтовых вод с учетом выпадающих атмосферных осадков;
- результаты разработки алгоритма расчета реальных польдерных систем, основанного на преобразовании параметров ПС, вычисленных в локальных системах координат, связанных с проводящими каналами, в общую систему координат ПС;
- результаты разработки алгоритмов повышенного порядка точности, численного решения дифференциальных уравнений необходимых для тестирования и выверки применимости экономичных разностных схем первого порядка точности, при проведении оптимизационных расчетов;
- разработка научно-технических рекомендаций по конструированию и эксплуатации польдерных систем в режиме оптимального управления уровнем грунтовых вод, обеспечивающих энергосбережение и высокую эффективность вегетации растений.
Научная новизна полученных результатов определяется проведенными комплексными исследованиями, в результате которых вместо применения отдельных математических моделей, разработанных для индивидуальных вариантов расчетов параметров ПС, применена созданная научно-обоснованная ИТНММ ПС и управления РУКС почвы мелиорированных земель, основанная на единой методологии моделирования ПС при задании естественных граничных условий: а) взаимодействия СПОК с ОМ; б) влагообмена между УГВ и корнеобитаемым слоем почвы. ИТНММ ПС обеспечивает сокращения время на решение конкретных задач моделирования и управления ПС при одновременном повышении качества полученных результатов (энергосбережение ресурсов, высокую эффективность вегетации сельскохозяйственных культур и минимизацию материальных и трудовых затрат при проектировании, строительстве и эксплуатации польдерных систем), в ходе которых:
- определен состав и структура задач, подлежащих решению при реализации ИТНММ ПС и управления РУКС почвы мелиорированных земель, по схеме - УГВ Нор - U Umax, при t min и технологических ограничений; разработана иерархическая структура ИТНММ ПС и управления УГВ и РУКС, основанная на применении математического аппарата теории гидродинамических (фильтрационных) течений жидкости и капиллярного потенциала влагопереноса, позволяющего однозначно описывать трехмерное пространственно-временное распределение воды и влаги;
- поставлена и решена задача создания научно-обоснованной интегральной трехмерной нестационарной математической модели ПС и управления РУКС почвы ОМ, путем управления уровнем грунтовых вод мелиорированных земель;
-построен новый эффективный алгоритм численного решения системы уравнений Сен-Венана вдоль ветвей направленного графа ПС, основанный на приведении к каноническому виду исходной системы дифференциальных уравнений, имеющих недивергентные члены, с учетом числа обхода каждого элемента СПОК и введении новых прогоночных соотношений, что в целом позволяет корректно задавать граничные условия, строить разностные схемы с учетом наклонов характеристик, а также автоматически выполнять законы сохранения потока воды в точках ветвления проводящих каналов;
- решены вопросы моделирования параметров СПОК ПС и реальных польдерных систем, в том числе с учетом рельефа местности;
- разработана и реализована математическая модель УГВ и РУКС почвы ОМ, основанная на двухмерном уравнении Буссинеска, включая дифференциальные уравнения переноса воды в ДТ и капиллярного потенциала влаги, при задании естественных граничных условий взаимодействия СПОК и ОМ, а также между УГВ и корнеобитаемым слоем; разработан алгоритм расчета потоков влаги от грунтовых вод в зависимости от испарения с поверхности почвы, транспирации для растений и выпадающих осадков;
- разработаны и реализованы алгоритмы численного моделирования сезонных потоков влаги от грунтовых вод в зависимости от поверхностного испарения и транспирации;
- разработаны алгоритмы численного решения систем нелинейных одномерных и двухмерных дифференциальных уравнений в частных производных гиперболического и параболического типов; выполнен анализ результатов численных расчетов параметров польдерных систем;
- разработана и реализована стратегия оптимального управления УГВ, обеспечивающая энергосбережение и оптимальную эффективность вегетации растений при РР УГВ; сформирован вектор целевых функций, характеризующий качества ПС, с компонентами: ИС УГВ, ВС УГВ и РР УГВ; разработаны функциональные ограничения параметров технологических процессов на основе учета факторов, нарушающих эффективный режим увлажнения в корнеобитаемом слое мелиорированных почв;
- выбраны и реализованы эффективные методы многокритериальной скалярной и векторной оптимизации;
- разработаны и реализованы алгоритмы оптимального управления режимом грунтовых вод с учетом выпадающих атмосферных осадков заданной интенсивности и длительности;
- разработаны научно-технические рекомендации по проектированию совершенных польдерных систем и их эксплуатации в режиме оптимального управления УГВ, обеспечивающие ресурсосбережение и высокую эффективность вегетации растений.
Практическая полезность исследования состоит в том, что в результате комплексных исследований на практике реализован математический аппарат теории гидродинамических и фильтрационных течений жидкости, а также капиллярного потенциала влагопереноса в зоне аэрации почвы, позволяющий однозначно описывать пространственно-временное распределение уровней и скоростей воды в СПОК, а также потоков влаги в ОМ, а с использованием разработанной стратегии управления УГВ устанавливать оптимальный РУКС почвы мелиорированных земель.
Применение ИТНММ ПС и управления РУКС почвы, путем управления УГВ мелиорированных земель, основанной на единой методологии моделирования ПС, учитывающей влагообмен между СПОК и ОМ, позволяет резко сократить время на решение конкретных задач моделирования ПС при одновременном повышении качества полученных результатов. Построенная стратегия оптимального управления РУКС обеспечивает энергосбережение и высокую эффективность вегетации сельскохозяйственных культур мелиорированных земель. В целом, применение разработанной ИТНММ ПС и стратегии управления УГВ обеспечивает минимизацию материальных и трудовых затрат при проектировании, строительстве и эксплуатации польдерных систем.
Полученные в работе методики и алгоритмы численных расчетов параметров СПОК ПС и реальных польдерных систем, в том числе с учетом рельефа местности, могут иметь самостоятельную практическую ценность, также как и разработанные алгоритмы численного решения системы уравнений Сен-Венана вдоль ветвей направленного графа ПС, автоматически учитывающие законы сохранения потока воды в проводящих каналах, основанные на приведении к каноническому виду исходной системы дифференциальных уравнений, имеющих недивергентные члены.
На основе модельных численных расчетов, проведенных по управлению РУКС почвы ОМ, выработаны научно-технические рекомендации по проектированию совершенных польдерных систем и их эксплуатации.
Реализация результатов работы. При непосредственном участии автора была разработана и реализована ИТНММ ПС и управление РУКС почвы мелиорированных земель по схеме - УГВ Нор - U Umax, при t min и технологических ограничений, построенная на положениях теории гидродинамических и фильтрационных течений жидкости, системного анализа конструкций и технологий польдерных систем, математической статистики, теории и практики управления и математического моделирования, сложных инженерно-технологических систем и реализованная на основе численных методов решения нелинейных дифференциальных уравнений, теории направленных графов, методов теории оптимального управления. Созданные комплексы алгоритмов, программ, технических рекомендаций по конструированию и эксплуатации польдерных систем используются при выполнении проектных работ по гранту РФФИ РАН № 06-01-00396 Оптимизация управления мелиоративными системами на основе математических моделей (руководитель проекта - Бобарыкин Н.Д., сроки разработки 2006-2008г.г.), в АО Калининградоблмелиорация, в ходе выполнения госбюджетной НИР в КГТУ по теме: Инвариантные методы расчета технологических процессов, а также в учебном процессе по дисциплинам: Математическое моделирование на ПЭВМ, Компьютерные технологии и Информатика.
Вся работа в целом, а также ее отдельные части могут быть использованы предприятиями агропромышленного комплекса и организациями, эксплуатирующими и разрабатывающими сети нефтегазопроводов с перекачивающими агрегатами и насосными станциями, сетями водоканалов и т.д.
Апробация работы. Основные результаты диссертационной работы на различных этапах ее выполнения докладывались и обсуждались на международных и межведомственных конференциях и совещаниях, в частности: Физике ионосферы (Москва, 1975 г.), Межведомственного семинара по моделированию ионосферы (Тбилиси, 1980), Правлении Государственного газового концерна Газпром (Москва, 1990г.), Межотраслевой конференции Творческий вклад молодежи в дело перестройки (Москва, 1990 г.), объединенных семинарах Московского НПО Нефтегазавтоматика и МИНГ им. И.М. Губкина (1990г.), МХТИ им. Д.И. Менделеева (Москва, 1991г.), Материалы международного симпозиума Новые технологии в нейрохирургии (Санкт-Петербург, ВМА, 2004г.), Международной научной конференции, приуроченной к 200-летию со дня рождения великого немецкого математика Карла Густава Якоби и 750-летию со дня основания г. Калининграда (Кенигсберга) (Калининград, КГУ, 2005г.), Международной научной конференции Инновации в науке и образовании - 2005г., Калининград, КГТУ, Институте математического моделировании РАН (Москва, 2005г.), Расширенном научно-техническом совете ОАО Белгорхимпром с участием Института математики НАН Беларуси, Белгосуниверситета, БеАН, Института проблем использования природных ресурсов и экологии НАН Беларуси, Белорусского НИИ геологоразведочного института (Минск, 2006г.), Институте математического моделировании РАН (Москва, 2007г.), Международной научной конференции: Российская наука и инженерная деятельность в социокультурном пространстве эксклавного региона: история, актуальные проблемы, перспективы развития, г. Калининград, КГТУ, 2007, и др.
Публикации. Основные научные результаты по теме диссертационной работы опубликованы в 46 печатной работе, в том числе: 1 монография (научное издание, 168с., 2004 г.). Автор имеет 18 научных трудов в изданиях, выпускаемых в РФ и рекомендуемых ВАКом для публикации основных результатов диссертации на соискание ученой степени доктора наук.
Структура диссертационной работы определяется общим замыслом и логикой проведения исследований.
Диссертация содержит введение, 7 глав, заключение и приложения, изложенные на 310 страницах компьютерного текста. В работу включены 102 рис., 6 табл., список литературы из 232 наименований и приложение, в котором представлены акты об использовании результатов работы.
СОДЕРЖАНИЕ РАБОТЫ
Введение содержит обоснование актуальности темы, формулировку цели и задач работы, основные положения, выносимые на защиту, и определяет содержание и методы выполнения работы. Изложены полученные автором основные результаты проведенных исследований, показана их научная новизна и практическая значимость.
В первой главе дан сравнительный анализ, различных математических методов поэлементного расчета параметров польдерных систем, который позволил выявить и систематизировать основные недостатки существующих моделей ПС и сформулировать задачи исследований. Обоснована актуальность разработки современных методов повышения плодородия переувлажненных и заболоченных земель на основе системы оптимального управления УГВ и, соответственно, РУКС почвы.
Отмечается, что метод водного баланса ограничивается необходимостью проведения комплексных экспериментальных исследований элементов водного баланса (уровни грунтовых вод, сток, влагозапасы и др.) переувлажненных земель, а также их водно-физических свойств. Также недостатком метода водного баланса является сложность оценки изменения экстремальных величин уровней и расходов влаги, что делает этот метод малопригодным для управления водным режимом польдерных систем.
Математическое моделирование мелиоративных систем, базирующееся на гидродинамических дифференциальных уравнениях, также имеют два основных недостатка. Во-первых, обособленность процессов переноса воды в проводящей сети открытых каналов и осушаемых земельных массивах друг от друга, и, соответственно, невозможности увязки такого фрагментарного моделирования отдельных процессов в единую математическую модель. Во-вторых, отсутствие системного подхода при моделировании польдерных систем, т.е. инвариантных алгоритмов, позволяющих проводить расчеты независимо от конфигурации ПС и числа проводящих каналов. Особенно, это сказывается для средних и больших мелиоративных систем, где число и изменчивость параметров, определяющих эти системы, резко возрастает.
Таким образом, постановка задачи математического моделирования и управления мелиоративными системами должна основываться на единой методологии и системном подходе и требует дальнейшего усовершенствования, в частности учета взаимодействия осушаемого массива со стоком каналов и рек, и должны быть направлены на создание высокоэффективных энергосберегающих технологий, обеспечивающих оптимальный РУКС почвы.
Во второй главе приведены численные методы решения обыкновенных, параболических и гиперболических дифференциальных уравнений, что обусловлено необходимостью выбора, а в ряде случаев - разработки методов численного решения системы дифференциальных уравнений, на которой строится интегральная трехмерная инвариантная математическая модель польдерных систем. Также, изучены вопросы корректной постановки граничных условий и повышения порядка точности аппроксимаций дифференциальных уравнений гиперболического и параболического типа, путем введения полуцелых временных слоев , с целью оценки применимости экономичных разностных схем, например, схем бегущего счета, имеющих первый порядок точности аппроксимаций производных.
В табл. 1 приведены основные конструкции разностных схем, как первого, так и второго порядка точности, а так же канонический вид системы гиперболических уравнений, полученной по методике предложенной автором диссертации (см. подраздел 2.1.3). При этом численное решение канонической системы уравнений определяется во всех точках характеристического треугольника, а краевые условия для первого уравнения системы задаются на левой, а для второго уравнения на правой границе.
Показано, что погрешность численного решения гиперболических скалярных уравнений при увеличении временного шага возрастает по разностным схемам, имеющих первый и второй порядка точности, и при = 1200 с не превышает 5%, а систем гиперболических уравнений не превышает 6.5%, что может оправдывать применения разностных схем первого порядка точности при численном решении дифференциальных уравнений гиперболического типа в частных производных.
Применимость таких разностных схем, особенно актуально с точки зрения экономии машинного времени на решение задач управления УГВ, реализуемых на основе оптимизационных решений при большом массиве значений варьируемых параметров.
В третьей главе представлен состав и структура интегральной трехмерной нестационарной математической модели ПС и управления УГВ и РУКС почвы мелиорированных земель базирующей на наличие для польдерных систем двух взаимодействующих контуров - сети проводящих открытых каналов и прелагающей к ней системы осушаемых земельных массивов.
Таблица 1.
Разностные схемы решения параболических и гиперболических уравнений
№ | Тип уравнений | Нач. условия | Гран. условия | Примечание | ||||||
1 | Параболическое уравнение второго порядка | |||||||||
Н(x; 0) = 3 м | N(0; t) = (t); | х0=0; хn=L1; Н - уровень воды; Ф- коэф. фильтрации, м2/с | ||||||||
Консервативная разностная схема порядка точности 0(+h2) | ||||||||||
i = 1, 2, . . . , n - 1; j = 0, 1, . . . , k-1. | ||||||||||
2 | Скалярное гиперболическое уравнение, V0>0 | |||||||||
N(x, 0)= 1028м-3 | N(0, t) = 1.21028 м-3 | N, V0 Цконцентрация и скорость потока частиц | ||||||||
Разностная схема бегущего счета первого порядка точности 0(+h) | ||||||||||
или , где = ; i = 1, 2, Е, n; j = 0, 1, Еk-1. | ||||||||||
Разностная схема второго порядка точности 0(2+h) | ||||||||||
или , i = 1, 2, Е, n; или , j = 0, 1, Еk-1 | ||||||||||
3 | Система гиперболических уравнений | |||||||||
; | N(x, 0)= 1028м-3 V(x, 0) = 0 | N(0, t) = 1.21028 м-3 V(xn, t) = 0.1 м/с | c = 4 м/c; R = 10-3 c-1 | |||||||
Канонический вид системы гиперболических уравнений | ||||||||||
; | ||||||||||
Разностная схема бегущего счета первого порядка точности 0(+h) | ||||||||||
; |
При построении ИТНММ ПС используются методы теории оптимального управления, системного анализа конструкций и технологий польдерных систем, математической статистики, теории и практики математического моделирование сложных инженерно-технологических систем, численные методы решения нелинейных дифференциальных уравнений гидродинамики и теории фильтрации жидкостей через пористые структуры, решение оптимизационных многокритериальных задач, теория направленных графов.
Формализованное математическое описание объекта. Для декомпозиции задачи моделирования польдерных систем вводится элементная база в виде множества где n-количество базовых элементов, включающего подмножества: М1 - насосных станций (в том числе и регулируемых); М2 - линейных каналов и прилегающих к ним осушаемых массивов; М3 - множество точек ветвления проводящих каналов; М4 - множество элементов, оснащенных измерительной аппаратурой; М5 = М51 М52; М51 - точки входа (атмосферные осадки, испарения, паводки и др.); М52 - точки выхода (УГВ, РУКС). На элементную базу накладывается структура направленного графа Г, соответствующая расчетной схеме ПС. Элементы множества М2 являются ребрами направленного графа Г, а элементы множеств М1 и М3 - его вершинами. Все вершины перенумеруем. Каждому ребру поставим в соответствие картет (четыре числа и тройка)
l: mn+1 in+1cor(a, L, S, Kf (in, ik, N)), (1)
где a, L, S - ширина, длина, площадь сечения канала; Kf - коэффициент фильтрации почвы; in, ik- номера вершин соответствующих началу и концу ребра; N - номер ребра между этими вершинами. Соответствие (1) является изоморфизмом, т.е. для любого mn+1 in+1 соответствует только один cor и наоборот. Таким образом, каждой расчетной схеме ПС соответствует только один граф Г и наоборот.
Каждому элементу ПС mk Mk, где k 1,2,Е,n поставим в соответствие вектор состояний Fkik = (hkik, ukik, Hkik, Ukik, Lkik), hkik, ukik - уровень и скорость воды в канале; Hkik - уровень грунтовых вод; Ukik - поток влаги в зоне аэрации почвы; Lkik - принимает значение 0 при нарушении технологических ограничений и 1 в противном случае. Параметры h, u, H, U для каждого базового элемента ПС могут быть разными, однако они не являются независимыми и связаны между собой уравнениями связи (математические модели базисных элементов ПС)
Nk ik(hkik, ukik, Hkik, Ukik) = 0. (2)
Математическое описание объекта также является формализованным описанием информационной базы данных программных комплексов моделирующих ПС. Сформулируем общую математическую постановку задачи моделирования технологического процесса и управления ПС на период времени t[t0, ty]. Формализованная постановка состоит в следующем. Информационное обеспечение расчетов ПС содержит:
- информацию о конфигурации объекта (т.е. cor для любого mkMkik Г), состояние насосных станций при t=t0 и при t=tr (регламентные по времени изменения состояний) (m1 i1 M1 t J0=[t0, tr]);
- информацию о предыстории процесса (т.е. фактические данные в дискретные моменты времени для идентификации параметров математических моделей, например, S, Kf) [(h*, u*, H*, U*, *,)( m4 i4 M4)], Q*( m1 i1 M1 Г), здесь звездочкой обозначены фактические значения параметров; Q*- параметры управления насосной станцией; *- параметры внешней среды;
- для статистического прогноза входных (выходных) параметров необходимо иметь фактические данные по более глубокой ретроспективе (или моделирующие их графики)(h*, u*, H*, U*, *,)( m5i5 M5 Г; t-j; j J1);
- кроме того, необходимо знать плановые показатели на выходе
(h*, u*, H*, U*,)( m5i5 M52 Г; t[t0, ty]); (3)
- информацию о состоянии технологического процесса (нарушение технологических ограничений) при
{Lkik}( mkik Mk Г; t = t0). (4)
По информационному обеспечению идентифицируется конфигурация ПС при t = t0, определяется начальное состояние параметров и элементов, состояние технологического процесса (включая аварийное), дается прогноз параметров входа-выхода на период управления t[t0, ty].
Пусть F = {Fkik} все возможные вектора состояний элементов ПС. Требуется найти частное подмножество (режим работы ПС) F0 F для t[t0, ty], чтобы выполнялись следующие условия:
отклонения от плановых показателей должно быть минимальным; если режимные параметры находятся за пределами технологических ограничений, то время выхода на нормальный режим работы ПС должно быть минимальным; параметры режима должны удовлетворять уравнениям связи, т.е. Nk ik(hkik, ukik, Hkik, Ukik) = 0 для t[t0, ty]; расчетные параметры должны адекватно отражать реальный процесс, т.е. для любого t [t0, ty] должен достигаться минимум суммы квадратичных отклонений расчетных значений от фактических (отсюда следует задача идентификации параметров модели ПС); должны выполняться условия качества, т.е. технологические параметры не должны выходить за пределы ограничений.
Выбор оптимального режима F0 F, t[t0, ty] проводится на множестве управлений {Q}(mn in Mn Г) с соблюдением условия минимума энергозатрат.
Сложность поставленной задачи моделирования и управления ПС заключается в том, что польдерные системы, состоящие из сети большого числа проводящих каналов и взаимодействующих с ними ОМ, имеющих разнообразные структуры, обусловлена необходимостью использования наукоемких технологий системного подхода к математическому моделированию и управлению ПС, включающего семантическую модель технологического процесса, декомпозицию польдерной системы на элементную базу на основе анализа ее конфигурации, теорию графа объекта, синтез разрешающей системы дифференциальных уравнений, методы интегрирования
Расчет оптимальной вегетации сельскохозяйственных культур
Рис. 1. Состав и структура задач, подлежащих решению при реализации интегральной трехмерной нестационарной математической модели ПС и управления УГВ и режимом увлажнения корнеобитаемого слоя почвы
моделирующей системы дифференциальных уравнений вдоль ветвей обхода графа ПС, стратегию оптимального управления УГВ и РУКС, синтез вектора целевой функции и функциональных ограничений, методы многокритериальной оптимизации и интерпретатор полученных результатов.
На рис. 1 представлен состав и структура задач подлежащих решению на стадии реализации ИТНММ ПС и управления РУКС почвы. В состав задач входит двенадцать задач подлежащих решению, а структура задач содержит два контура управления - управления режимом увлажнения корнеобитаемым слоем почвы (блок 1) и управления УГВ (блок 5), связанных критерием управления уровнем грунтовых вод УГВ Нор, при t min.
Остальные задачи направлены на разработку интегральной трехмерной нестационарной математической модели ПС, состав и структура которой приведена на рис.2. При математическом моделировании польдерные системы представляются в виде множеств пассивных, активных элементов и точек соединений проводящих каналов, на которые накладывается структура направленного графа, соответствующая расчетной схеме ПС. Каждый элемент множеств получает свое математическое описание в рамках теории гидродинамики течений (фильтрации) жидкости. Синтез разрешающей системы моделирующих уравнений, описывающих динамику воды в сети проводящих открытых каналов с соответствующими начальными условиями, осуществляется на ребрах, а формирование граничных условий, выражающих законы сохранения расхода воды, в вершинах графа ПС, которые являются условиями сшивания системы моделирующих уравнений вдоль ветвей обхода графа ПС. Производительность активных элементов - насосных станций, при их наличии в ветвях графа ПС, используются в качестве соответствующих граничных условий. Моделирование режима УГВ и РУКС в ОМ, прилегающих к сети проводящих открытых каналов, основывается на численном решение двухмерного уравнения Буссинеска и уравнения для капиллярного потенциала и задании естественных граничных условий взаимодействия СПОК и ОМ. При наличии дренажных систем к указанной системе дифференциальных уравнений подключается и дифференциальное уравнение переноса воды в дренажных трубах.
Таким образом, задача моделирования режима увлажнения корнеобитаемого слоя в ОМ сводится к решению трехмерной нестационарной задачи, в которой пространственно - временное распределение УГВ вычисляются в горизонтальной, а потоков влаги - в вертикальных плоскостях.
Состав и иерархическая структура интегральной трехмерной нестационарной математической модели ПС и управления РУКС, соответствующая решаемым задачам на стадии реализации этой модели, приведена на рис. 2.
Синтез разрешающей системы уравнений удобно проводить на основе описаний отдельных каналов, входящих в каждый контур обхода графа ПС. Условием сшивания уравнений в точках соединения каналов (точка 2 на рис. 3), должно являться условие сохранения потока Q прихода и расхода воды в этих узлах.
ДЕКОМПОЗИЦИЯ ПОЛЬДЕРНОЙ СИСТЕМЫ НА ЭЛЕМЕНТНУЮ БАЗУ
НА ОСНОВЕ АНАЛИЗА КОНФИГУРАЦИИ ПС
РАСЧЕТ ЧИСЛА КОНТУРОВ ОБХОДА ГРАФА ПС
Контур обхода Контур обхода Контур обхода
L=1 L=i L=K
СОСТАВ L-ГО КОНТУРА ОБХОДА ГРАФА ПС
Открытые каналы с Открытые каналы с Открытые каналы с
номерами, для L=1 номерами, для L= i номерами, для L= K
СИНТЕЗ РАЗРЕШАЮЩЕЙ СИСТЕМЫ УРАВНЕНИЙ
СТРАТЕГИЯ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ УГВ И РУКС
Рис. 2. Иерархическая структура интегральной трехмерной нестационарной математической модели ПС и управления УГВ и РУКС
Совокупная система уравнений интегрируется вдоль контуров обхода, а граничные условия задаются в виде расхода (прихода) воды естественным образом. Интегрирование вдоль контуров обходов циклов графа ПС проводится только для решения уравнений Сен-Венана, описывающих уровневый режим в каналах, в то время как остальные дифференциальные уравнения интегрируются с привязкой к каждому открытому каналу.
4 x3
Y y2 y4
y1 y3 III
1 I 2 II 3 x2
x1
IV 5
x4
0 X
Рис. 3. Расчетная схема ПС, состоящая из четырех проводящих каналов
Таблица 2
№ | Cоедин. каналов | Длина, м |
1 | 1-2 | 1200 |
2 | 2-3 | 1800 |
3 | 2-4 | 2000 |
4 | 2-5 | 1600 |
1 - ЫЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ: 1, 2
И ВЕРШИН: 1, 2, 3;
2 - ОЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ: 1, 3
И ВЕРШИН: 1, 2, 4;
3 - ИЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ: 1, 4
И ВЕРШИН: 1, 2, 5.
Уровни и скорости движения воды в проводящих открытых каналах описываются нестационарной системой нелинейных уравнений Сен-Венана в частных производных гиперболического типа:
(5)
где h, u -уровень и скорость движения воды в канале, м и м/с; uq- скорость бокового притока, м/с; q - боковой приток, м2/с; g - ускорение свободного падения, м/с2; J0 - продольный угол наклона дна: Jf - уклон трения; Q - расход воды в канале, м3/с; F - площадь живого сечения, м2.
Уровень грунтовых вод описывается нестационарным двумерным уравнением Буссинеска в частных производных параболического типа:
, (6)
где Н - уровень грунтовых вод, м; - коэффициент водоотдачи; Kf - коэффициент фильтрации, м/с; - функция источника (стока) влаги, м/с.
Динамика движения воды в дренах описывается дифференциальным уравнением в частных производных гиперболического типа следующего вида:
, (7)
где Qдр - расход воды в дренах, м3/с; , d - площадь сечения и диаметр ДТ, м2 и м; Н0 - начальное значение УГВ, м; р - глубина закладки ДТ, м; - коэффициент гидравлического сопротивления.
Качество польдерных систем определяется интенсивностью снижения и равномерностью распределения уровня грунтовых вод в осушаемом массиве. Факторы, влияющие на интенсивность равномерного снижения уровня грунтовых вод R, могут быть описаны следующей многомерной функцией:
R = f (F, Qн, Ho, t, a, d, р, L1, L2, Кf, hO, hB, Ф), (8)
Время, необходимое для снижения уровня грунтовых вод до интервала определенных значений, необходимого для вегетации растений с учетом функционала (6), определяется как:
t = Т(F, Qн, Ho, a, d, р, L1, L2, Кf, hO, hB, Ф). (9)
На все параметры, входящие в функциональные зависимости (8) - (9), накладываются функциональные ограничений параметров технологического процесса на основе учета факторов, дестабилизирующих эффективный режим увлажнения в корнеобитаемом слое мелиорированных почв.
Таким образом, в рамках интегральной трехмерной нестационарной математической модели ПС и управления РУКС система дифференциальных уравнений в частных производных (5) - (7) описывает пространственно - временные распределения уровня воды в СПОК и УГВ в горизонтальных плоскостях.
К этой системе уравнений необходимо добавить дифференциальное уравнение капиллярного переноса потенциала влаги F в вертикальной плоскости, учитывающего транспирацию воды корневой системой растений G:
, (10)
где λ∂∂ - поток влаги м3/с3; x, y и - коэффициенты влагопроводности почвы вдоль горизонтальной координатной плоскости x, y и вертикальной координатной оси z, м2/с; F - капиллярный потенциал почвенной влаги, м2/с2; z Цвертикальная координатная ось с началом координат на поверхности грунтовых вод и направлена вверх; = 0.02G - коэффициент скорости потери влаги в почве за счет транспирации воды корневой системой растений G.
Алгоритм численного решения трехмерных дифференциальных уравнений параболического типа является весьма трудоемким по затратам машинного времени. Однако, стратегия управления уровнем грунтовых вод, направленная на установление равномерного распределения УГВ в горизонтальной координатной плоскости Х, У, позволяет не учитывать горизонтальные переносы влаги от поверхности грунтовых вод, вызванные практически нулевыми горизонтальные градиентами влажности почвы. Указанное благоприятное обстоятельство РР УГВ, позволяет рассчитывать только вертикальный перенос капиллярного потенциала влаги от поверхности грунтовых вод, описываемый одномерным дифференциальным уравнением:
. (11)
Начальные условия. В начальный момент времени при t0 =0 для системы нестационарных уравнений в частных производных (5) - (7) и (11) начальные условия задаются в виде:
h(x, t0) = H(x,y,t0) = 3; u(x, t0) = Q (x,y,t0) = 0. (12)
В качестве начальных условий для функции капиллярного потенциала влаги F, описывающего уравнением (11), задавалась линейная функция таким образом, что на поверхности почвы (z = H) F = 0, а на поверхности грунтовых вод (z = 0), значение функции F соответствовало значению потока влаги Umax.
Граничные условия. Для системы уравнений Сен-Венана (5), описывающих уровневый режим в СПОК на концах контуров ПС, задаются значения расхода воды QH(t) = QН, равные производительности насосов QН, если насосы (насос) отсутствуют, то естественно, QH(t) = 0.
Граничные условия для дифференциального уравнения Буссинеска (6), описывающего уровневый режим грунтовых вод, задаются для каждого канала на четырех сторонах прямоугольника осушаемого массива, прилегающего к данному каналу, следующем образом:
∂∂; ∂∂; ∂∂ , (13)
где L1 - длина открытого канала и осушаемого массива, м; L2 - ширина осушаемого массива, м; LHK-дополнительное фильтрационное сопротивление на контуре взаимодействия руслового потока с грунтовым потоком, м.
Для дифференциального уравнения (10), определяющего расход воды в дренах, краевые условия задаются только на левой границе для каждого ОМ в виде (см. подпункт 5.2.1), при УГВ выше уровня воды в проводящих каналах (процесс осушения мелиорированых земель) > 0:
, (14)
а при < 0
. (15)
На поверхности почвы (z = H) граничное условие для капиллярного потенциала влаги F, задавалось в зависимости от значения потока суммарного испарения E и интенсивности осадков X. В данных расчетах поток Е принимался равным нулю и значение производной капиллярного потенциала влаги F по координате z, через поток влаги U определялся следующем образом:
, (16)
где Е- поток суммарного испарения, м3/с; Х - поток атмосферных осадков, м3/с .
Выражение для вычисления значения капиллярного потенциала влаги на поверхности грунтовых вод (z = 0) через частную производную ∂∂ определяется как
. (17)
Таким образом, приведенная выше система нелинейных одномерных и двумерных уравнений в частных производных, совместно с начальными и граничными условиями, представляет собой замкнутую систему моделирующих уравнений, позволяющих рассчитывать параметры ПС.
Интерпретатор полученных результатов. Так, как интегрирование моделирующей системы дифференциальных уравнений проводится в локальных системах координат zl (ЛСК l-го канала), в которых ось хl направлена вдоль проводящего канала, то при анализе численных результатов расчетов, необходимо значения рассчитанных параметров представлять в общей системе координат. Вектор преобразований координат zl из ЛСК в общую систему координат Zi (ОСК) имеет вид:
Zl = Z0l + C(l)DLl , (18)
где Z0l Ц вектор координат начала l-го проводящего канала в ОСК; C(l)-матрица направляющих косинусов; DLl-вектор длин l-го проводящего канала в ЛСК; l = 0, Е, k; k - число проводящих каналов ПС.
Вектор преобразования координат (19) в матричной форме записывается следующим образом:
, (19)
где dl - длина l-го проводящего канала в ЛСК.
Таким образом, алгоритм вектора преобразований координат zl ОМ и проводящих каналов из ЛСК в общую систему координат Zl (ОСК), полностью определен соотношениями (18) - (19), а, следовательно, все искомые функции, вычисленные в узлах разностной сетки осушаемого массива, могут быть интерпретированы в ОСК.
В четвертой главе решена составная часть задачи моделирования и управления ПС - задача моделирования параметров сетей проводящих каналов, в том числе и с учетом рельефа местности реальных ПС. Приведен разработанный эффективный алгоритм численного решения системы уравнений Сен-Венана вдоль ветвей направленного графа ПС, основанный на приведении к каноническому виду исходной системы уравнений, имеющих недивергентные члены, с учетом числа обхода каждого элемента СПОК и введении новых прогоночных соотношений, что в целом позволяет корректно задавать граничные условия, строить разностные схемы с учетом характеристик, а также автоматически учитывать законы сохранения потока воды в СПОК ПС.
Численное решение системы уравнений Сен-Венана, имеющих недивергентные члены, не может быть непосредственно построено из-за двух проблем. Запишем систему уравнений Сен-Венана (8) в следующей форме:
(20)
где ; ; ; а = 4 м - ширина канала.
Первая проблема, связанная с недивергентным видом уравнений Сен-Венана и задания граничных условий в зависимости от наклонов характеристик решается сравнительно просто. Система дифференциальных уравнений (20), с помощью некоторых преобразований, приводится к каноническому виду (см. подраздел 2.1.3):
(21)
где ; -скорость звука в воде; ; ; .
При реалистическом предположении u < c уравнения (21) имеют два семейства характеристик - с положительным тангенсом угла наклона характеристик (первое уравнение dx/dt = u + c > 0) и с отрицательным (второе уравнение dx/dt = u - c < 0) и, следовательно, для первого уравнения характеристики приходят на правую границу, а для второго уравнения - на левую границу, а соответственно, численное решение определяется во всех точках характеристического треугольника, образованного характеристиками.
Вторая проблема, связанная с выполнением закона сохранения потоков в узлах расчетной схемы ПС, значительно сложней первой, рассмотренной выше. Оказалось необходимым изменить общепринятый алгоритм численного решения этой системы уравнений, основанный на методе матричной прогонки.
Это связано с необходимостью выполнения равенства входящих и выходящих из узлов потоков Q (закон сохранения потока в узлах проводящей сети), которое обеспечивается путем подсчета числа обхода каждого открытого канала и постулирования формулы, определяющей значения скорости движения воды в канале и некоторых других концепций. Рассмотрим более подробно такой алгоритм численного решения уравнений Сен-Венана, основанный на разностных схемах бегущего счета:
i=1,2,ЕNх Ц1; j= 0,1,ЕK - 1; =; =, (22)
где Nх и K - число пространственных узлов по оси Х и временных слоев на разностной сетке; , h - шаги интегрирования по времени t и оси Х.
Система разностных уравнений (22) записывается в явном трехточечном виде, следующем образом:
; (23)
где
;;;
=; .
Вектор искомых функций Тij+1 определяется следующим рекуррентным соотношением:
; i = n - 1, Е, 0; j =0, Е, k Ц1. (24)
Прогоночная матрица Ei и вектор Wi определяются в прямой прогонке (индекс i возрастает) по следующим формулам:
(25)
Обратная матрица в формулах (25) имеет вид:
; i = D11i⋅ D22i - D12i⋅ D21i;
D21i = B21i; D22i=B22i; D11i =B11i- (C11i⋅ E11i+ C12i⋅ E21i);
D12i = B12i - (C11i⋅ E12i +C12i ⋅ E22i).
Компоненты матрицы прогоночных коэффициентов Еi+1 определяются как:
E11i+1 = - A21i⋅ D12i/i; E12i+1 = - A22i⋅D12i/i;
E21i+1 = A21i⋅ D11i/i; E22i+1= = A22i⋅D11i/i,
а компоненты вектора прогоночных векторов Wi+1:
W1i+1 = (D22i ⋅W1i - D12i⋅W2i)/ i;
W2i+1 = (D11i⋅ W2i - D12i⋅ W1i)/ i, i = 0,Е, Nx - 1.
Для начала расчета прогоночных коэффициентов (25) необходимо задать их начальные значения в точке i = 0 с учетом наклона характеристик системы уравнений (21) и значений расхода воды на левой границе контура обхода графа ПС (для случая, когда в начале контура обхода установлена насосная станция). Так как на левую границу контура обхода ПС приходят характеристики второго уравнения, что следует из записанных в каноническом виде уравнений (21), то к граничному условию присоединяется это уравнение, и система разрешаемых уравнений на левой границе имеет вид:
- F0,j+1 u0,j+1 = - QT ; (26)
A210 F1,j+1 +A220 u1,j+1 - B210 F0,j+1 - B220 u0,j+1 = F20 .
Или в векторном виде:
B0 T0,j+1 = A0 T1,j+1 + F0 . (27)
Рекуррентное соотношение (24) для точки i = 0 в векторной форме записывается в следующем виде:
T0,j+1 = E1 T1,j+1 + W1 . (28)
Умножая слева соотношение (28) на обратную матрицу В0-1 и сравнивая полученное выражение с формулой (26), определяем начальные значения прогоночных коэффициентов E1 и W1:
E111 = A210 / B210; E121 = A220 / B210; E211 = 0; E221 = 0; (29)
W11 = - QT B220 / (F0,j+1 B210) + F20 / B210; W21 = QT / F0,j+1.
Искомые функции Fi,j+1,ui,j+1 в узлах разностной сетки рассчитываются обратной прогонкой (индекс i - уменьшается) на основе скалярных выражений:
Fi,j+1 = E11i+1 Fi+1 + E12i+1 ui+1,j+1 + W1i+1; i = Nx - 1,Е, 0;
ui,j+1 = E21i+1 Fi+1 + E22i+1 ui+1,j+1 + W2i+1. (30)
С целью учета сохранения расхода воды и его распределения по каналам в узлах расчетной схемы ПС, введем функцию скорости движения воды в русле канала в следующем виде:
ui+1,j+1 = (F1i - D11i Fi+1,j+1) / D12i . (31)
Тогда прогоночные соотношения (30), после подстановки формулы (36), преобразуются в рекуррентное выражение:
Fi, j+1 = (32)
ui, j+1 = ,
где KF1i+1 = E11i+1 - E12i+1 D11i /D12i; KF2i+1 = W1i+1 + E12i+1 F1i /D12i; KF3i+1 = F1i / D12i; KF4i+1 = D11i /D12i; s - количество обходов открытых каналов. Далее по тексту N = Nx.
При наличии откачивающих насосных станций на правых концах контуров обходов графа ПС задаются значения расхода воды QH. В результате решения квадратного уравнения определяется выражение для расчета значений площади живого сечения FN,j+1 и скорости движения воды в русле каналов uN,j+1, которое образуется при умножении второго уравнения системы уравнений (31) на значения FN,j+1:
uN,j+1 = KF3N / 2 - (KF3N2 / 4 - KF4N QH)0.5;
FN,j+1 = (KF3N - uN,j+1) / KF4N . (33)
Возможность введения новой функции (21), определяющей значения скорости движения воды в руслах каналов, обуславливается семейством характеристик первого канонического уравнения системы уравнений (21) dx / dt = u + c, которые приходят на правую границу при переходе на следующий временной слой. Поэтому, при задании граничных условий на правой границе контура должно использоваться это уравнение и значение расхода воды, аналогично методике задания граничных условий на левой границе (см. систему уравнений (26)). Проведя аналогичные расчеты расчетам граничных параметров на левой границе, можно получить формулу (31). Имея рассчитанные значения FN-1,j+1 и uN-1,j+1 и используя их в качестве УочередныхФ граничных условий, рассчитываются параметры в следующей точке N-2 и т. д.
Суммирование новых прогоночных коэффициентов KF по количеству обходов каждого канала в процессе обхода всех контуров графа ПС необходимо для учета разделения расхода воды Q при разветвлении одного канала на несколько и наоборот (закон сохранения расхода воды в узлах расчетной схемы ПС).
Таким образом, приведенный алгоритм численного решения системы дифференциальных уравнений Сен-Венана (5) автоматически просчитывает число обхода каждого проводящего канала, используя граф СПОК ПС, формируя реальный суммарный поток воды Q, который естественным образом разделяется в точках ветвления каналов.
Моделирование сетей проводящих каналов ПС. Для водо-изолированных и горизонтально расположенных проводящих каналов (приток воды в открытые каналы равен нулю q = 0), имеющих прямоугольное сечение, система дифференциальных уравнений Сен-Венана (8), записывается в следующем виде:
(34)
Начальные условия для системы уравнений (34) задаются следующим образом:
F(x, 0) = 12 м2; u(x, 0) = 0 м/с. (35)
Граничные условия задаются в виде потока воды Q, используя краевые условия 1Цго рода. На левой границе контуров обходов графа ПС в начале первого проводящего канала, при наличии насосной станции, задается ее производительность
Q(0, t) = 1.0 м3/с, (36)
а на правой границе контуров обходов графа ПС c учетом непроницаемости торцевых стенок открытых каналов
Q(xn, t) = 0.0 м3/с. (37)
Интегрирование системы дифференциальных уравнений СенЦВенана (34) вдоль контуров обходов графа Г польдерных систем осуществляется на основе алгоритма численного решения систем гиперболических уравнений (22) - (33).
Результаты численных расчетов параметров сети проводящих каналов реальной польдерной системы, действующей на территории Калининградской области, состоящей из тринадцати проводящих каналов различной длины и шести ветвей обхода графа ПС (рис. 4 и табл. 3) приведены на рис. 5.
Проводящие каналы имеют разную длину (см. табл. 3), но одинаковую ширину, равную а = 4 м, при этом шаги интегрирования по координате h = 15 м и по времени = 120 с задаются одинаковыми для каждого канала. Это связано с необходимостью, с одинаковой точностью аппроксимировать систему дифференциальных уравнений конечно-разностными уравнениями, описывающими динамику воды в проводящих каналах различной длины. Действительно, точность аппроксимации дифференциальных уравнений полностью определяется шагами интегрирования по координате h и времени и выбранной разностной схемой. Сравнительно небольшое значение шага интегрирования по времени определяется не условием счетной устойчивости алгоритма, а характерными временами процесса переноса воды в открытых каналах в момент запуска насосной станции.
У
7
VI
6 VIII
9
V
5 VII 8
IV IX 14
4 10
III XIII
3 X 11 XII 12
II
2
I XI
1
0 13
X
Рис. 4. Расчетная схема польдерной системы, состоящей из тринадцати проводящих каналов различной длины и шести контуров обходов графа ПС
Таблица 3.
Состав и структура ветвей направленного графа ПС
№ | Cоедин. каналов | Длина, м |
1 | 1-2 | 495 |
2 | 2-3 | 1000 |
3 | 3-4 | 1000 |
4 | 4-5 | 1000 |
5 | 5-6 | 1000 |
6 | 6-7 | 2500 |
7 | 5-8 | 1000 |
8 | 6-9 | 3900 |
9 | 4-10 | 3000 |
10 | 3-11 | 2500 |
11 | 11-13 | 2000 |
12 | 11-12 | 1500 |
13 | 12-14 | 1000 |
1 - ЫЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ:
1, 2, 3, 4, 5, 6 И ВЕРШИН: 1, 2, 3, 4, 5, 6, 7;
2 - ОЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ:
1, 2, 10, 11 И ВЕРШИН: 1, 2, 3, 11, 13;
3 - ИЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ:
1, 2, 3, 9 И ВЕРШИН: 1, 2, 3, 4, 10;
4 - ЫЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ:
1, 2, 3, 4, 7 И ВЕРШИН: 1, 2, 3, 4, 5, 8;
5 - ЫЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ:
1, 2, 3, 4, 5, 8 И ВЕРШИН: 1, 2, 3, 4, 5, 6, 9;
6 - ОЙ КОНТУР ОБХОДА ГРАФА СОСТОИТ ИЗ ДУГ:
1, 2, 10, 12, 13 И ВЕРШИН: 1, 2, 3, 11, 12, 14.
Действительно, как показывают численные эксперименты, подтверждающие экспериментальные данные, переходной процесс продолжается около одного часа в зависимости от площади живого сечения и длины проводящих каналов (см. рис. 5, b-d).
На рис. 5 приведены результаты численных расчетов пространственно-временных распределений площадей живого сечения F и потоков воды Q для открытых проводящих каналов ПС (время t приводится в часах, расстояния - в км). При этом использовались следующие обозначения: F1i1, F2i2, F3i3, F4i1, F5i2, F6i3 и Q1i1, Q2i2, Q3i3, Q4i1, Q5i2, Q6i3 - пространственные распределения площадей живого сечения F и потоков воды Q в каналах вдоль 1-го, 2-го, 3-го 4-го, 5-го и 6-го контуров обходов графа польдерной системы для момента времени t, равного одному часу; (Fj)1,0, (Fj)5,67, (Fj)6,167, (Fj)7,67, (Fj)10,167, (Fj)11,80, (Fj)13,40 и (Qj)1,0, (Qj)5,67, (Qj)6,80, (Qj)7,35, (Qj)4,65, (Qj)11,80, (Qj)13,40 - временные распределения площадей живого сечения F и потоков воды Q в каналах, для номеров пространственных узлов i, равных 0, 35, 40, 67, 80, 167, что соответствует значениям координат xi - 0 км, 0.525 км, 0.6 км, 1.0 км, 1.2 км, 2.5 км для 1-го, 5-го, 6-го, 7-го, 10-го, 11-го и 13-го проводящего канала.
Включение насоса, установленного в начале первого канала, порождает волну. Волновой процесс имеет затухающий характер со временем затухания порядка одного часа (см. рис. 5 b-d). По результатам численных расчетов, приведенных на рис. 5 b, амплитуда колебаний уровня воды Н (Н= F/a) относительно меньше, чем амплитуды колебания скорости воды в проводящих каналах (здесь и далее - (Qj)5,67 - значение потока воды Q на j-временном слое 5-го канала в 67-ом пространственном узле). Сдвиг фаз волновых процессов для различных проводящих каналов польдерной системы обуславливается временем добегания волны до соответствующего канала.
Отметим, что в точках ветвления проводящих каналов потоки воды Q делятся пропорционально длинам расходящихся каналов (см. рис. 5 с, Q1i1- распределение потока воды Q вдоль 1-го контура обхода графа Г ПС), так как в более коротких проводящих каналах быстрее падает уровень воды и, соответственно, значения потоков имеют меньшие значения по абсолютной величине (знак минус указывает, что поток воды Q направлен в противоположную направлению оси Х сторону).
Рис. 5. Пространственно-временные распределения площадей живого
сечения открытых каналов F, потоков Q и скоростей воды u для тринадцати проводящих каналов различной длины и шести контуров обходов графа ПС
Расчет параметров СПОК ПС с учетом рельефа местности. Для водо-изолированных проводящих каналов (приток воды в каналы равен нулю, q = 0), имеющих прямоугольное сечение (F = a h, где а - ширина проводящего канала), система уравнений Сен-Венана (8) записывается в следующем виде:
(38)
Для двух соединенных вместе проводящих каналов длиной l1 и l2, вертикальное сечение с обозначением начального h0 и установившегося уровня воды h1, h2 (ломаная и сплошная прямая, соответственно), приведено на рис. 6. При этом продольный угол наклона дна первого проводящего канала равен J01 = 0, а второго - J02 = Hyk2 / l2 (Hyk2 и l2 Ц высота уклона и длина второго проводящего канала).Как следует из геометрических построений, приведенных на рис. 6, задачу математического моделирования распределения уровней воды в сети проводящих каналов необходимо решать в два этапа в общей системе координат Н0Х. На первом этапе, на основе численного интегрирования системы дифференциальных уравнений Сен-Венана (38), вычисляется установившийся уровень воды вдоль контуров обходов графа польдерной системы, исходя из начальных значений уровней воды h0, заданных в локальных системах координат (ЛСК) с осями х1 и х2, направленными вдоль проводящих каналов и преобразования этих уровней в общею систему координат (ОСК) Н0Х (пунктирная кривая на рис. 6). На втором этапе, с учетом углов наклонов J01 и J01 координатных осей х1 и х2 ЛСК, связанных с первым и вторым проводящим каналом к координатной оси Х ОСК, значения уровней воды h, вычисленные вдоль контуров обходов графа ПС, преобразуются из ОСК в ЛСК, связанными с проводящими каналами, в которых значения уровней воды равны h1 и h2.
H
h0 h1 h2= h1 - J02Х2
J02 Hyk2
l2
l1 0 Х
Рис. 6. Установившийся уровень воды в двух взаимодействующих проводящих каналах (сплошная кривая), имеющих разные углы наклонов дна J0.
Интегрирование системы дифференциальных уравнений СенЦВенана (38) вдоль контуров обходов графа Г польдерных систем осуществляется на основе алгоритма численного решения систем гиперболических уравнений (22) - (33).
Начальные условия для системы уравнений (38) задаются в виде:
h1(x, 0) = 3 м; h2(x, 0) = 3 + J02x2 м; u(x, 0) = 0 м/с. (39)
Граничные условия задаются в виде потока воды Q, используя краевые условия 1Цго рода. На левой границе контуров обходов графа ПС, при наличии насосной станции, задается ее производительность:
Q(0, t) = 1.0 м3/с, (40)
а на правой границе контуров обходов графа ПС c учетом непроницаемости торцевых стенок открытых каналов:
Q(xn, t) = 0.0 м3/с. (41)
Результаты численных расчетов. Рассмотрим польдерную систему, состоящую из четырех открытых проводящих каналов одинаковой длины и ширины (а = 4 м), но с разными продольными углами наклона дна J0, имеющую три контура обходов графа (см. рис. 3 и табл. 2). Значения шагов интегрирования по координате х (h = 12 м) и времени t ( = 120 с) задавались равными для каждого проводящего канала. По нелинейности и связанности системы дифференциальных уравнений Сен-Венана (38), проводились три итерации по всем контурам обхода графа ПС.
На рис. 7 приведены результаты численных расчетов пространственно - временных распределений площадей живого сечения F и потоков воды Q для открытых проводящих каналов ПС (время t приводится в часах, расстояния - в км). При этом использовались следующие обозначения: h1i1, h2i2, h3i3 и Q1i1, Q2i2, Q3i3 - пространственные распределения уровней h и потоков воды Q в каналах вдоль 1-го, 2-го и 3-го контуров обходов графа польдерной системы для момента времени t, равного двум часам; (Fj)1,10, (Fj)2,250, (Fj)3,250, (Fj)4,250 и (Qj)1,10, (Qj)2,10, (Qj)3,10, (Qj)4,10 - временные распределения площадей живого сечения F и потоков воды Q в каналах, для номеров пространственных узлов i, равных 10, 250, что соответствует значениям координат xi - 0.12 км, 3 км для 1-го, 2-го, 3-го и 4-го проводящего канала.
Рис. 7. Пространственно-временные распределения площадей живого сечения F, уровня h и потоков воды Q для открытых проводящих каналов ПС, для случая а)
В результате процесса установления уровней воды h в двух соединенных вместе проводящих каналах (время установления около двух часов, см. рис. 7 b-d), имеющих разные продольные углы наклонов дна J0, установившиеся уровни воды h в ОСК (рис. 7 а), имеют разные линейные распределения в зависимости от значений продольных углов наклонов дна каналов J0. Наиболее интенсивные перераспределения уровней воды h в проводящих каналах происходит в течение первого часа процесса и носит волнообразный характер, о чем и свидетельствуют временные распределения площадей живого сечения F и потоков воды Q, приведенные на рис. 7 b-d. Установившиеся значения площадей живого сечения проводящих каналов F, выбранные в точках - х250 = 3 км, обратно пропорциональны продольным углам наклонов дна J0 (см. рис. 10 а). Абсолютные значения потоков воды Q для открытых проводящих каналов ПС при времени t 2 часов стремятся к нулю, что следует из рис. 10 с-d. Пространственные распределения уровней h и потоков воды Q имеют линейную зависимость, что характерно и для горизонтально расположенных проводящих каналов.
В пятой главе сформулирована и реализована задача двухмерного моделирования динамики уровневого режима грунтовых вод польдерных систем в горизонтальной плоскости, являющая составной частью ИТНММ ПС. Особое внимание при этом обращалось на разработку экономичных численных методов решения двумерного дифференциального уравнения Буссинеска и переноса воды в ДТ, увязки граничных условий на контуре взаимодействия руслов сети проводящих каналов и осушаемых массивов с учетом стока воды по ДТ, а также построению инвариантного алгоритма, синтезирующего результаты расчетов УГВ, полученные в ЛСК и преобразования их в ОСК.
Расчет УГВ с учетом дренажа и заданием граничных условий в области сопряжения проводящей сети с осушаемым массивом.
Граничные условия. Для вывода формулы, позволяющей рассчитывать значения потока воды в ДТ на левой границе у = 0, рассмотрим построения, приведенные на рис. 8.
УГВ
h ghp gHp Нр H H0 р
d
zp=H0 - p
Рис. 8. Схема системы проводящего канала и прилегающего к нему ОМ для вывода формулы значения потока воды в ДТ на левой границе y= 0 ОМ
Запишем уравнение Бернулли в вертикальных плоскостях Н0, d и h, пренебрегая членами уравнения с квадратами скоростей подъема (спуска) воды в канале и ОМ по сравнению с квадратом скорости воды в ДТ, следующим образом:
, (42)
где - плотность воды, кг/м3; V - скорость воды при выходе из дренажной трубы, м/с; р - глубина закладки ДТ, м; hp= h - zp; Hp= H - zp; Нр - расстояние между УГВ и горизонтом закладки дрен.
Выражая явно скорость воды при выходе из дренажной трубы V и, учитывая, что поток воды в ДТ, равен Q = V, получим соотношение для вычисления краевого условия на левой границе для дифференциального уравнения (6) в следующей форме, при УГВ выше уровня воды в проводящих каналах (процесс осушения мелиорированых земель)
. (43)
Уравнение движение сплошной среды по трубопроводу относительно потока жидкости Q записывается в виде:
, (44)
где Р - давление движущей среды; F - внешние массовые силы.
С учетом построений и обозначений, приведенных на рис. 11, дифференциальное уравнение (44), для переноса воды в ДТ, запишется в следующей форме:
. (45)
Алгоритм численного решения уравнения переноса воды в дренах. Запишем дифференциальное уравнение (45) в следующем виде:
, (46)
где ; .
C учетом знака Q > 0 разностная схема бегущего счета для дифференциального уравнения (46) запишется следующим образом:
, (47)
где m = 1, Е, M; j = 0,Е, k-1.
Умножая левую и правую часть разностного уравнения (47) на , приводя подобные члены, явно выразим искомую функцию Qmj+1
, где ; (48)
m = 1, Е, M; j = 0,Е, k-1,
Такая методика построения численного решения двухмерного уравнения Буссинеска подразумевает проведение итерационного процесса до получения необходимой точности.
Таким образом, присоединяя краевое условие на левой границе (43) к разностному уравнению (47), построенному по схемам бегущего счета, получим экономичный с вычислительной точки зрения алгоритм численного решения дифференциального уравнения (46), описывающего пространственно-временное распределение потока воды вдоль дренажных труб.
Алгоритм численного решения двумерного дифференциального уравнения Буссинеска (6) реализовался на основе двухшаговой консервативной разностной схемы, с введением полуцелых временных слоев, обеспечивающей второй порядок точности по обеим переменным, т.е. 0(2+h2) и счетную устойчивость. Следуя приведенному алгоритму в гл. 2, дифференциальное уравнение (6) запишется в виде следующих двух разностных уравнений:
(49)
где hy - шаг интегрирования по оси у; i = 1, 2, Е, n - 1 - номер пространственного узла по оси х; m = 1,2,Е,M-1 - номер пространственного узла по оси у; j = 0,Е,k-1 - номер временного слоя; hy - шаг по оси у.
Алгоритм численного решения двумерного уравнения (6) на разностной сетке строился на основе Т-образных разностных шаблонов (см. гл. 2, рис. 2.7 - 2.8, шаблоны выделены жирными линиями).
Моделирование УГВ реальных польдерных систем. Особенности расчета УГВ реальных польдерных систем связаны с необходимостью моделировать уровневый режим грунтовых вод ОМ в общей системе координат, математически увязывая модели расчетов УГВ для одиночных ОМ, выполненных в локальных системах координат с осью Х, направленной вдоль проводящего канала (см. рис. 3). На рис. 9 приведен фрагмент польдерной системы, состоящий из l-1, l и l+1 проводящих каналов и прилежащих к ним ОМ с указанием их длины: dl-1, dl и dl+1 и ширины ОМ: L2l-1, L2l и L2l+1.
Y yl+1
yl-1
L2l+1 dl+1
xl
L2l-1 Xl+1
dl
xl-1 L2l yl
X0l-1,Y0l-1 dl-1 X0l,Y0l
0 X
Рис. 9. Фрагмент ПС, состоящий из l-1-го, l-го и l+1-го проводящих каналов, с указанием ЛСК и дренажных труб (выделенные прямые)
Для фрагмента польдерной системы, изображенного на рис. 9 синтезируем общую матрицу преобразований координат С(L) из ЛСКL (L= l - 1, l, l +1) в ОСК по формулам (18) и (19), исходя из преобразований координат в следующем виде:
, где С(L)= (50)
С(L) - матрица преобразований координат из ЛСКL в ОСК, сформирована путем записи на главной диагонали матриц преобразований для каждого канала, при этом остальные матричные элементы равны нулю.
Необходимо отметить, что приведенный выше алгоритм синтеза общей матрицы преобразований координат С(L) из ЛСКL (L= l - 1, l, l +1) в ОСК, может легко быть распространен на ПС с произвольным числом проводящих каналов и ОМ, что позволяет говорить об инвариантности данного алгоритма.
Результаты численных расчетов. В качестве примера расчета основных параметров реальных ПС, включая расход воды в дренах, задавалась польдерная система, состоящая из трех последовательно соединенных проводящих каналов: длиной dl-1 = 300 м (Nxl-1 = 60), dl = 100 м (Nxl = 50), dl+1 = 300 м (Nxl+1 = 60) и шириной прилегающих к ним ОМ: L2l-1 = 150 м (Nyl-1 = 15), L2l = 150 м (Nyl = 15), L2l-1 = 100 м (Nyl-1 = 15,). Шаг интегрирования по времени составлял 3600 с и в момент времени, когда уровень воды в канале достигал 2.2 м насос выключался. Далее повышался до 2.4 м и насосная станция включалась (решались совместно система дифференциальных уравнений Сен-Венана для проводящего канала, двухмерное дифференциальное уравнение Буссинеска и дифференциальное уравнение переноса воды в дренах (5.10)), моделировалась закладка дрен диаметром d = 0.15 м на глубине р = 2.5 м, перпендикулярно руслу канала. При этом междренное расстояние составляло 25 м, и дрены, закладывались по всей длине осушаемого массива (канала).
На рис. 10 приведены результаты численных расчетов пространственных зависимостей УГВ в общей системе координат Х0Y вдоль осей Х - (а), (с) и Y - (b), (d).
Рис. 10. Пространственные распределения УГВ вычисленные в ОСК
вдоль осей Х - (а), (с) и Y - (b), (d)
Верхних два рисунка (см. рис. 10 а - b) подразумевают нижнюю часть рис. 9, до 3-го канала включительно и означают следующее: распределение УГВ вдоль оси координат хl-1 на расстояние у2 = 20 м (а) и вдоль - уl на расстоянии х5 = 25 м. Для верхней части рис. 9 приведено распределения УГВ вдоль оси Х ОСК, учитывающее влияние только 1-го и 3-го канала, отстоящее от нее на расстоянии Y = 110 м (с). На рис. 10 d указано распределения УГВ вдоль оси У ОСК, учитывающее влияние только 2-го и 3-го канала, отстоящее на расстоянии Х = 310 м. Необходимо отметить, что численные расчеты УГВ в ОСК дают адекватное представление о пространственно-временном распределении УГВ для всего ОМ, позволяющее оценивать качества ПС.
В шестой главе завершено решение задачи формирования и реализации ИТНММ ПС и управления РУКС почвы мелиорированных почв, подключением к горизонтальным процессам переноса влаги вертикального переноса влаги, моделирование которого основано на численном решении дифференциального уравнения капиллярного потенциала влаги (11) в зоне аэрации почвы. Заменяем исходное дифференциальное уравнение в частных производных параболического типа (11) конечно - разностным во внутренних узлах по консервативной схеме
τ⋅⋅λλ⋅λλ⋅λ i = 1, 2, . . . , N - 1, j = 0, 1, . . . , k-1. (51)
Эти разностные уравнения приводятся к явному каноническому трехдиагональному виду:
Ai Fi+1,j+1 - Bi Fi,j+1 + Ci Fi-1,j+1 = - Di , (52)
где λλ⋅ ; λλ⋅λ; τ ; τ.
Отметим, что критерий устойчивости счета методом прогонки к ошибкам округления выполнен при любых значениях τ и h, так как
⋅τ> , (53)
а значит полученная консервативная дивергентная разностная схема (51), аппроксимирующая исходное дифференциальное уравнение (11), безусловно - устойчива и аппроксимирует его с первым порядком точности по времени и со вторым по координате 0(τ, h2). Решение уравнения для капиллярного потенциала (11), после его приведения к трехдиагональному виду (52), реализуется на основе метода прогонки.
Управление режимом увлажнения корнеобитаемого слоя. Рассмотрим алгоритм расчета оптимального уровня грунтовых вод, например, для условий глубокозалежных болот (многолетние травы, U0 =86.5 мм). Оптимизируя целевую функцию потока влаги U(H, h) (см. раздел 6.3, (6.6)), при управление потоком влаги Umax max и Umax U0 и функциональных ограничениях 5 ≤ h ≤ 50 и 20 ≤ Н ≤ 60, вычисляем оптимальное значение уровня грунтовых вод Ноп = 20 см, высоту корнеобитаемого слоя почвы hоп = 14.9 см и максимально возможный поток влаги в зону аэрации Umax= 86.44 мм (см. рис. 11). По вычисленному значению максимально возможного потока влаги Umax от поверхности грунтовых вод, определяется краевое условие для капиллярного потенциала влаги F, при z = 0. На рис. 12 приведены значения капиллярного потенциала влаги F и потоков влаги U для трех моментов времени 1.5сут, 5сут и 10сут, рассчитанные на основе численного решение дифференциального уравнения (6.4) с начальными и граничными условиями (6.5) и (6.7) и оптимальном уровне грунтовых вод Ноп.
Рис.11. Варьирования потоков влаги в зону аэрации Ui,j для трех значений
высот корнеобитаемого слоя почвы h0= 5 cм, h10= 14 cм, hN= 50 cм
и уровней грунтовых вод Н0=20 см, Н10 = 28 см, НN = 60 cм
Необходимо отметить, что при постоянном значении уровня грунтовых вод Н и подпитывающего потока влаги Umax, значения капиллярного потенциала F со временем возрастают.
Рис. 12. Рассчитанные зависимости оптимальных значений капиллярного
потенциала F и потоков влаги U от высоты z для трех моментов времени
Функция капиллярного потенциала F от высоты z носит ярко выраженную экспоненциальную зависимость.
Расчет потока влаги от УГВ с учетом испарения, транспирации и выпадающих осадков. Испарение зависит в основном от метеорологических условий, солнечной радиации, условий водного режима и водно-физических свойств почвы. А транспирация корневой системой растений почвенной влаги определяется их биологическими свойствами, обусловленными взаимодействием среды, метеорологических условий, почвы и агротехники.
Метод теплового баланса основан на расчете испарения по затратам энергии на этот процесс. Зная все составляющие теплового баланса испарения водяного пара у поверхности почвы, величину испарения можно получить по его остаточному члену. Уравнение теплового баланса процесса испарения водяного пара имеет вид:
где Qпоч - поток тепла в почву; L - скрытая теплота испарения; E - испарение; Р - турбулентный отток тепла в атмосферу; Rб - радиационный баланс поверхности почвы.
На рис. 13 приведены расчеты пространственно-сезонные зависимостей капиллярного потенциала влаги F (а и c) и его потока U (b и d) учитывались процессы сезонной транспирации G.
Рис. 13. Пространственно-сезонные зависимости капиллярного потенциала
влаги F (а и c) и его потока U (b и d), поверхностное испарения Е
и транспирация G задается в виде сезонных зависимостей (e, f).
Следует отметить, что учет процесса сезонной транспирации G приводит к уменьшению значений капиллярного потенциала влаги F практически в 1.8 раза (см. рис. 13 а, е). В это же время, временные зависимости капиллярного потенциала влаги F и его потока U имеют изрезанный характер, повторяющий сезонный ход транспирации G (рис. 13 с и d). Влияние процесса испарения Е на пространственно-сезонные зависимости капиллярного потенциала влаги F (а и c) и его потока U (b и d), сказывается значительно слабее, чем процесса транспирация G. Увеличение капиллярного потенциала влаги F со временем определяется только грунтовыми водами, при этом пространственное распределение потока влаги U имеет ярко выраженную убывающую к поверхности почвы линейную зависимость
Приведенные численные результаты свидетельствуют об эффективности и работоспособности алгоритма расчета сезонных потоков влаги от грунтовых вод с учетом испарения с поверхности почвы, транспирации растений и выпадающих осадков.
В седьмой главе приведены результаты по расчету оптимальных значений параметров при проектирование совершенных конструкций ПС и эксплуатации уже существующих, таких как: производительность насосной станции, диаметр и глубина закладки дрен, междренное расстояние и некоторых других, свидетельствуют о применимости ИТНММ ПС и управления РУКС почвы ОМ при проектировании совершенных ПС и эксплуатации уже существующих.
В качестве примера расчета оптимальных значений варьируемых параметров совершенных конструкций ПС в процессе управления уровнем грунтовых вод УГВ Нор в осушаемых массивах, моделировалась польдерная система, состоящая из проводящего канала, длиной 200 м (число пространственных узлов вдоль оси х составляло Nх = 10) и прилегающего к нему осушаемого массива шириной 200 м (Ny = 10, длина ОМ соответствовала длине проводящего канала, шаг интегрирования по времени составлял =300 с).
Стратегия оптимального управления уровнем грунтовых вод, обеспечивающая энергосбережение и эффективность вегетации растений мелиорированных земель сводится к следующему: ИС УГВ max, BC УГВ min и РР УГВ min, при t min и учетом технологических ограничений.
Рис. 14. Временная разверстка уровней воды в канале (сплошная кривая)
и в осушаемом массиве для трех значений: у1=20м, у5=100м, у10=200м (а) и режим работы насосной станции (b)
На рис. 14 приведена временная разверстка процесса управления уровневым режимом грунтовых вод (a) и насосной станции (b), здесь и далее время в час. Отметим, что процесс управления уровневым режимом грунтовых вод УГВ Нор длительностью 1.25 ч, обеспечивает приведения уровней грунтовых вод по всему осушаемому массиву и воды в канале к заданному значению отметки УГВ в 2.8 м (см. рис. 14 а). При этом в первые двадцать минут насосная станция, как это следует из рис. 14 (b), обеспечивает максимально возможную производительность откачки воды из проводящего канала, равную QH = 0.4 м3/с (контуром управления является осушаемый массив т.к. все критерии качества ПС связанны с УГВ ОМ, а не проводящим каналом), что приводит к резкому спаду уровня воды в канале, соответственно и к спаду и выравниванию УГВ в ОМ. За тем, в ходе управления уровнем грунтовых вод УГВ Нор на заданную отметку, составляющую 2.8 м, насосная станция выключалась и кривые спада УГВ сближались с отметке в 2.8 м. Максимально возможная производительность откачки воды из проводящего канала, равная QH = 0.4 м3/с определяется стратегией оптимального управления УГВ УГВ Нор, для которой интенсивность снижения уровня грунтовых вод R (8) стремится к максимуму, а уровень воды в проводящем канале постоянно пополняется стоком воды из ОМ.
Из результатов по управлению уровнем грунтовых вод УГВ Нор, действующих польдерных систем в штатных и нештатных ситуаций следует, что поставленные цели стратегией управления уровнем грунтовых вод УГВ соответствующие управлению режимом увлажнения корнеобитаемого слоя почвы - поток влаги U Umax выполнены в полном объеме и с высокой точностью (за минимальное время управление УГВ устанавливается заданное значение уровня Нор, при его равномерном распределение по всему ОМ).
ЗАКЛЮЧЕНИЕ
1. Выполненный анализ современного состояния теории и практики актуальной задачи - повышения плодородия мелиорированных земель сделал возможным сформулировать цели и задачи научных исследований создания научно-обоснованной интегральной трехмерной нестационарной математической модели ПС, обеспечивающей моделирование уровневого режима грунтовых вод и на основе управления им рассчитывать оптимальный режим увлажнения корнеобитаемого слоя почвы по схеме - УГВ Нор - U Umax, при t min в области допустимых значений технологических параметров.
2. Впервые поставлена и разработана научно-обоснованная интегральная трехмерная нестационарная математическая модель ПС управления режимом увлажнения корнеобитаемого слоя почвы мелиорированных земель, основанная на применении единой методологии моделирования ПС и оптимального управления режимом грунтовых вод, учитывающая взаимодействие между сетью проводящих открытых каналов и осушаемых массивах, позволяющая сократить время на решение конкретных задач максимально возможного увеличения плодородия переувлажненных земель для каждого вида сельскохозяйственных культур. При этом проводится декомпозиция исходной постановки задачи моделирования ПС на подзадачи с последующим синтезом системы разрешающих дифференциальных уравнений, осуществляемым на ребрах направленного графа, с заданием граничных условий в вершинах графа.
3. Созданы алгоритмы и программное обеспечение численного решения разрешающей системы дифференциальных уравнений Сен-Венана, описывающих динамику воды в СПОК, для которого автоматически выполняются законы сохранения потоков воды в точках ветвления проводящих каналов; реализация алгоритмов расчета СПОК и реальных польдерных систем, в том числе с учетом рельефа местности.
4. Впервые поставлены и решены задачи моделирования: динамики УГВ и управления РУКС почвы ПС, основанные на двухмерном уравнении Буссинеска, дифференциальном уравнении переноса воды в дренажных трубах и трехмерном уравнении переноса капиллярного потенциала влаги, при задании естественных граничных условий: а) взаимодействия СПОК с ОМ; б) влагообмена между УГВ и корнеобитаемым слоем почвы; потоков влаги от грунтовых вод с учетом поверхностного испарения, транспирации растений и выпадающих атмосферных осадков, основанных на численном решении дифференциальных уравнений капиллярного потенциала и теплового баланса почвогрунтовой системы; сезонных потоков влаги от грунтовых вод с учетом поверхностного испарения и транспирации влаги.
5. Впервые создан алгоритм расчета оптимальной вегетации сельскохозяйственных культур; формирование целевой функции для потока влаги U(Н, h) и функциональных ограничения параметров, дестабилизирующих вегетацию растений; построена стратегия управления РУКС; алгоритм расчета реальных польдерных систем, основанного на преобразовании параметров ПС, вычисленных в локальных системах координат, связанных с проводящими каналами, в общую систему координат ПС.
6. Представлены результаты разработки алгоритмов повышенного порядка точности, численного решения дифференциальных уравнений необходимых для тестирования и выверки применимости экономичных разностных схем первого порядка точности, при проведении оптимизационных расчетов;
7. Практическое значение результатов работы определяется тем, что они нашли применение в комплексной мелиорации заболоченных и переувлажненных земель, без которой не могут быть решены вопросы продовольственной безопасности страны и направлены на повышения плодородия земель сельскохозяйственного назначения.
8. Научная и производственная значимость настоящей работы подтверждены опытно-промышленными результатами эксплуатации польдерных систем, положительными отзывами при апробации результатов монографии.
Основные результаты диссертации опубликованы в работах:
1. Бобарыкин Н. Д., Латышев К. С. Структура дневных потоков ионов в задаче с верхней границей. Всесоюзная конференция по физике ионосферы. Тезисы докладов, часть II, М., 1976, с. 34.
2. Бобарыкин Н. Д., Латышев К. С. Высотная структура скоростей и потоков ионов с учетом силы инерции и связанные с ней особенности численного решения моделирующих уравнений: Сб. науч. трудов Диагностика и моделирование ионосферных возмущений, 1978.- М.: Наука, с. 23-32.
3. Латышев К. С., Бобарыкин Н. Д., Медведев В. В. Разностные методы решения системы одномерных газодинамических уравнений в задачах моделирования ионосферы //Ионосферные исследования.-1979.-№28.- М.: Сов. Радио, с. 37-49.
4. Бобарыкин Н.Д., Латышев К.С. Особенности построения численных алгоритмов в задачах ионосферного моделирования. Тезисы докладов V-го Междуведомственного семинара по моделированию ионосферы, Тбилиси, 1980, с. 6.
5. Бобарыкин Н. Д., Латышев К. С., Осипов Н. К. Нестационарный полярный ветер - причины и следствия //Геомагнетизм и аэрономия.-1981.-Т. 21.- №4.- С. 698 - 703.
6. Бобарыкин Н.Д., Латышев К.С. Расчет ионосферных параметров вдоль геомагнитной силовой линии с учетом инерционных членов. Сб. Результаты обработки геофизических данных в МЦД Б-2, 1982.- М., Мировой центр данных АН СССР, с. 23-33.
7. Бобарыкин Н.Д., Латышев К.С. О роли вертикальных переносов в формировании плазмопаузы на больших высотах. Сб. Результаты обработки геофизических данных в МЦД Б-2, 1982.- М., Мировой центр данных АН СССР, с. 54-62.
8. Бобарыкин Н. Д., Латышев К. С., Осипов Н. К. Температурный режим и характерные времена нестационарности полярного ветра //Геомагнетизм и аэрономия.-1984.- №4,- С. 23 Ц30.
9. Бобарыкин Н. Д. Математическое моделирование технологических процессов в тренажерах установок газоперерабатывающих предприятий на базе персональных компьютеров. Диссертация в форме научного доклада на соискание ученой степени кандидата тех. наук - 05.13.16, М., МХТИ им. Д.И. Менделеева, 1991, с. 20.
10. Бобарыкин Н. Д., Бобарыкина Е.Н. Инвариантный метод расчета кинематических характеристик механизмов на ПЭВМ: Сб. науч. трудов КГТУ.-1998.- С. 41-43.
11. Бобарыкин Н. Д., Рябой В.Е. О применении метода комплексного математического моделирования динамики уровневого режима польдерных систем в процессе изучения сельскохозяйственных дисциплин: Сб. науч. трудов КГТУ.-1998.- с. 38-40.
12. Бобарыкин Н. Д. Состав и структура инвариантной нестационарной математической модели польдерных систем, включая алгоритмы оптимального управления режимом грунтовых вод осушаемого массива: Сб. науч. трудов КГТУ Математическое моделирование и численные методы решения интегрально-дифференциальных уравнений, 2003, с. 4-18.
13. Бобарыкин Н. Д. Расчет влагообмена с грунтовыми водами на основе дифференциального уравнения для капиллярного потенциала: Сб. науч. трудов КГТУ Математическое моделирование и численные методы решения интегрально-дифференциальных уравнений, 2003, с. 19-30.
14. Бобарыкин Н. Д. Расчет уровневого и скоростного режимов движения воды в каналах польдерных систем на основе решения дифференциальных уравнений Сен-Венана: Сб. науч. трудов КГТУ Математическое моделирование и численные методы решения интегрально-дифференциальных уравнений, 2003, с. 31-41.
15. Бобарыкин Н. Д. Расчетов уровневого режима грунтовых вод осушаемого массива, примыкающего к проводящему каналу на основе решения двумерного дифференциального уравнения Буссинеска: Сб. науч. трудов КГТУ Математическое моделирование и численные методы решения интегрально-дифференциальных уравнений, 2003, с. 42-45.
16. Бобарыкин Н. Д. Расчет уровневого режима грунтовых вод осушаемого массива с учетом дренажа в области сопряжения проводящей сети с осушаемым массивом: Сб. науч. трудов КГТУ Математическое моделирование и численные методы решения интегрально-дифференциальных уравнений, 2003, с.46-52.
17. Бобарыкин Н. Д. Оптимальное управление водным режимом осушаемого массива и выбор вектора целевых функций и ограничений для польдерных систем: Сб. науч. трудов КГТУ Математическое моделирование и численные методы решения интегрально-дифференциальных уравнений, 2003, с. 53-60.
18. Бобарыкин Н. Д. Оптимальное управление режимом грунтовых вод с учетом выпадающих атмосферных осадков: Сб. науч. трудов КГТУ Математическое моделирование и численные методы решения интегрально-дифференциальных уравнений, 2003, с. 61-65.
19. Рябой В. Е., Бобарыкин Н. Д. Структура, состав и основные задачи гидрологических расчетов польдерных систем: Сб. науч. трудов КГТУ Математическое моделирование и численные методы решения интегрально-дифференциальных уравнений, 2003, с. 91-105.
20. Бобарыкин Н. Д. Оптимальное управление режимом грунтовых вод на основе инвариантной нестационарной математической модели польдерных систем (монография, научное издание). - Калининград: КГТУ, 2004.-16ышев К.С. Моделирование уровня грунтовых вод реальных польдерных систем //Математическое моделирование. РАН.-2007, т. 19, №
43. Бобарыкин Н.Д., Графова Е.Н., Латышев К.С. Расчет уровневого режима грунтовых вод в осушаемых массивах польдерных систем с учетом дренажа //Математическое моделирование. РАН.-2007, т. 19, №
44. Бобарыкин Н.Д., Графова Е.Н., Латышев К.С. Численные расчеты оптимального режима увлажнения корнеобитаемого слоя почвы от поверхности грунтовых вод //Математическое моделирование. РАН.-2007, т. 19, №
45. Бобарыкин Н.Д., Графова Е.Н., Смертин В.М., Латышев К.С. Об использовании разностных схем первого порядка точности, аппроксимирующих дифференциальные уравнения в задачах моделирование польдерных систем. Тезисы докладов международной научной конференции: Российская наука и инженерная деятельность в социокультурном пространстве эксклавного региона: история, актуальные проблемы, перспективы развития, г. Калининград, КГТУ, 2007, с. 51-52.
46. Бобарыкин Н.Д. Математическое моделирование и управление режимом увлажнения корнеобитаемого слоя почвы от поверхности грунтовых вод. Тезисы докладов международной научной конференции: Российская наука и инженерная деятельность в социокультурном пространстве эксклавного региона: история, актуальные проблемы, перспективы развития, г. Калининград, КГТУ, 2007, с. 18-19.
Н.Д. Бобарыкин
Авторефераты по всем темам >> Авторефераты по техническим специальностям