Брэббия К., Теллес Ж., Вроубел Л. - Методы граничных элементов (1050609), страница 29
Текст из файла (страница 29)
(4.75) Все поверхностные граничные элементы, описанные в гл. 3, могут применяться и в данном случае, включая несогласующиеся элементы. В частности, если требуется интегрировать по области, следует применять и элементы в форме трехмерных ячеек, рассматривавшиеся в равд.
3.5. 4.8. Осесимметричные задачи Полагая, что все граничные и внутренние значения являются осесимметричнымн, уравнение (4.34) можно записать в цилиндрической системе координат: вл 2Л сД) и Д, 1„)+ й ~ Г]и(х, 1) ~ до(4, х, 1„, 1)вЮ(х) )с(х)в[Г(х)е[[= еа г О ер 2л = й] ] в[(х, 1) ~и'Д, х, Ею 1)в[6(х)Я(х)вгГ(х)еИ+ ев г о 2л + ~ и,(х, 1,) ~ и'Д, х, 1р, 12)еЮ(х) 12(х)Ж(х), (4.76) о где à — часть контура образующей тела вращения, Й вЂ” часть ограниченной этим контуром плоской области, получающаяся при пересечении поверхности И и контура Г с полуплоскостью Я'Л (см. рнс. 2.37). Записывая трехмерное фундаментальное решение (4.72) в цилиндрической системе координат и интегрируя его по окружности радиуса ег (х), лежащей в плоскости Л = г, (х), получим й*Д, х, 1„, 1)-= ~ и'Д, х, 1„, 1)еЮ(х) = о 2л о (4.77) где 2 = Яо Щ + )со (х) + [Я ($) — Е(х)]2.
Это выражение для осесимметричного фундаментального решения можно представить в виде [30] 1 — )2 (гь) о (х), г, — модифициРованнаЯ фУнкциЯ БесселЯ первого рода нулевого порядка. Производную фундаментального решения по нормали к граничному контуру можно найти, продифференцировав выражение (4.78), что дает 8 на(/ )вег ехр ( — м, ) [[ «х) ( — ', 7 )~ Д) 12 ( ~у, )1 пн(х) — [3Д) ~(х)] 12( ~~, ) пг (х)~, (4.79) 1вт 186 Глава 4 Задачи теории теллолрооодноеош где /, — модифицированная функция Бесселя первого рода н первого порядка. Йз приведенных выражений можно видеть, что при /с ($) — ч- 0 имеем 1 — О, /О (1/2йт) — 1, !с (1/2йт) — 0 и, следовательно, кольцевой источник переходит в точечный источник с интенсивностью 2п на оси вращения.
Подставляя выражения (4.78) н (4.79) в уравнение (4.76), получим граничное интегральное уравнение для зависящих от времени функций: ср с(е)и($, /о)+й~ ~и(х, 1)с]О(й, х, 1;, 1)д(х)с[Г(х) И = с,р ср = ] ] у(х 1)й*($, х, 1р, 1)]с(х)с(Г(х)с(1+ сч г + ~ ио(х, 1а) й~(й, х, /а, 1~) /4 (х) сй(х). (4.80 ) Для решения уравнения (4.80) можно использовать численные процедуры, которые обсуждались ранее. Для простоты далее будут использоваться функции и и с/ кусочно-постоянного типа.
После введения дискретного представления поверхности и внутренней части рассматриваемой области с помощью соответственно граничных элементов и ячеек получается аналогичное (4.39) уравнение, матрицы которого определяются выражением (4.41). Для того чтобы выполнить интегрирование по времени аналитически, требуется перейти к новым переменным. Вводя обозначения с = 1/в, у = 2/(4/ст), а = 2/(4й Л/с), (4 .81) для интеграла от функции ио получим сс о/ о у. 1 й* сй' =, ] !Д2су) у — 'сое — о йу.
(4.82) Функцию Бесселя /о можно разложить в ряд вида [10 ] (4,83) л=о и после подстановки этого выражения в формулу (4.82) и почленного интегрирования найдем Сч 1 е2" Р с[1= ~~, ~~)' ' ~Г(2п+ ~~, с,) л=о — Г(2п+ —,', ас)]. (4.84) Для интеграла от функции с]о имеем с]асИ = — [/4(х)па(х) — [ЕЯ) — Е(х)]пх(х)] Х )!/2 с с-с а/ ос х ~ /О(2су)у'сое оау — /с($)пв(х) ~ /,(2су)у'с'е-вйу . (4.85) ас ас с Представляя функцию Бесселя /, в виде [10] %'~ (Оу)ол+' /,(2су) =- т -о (4.86) получим сс ,1 [[/[ (Х) П„(Х) — [Е (й) — Е (Х)] Пх (Х)] лС си (ив)'с2 ~ ~С-2 х Ъ~ъ о2 ~Г(' с-[- 3, ас ) — Г(2п+ —, ас/1— л О 2л+ С вЂ” й($) па(х) ~ (л1), (а+1) ~Г(2п+ 2 1 ас 2/ л=о Г (2п+ —, асч~~ . (4,87) Все неполные гамма-функции, фигурирующие в приведенных выше выражениях, можно выразить через функции Г (1/2, а) с помощью рекуррентной формулы [18] Г(п+ 1 а) пГ(п а) + але-а (4.88) В соответствии с введенными обозначениями (4.81) величина с изменяется от 0 (при /с ($) = О, /с (х) = 0 или в — ~ оо) до 0,5(при $ = х).
Ряды в выражениях (4.84) и (4.87) сходятся очень быстро для малых значений с и медленно при с- 0,5. При с = 0,5 эти ряды расходятся из-за наличия сингулярности для $ = х. Таким образом, видим, что с вычислительной точки зрения не следует использовать ряды (4.83) и (4,86) для значений с, лежащих в окрестности с = 0,5. Чтобы преодолеть эту трудность, можно воспользоваться аснмптотическими представлениями функций Бесселя, которые будут !вв Глава 4 Зидаоа теории твллолроводлооит Оу 94 ~ (а) = (2п — 1)'(2а — 3)'...
1, Ке (п) = ( — 1) [4 — (2п — 1)2] [4 — (2а — 3)']... [4 — 1]. Тогда интегралы по времени принимают вид у! 1'-,, ° 1 йе ([1 =,, Е2 (Вг„т) — Е2 (Вг) + (4.89) (4,90) (4.91) 49 .[.~ ' „[Г( —, 94 ) — Г( — а В))), (4.92) л! (!6о) 2! [ [*В= — ' „[- 4 Ь-'4- — Е-'Ч((В(2)-В())..(*Н. 2! 2 + [Е($) — 2(х)] пх (х)]+ [Я (х) а„(х) — [Е($) — Е(х)] ах(х)] х 94 х р — ' — „[Г(1 — п, В! 2) — Г(1 — п, Вг)] — Я(5)па(х) х с~ 12 (л) ви Л !(в)" х ~ ' " „]Г(1 — п, В! 2) — Г(1 — п, В )], (4.93) й) 72(л)аи л1(!в )" где Ь = 1 — 2с, В = аЬ. Неполные гамма-функции можно теперь выразить через Г (О, В) с помощью рекуррентных формул [10]: 1 Г Г( — п, В) = — — ~Г(1 — п, В) — — „~, Г(0, В) = Ет(В).
(4.94) Когда с принимает значения, близкие к 0,5, а значения у малы на части интервала интегрирования (ау „ау), представления (4.89) справедливы для больших значений их аргументов. Так, для боль- ших значений у можно записать [10]: и (4.90) нельзя применять непосредственно.
С другой стороны, выражение (4.82) можно записать в виде В' ~ йе([1 = ~ 29(2су)у-)/ое — о([у+ 2й (иу)[а Е-2 В! .). [ [,(2 у)у- У вЂ” 49), (4.99) В' ко[да значения у достаточно велики на интервале (а', ау). Теперь первый интеграл в правой части этого равенства можно вычислить по формуле (4.84), а для вы- 4 и числения второго интеграла использовать разложение в СГ ряд (4.89). Аналогичный под- '~~, и ход можно использовать и при вычислении интеграла сФ по времени от функции (]*, Заключительным шагом численного решения граничного интегрального уравнеl. с ния (4.80) является вычисление пространственных интегралов. Коэффициенты НЫ и Рис.
4 !4. Геометрия иииейиого Влемеита. Ви (2' Ф /) окончательной системы уравнений (аналогичной системе (4.40)) можно вычислить, используя шеститочечную квадратурную формулу Гаусса. Однако с диагональными коэффициентами Нп н (уп следует обращаться более осторожно, поскольку при их нахождении требуется вычислять сингулярные интегралы. Коэффициенты 0„представляют собой интегралы с логарифмической особенностью. Разложением интегральной показательной функции в формуле (4.92) можно изолировать слагаемое с логарифмической особенностью и проинтегрировать его аналитически [26].
Остальные слагаемые являются несингулярными, поэтому в расчетах можно использовагь стандартную квадратурную формулу Гаусса. Коэффициенты Н„содержат сумму логарифмической и вида 1/Ь особенностей. Первая из них является интегрируемой, а вторая интегрируется только в смысле главного значения. Тем не менее для случая постоянных или линейных элементов (рис.
4.14) можно записать па(х) =сова, их(х)=в1па, 42($) — !4(х) = — 21 2 в1па, Я $) — Е (х) = 21 — сова. (4.96) Задачи теории 1пепаапроеадиасти Глава 4 Таким образом, первое слагаемое в правой чаоти выражения (4.93)„ которое имеет особенность вида 1Ъ, будет тождественно равно нулю. Выделяя первые слагаемые каждого ряда в выражении (4.93), содержащие логарифмические особенности, можно вычислить их аналитически 126), а остальные расчеты проводить с помощью стандартной квадратурной формулы Гаусса. Свободные коэффициенты с, обусловлены появлением скачка в интеграле от функции 4ч в точках, принадлежащих границе Г, В рассматриваемом случае можно показать, что значения с ($) совпадают с найденными для двумерных задач (см.
равенство (4.65)). Пример 4.6. Сначала рассмотрим сплошной цилиндр с единичными начальными условиями и граничными условиями вида О,В О,В 0,4 Рнс. 4.!6. Зависимость и от раднуса на поверхностях Я вЂ” ~1. Сплошная кривая соответствует аналнтнческому решению, точки — решению по методу граянчных элементов.
0,2 и= О при ах =а, д = 2и при Я = ~1. ХОО 0,50 я 0,21 ДО Схема дискретного представления показана на рнс, 4.14. Отметим, что благодаря симметрии относительно оси )х необходимо рассматривать лишь половину поперечного сечения. Численные исследования проводились для поперечного сечения размерами а = 1, 1 = 1; для простоты коэффициент й, характеризующий физические свойства материала, также полагался равным единице.
Полученные здесь решения сопоставляются с аналитическим решением [14) на рис, 4.15 и 4.16, где видно хорошее их согласование. При расчетах использовался шаг по времени Ы = 0,025. Пример 4.7. В этом примере исследуется задача теплопроводностн для эллипсоида, температура которого в начальный момент времени равнялась нулю, а температура поверхности равнялась единице. Параметрические уравнения точек эллипсоида, лежащих в плоскости ЙЯ, имеют вид 0,0 Од Рнс.