Голдаев С.В., Ляликов Б.А. Основы математического моделирования в теплотехнике

Подождите немного. Документ загружается.

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

какая задача близка к рассматриваемой. Поэтому целесообразнее использо-

вать безразмерные переменные. Кроме того, на основании «Пи – теоремы»

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

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

рая теорема подобия) [20], [22], [28-31].

В последующих разделах пособия на конкретных примерах будет неод-

нократно применяться описанный подход.

Большинство мето

дов являются итерационными. Если последователь-

ность приближений некоторым образом стремится к пределу, то процесс на-

зывается сходящимся, а если нет, – то расходящимся. Сходимость метода

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

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

ния, зависит от начального приближения и требуемой точности. Чем ближ

начальное приближение к истинному решению, тем быстрее оно будет дос-

тигнуто.

Прежде чем выбрать тот или иной метод, необходимо проанализировать

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

цию, или систему уравнений аналитическому исследованию, в результате ко-

торого выявляется возможность использования данного метода. При этом

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

ветствующим образом преобразована, с тем, что

бы можно было эффективно

применить численный метод.

После выбора метода решения желательно составить четкое описание по-

следовательности вычислительных и логических действий, обеспечивающих

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

требованиями к форме и содержанию записи алгоритма являют

ся его нагляд-

ность, компактность и выразительность. В практике математического обес-

печения широкое распространение получил графический способ описания

алгоритма. Этот способ основан на представлении отдельных элементов ал-

горитма графическими символами, а всего алгоритма – в виде блок-схемы.

Внутри графических символов записываются словесно или посредством сим-

волов произвольные действия. Степень детализации алгоритма в виде блок –

схемы за

висит от его сложности, от степени использования стандартных ал-

горитмов [21], [24-26].

1.4. Отладка программы

Все ошибки, возникающие в процессе подготовки и решения задачи,

можно свести к следующим: ошибки в алгоритме; ошибки программирова-

ния; ошибки при подготовке информации; ошибки ЭВМ.

Ошибки в алгоритме могут быть как следствием неточности математиче-

ского описания зад

ачи, так и следствием неправильно выбранного численно-

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

тельность или расходимость итерационного вычислительного процесса, ко-

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

гда решение либо принципиально не может быть получено данным методом,

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

выявляются либо после получения результата, либо непосредственно в про-

цессе счета, когда ЭВМ не выходит на окончание расчетов вследствие «за-

цикливания» в некоторых частях программы. Для выявления такого рода

ошибок осуществляется контроль по физическому смыслу получаемых ре-

зультатов, либо производится пересчет с мень

шим шагом [13], [21], [24].

Ошибки программирования обусловлены тем, что алгоритм решения не-

правильно записан на соответствующем языке. Нарушения синтаксиса языка

выявляются на этапе трансляции, и информация о них выдается в виде диаг-

ностических сообщений [32]. Определенные трудности возникают в больших

программах при допущении механических ошибок, т.е. пропуске символов,

замене операций, неправильной расст

ановке скобок и т.п. Такие ошибки не

воспринимаются как синтаксические, но они приводят к неверному результа-

ту. Чтобы этого избежать, следует сокращать длину сложных выражений,

используя процедуру разбиения их на ряд более простых; в строке записы-

вать один оператор, чаще применять ко

мментарии [24], [25], [26], [32].

1.5. Установление адекватности модели

Математическая модель объекта является лишь определенным в рамках

принятых допущений аналогом. Поэтому значения переменных, получаемых

на модели и в объекте, различаются. Этап установления адекватности (соот-

ветствия) модели объекту является заключительным в последовательности

этапов, выполняемых при ее разработке. Результаты, получаемые на этом

этапе, могут привести к существенным изменениям решения задачи, начиная

с пост

ановки. [13], [21-24].

Прежде чем преступить к проверке и установлению адекватности, необ-

ходимо выработать критерий, который позволил бы сделать заключение о

соответствии модели и объекта. При моделировании различают качественное

и количественное соответствие. Можно найти средневзвешенные отклонения

между расчетными и экспериментальными значениями. Полученные число-

вые значения не отражают степени со

ответствия модели и процесса, а позво-

ляют сделать лишь заключение о характере поведения модели и процесса, ее

качественном соответствии.

Критерий адекватности является функцией некоторой совокупности па-

раметров, входящих в уравнения модели. Изменяя значения этих параметров

и находя минимум критерия адекватности, делают заключение о количест-

венном соответствии модели и процесса. После этого появляется возмож-

ность экстраполяции модели за пределы условий определения параметров.

Одним из наиболее распространенных критериев адекв

атности является

сумма квадратов отклонений экспериментальных и полученных на модели

значений переменных.

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

1.6. Применение математической модели

Математическая модель после установления адекватности поступает в

эксплуатацию как прикладная программа и может использоваться в двух ва-

риантах: автономно или в совокупности с другими программами [5], [11-13].

Прикладная программа обычно оформляется в соответствии с определенны-

ми требованиями так, чтобы обеспечить возможность ее хранения и авто-

номного или системного использования; такой программный продукт назы-

вают мод

улем. Для указания характеристик каждый модуль должен сопро-

вождаться своеобразным паспортом, в котором может содержаться следую-

щая информация: 1) описание задачи, математическая формулировка с пе-

речнем принятых допущений и описание алгоритма решения; 2) название

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

входных и выходн

ых параметров; 5) описание схем реализации для многоце-

левых модулей, с указанием входа и выхода для каждой схемы; 6) указание

оператора ввода-вывода с определением вводимых и выводимых перемен-

ных; 7) указание характеристик по быстродействию, объему занимаемой па-

мяти; 8) указание ресурсов ЭВМ для выполнения модуля; 9) описание ис-

ключительных ситуаций и рекомендации по их преодолению; 10) список

други

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

ние контрольного (тестового) примера, исходных данных и результатов рас-

чета.

Для удобства обработки результатов используются подпрограммы по-

строения графиков, таблиц и т.д.

2. АППРОКСИМАЦИЯ ФУНКЦИЙ

2.1. Интерполирование функций

В инженерной практике получил распространение табличный способ за-

дания функций, при котором для конечного множеств

а значений аргумента

xxxx ,...,,

= известны, полученные экспериментально соответствующие

значения функции

yyyxfy ,...,,:)(

= . Они являются результатами наблю-

дений или измерений. Аналитическое выражение функции в таком случае

неизвестно, что не позволяет вычислить ее значения в промежуточных точ-

ках ),...,1,0(, nixx

== . Для отыскания этих значений строят аппроксими-

рующую (приближающую) функцию ϕ(x), расчеты по которой либо совпа-

дают, либо в определенном смысле приближаются к экспериментально по-

лученным величинам. Построение функции называется интерполированием

[20], [22], [23].

Интерполирование применяют и в случае, когда аналитический вид

функции ϕ(х) известен, но сложен и требует большого объема вычислений

для определения отдельных з

начений функции. Над такой функцией, кроме

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

того, трудно выполнять математические операции (например, дифференци-

рование или интегрирование). Замена f(х) приближенной функцией ϕ(х) по-

зволяет упростить вычисления. Для этого по известному выражению f(х) вы-

числяют небольшую таблицу ее значений и по ним строят простую интерпо-

лирующую функцию ϕ(х). Последнюю используют для нахождения прибли-

женных значений f(х) в других точках.

Широкое распрост

ранение получила алгебраическая интерполяция, при

которой в качестве интерполирующих функций ϕ(х) применяют алгебраиче-

ские многочлены. Это объясняется тем, что полиномы вычислять несложно,

кроме того, достаточно просто выполнять над ними математические опера-

ции. Алгебраическое интерполирование функции y = f(х), заданной своими

значениями

)(),...,(),(

10 n

xfxfxf

в точках

xxx ,...,,

отрезка [a, b], состоит в

приближенной замене этой функции на данном отрезке многочленом P

(x)

степени n, который в точках

x , принимает те же значения, что и функция

f(х): )()(

iin

xfxP = , i = 0, 1, ..., n. Точки

x (i = 0, 1, ..., n) называют узлами

интерполяции.

Например, при линейной интерполяции заданные точки (x

, y

) соединя-

ются прямолинейными отрезками, и функция f(x) приближается ломанной с

вершинами в данных точках. Уравнения каждого отрезка ломанной в общем

случае разные. Поскольку имеется n интервалов (x

i-1

, x

), то для каждого из

них в качестве уравнения интерполяционного многочлена используется

уравнение прямой, проходящей через две точки (x

i-1

, y

i-1

) и (x

, y

), в виде

−

Отсюда

()()

.,/

1111

−−−−

−

−=−−=

≤

iiiiiiiii

iiii

xaybxxyya

xxxbxay

(2.1)

Таким образом, при использовании линейной интерполяции сначала оп-

ределяется интервал, в который попадает значение аргумента x, а затем он

подставляется в формулу (2.1), и находится приближенное значение функции

в этой точке.

Пример 2.1. Водяной пар при давлении p

= 1,8 Мпа и температуре

=523 К дросселируется до давления p

= 1,0 МПа. Определить температуру

Решение. Основываясь на том, что процесс дросселирования происходит

при постоянной энтальпии [34], имеем h

= h

. В исходном состоянии пар на-

ходится в перегретом состоянии, и из таблиц [35] выписываем: h

= 2911,1

кДж/кг. Данная величина энтальпии при давлении p

= 1,0 МПа находится

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

между следующими величинами энтальпий, которым соответствуют значе-

ния температур: 15,513,5,2920;15,503,9,2897

)2(

)1(

==== ThTh

Для определения искомой температуры T

воспользуемся формулой (2.1).

В рассматриваемом примере

() ()

()

() ()

()

() ()

.1,2378;26,2/

;

=−==−−=

aThbTThha

baTh

Поэтому

()

.95,508/

KabhT

−=

При построении единого для всего отрезка [x

] интерполяционного

многочлена, его график должен проходить через все заданные точки. Одна из

таких кривых строится с помощью интерполяционного многочлена Лагран-

жа, который имеет вид

(

)

(

)

(

)

()()()

()

−−−

−

02010

...

)( xf

xxxxxx

()()

(

)

()()()

()

−−−

−

12101

...

xxxxxx

()()

(

)

()()( )

()

...

10100

110

xxxxxx

−

−−−

−

(2.2)

Из формулы (2.2) можно получить более простые варианты интерполяций

(линейную или квадратичную).

Данный алгоритм реализован в виде программы для ЭВМ (П1.2). Обра-

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

Пример 2.2. При определении к.п.д. цикла Ренкина необходимо найти эн-

тальпию отработавшего влажного пара (рис.2.1) с параметрами (P

= 0,04 бар;

x = 0,718 ; t′

= 28,4 0C) на выходе из турбины по известной формуле

[33]:

= h′

+ x r,

(2.3)

где h′

= f (t′

), r = f(t′

) соответственно энтальпия воды и теплота паро-

образования в состоянии насыщения на выходе из конденсатора.

Для определения h2x необходимо проинтерполировать h′

и r, используя

справочные данные [34].

Для решения задачи в память ЭВМ из входного файла вводятся таблич-

ные массивы: h′, t′

r [34], выполняется обращение к процедурам интерполя-

ции (П.1.2. или П.1.3.), производится вычисление по формуле (2.3) и вывод

результата в выходной файл. Фрагмент программы, файлы с исходными дан-

ными и с результатами к примеру 2.2. представлены ниже.

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

Рис. 2.1. Цикл Ренкина в TS - диаграмме

program RTP;

uses Uchismet;

var h2x, h3s1, x, r, t2s1 : real; tsm, hs1m, rm : mac22; n : integer;

begin

assign(finp,’c:\tp55\st\G6A54\dint1.dat’);

reset(finp);

assign(fout,’c:\tp55\st\G6A54\rint1.dat’);

rewrite(fout);

readln(finp,n); {Ввод количества точек в базовой таблице интерполя-

ции}

readln(finp,t2s1,x); {Ввод значений t2s1,x }

for i:=1 to 17 do read(finp,tsm[i]); {Ввод массива температур насыщения }

for i:=1 to 17 do read(finp,rm[i]); {Ввод массива теплоты парообразования }

for i:=1 to 17 do read(finp,hs1m[i]); {Ввод массива энтальпии насыще-

ния }

inter(n,tsm,hs1m,t2s1,h3s1); {Обращение к процедуре интерполяции

методом Лагранжа}

linter(n,tsm,rm,t2s1,r); {Обращение к процедуре линейной интерполя-

ции}

H2x:=H3s1+X*r; {

Вычисление энтальпии вла

жного пара}

writeln(fout, {Вывод результатов расчетов}

′Для влажного пара с параметрами t′

= 28,6 бар; x = 0,718 - энталь-

пия = ′,

h2x:7:1, ′ кДж/кг′);

close(finp);

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

close(fout);

end.

Файл с исходными данными

28.6 0.718

20 22 24 30 35 40 45 50 55 60 100 150 200 250 300 350 374

2455.3 2452.8 2442.8 2431.9 2420.0 2406.0 2395.8 2384.1 2371.0 2359.0

2258.41 2115.95 1942.065 1716.22 1405.33 893.73 139.46

84 92 101 125.7 146.7 167.6 188.6 209.5 230.5 251.4 419 628.5 838

1047 1257 1466.5 2033.4

*** Файл с результатами расчета ***

Результаты интерполирования :

для t2s1 = 28.6 C - H3s1 = 119.9 кДж/кг;

для t2s1 = 28.6 C - r = 2433.5 кДж/кг.

Энтальпия:

для параметров влажного пара t2s1 = 28.6 C, x = 0.718 : h2x = 1867.2 кДж/кг

2.2. Подбор эмпирических формул

В случае обработки опытных данных, полученных в результате наблюде-

ний или измерений, нужно иметь в виду ошибки этих данных. Они могут

быть вызваны несовершенством измерительного прибора, субъективными

причинами, случайными различными факторами и т.д. Ошибки эксперимен-

тальных данных можно условно разбить на три категории по их происхожде-

нию и величине: систематические, случ

айные и грубые [20].

Систематические ошибки обычно дают отклонение в одну сторону от ис-

тинного значения измеряемой величины. Они могут быть постоянными или

закономерно изменяться при повторении опыта, а их причина и характер из-

вестны. Систематические ошибки могут быть вызваны условиями экспери-

мента (влажностью, температурой среды и др.), дефектом измерительного

прибора, его плохой регулировкой и т.д. Эти ошибки можно устранить на-

ладкой аппаратуры или введ

ением соответствующих поправок.

Случайные ошибки определяются большим числом факторов, которые не

могут быть устранены либо достаточно точно учтены при измерениях или

при обработке результатов. Они имеют случайный (несcимметричный) ха-

рактер, дают отклонения от средней величины в ту и другую сторону при по-

вторении измерений и не могут быть устранен

ы в эксперименте, как бы тща-

тельно он не проводился. Статистическая обработка экспериментальных

данных позволяет определить величину случайной

ошибки и довести ее до

некоторого приемлемого значения повторением измерений достаточное чис-

ло раз.

Грубые ошибки явно искажают результат измерения; они чрезвычайно

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

падают при повторении опыта. Измерения с такими ошибками отбрасывают-

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

ся и в расчет при окончательной обработке результатов измерений не прини-

маются.

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

ошибки. Они, вообще говоря, могут быть уменьшены до сколь угодно малой

величины путем многократного повторения опыта. Однако это не всегда це-

лесообразно, поскольку могут потребоваться большие материальные или

временные ресурсы. Значительно дешевле и быстрее мо

жно в ряде случаев

получить уточненные данные хорошей математической обработкой имею-

щихся результатов измерений.

Пусть при изучении неизвестной функциональной зависимости между

параметрами y и x в результате серии экспериментов произведен ряд измере-

ний этих величин и получена таблица значений (табл.2.1).

Таблица 2.1

… x

… y

Задача состоит в том, чтобы найти приближенную зависимость

(

)

,xfy

(2.4)

значения которой при x = x

(i = 0, 1, …, n) мало отличаются от опытных дан-

ных y

. Подобная приближенная функциональная зависимость (2.4), полу-

ченная на основании экспериментальных данных, называется эмпирической

формулой.

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

щего вида этой формулы и определения наилучших значений содержащихся

в ней параметров. Общий вид формулы иногда известен из физических сооб-

ражений или после применения теории подобия [28-30].

Если характер зависимости неизв

естен, то вид эмпирической формулы

может быть произвольным. Предпочтение обычно отдается наиболее про-

стым формулам, обладающим достаточной точностью. Они первоначально

выбираются из геометрических соображений: опытные данные наносятся на

график, и примерно угадывается общий вид зависимости путем сравнения

полученной кривой с графиками известных функций (многочлена, показа-

тельной или логарифмической функций и т.п.), которые имеются в справоч-

ной литературе [36

], [37].

Предположим, что тип эмпирической формулы выбран, и он имеет вид:

(

)

,...,,,,

10 n

aaaxgy =

(2.5)

где g – известная функция, a

, a

, …, a

– неизвестные параметры. Задача со-

стоит в том, чтобы определить такие значения этих параметров

, при которых

эмпирическая формула обеспечивает хорошее приближение данной функции,

значения которой в точках x

равны y

Разность между этими значениями (отклонения) обозначим через ε

. Тогда

()

....,,,,,

210 ini

yaaaaxg

−

(2.6)

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

Задача нахождения наилучших значений параметров a

, a

,…,a

заключа-

ется в некоторой минимизации отклонений ε

. Существует ряд способов ре-

шения такой задачи (методы выбранных точек, средних, наименьших квадра-

тов и др.). Первые два из перечисленных являются сравнительно простыми,

однако в ряде случаев получаемые с их помощью аппроксимации не облада-

ют достаточной точностью. Опишем процедуру использования метода наи-

меньших квадратов (МНК), с помощью которого, например, были получены

многие критериальные зави

симости для безразмерных коэффициентов теп-

лоотдачи при свободной и вынужденной конвекции газов и жидкостей [28-

31], интерполяционные формулы, обобщающие результаты вычислительного

эксперимента [14], [17], [35], [36].

Параметры a

, a

,…, a

эмпирической формулы (2.4) находятся из условия

минимума функции S, представляющей собой сумму квадратов отклонений

(2.7) для всех точек x

, x

,…, x

()

[]

....,,,,

∑∑

−==

ini

yaaaxgS

(2.7)

В качестве эмпирической функции на практике широко используется по-

лином степени m<n:

....),(

210

xaxaxaaax

++++= .

(2.8)

Для отыскания коэффициентов

aaa ,...,,

полинома используется необ-

ходимое условие минимума функций многих переменных [37]. Оно состоит

в равенстве нулю частных производных по каждому параметру

a :

()

∑

=−=

∂

=−=

∂

=−=

∂

iim

xyp

.02

............

,02

(2.9)

В результате получается так называемая нормальная система m+1 линей-

ных алгебраических уравнений, решение которой дает искомые m+1 значе-

ний коэффициентов многочлена

aaa ,...,,

⎪

⎩

⎪

⎨

⎧

=+++

,...

,.............................................

,....

,...

1100

11111010

00101000

mmmnmm

cababab

(2.10)

где

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

....,,1,0,,,

mlkyxcxb

ikl

===

∑∑

(2.11)

Как видно, порядок системы (2.10) равен степени аппроксимирующего

полинома (2.8), поэтому ее размерность составляет несколько единиц. Реше-

ние таких систем обычно получают прямыми (точными) методами, которые

просты и наиболее универсальны [18], [20], [22,23], [32].

Для решения системы, содержащей два-три уравнения, можно использо-

вать правило Крамера, выражающее неизвестные через отношения определи-

телей [20], [32], [36].

Поскольку вычисление определителей высокого порядка является трудо-

емкой процедурой, то в случае большего количества уравнений удобнее

применять метод Гаусса, называемый еще методом последовательного ис-

ключения неизвестных [18], [20], [22], [23], [32]. Согласно этому алгоритму,

вначале с помощью первого уравнения системы (2.10) исключается a

из

всех остальных ее уравнений. Затем с помощью второго уравнения исключа-

ется a

из третьего и всех последующих уравнений. Этот процесс называется

прямым ходом, и он продолжается до тех пор, пока в левой части последнего

m – го уравнения не останется лишь один член с неизвестным a

. Обратный

ход состоит в последовательном вычислении искомых неизвестных: из реше-

ния последнего уравнения находится a

. Далее, используя это значение, из

предыдущего уравнения вычисляется a

m-1

и т.д. до a

В качестве иллюстрации применения описанной выше методики подбора

эмпирических формул приведем решения нескольких задач, которые форму-

лируются следующим образом [38].

Пример 2.3. Зависимость мольной теплоемкости кислорода от темпера-

туры характеризуется следующими экспериментальными данными:

T, K

273 573 873 1273 1573 1973

кДж/(кмоль⋅К)

29,278 31,836 34,206 35,919 36,756 37,720

Воспользовавшись МНК, выразить зависимость μcp = f(T) в виде

210

TaTaac

++=

(2.12)

Проверить точность полученной аппроксимации путем сравнения с таб-

личными значениями [35].

В данном случае имеем m = 2, n = 5, а система уравнений (2.10) запишет-

ся следующим образом: