CC BY
69
Применение сопряженных уравнений в задачах верификации и валидации расчетов
Алексеев А.К.*, Бондарев А.Е.** *РКК ««Энергия», Россия, Королев ** Институт прикладной математики им.М.В.Келдыша РАН, Россия, Москва e-mail: Aleksey.k.alekseev@gmail.com bondakeldysh.ru
1. Введение
Сопряженные уравнения в настоящее время широко применяются при решении оптимизационных и обратных задач в самых разных областях приложения вычислительной математики. Использование сопряженных уравнений позволяет определять чувствительность искомых функционалов в различных задачах к различного рода возмущениям. Развернутая теория сопряженных уравнений и различных способов их применения в практических расчетах представлена в работах [1,2]. В задачах вычислительной механики жидкости и газа сопряженные уравнения применяются при решении задач, связанных с необходимостью идентификации и управления изучаемым течением. Также теория сопряженных уравнений может быть применена в задачах верификации и валидации численных результатов и математических моделей в целом.
За последнее десятилетие в вычислительной математике и математическом моделировании сложилась новая самостоятельная дисциплина - верификация и валидация (Verification&Validation -V&V). Основная цель данной дисциплины -обеспечить уверенность исследователя и пользователя во всей цепочке исследования: физическая модель - математическая модель - численный алгоритм - программная реализация алгоритма - организация расчета - обработка и представление результатов. Установление новых стандартов надежности программ и результатов расчетов привело к формализации понятий верификации (подтверждения правильности численной модели) и валидации (подтверждения соответствия физической модели задаче реального мира) в вычислительной аэрогазодинамике [35]. В рамках V&V стандартные методы, включающие асимптотическую сходимость по сетке и проверку на тестовых задачах, дополняются методами, основанными на апостериорной оценке численной погрешности расчета [6-8], оценке погрешности выбора модели [9], В подобных задачах визуализация играет ключевую роль, например, при адаптации сетки [7] или анализе оптимального расположения датчиков [10]. Как отмечается в работе [11], именно синтез математического аппарата с инструментами визуализации позволяет создавать надежные и востребованные средства V&V.
Использование поля сопряженных параметров в работе [6] позволяет рассчитать и визуализировать вклады ошибки исходных данных, ошибки физической модели и погрешности дискретизации в погрешность целевого функционала. Это может использоваться для уменьшения погрешности и для получения решения с необходимой гарантированной точностью. Данная работа является логическим продолжением работы [6] в области применения сопряженных уравнений и визуального представления сопряженных параметров в качестве инструментов V&V.
Данная работа ставит своей целью показать возможность применения сопряженных уравнений к оценке погрешности расчета (верификация) и погрешности
физической модели (валидация). С применением сопряженных уравнений проводится оценка и прямая визуализация погрешности для различных случаев. Визуальное представление позволяет определить местоположение необходимых сеточных сгущений с целью уменьшения погрешности, а также непосредственно показать погрешность, вносимую влиянием тех или иных членов системы уравнений математической модели.
2. Задачи верификации
Рассмотрим существенные новые элементы верификации, а именно апостериорную оценку погрешности расчета (учитывающую регулярность данного решения и не содержащую неопределенных констант).
Апостериорная оценка погрешности расчета
Начнем с апостериорной оценки погрешности в ее варианте для конечно-разностных схем [6-8].
М/ = w в Ос Я" ,
/(дО) = /в(х)е ¿2(ЭО) . (1)
(управлений на границе нет).
Пусть численное решение получено некоторым конечно-разностным методом Мк/к = w. (2)
В результате решения получим сеточную функцию /. Предположим
существование функции / е Ик+" (О) совпадающей с сеточной в узлах. Конечные разности в N к /к разложим в ряд Тейлора, что дает дифференциальное приближение схемы [12].
М/ + 8к (/) = w , или М/ = w + А w . (3)
Здесь 8к (/) ошибка аппроксимации, содержащая члены ряда Тейлора и служащая возмущающим источником А w = -8к (/). Соответствующее уравнение
N г А/ = А w = -8к (/), О с Я", А/(дО) = 0. (4)
Согласно (16)
Ае = - (А w, V ) ¿2( а) = | 8к (/) V й О
а
Сопряженные уравнения:
М* V + е/ = 0, в О, V = 0, на дО .
(5)
(6)
В качестве примера рассмотрим расчет и визуализацию погрешности на примере двумерного сверхзвукового вязкого течения, рис. 1.
у
А
в
Цу) *
*
Б
С
Рис. 1
Используем недивергентную конечно-разностную аппроксимацию
параболизированных уравнений методом вдоль X.
д (pU ) + ЦррП =
U
U
дX д_Ц_ дX д_У_ дX
де
+ V
+ V
де
дY
д_и_
д Н
дУ_
дY
= 0 1
Р 1
Навье-Стокса. Расчет проводится маршевым
дP
дX дP
и— + V— + (к-1)е
дХ дН
д Н дУ ^
V дК д¥
1
д 2 и
Яе р 4
д Н д 2 V
= 0
3 р Яе д н
=0
1 К д2е 14
ди^2
р ЯеРг д¥2 р 3Яе
д¥
= 0
(7)
(8)
(9)
(10)
2
P= pRT; e= С T= Щк-1)Т; (Х,Г)е О=(0<Х< Xmax; 0<Г<1); На входной границе (А (Х=0), Рис. 1): е(0,Г)=е^ (Г); Р02)=р„ Г)); и(0,Г)=и^ (Г); У(0,Г) = У^ (Г); на боковых границах В, Б (Г=0, Г=1) выполнены условия с$/сГ=0. В качестве целевого функционала, погрешность которого мы будем рассчитывать, выберем температуру или плотность потока в некоторой точке.
Обозначим Т (ХеЛ , ) как е(8).
е(8/) = Г Т (X, Н )8(Н - )8(X - Хе*( )йхйу
О ( 11)
Нам нужно рассчитать градиент целевого функционала по отношению к погрешности исходных данных, погрешности модели и погрешности численной аппроксимации. Соответствующие сопряженные переменные рассчитываются с помощью следующей системы уравнений.
дЧ
дЧ
и-р+ (к- 1) д (Ч vе 7 Р) + (к- 1) д (Ч ие 7 Р)
дХ
К-1 ( де_
дН
дН
дН
^ +
дХ
и
Г
1 дР
1 д2и
р2 дХ Р2 Яе дН:
дХ дР
4 д^
дН 3Яе дН2
У» -
'7
К д2е + (и 2 ^
ЯеРгдН2 3Яе1 Н
V -8(X - Хе*)8Н ) = 0
(12)
Источник в уравнении для ^р соответствует расположению оцениваемой точки.
и
дг
дХ
и_ +д^*иХ± +р
дН
дН
р Яе
и
д¥р дХ д
' д¥_ V дХ
V V +
дН
8 ди
3Яе дН
де дХ
дХ
= 0
(13)
д (и V V ) + V дЩ V
+ р
дХ д*¥
дН дН
дН
Р т.
ди
V дН
4
V и
д 2 (
де ш
дН е
3Яе дН
= 0
(14 )
2
1
2
е
2
д (и У е ) д (V У е ) к- 1 ( др
дх
д¥
д¥
У V
др_ дх
У
и
- (к-1)
ди дV
дх дУ
У е +
, „ ч дУ V , ч ддУ и + (к- 1)-^ + (к - 1)-и
к
12 С ш Л т е
дУ дх Яе Рг дУ2 ^ р ) (15)
Величины (Ур, У, Уе) являются сопряженными аналогами плотности, компонент скорости и энергии.
С (Х-Хтах):
У
x=х„
и ,V ,е
= 0;
(16)
Начальные условия
и¥р +Уие 1р-д(У-Уы) = 0;
Выражение для Ур соответствует расположению оцениваемой точки на границе
Хтах-
Граничные условия Е,В (7=0; 7=7):
¿у
^ = 0;
дУ (17)
Сопряженная задача решается в обратном направлении по X. При решении сопряженной задачи используется тот же конечно-разностный алгоритм, что и для решения прямой задачи. Параметры течения, входящие в коэффициенты, берутся из решения прямой задачи.
Погрешность дискретизации в каждой расчетной точке можно представить в качестве ряда Тейлора, соответственно конечно-разностное представление задачи оказывается эквивалентно исходной форме задачи (в виде ЧДУ) с дополнительными источниковыми членами.
Здесь использована разностная схема первого порядка точности, позволяющая получить наименее громоздкие выражения для оценок ошибки и достаточно большую погрешность аппроксимации в сравнении с другими источниками ошибки (такими, как разрядность компьютера).
Для оценки точности конечно-разностного решения рассмотрен подход, основанный на использовании оценки погрешности с помощью ряда Тейлора в форме Лагранжа и сопряженных уравнений в непрерывной постановке.
Вариация целевого функционала в зависимости от ошибки дискретизации имеет
вид:
де
& = — &у = Ур,ур + Уи ,У#и у + Уv + Уе^е
диу
(18)
В качестве примера приведем оценку погрешности аппроксимации производных вдоль маршевой координаты. Для ее оценки разложим разностные
выражения в окрестности и%, используя ряд Тейлора в форме Лагранжа.
ипки"к+ - икп = и ди+1Щ
к кх дх 2 к
к
х ,п
д 2и (хп + ак%кх ,п, Ук) дх 2
(19)
Соответствующая компонента погрешности целевого функционала Аре^ имеет
вид
Ае(3и) = 1 \
а
к
д 2и (хп +аккх, п, Ук)
х, п
дх2
Уи^хйу
20)
(
№, N
1 1 , 1'
Ее дискретную форму Ае(д) = — 2
к,
Э 2П (хп + а1кх, Ук)
к=1, п=2\,
первом порядке точности акпкх,п можно представить как
Э*2
Ши, кку, ккх, п в
1
т, N1
Ае(3и) = - 2
к=1, п=2
к
х, п
V
Э 2и_(хй1ук)
Эх 2
+ аккх, п
3
2 Э и (хп»Ук )
Эх3
Шп, кку,ккх, п
(21)
Первая часть этой суммы может быть использована для уточнения функционала, она является аналогом систематической ошибки в лабораторном эксперименте.
Уточнение функционала
Л т,N1 Э 2и (* х ) (22)
Ьрхц = Ае(ЗТ) = -2 2 Э и( 2хк) ^,кку,ккх2,п
2 к=1,п=2 Эх
Неустранимая погрешность создается второй частью
Nx,Nt Э 3,
3 ,кку к. Она является аналогом случайной ошибки в
1 2 гп^3 Э3и(*п,хк)
к=1,п=2
2
аккх, п
Эх3
лабораторном эксперименте. Для нее можно получить верхнюю оценку ошибки функционала (после уточнения):
1
т, N
- 2 а!ккх,п
3
3 Э и(хп, Ук )
т, N
Эх3
2
к к3 Э и(хп > Ук ) Шп
ку,ккх,п Ши ,к
к=1,п=2 к=1,п=2
Верхняя оценка ошибки функционала
3 Э3и(хп, Ук )
Эх3
. Nx, N
Р? = 1 2
2 к=1, п=2
ку,ккх, п-
Эх3
п
Ши ,к
= АРх,1 (23) (24)
Соотношения (22) и (24) позволяют получить следующую оценку погрешности численного решения
Р - АрСх°1 - Рехас* < Д^хТ
Рассматриваемый подход позволяет определить устранимую (систематическую) погрешность расчета и верхнюю оценку случайной ошибки. В этом смысле он эквивалентен лабораторному эксперименту.
В качестве иллюстрации рассчитана погрешность плотности за скачками уплотнения (по линии симметрии течения на границе вытекания). Соответствующее поле плотности представлено на рис. 2. Поле сопряженной плотности изображено на рис. 3. Вторая производная плотности представлена на рис. 4, на равномерной сетке это поле соответствует полю максимальных локальных погрешностей аппроксимации в уравнении неразрывности. Пространственная плотность систематической
1 Э 2и (*п, хк)
погрешности
Эх2
Шип,кку,ккх п (22) представлена на рис. 5, а плотность
1
неустранимой погрешности —
к к3 Э и(хп> Ук) Шп
пу,кпх, п :~"3 Ши ,к
Эх3
(24) на рис. 6,7.
2
Рис. 2. Плотность
Рис. 3. Сопряженная плотность
Рис. 4. Вторая производная плотности
Рис. 5. Плотность устранимой погрешности (22)
Следует отметить, что стандартные методы повышения точности, основанные на сгущении сетки в областях больших градиентов, приведут к сгущению сетки в соответствии с рис. 2, в то время как зоны реальной погрешности расчета соответствуют рис. 6 и 7 (для разных точек оценки, к=50, к=70).
1 00
60
>
40
20 -
50 100 X 150 200
150 200
Рис. 6. Плотность неустранимой погрешности, оценивается плотность течения в точке к=50
Рис. 7. Плотность неустранимой погрешности, оценивается плотность течения в точке к=70
Таким образом, зоны максимального вклада ошибки дискретизации в погрешность искомой величины и, соответственно, зоны сгущения сетки существенно зависят не только от поля течения, но и от выбора оцениваемого параметра, поэтому не могут быть определены без решения сопряженной задачи.
Следует отметить, что при наличии значительного числа ценных функционалов ег, г = 1... Ме, и разных наборов параметров, содержащих ошибку /■, у = 1...М/ ,
возникает необходимость решения 2 • Ме- N / практически совпадающих задач. Это
1 00
80
60
40
20
5 и
1 и 0 150
X
200
X
обстоятельство делает данный подход к определению погрешности перспективным для реализации параллельных вычислений с точки зрения концепции многозадачного параллелизма.
3. Задачи валидации
Еще одним источником погрешности расчета является отклонение используемой физической модели от точной. Доказательство соответствия этих моделей является предметом валидации.
В том случае, когда отклонение сводится к наличию или отсутствию в системе уравнений некоторых членов, не приводящих к радикальной перестройке течения, использование сопряженных уравнений позволяет количественно определить вносимую ими погрешность. Ее малая величина может являться обоснованием верности применяемой физической модели и, соответственно, существенным элементом валидации. В качестве примера рассмотрим влияние вязких членов в уравнениях (12-15) на оцениваемую величину. Пусть возмущенные величины описываются параболизованными уравнениями Навье-Стокса
и
д и
дХ
+ V
1
д и 1 д Р
д Н р д Х Яе р д Н
и
=0
а невозмущенные -уравнениями Эйлера
и
д и
+ V
и
1 д Р
= 0
д Х д Н р д Х
Вариация целевого функционала (плотности в некоторой точке) принимает вид
8е = Г А р8 (х - хе ) + Г
и дА^ + V дА и
1
д2 и
дХ дН
Влияние вязкости на целевой функционал имеет вид
Яе р д Н
и
8е
= -±- Г
Ке О
1д 2 и щи + 4
д 2 V
р д У2
3р д Н2
Шv +
К
рРг д Н2 е 3р
' ди^2
Н
Х¥е
йХй
Н(25)
//
При этом соответствующие сопряженные уравнения в отличие от системы (1215) не содержат вязких членов. Пространственная плотность вклада вязких членов в погрешность расчета целевого фу| 1кционала представлена на рис. 8.
2
2
2
100
20
X
Рис. 8. Плотность погрешности, вызванной влиянием вязких членов Таким образом, сопряженные параметры дают возможность оценки влияния вариации физической модели на целевой функционал. Рис. 8 показывает те зоны, в которых влияние рассматриваемой вариации наиболее велико. В рамках рассматриваемой относительно простой задачи эти зоны интуитивно понятны и без данной визуализации - они соответствуют пересечению скачков уплотнения и
характеристик, приходящих в точку измерения. В более реалистичных задачах такого сорта интуитивные оценки затруднены или невозможны, поэтому поля сопряженных параметров могут служить весьма важной эвристической информацией. Визуальное представление этих полей позволяет определять в расчетной области зоны наибольшего влияния вариации физической модели.
Заключение
Представленный в настоящей работе подход, согласно классификации [13] может быть отнесен к методам анализа данных, так как основан на синтезе математического аппарата теории сопряженных операторов и уравнений и методов визуального представления. Следует отметить перспективность рассмотренного в данной работе подхода для практического решения актуальных задач верификации и валидации (V&V).
Визуализация поля градиента целевого функционала и связанных с ним полей сопряженных параметров позволяет не только качественно и количественно оценить вклад погрешности аппроксимации, но и выработать рекомендации по выбору расчетных сеток и необходимого числа итераций.
В задачах валидации применение данного подхода к изучению влияния членов уравнений, описывающих малые (отбрасываемые) эффекты на градиент целевого функционала также предоставляет количественную и качественную информацию, которая может быть визуально представлена.
Важным фактором является то, что задачи, решаемые с помощью подобного подхода, при наличии значительного числа ценных функционалов, и разных наборов управляющих параметров, либо параметров, содержащих ошибку, приводят к необходимости решения большого количества практически совпадающих однотипных задач с различными наборами данных. Это обстоятельство делает рассматриваемый в данной работе подход перспективным с точки зрения многозадачного параллелизма для применения параллельных вычислений.
Литература
[1] Марчук Г.И. Сопряженные уравнения и анализ сложных систем, М., Наука, 1992, 336 с.
[2] Марчук Г.И. Сопряженные уравнения и их применения / Труды Института математики и механики УрО РАН, 2006, т. 12, N 1, С.184-195.
[3] Guide for verification and validation of computational fluid dynamic simulations / AIAA Guide, AIAA G-077-1998.
[4] Roache P. J. Verification and Validation in Computational Science and Engineering. Hermosa Publishers, Albuquerque, New Mexico, 1998.
[5] Oberkampf W. L., Blottner F. G. Issues in Computational Fluid Dynamics Code Verification and Validation / AIAA Journal, v.36, N 5, 1998, pp.687-695.
[6] Алексеев А.К., Бондарев А.Е. Визуализация переноса погрешности при расчете поля течения // «Научная визуализация в прикладных задачах», Сб. науч. тр., Москва, МГУ, 2003, с.4-13.
[7] Venditti D., Darmofal D. Grid Adaptation for Functional Outputs: Application to Two-Dimensional Inviscid Flow / J. Comput. Phys., 176 (2002), pp.40-69.
[8] Алексеев А.К. Апостериорная оценка погрешности конечно-разностного решения с помощью сопряженных уравнений и дифференциального представления / ЖВМ и МФ, т. 45, N7, с. 1213-1225, 2005.
[9] Braack M., Ern A. A Posteriori Control of Modeling Errors and Discretization Errors / Multiscale Model. Simul., V. 1, N. 2, , 2003, pp. 221-238.
[10] Alekseev A., Navon M. Criteria of optimality for sensors' location based on adjoint transformation of observation data interpolation error / Int. J. Numer. Meth. Fluids, V. 62, Issue 1, 2010, pp. 74-89.
[11] Kirby R. M., Silva C. T. The Need for VerifiableVisualization // IEEE Computer Graphics and Applications, Volume 28, Number 5, 2008, pp.78-83.
[12] Шокин Ю.И., Яненко Н.Н. Метод дифференциального приближения, Наука, 1985.
[13] Бондарев А.Е., Галактионов В.А., Чечеткин В.М. Анализ развития концепций и методов визуального представления данных в задачах вычислительной физики / Журнал вычислительной математики и математической физики, 2011, Т. 51, N 4, С. 669-683.
