Построение ортогональной системы координат для выпуклого тела с вырезом произвольной формы Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ НАГИ

Т о м XX 1 989 №4

УДК 536.24

ПОСТРОЕНИЕ ОРТОГОНАЛЬНОЙ СИСТЕМЫ КООРДИНАТ ДЛЯ ВЫПУКЛОГО ТЕЛА С ВЫРЕЗОМ ПРОИЗВОЛЬНОЙ ФОРМЫ

И. В. Петухов, С. М. Решетько

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

При решении сопряженных задач конвективного теплообмена приходится решать собственную тепловую задачу для обтекаемого тела. Алгоритм решения такой задачи значительно проще в ортогональной системе координат. Кроме того, имеются сведения [1], что в неортогональной системе координат ухудшается сходимость итерационного процесса.

В общей постановке построение ортогональной системы координат связано с численным решением уравнения Лапласа [2]. Предлагаемый метод значительно проще, но накладывает более жесткие ограничения на форму тела. Однако эти ограничения не являются существенными при рассмотрении тепловой задачи для передней части плоских или осесимметричных тел.

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

1. Построение семейства координатных кривых т) = const. Рассмотрим криволинейную ортогональную систему координат rj в области плоского тела. Пусть кривые семейства £=const заданы уравнением

/(*, У. &) = 0, (1.1)

где f — известная аналитическая функция своих аргументов. Соотношения, необходимые для построения кривых семейства г) = const, можно получить из выражений

Л*ч+ЛУч = 0. Л*е+/у)!5 + /е=0 и из условия ортогональности кривых |=const, Г] = const

/У^-Л^ = о.

Эти соотношения имеют вид

VDG^fzl, D =fl + /у , (1.2)

D* Я6 - - Л /£ Я, л ,=/,2 fyy - 2/, fyfxy+ fy /,, , (1.3)

Dxt = —fxU яуе=-ЛЛ> (!-4>

/5жч=+/,Я, К5 у,--=+/, Я, (1.5)

где G = ']/.x|+jyj;, — параметры Ламэ. Параметр Я и

декартовы координаты *. у кривых tj = const определяются из решения дифференциальных уравнений (1.3), (1.4) при начальных значениях Я = Я0(т]), x = x0(t\), у = у0(?)), где Я0 == V+ у1ч, заданных на кривой 5 = £0. Знаки в (1.5) легко находятся в частной задаче.

Ниже уравнение (1.1) для кривой £=const будем задавать в прямоугольных координатах о, т сопровождающей системы

а — — х sin В у cos б, х == X cos 0 -f у sin 6, (1.6)

где 0=0(1) —угол между касательной к кривой 1 на рис. 1 и осью ,х в точке Р, через которую проходит кривая 1 = const, т. е. в виде

/■(о, £) = 0. (1.7)

Принимая f(x, у, %)—f (а, т, £), получим

Л — А + (—Л х + /с °) ве» | ^ щ

—/«i sin 0 +/С cos 0, fy = cos0 + fz sin 0 . |

Выражения для ВиА примут вид

D — fa 4- f? , А = flu -2f'f.U + • (1 -9)

Вместо (1.4) будем иметь уравнения

D{4 + x%0 = = (1.10)

при начальных условиях о = о0 (г(), -с = т0 (ttj) на кривой £ = £0- После решения уравнений (1.10) координаты х, у определяются из (1.6), а параметры О, Нм производные х%, у$, хп, уп, которые необходимы для решения физической задачи, —из (1.2) — (1.5) и (1.8), (1-9). Уравнение (1.3) решается при начальном условии Я0=У +

2. Ограничивающие кривые. Будем предполагать, что ограничивающие кривые / и 2 на рис. 1 имеют непрерывную касательную и ку-сочно-непрерывную неотрицательную кривизну, т. е. не являются вогнутыми. Предположим также, что существует прямая х=0 (необязательно ось симметрии), ортогональная к обеим кривым. Отрезок О этой прямой является локальным экстремумом толщины тела.

Уравнения ограничивающих кривых будем рассматривать в параметрической форме. Ниже параметры кривых 1 и 2 обозначены р и значения величин на этих кривых и на оси л: = 0 — нижними индексами 1, 2 и 0; производные по д, р и 0 — теми же буквами в роли нижних индексов; производные по | от функций, зависящих только от I, — штрихом. Если функция £(1) рассматривается в виде g{q, р, 0)

ИЛИ g*(g,p), где g* = g И gq = gq, то

g' = gq4'-\-gpP' + g^'-gpP' + gq<j' , g*p = gp + g9 V

Параметры р и 0 будем считать известными функциями от £ и р. В соответствии с (1.6)

°1 = °1 {р, 6) . 'Ч = *1 (Р, 6) , = °2 (Я> 6) . 4 = *2 (Я, 6) •

Приведем некоторые соотношения на ограничивающих кривых. Кривая 1. Уравнение кривой х1~х1 (р), У,=У! (р). Условия на

оси л: = 0: xio = xlpfi0 = O, yip0 = 0. Длина кривой: slp=Vxip-hyfP; xlp = sl pcos 0, yXp — Sj p sin 0. Угол 0 = ц между нормалью к кривой и осью — у: tg9 = у1р1х1р. Производные от о* = — хг sin 0 + cos 0 и tj = cos б -|- У1 sin 0:

Рис. 1

(2.1)

^1 р 0, Ol 0 --------1 Xj, Tj р _- Sj р, g : (

Ti Р = ■т! ррР'> ■'l рр ~ рр COS 0 Л' У\рр sin I

(2.2)

*>; ) е- >

Кривизна = dbjdSi кривой:

si р — °i рр> °i рр — ~ ^-1 рр s^n ® рр cos ®; | (2.3)

Кривая 2. Уравнение кривой: x2 = x2(q), УЪ=У2(Я)- Условия на оси л: = 0: x20=x2qqQ=O, у2д0=0. Длина s2 кривой: s2 q=Vxlq-\-ylq', x2q=s2q C0SX2. У2 q = S2 q sin /г- Угол x2 между нормалью к кривой и осью —у. tgX2 =^29/^2в- Угол Фг между нормалью к кривой и

осью о: *]>2 = 0— Х2. Производные от о2 = — х2 sin 6 + у2 cos б и х2 =

— х2 cos б -\-уъ sin 0:

о2(? = — JC29sin0 + Уг„ cos6 = — s2?sin^, cs2e = —^ 1 ^.4)

? = ? cos 0 + у2 ? sin 0 = s2 9 cos'!»2, 126=о,. I

Кривизна k2 = dii/dso кривой:

S2? = 02 qq cos ф2 + т2 qq sin ф2, S2? k2 = ^2q °2 qq ®2 q T2 > (2.5>

02 qq = - *2 ?<7 sin 0 + y2 qq cos 6, x2 qq = *2 и cos 0 + j/2 qq sin 6,

X2q = S2qk2, X2-ХидЯ'- (2-6>

Ниже будем считать, что

—те/2 < фг < тс/2. (2.7)

Из условия невогнутости кривых /, 2 и из (2.2), (2.4) следует

ki^-0, k2>0, т,р ]> 0, т2в^>0. (2.8)

Приведем примеры выбора параметра £ = р, q, когда кривой £=1,2 или ее передней частью является

а) эллипс + (,уг —Уи)2!^ = 1:

хг = a, sin t, yf = + bt cos t;

б) парабола kl0 x] = 2 (у, —yi0):

•*« = *, J;/=3'/o +y £/»

в) гипербола jcf = tg2 ^ [ (Уг—^iJ2 —(Уго—^а)2 ] c асимптотой *«tg <р, (У — УtaY-

Xi-= tg <Pt (y 10 - yio) sh .y, = .yte + (у, 0 — y/e) Ch *.

3. Семейство эллипсов — гипербол. Для координатных кривых семейства 1 = const примем

хз°г _ 2т°’ + Хт02, о° = о — б„ т°==т — т,. (3.1)

где х==х($), Х = Х(£). Уравнение (3.1) описывает эллипсы и гиперболы при разных и одинаковых знаках у х, X с вершинами на

0° = 02 О , х° = х2 х

кривой 1, ортогональные к кривой /; х— кривизна кривой в вершине.

Обозначим х„=^ ф = ^т/йа. Имеем ха2 = tg ф2,

(1 + ^г) Ха2=Х02, О2>0, 1+ЛХ2>0. (3.2)

Последнее неравенство следует из (2.7) и ограничивает |х°| для эллипсов. Условие ортогональности кривых (3.1) и '2 имеет вид [см.

(2.4)]

^2?+ О2? = 0 . (3.3)

Полагая в (1.7) 2/- = хо°2— 2х° — Хх02, получим для (1.8,) (1.9)

А=-у (*' с08 -)/ т°2) - /; х° в'-/-(х1р р' -0° ео, | ^

‘/а=хо°, /;=—(! +Хх°), Л = х . ]

Для построения кривых (3.1) при условии (3.3) достаточно задать один из параметров х, X. Основная трудность заключается в правильном выборе этого параметра.

Представим уравнение (3.1) в безразмерном виде

;°2 = 2х° + X ^°2, X = х — 2 ,

(3.5)

X —х2Х, о2х = ех, е — х2/о2 .

Обозначая ч = (йх°/(1<з0)2, будем иметь

^2 = ^, (х----- 1) (V- 1)=1.

Из (3.1), (3.2) следует, что знаки * и г,2 совпадают со знаком х2, т. е. *>0, V > 0. Последние неравенства возможны лишь при условии

х > 1 , у> 1 . (3.6)

Таким образом, семейство (3.5) в области 0<о°<1, 0<х°<;1 меняется от окружности О02 + (х°—I)2 = 1 при х -> 1 ДО прямой х° = о° при х со. Эти две предельные линии исключаются. Если в^О, то при переходе от плоскости 0о х° к плоскости О0 х° эллипсы Х<0 переходят в эллипсы, а гиперболы Х>0 — в гиперболы. При £ = 0 любая кривая (3.5) переходит в прямую х° = 0. В этом случае X = оо, поэтому параметр X вычислять не следует. Параметр х ограничен, если х ограничен. Ниже параметр х>1 будем считать заданным.

Приведем в новых обозначениях условие ортогональности (3.3) '•"^2 'С29 + °2О2?==0 , у(х—1)=х, (3.7)

уравнение для о из (1.10) и выражения (3.4)

;(а°-^) в' + в°*(/5/р—/б2/0г)=-о (3.8)

I +кх° = У \ +\*а°\ (3.9)

о=/?+/?, £=ГеИ°, /; = -(\+и°),

Л= у Ъ (в02 - *°2) *' - (х° /; ^ , + /, а2, д') -

— (1 — х°) А^рР' + (°° — <с°) Кл + вг/т) в' , (ЗЛО)

/51 = ^1РР',

/£2= [ (1 +х) т29~в;о,9] ^'=(х - 1 + е2 ух) х2? , (3.11)

а также уравнение (1.3)

02Яе = -х/5Я. (3.12)

Из (3.2) следует Б>0.

Уравнения (3.8), (3.12) решаются при начальных условиях

— —в о — ° "1 / і . g и и

O^oJy)), Н0 = 020 вот) у 1 + 1 + Ъъ4 •

Если линия ? = $о описывается прямой, то H0 — al0aln. За координату г) удобно принять оо = т), 0<к)<1. Координаты, х, у вычисляются из (1.6), где

о = ot + 02 0° , X = Tj -f х2 Х° .

В выражениях (1.5) следует брать нижние знаки.

4. Процедура решения. Будем считать, что р(£)—известная функция, причем р'>0. В частности, за координату £ можно принять \=р—Ро или ? = при s,0 = 0. Предположим также, что х =

— *{q, р, 0) = х* (q, р) — заданная функция своих аргументов.

Для построения кривых S = const необходимо определить о2, х2, т. е. q (£), а для кривых т] = const — еще и q'{%). Величины q и q' определяются из условия (3.7). Обозначим

Ф = ^2'2г + О0°2(7 •

Дифференцируя (3.7) по S, получим с учетом (2.1)

Ф?<7' + ф;/>' = 0,

®9 = V(X2X2?? + X29) + 0202??+ 029 + V?X2X2(J,

Фр= - (vx,p— V^Xs) Xt4+ [(у— 1) (02X2? + t202J + V9X2X2J Єр

где (х — 1) V,= — (v — l) x, ; t = q, p, 6.

Здесь использованы выражения (см. (2.2), (2.4))

а°2 —a2qq'— xlW, х2' = х2<? q' — xlppr + о°2 0' .

(4.1)

Решение уравнения (3.7) можно свести к нахождению корня уравнения х1 = Т(д), где Т (д) определяется из соотношения

v (т2 — Т) хїг + о2 о2<? = 0 .

Здесь р и 6 фиксированы. Дифференцируя по д и подставляя в выражение для Тд значение 7’=х1, получим

где дп — значение д из предыдущего приближения. За д0 можно взять д, вычисленное для предыдущей кривой I = const. Производная д’ определяется из (4.1).

Функции q, р и в являются непрерывными. Производная q' претерпевает разрыв на кривой сопряжения Г, т. е. на кривой |=const, которая проходит через точку разрыва кривизны одной из кривых 1, 2. На кривой сопряжения значения q' = q'~ и q/+ определяются до и после разрыва. Разрыв q' влечет за собой разрыв /£, Н^, Xt, у^ и G. Кривая Г должна принадлежать к кривым координатной сетки.

Пусть кривая 2 в малой окрестности оси х = 0 имеет форму клина с вершиной г/20 (кривая 2' на рис. 1). Тогда часть кривых семейства

1 = const, начиная с прямой 0, исходят из одной точки у2о- Для этих кривых <7 = <7о, q' = 0 и уравнения (3.7), (4.1) можно игнорировать. Последней из них является кривая, ортогональная к 2' в вершине у20. Эта кривая также является кривой сопряжения, поскольку q'~=0, а 0. Соответствующее значение | и величину q'+ можно определить из (3.7), (4,1) при q = q0.

Кривые семейства т] = const строятся на основе решения уравнений (3.8). При численном решении этих уравнений, а также уравнений (3.12) методом Рунге — Кутта следует учитывать разрыв производных по £ на Г. Отметим также, что (3.9) не годится для вычисления т° при малых К или А,=0. Для вычисления корня т° уравнения (3.5) при заданном о° лучше всего использовать метод касательных

Здесь за нулевое приближение ^=1 — ■/1 — о°2 берется решение для окружности X = — 1. Нетрудно показать, что хп+1^хп, о°<о, т. е. процедура корректна при любом Х>—1.

5. Ограничительные условия для координатного параметра %. Обозначим: (х—I)2Ф„ = Ф?х2?, (х — I)2Ф£ = — Фрт1рт2(?. Используя в (4.1) соотношение

ф? = фч — ^ Vх2?» ф<? = (х — 1 + Iх)х + [(х — I)2 + Н К . (5-2)

^Х2 °2 °2qq °2 ^2 ^2qlX2q >

которое следует из (2.5), (3.7), получим с учетом (2.3)

(5.1)

фр = фр' фр = (х — 1)(х- ло + ^1, (5.3)

где ^,=0°^!, А2 = 02А252?Л*г. }1 = Х°Х=:е°2Ч.

В точках пересечения кривых 1=сопэ1 параметр Ламэ й из (1.2), а значит, и /% обращаются в нуль. Из (3.11), (3.6) и (2.8) следует, что при р'>0 условиями непересекаемости являются

/Е>о, д'>о.

Для определения д' необходимо: Ф?=^0. Из (5.2) видно, что Ф? = Ф9>0 на прямой ,1о = 0 или при х = ч(р). Учитывая (2.8) и

(5.1), примем Ф9>0, а значит, и Фр>0,

ФЯ — т2х<?/т29>0. (5.4)

Ф7+^р1^р>0. (5.5)

Случай я = Яо, ^' = 0 можно рассматривать как предельный при О, &2-»со. Используя в (5.1) и (5.2) выражение °°2х'2 —

= ^2 х2д Я'> которое следует из (2.6), получим

[(х — 1)2 + Н02Х2=Фр'Ч/,/> (5-6)

где /2 — угол наклона кривых в точке у2о- Уравнение (5.6) выражает соотношение X2==®, — Фг- Его можно использовать для определения р', если принять $ = Хг- Из (5.6) также следует условие

(5.5), так как Х2>0.

Рассмотрим теперь условие, при котором /е>0. Из (5.1), (5.3),

(2.1) имеем

Х2Х' = Х2« Я' Фр Т1 рР'•

Подставляя это выражение в (3.10), получим с учетом (3.7) и (3.5)

и = х° К1 ?) ? ^яЯг + (1 - х°)2 Ер ^рР',

X = X2 4- {XX (х + 2м) + [(х — I)2 + [Ах] Еъ , Рд2 = X — 1 + (XV,

х2 Ер = X2 — х&, 5° + -4-(х — 1 + ц) 5° 2,

(5.7)

где 5° = ■ = X ----5° —0, 5°=1, 0<5°< 1.

о° + т° 1— т° 1 2

При <7 = <70, д'— 0 выражение (5.7) принимает вид

/£=т°(1 -х°) [(х_ 1)2/; + йо;х; + (1 -*уЕ*р^рР>.

Напомним, что ц = т°х>0. Так как /^,>0, /7?2>0, Хг>0, то /5>0 всюду при Е*р> 0 и >0; /? = 0 только в точке _у20 при д' = 0.

3—«Ученые записки» № 4 ■ 33

Неравенство Е*р > 0 выполняется при условии

2 xs — (х-f-1—при )

х — 1 + (г — &,/2>0 при jj-> 1. J

(5.8)

Условие (5.8) можно ослабить, заменив более простым

Х> V

(5.9)

Условие (5.8) или (5.9) является достаточным для непересекаемости кривых £ = сопз^ Оба условия являются ограничительными лишь при &!>1. Если выполнено (5.9), то условие (5.5) удовлетворя: ется автоматически на прямой т° = 0 или при х = х(<?), а оба условия

(5.4), (5.5)—при х = сопб1.

Из (3.12) имеем для параметра Ламэ Я

т. е. кривые семейства т} = const тоже нигде не пересекаются (НфО).

Можно показать, что условиям (5.4), (5.5) удовлетворяет следующее выражение для координатного параметра х:

если ш и а подобраны на основе (3.6), (5.9). Если боковая граница О или 3 на рис. 1 является прямой, то она не накладывает ограничений на х. Если же кривая 3 задана и т°3^0, то ш и а связаны соотношением (о + ао°3 = х3. При рассмотрении тепловой задачи для тела при наличии оси симметрии лг = 0 форма боковой границы 3 не существенна. В этом случае можно принять одно из следующих выражений для х

Второе выражение эффективнее при быстром уменьшении о°(£).

6. Примеры построения координатной сетки. В качестве иллюстрации приведем результаты построения координатной сетки |, т] для затупленного клина (или для разреза затупленного конуса) с гиперболическим вырезом. Внешняя граница 1 образуется окружностью радиуса г1

сопряженной с прямолинейной образующей клина хх = tg <р (у! — у1а), внутренняя граница 2— гиперболой

асимптота которой х = ^у(у—_у2а) параллельна образующей клина.

Рассматриваемая конфигурация зависит от трех параметров: угла полураствора клина ф, относительной толщины носовой части А

х = со —]— аа°, со = const, а = const,

х = 1 + (£, о°),

— 1 4" шах °2 ‘

(5.10)

А + (у і — у и)2 = г

.2

1 »

х\ = tg2 <Р [(у2 - У2ау — (_у20 — у2аП

Уго ~ .У 10 Г\

и относительного расстояния между прямолинейной образующей клина и асимптотой гиперболы б (

Приведенные ниже иллюстрации. отвечают следующим комбинациям параметров:

На рис. 2 и 3 Приведены координатные сетки rj, построенные при выборе к по (6.2). Буквой Г обозначены кривые сопряжения g = const, которые проходят через точку сопряжения на кривой 1.

Случай ф = 0 интересен как предельный. При ф->0 гипербола 2 вырождается в прямую (вырез) на оси симметрии (предельная конфигурация не зависит от б). Часть кривых семейства 1 = const исходит из одной точки — вершины выреза у = уго- Последней такой кривой является вертикальная прямая сопряжения Г на рис. 2 при Д=1 и Г' на рис. 3 при Л = 3. Линия сопряжения Г' появляется лишь при ф = 0. Следовательно, при малых ф малая окрестность Г' является областью резкого изменения величин д', /£, Н^, х^, и G.

Выражения (5.10) принимают вид

;=1 + д;

* = 1 +- °з/Г1 •

(6.1)

(6.2)

8 = 0,01; Д= 1,3; ? = 0; 2,5°; 45°.

Г

Г

а) ср = 0 ; /5 = 1,0

S) <р=2°5- ; А =1,0

Рис. 2

Г

a) ip = 0 ■,< i=S,0

Рис. 3

S) <р - 2°5 ■; й=3,0

Ф-0; & = 1,0

На рис. 4 изображены относительные величины Hi = Hl/H 10 параметра Ламэ на внешней границе 1. Штриховые и сплошные кривые получены при выборе и по (6.1) и (6.2) соответственно. Видно, что выбор х по (6.2) приводит к более медленному убыванию Ну с ростом £. Этот эффект наиболее значителен при Д = 3 и ср=0; 2,5°_(функции Яt для этих двух случаев совпадают); при больших Д выбор х по (6.2) приводит к существенно более равномерной сетке, особенно при небольших ф.

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

ЛИТЕРАТУРА

1. Rieger Н., Projahn U. and Beer Н. Fast iterative solution of poisson equation with neumann boundary conditions in nonorthogonal curvilinear coordinate systems by a multiple grid method. — Numer. Heat Transfer, 1983, vol. 6, N 1.

2. Papantonis D. E. and Anthanassiadis N. A. A numerical procedure for the generation of orthogonal body-fitted coordinate systems with direct deermination of grid points on the boundary. — Int. J. Num.

Meth. Fluids, 1985, vol. 5, N 3.

Рукопись поступила 31/V 1988 г.