CC BY
71
УДК 001.57/519.876.5/004.942
МОДЕЛИРОВАНИЕ ФОРМИРОВАНИЯ СЛОИСТОЙ СТРУКТУРЫ И ПОРИСТОСТИ ПЛАЗМЕННЫХ ПОРОШКОВЫХ ПОКРЫТИЙ С УЧЕТОМ ИЗМЕНЯЕМОЙ ТОПОЛОГИИ ПОВЕРХНОСТИ ПРИ НАПЫЛЕНИИ
В.А. Бледнов, В.И. Иордан*, О.П. Солоненко
Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, г. Новосибирск *Алтайский государственный университет, г. Барнаул E-mail: nightmail@list.ru; jordan@phys.asu.ru
Разработаны и реализованы в среде Delphi 7.0 алгоритмы в составе программного комплекса для моделирования процесса укладки сплэтов (растекшихся и затвердевших на поверхности основы капель расплава) при напылении покрытия сучетом динамически изменяемой топологии его поверхности, а также формирования его слоистой структуры. Для прогнозирования сценариев формирования сплэтов и оперативной оценки их размеров используются экспериментально апробированные теоретические решения. Проведено тестирование программного комплекса сучетом практики напыления.
Ключевые слова:
Моделирование, алгоритм, программный комплекс, капля расплава, кластер покрытия, плазменное напыление, слоистая структура, пористость, топология поверхности, сплэт.
Key words:
Modeling, algorithm, bundledsoftware, droplet of melt, coating cluster, plasma spraying, lamellarstructure, porosity, surface topology, splat.
тура и размер частицы, температура основы и топология ее поверхности.
Для оценки необходимых вычислительных затрат предположим, что необходимо провести моделирование формирования фрагмента покрытия с поперечными размерами Ьх и Ьу, толщиной Н, а напыляемый порошок считается монодисперсным с размером частиц Бг Таким образом, если, в первом приближении, пренебречь объемом пор в покрытии, то оценка необходимого количества сплэтов для полного заполнения объема покрытия У^хЬуН, представится как Мр=6¥/(лБрг). Т. е., если фрагмент покрытия имеет площадь основания 5=1 см2 (ЬХ=ЬУ=1 см) и толщину Н=100 мкм, то при размере напыляемых частиц Бр=10, 25, 50 и 100 мкм необходимо выполнить моделирование формирования сплэтов в количестве А^1,9.107, 1,2.106, 1,5.105и2,0.104, соответственно.
Детальное моделирование процесса формирования одного сплэта требует расчета трехмерной нестационарной краевой задачи со свободной границей для уравнений Навье-Стокса совместно с уравнениями сопряженного конвективно-кон-дуктивного теплообмена и фазовых превращений в растекающейся частице, а в ряде случаев - ив подложке. Если допустить, что для такого расчета требуется 1 с, что является существенно заниженной, как минимум на один-два порядка, величиной (даже при моделировании процессов формирования сплэта на гладкой подложке), тогда время вычислений составит, соответственно, 4,^220 дн., 14 дн., 2 дня и 5,5 ч. Следовательно, в настоящее время такой подход носит скорее теоретический, нежели практический интерес, поскольку моделирование соударения десятков и сотен тысяч капель с подложкой и напыляемым покрытием требует колоссальных вычислительных затрат.
Для прогнозирования структуры и свойств напыляемого покрытия в технологической цепи «плазмотрон - струя - покрытие» важную роль играет моделирование его слоистой структуры, формируемой из полностью расплавленных частиц, рис. 1.
Рис. 1. Схематическое представление слоистой структуры плазменного покрытия, сформированного из полностью расплавленных частиц. Граница между: 1) покрытием и основой, 2) напыленными слоями, 3) сплэтами
Как известно [1], при плазменном напылении, даже, при максимально возможном расходе порошкового материала, частицы которого подвергаются полному плавлению, покрытие формируется путем послойной укладки отдельных сплэтов - растекшихся и затвердевших на основе (подложке или предварительно напыленном слое) капель расплава. Следовательно, достоверность результатов моделирования структуры напыляемого покрытия во многом определяется точностью расчета укладки отдельных сплэтов при заданных значениях ключевых физических параметрах (КФП) соударения капли расплава с основой: скорость, темпера-
Альтернативой ему, является использование теоретических или эмпирических зависимостей, характеризующих диаметр сплэтов в заданном диапазоне КФП, что существенно (на несколько порядков) ускоряет вычислительную процедуру моделирования слоистой структуры покрытий [2-4] и др.
В настоящей статье представлены результаты, связанные с разработкой вычислительного алгоритма и программного комплекса для моделирования процесса укладки сплэтов при напылении покрытия с учетом динамически изменяемой топологии его поверхности, а также формирования пористости и слоистой структуры покрытия с использованием экспериментально апробированных теоретических решений [5], позволяющих предсказывать сценарий формирования сплэта и с приемлемой для практики точностью проводить оперативную («мгновенную») оценку его толщины и диаметра при заданных значениях КФП и теплофизиче-ских свойствах материалов напыляемого порошка и основы.
Процесс напыления характеризуется определенной степенью «стохастичности», подтверждающейся в ходе наблюдений наличием у частиц разброса их КФП. А именно, с выхода плазматрона нагретые частицы потока в виде капель расплава, летят с различными значениями скорости ир в направлениях, близких к направлению нормали к подложке в пределах небольшого телесного угла, при соударении с ней поочередно принимают форму сплэтов, образуя «пятно» напыления, состоящее из множества слоев сплэтов, расположенных друг на друге. Последовательным смещением подложки относительно оси плазменной струи формируется защитное покрытие на подложке.
Процесс моделирования и анализа структурных характеристик пористых покрытий промышленного образца, размеры которых могут достигать несколько десятков сантиметров и более, требует больших временных вычислительных затрат на однопроцессорной машине. Поэтому в процессе моделирования можно ограничиться «кластером покрытия (КП)» - частью объема всего покрытия в виде прямоугольного параллелепипеда (нижнее основание - плоская грань в виде прямоугольника или квадрата), верхнее основание которого в результате процесса напыления оказывается шероховатой поверхностью, рис. 2, а. В объеме КП наблюдается слоистая структура с наличием пор. Размеры нижнего основания КП (на подложке) необходимо задавать порядка нескольких миллиметров для того, чтобы вдоль каждой размерности X и У умещалось хотя бы 20-30 сплэтов (их диаметры могут составлять примерно от 10 и до 100...200 мкм). В противном случае структура отдельно моделируемого КП не будет адекватна структуре аналогичного КП в составе полностью моделируемого покрытия.
В процессе моделирования КП используется дискретная прямоугольная сетка «узловых» точек в плоскости ХУ (Дх и йу - шаги дискретизации, со-
ответственно по осям X и У). Вычисления в алгоритме укладки сплэтов производятся в вещественном формате, но результаты вычислений X, У-коор-динат округляются до целых ^-координаты остаются вещественными). Размеры основания КП на подложке задаются в относительных единицах целыми числами Лх, Лу (реальные размеры равны, соответственно, ЛДх и ЛД). Объем КП равен суммарному объему «элементов разбиения (ЭР)», рис. 2, б, каждый из которых имеет площадь основания Д5=ДхДу, а их количество в объеме КП равно Л=МхЛг Тонкий вертикальный слой толщиной в одну дискретную величину Дх или йу определяет понятие «шлиф».
Рис. 2. Модели структур: а) КП; б) его произвольного элемента разбиения
Учет стохастичности процесса напыления при моделировании реализуется заданием соответствующих гистограмм распределений параметров напыляемых частиц. Т. е. в соответствии с заданными законами распределений датчиками «псевдослучайных чисел» генерируются: xp, yp - целочисленные (обезразмеренные) координаты капли расплава в «лобовой точке» столкновения с напыляемой поверхностью; up - ее скорость в направлении нормали к подложке; Tp - температура частицы; Dp - диаметр частицы. Кроме того, в алгоритме укладки сплэтов на поверхность учитывается «расход» частиц Np, а также ряд теплофизических свойств материала частиц. Последовательный характер укладки сплэтов («splat by splat») в алгоритме позволяет гибко оперировать и дополнять структуру данных, сохраняющую результирующую информацию о структуре последовательно напыляемого покрытия.
Процесс растекания каждой капли на шероховатой поверхности очень сложный, рис. 3, а, поэтому в процессе моделирования принята простая модель сплэта, рис. 3, б, в виде цилиндра с диаметром Ds=2Rs и много меньшей высотой hs.
Для укладки первых слоев сплэтов необходимо знать характеристики материала подложки и ее температуру Ть, так как от этого существенно зависит сценарий растекания и затвердевания капли, а значит и расчетные формулы для размеров сплэта [5]. Параметры, характеризующие свойства напыляемых частиц и подложки, считываются
из базы «справочных» данных программного комплекса. Важно учитывать, что в процессе формирования покрытия его толщина возрастает и влиянием материала подложки на процесс растекания капли в дальнейшем можно пренебречь. С учетом выше сказанного, в программном комплексе специальная функция вычисляет значения параметров сплэта к, и В, [5]. Алгоритм укладки каждого сплэ-та на напыляемую поверхность с учетом ее топологии состоит из двух этапов.
На первом этапе алгоритма вычисляются в вещественном формате «оценки» Z-координат так называемого «опорного» массива вершин у= z¡J); ¡=хр-Я1,...,хр,...,хр+Я-/=ур-Л„...,ур,...,ур+Л5}, посредством которого на втором этапе методом «сплайн-аппроксимации» строится гладкая поверхность Z=F (X, У) нижнего основания укладываемого на напыляемую поверхность сплэта. Т. е., при формировании «опорного» массива вершин в рассмотрение попадают те точки поверхности, для которых целочисленные X, У-координаты (проекция точек поверхности на плоскость ХУ - на плоскость подложки) образуют квадрат со стороной В, и центр квадрата совмещен с проекцией лобовой точкой. Z-коорди-наты опорных вершин, проекции которых на плоскость ХУ не «захватываются» вписанным в квадрат кругом (проекцией сплэта), приравниваются Z-координатам соответствующих точек напыляемой поверхности.
В процедуре оценки Z-координат опорных вершин нижнего основания сплэта принято вполне обоснованное допущение: центральная часть сплэта - ядро сплэта (на рис. 3, б, выделено черным цветом), также представляющее из себя цилиндр диаметром В0=2Я0=кВВр, немногим больше диаметра Вр исходной капли (1<кв< 1,3), своим нижнем основанием полностью повторяет контактную часть напыляемой поверхности, так как жидкая капля расплава при столкновении с поверхностью имеет большое напорное давление. Поэтому Z-ко-ординаты опорных вершин, на которые «опирается» ядро сплэта, приравниваются Z-координатам соответствующих точек напыляемой поверхности.
Периферийная часть сплэта в виде кольца, ограниченного двумя окружностями с диаметрами В0 и В„ образуется за счет растекания капли за пределы ядра сплэта вдоль шероховатой поверхности, которая в зоне кольца сплэта может иметь углубления, в некоторые из которых жидкость капли может и не затечь, образуя поры. Образование пор в углублениях кольца сплэта объясняется превосходством высокой радиальной над направленной вниз нормальной составляющей скорости растекания (жидкость пролетает по инерции).
Рассмотрим более подробно процедуру «оценивания» Z-координат опорных вершин для периферийной части сплэта (кольца сплэта), для которой будем пользоваться «шаблоном узловых точек (ШУТ)» проекции сплэта на плоскость ХУ, рис. 3, в, с учетом исходной дискретизации X, У-координат. Необходимо отметить, что ШУТ с координатами его центра (0,0) реализуется «одноразовым» вызовом отдельной программной функции с учетом на максимально возможный диаметр сплэта еще до выполнения процедуры оценивания Z-коорди-нат опорных вершин для кольца сплэта и хранится в памяти программы до конца моделирования. Для каждого конкретного диаметра сплэта из памяти извлекается лишь только необходимая часть данных ШУТ и производится сдвиг этих координат, соответственно, на величины хр и ур (таким образом осуществляется «привязка» ШУТ к корректным значениям координат сплэта).
Процедура оценивания Z-координат опорных вершин для кольца сплэта учитывает физический аспект растекания капли в виде «цилиндрической» волны. Т. е. на рис. 3, в, отображен порядок последовательного оценивания Z-координат опорных вершин обходом по часовой стрелке ломаной кривой ШУТ, являющейся дискретным представлением дуги окружности - текущего состояния границы фронта этой волны. Стрелки, начинающиеся от границы ядра сплэта радиуса Я0 и заканчивающиеся в текущей точке ШУТ для оценивания ее Z-координаты, поясняют тот факт, что обход предыдущих замкнутых ломанных кривых уже произведен.
Непосредственно оценивание текущей Z-коор-динаты опорной вершины можно пояснить рис. 4, который отображает в виде «жирной» кривой линии профиль напыляемой поверхности вдоль выделенного радиального направления, показанного на рис. 3, в. На рис. 4 точка А соответствует «лобовой» точке столкновения с координатами (хр,ур,£р), с которой совмещен центр сплэта. Как уже отмечалось выше, независимо от топологии шероховатости в зоне напыляемой поверхности, которую покрывает ядро сплэта (на рис. 4 от точки А до В), за счет большого напорного давления капли образуется затвердевший слой сплэта без пор. За пределами ядра сплэта в выделенном радиальном направлении возможны различные сочетания его участков (рис. 4), помеченных значениями признака К (К=0, К=1, К=2, К=3).
Рис. 4. Образец профиля поверхности вдоль выделенного радиального направления
На рис. 4 рассмотрен вариант сочетания участков, в котором Z-координаты точек профиля поверхности непосредственно за ядром сплэта не превосходят Z-координату точки В на величину, большей Н0=ккН, (коэффициент кь немногим больше 1). Такие и подобные этому участки помечаются в алгоритме признаком К=1, т. е. на таких участках профиля допускается натекание жидкости растекающейся капли. Если же превышение Z-коор-динаты точки В оказалось бы на величину, большей к0=кнк, (величина кь с увеличением размера препятствия немного увеличивается), тогда этот участок был бы помечен признаком К=3, и на этот участок препятствия натекание жидкости капли не допускалось бы (на рис. 4 показан подобный участок).
За препятствием жидкость, обтекая его с боковых направлений, в данном радиальном направлении может затечь (на рис. 4 участок помечен признаком К=2) только в случае, когда профиль поверхности оказывается ниже расчетной траектории (на рис. 4 помечена пунктиром) вязкого растекания капли, определяемой с помощью радиальной (ыг=Рп) и нормальной (иг=-^2) составляющих скорости растекания капли из дифференциального уравнения (производная от функции профиля в радиальном направлении равна тангенсу угла затекания жидкости, равного отношению составляющих скорости растекания с учетом поправочного параметра /)
¿7
а' = и =
пт и
-и
(1)
Решением (1) с учетом начального условия 1(гг)=1г является зависимость
7(т) = 7С + Я0[(г/Я0)(с /Я0уи]. (2)
Точки расчетной траектории (2) растекания капли, которые оказываются выше точек профиля поверхности, на рис. 4 определяют участок с признаком К=0 и именно они сохраняются в массиве опорных вершин. Участок с признаком К=0 может следовать только за участком ядра сплэта или только, как показано на рис. 4, за участком К=1 (и не может следовать за участками с признаками К=3 и К=2).
Таким образом, в массив опорных вершин {/,}, соответствующих точкам нижнего основания сплэ-та, вместо точек профиля каждого радиального направления из участков типа К=0 войдут точки расчетных траекторий этих участков, а также все точки профиля оставшихся участков каждого из этих направлений. Как уже отмечалось выше, кроме опорных вершин, на которые «опирается» сплэт, в массив опорных вершин {/,} войдут и точки поверхности за пределами сплэта, проекции которых на плоскость вместе с проекцией сплэта образует квадрат, центрированный в лобовой точке (иначе говоря, проекция кругового сплэта вписана в квадрат).
Сформированный для второго этапа алгоритма укладки сплэта массив опорных вершин {/,}, проекция которых на плоскость ХУ образует квадрат, центрированный в лобовой точке, позволяет аппроксимировать гладкую поверхность с использованием «составных рациональных В-сплайновых поверхностей» в параметрической форме [6]. Однако в результате такой сплайн-аппроксимации некоторые расчетные Z-координаты точек В-сплайновой поверхности могут оказаться ниже соответствующих Z-координат точек напыляемой поверхности, на которую укладывается сплэт. В этом случае, исходя из физического смысла, в соответствии с алгоритмом такие «заниженные» расчетные Z-координаты точек В-сплайновой поверхности необходимо заменять на Z-координаты соответствующих точек напыляемой поверхности, на которую укладывается сплэт.
Как уже говорилось выше, в любом случае Z-координаты точек нижнего основания сплэта из зоны ядра должны быть также восстановлены прежними значениями, но тогда возможно возникновение больших перепадов (разрывов) нарас-четной поверхности нижнего основания сплэта в точках круговой границы между ядром и периферийной частью сплэта. Поэтому, чтобы исключить возможность таких перепадов (разрывов), на втором этапе алгоритма используется одно из свойств В-сплайновых поверхностей [6], заключающееся в следующем. Если продублировать одну из точек, хотя бы 2 раза (сделать ее «повторной» с кратностью 3), тогда в результате сплайн-аппроксимации расчетная В-сплайновая поверхность практически проходит через эту опорную повторную вершину. Следуя этому свойству, в алгоритме обеспечивается кратность повторности вершин круговой границы ядра сплэта не ниже 3 (обе размерности опорного массива вершин возрастают на 4).
Таким образом, по расширенному опорному массиву вершин рассчитывается гладкая поверхность ¿=¥(Х, У) (точнее говоря, массив сглаженных вершин) с использованием «составных рациональных 5-сплайновых поверхностей» в параметрической форме [6] согласно выражениям Г V /-.. .лл
Я. 1(и, V) =
X,,(и,V)
\, (и, V)
V 1
(и, V)
3 3
-1+г I -1+1
к=0 I=0_
3 3
""щл (и )п1 <У)
к=0 I=0 1=0,1,2.
0 < и, V < 1;
• , т;;=0,1,2,..., п; т=п=Б+4;
3 3
""пк(и )п1(v) =1
(3)
к=0 I=0
Функциональные коэффициенты п^ (и) определяются по формулам п0(и)=(1-и)3/6, п1(и)=(3и3—6и2+4)/6, п2(и)=(—3и3+3и2+3и+1)/6, п3(и)=и3/6 и подчиняются условию «нормировки» (3). При этом кубические многочлены п(и) рассчитываются по аналогичным формулам с учетом замены параметра и на у. Коэффициенты определяют вес соответствующей вершины. Чем выше значение какого-либо из этих коэффициентов по отношению к другим, тем ближе поверхность подходит к соответствующей вершине. Таким образом, эти весовые коэффициенты позволяют «регулировать» форму нижнего основания сплэта и, тем самым, «регулировать» форму и «эффективные» размеры пор, образующихся под укладываемым на напыляемую поверхность сплэтом.
В результате аппроксимации на 5-сплайновой поверхности расчетные значения координат вершин, соответствующие круговой границе нижнего основания ядра сплэта, практически будут совпадать с их исходными координатами, так как кратность повторений вершин не менее трех. Следовательно, на границе ядра и периферийной части нижнего основания сплэта обеспечивается условие «сшивания» этих 2-х частей.
Затем из массива пересчитанных в результате сплайн-аппроксимации вершин отбрасываются предварительно добавленные 4 строки и 4 столбца повторных вершин, а также и те вершины поверхности, которые не входят в зону кругового сплэта. Далее, как уже отмечалось выше, восстанавливаются исходные значения координат вершин, соот-
ветствующих зоне ядра сплэта, и исходные значения координат вершин из периферийной части сплэта, «помеченных» на первом этапе алгоритма признаками К=1, К=2, и К=3. На вершины с признаком К=3 из-за их высоких значений Z-коорди-нат жидкость капли натекать не может. Поэтому величина эффективной площади ^ поверхности нижнего основания сплэта, на которую затем равномерно распределяется весь объем сплэта, уменьшается на величину площади, которую составляют точки с признаком К=3.
Эффективную площадь ^ поверхности нижнего основания сплэта образуют точки зоны ядра сплэта, точки расчетных траекторий (2) из участков типа К=0, точки профилей поверхности из участков типа К=1 и К=2 всех радиальных направлений. Т. е. исходную расчетную толщину к сплэта необходимо скорректировать (увеличить) с учетом равенства объемов укладываемого сплэта и исходной капли как к=к(пК?/8/). Затем в соответствующей структуре данных фиксируются координаты точек верхнего основания, параллельного нижнему основанию укладываемого сплэта с высотой к,. Тем самым завершается второй этап алгоритма моделирования укладки на напыляемую поверхность очередного сплэта. Алгоритм укладки на напыляемую поверхность очередного сплэта циклически повторяется до тех пор, пока не будет «исчерпано» заданное значение расхода частиц Затем, как показано на блок-схеме, рис. 5, рассчитываются структурные характеристики КП (шероховатость, адгезия и когезия).
В качестве иллюстрации работы программного комплекса моделирования процесса формирования слоистой структуры плазменных порошковых покрытий и его тестирования для тестовой медной подложки форматом 2x2 мм и напыляемых на нее частиц никеля в количестве 104 и размерами Хр=20...40 мкм получены следующие значения характеристик КП, рис. 6, 7. Общая пористость КП (отношение объема пор к общему объему КП) равна Рет=18,5 %, диапазон распределения пористости по 4 слоям шлифа составляет 16.19%, рис. 7, а, шероховатость (перепад высоты между максимально высокой и минимально низкой точками поверхности КП) АИШ=60 мкм, средняя высота КП Нш=433 мкм, адгезионная прочность (отношение покрытой сплэтами площади подложки КП к площади самой подложки) равна 68 %, когезионная прочность в сечениях КП (отношение суммарной контактной площади сплэтов из двух соседних слоев к площади сечения этих слоев) равна 75.90 %.
Установка
параметров
процесса
напыления
покрытия
Инициализация ШУТ сплэтов и структур данных для кластера покрытия
Расчет толщины и диаметра очередного сплэта
Определение опорного массива вершин с помощью ШУТ сплэта
Сплайн-аппроксимация формы сплэта и ее сохранение в структуре данных кластера покрытия
Расчет шероховатости адгезии и когезии для кластера покрытия
Рис. 5. Блок-схема процесса моделирования формирования слоистой структуры КП
Рис. 6. Внешний вид программного интерфейса с результатом моделирования КП
• i . ч- '
' s, У
Гг *
v ( jf - _ /Ц
L v- - ,. * > -v J
V 1
и:
У. L. Л
ч
а бе
Рис. 7. Результаты моделирования: а) пример изображения шлифа с гистограммой пористости по слоям; б) изображение первого монослоя КП (вид сверху); в) изображение горизонтального разреза КП (вид сверху) на высоте 42 мкм от подложки
Выводы
С использованием среды программирования Delphi 7.0 и алгоритма укладки сплэтов на поверхность с произвольной шероховатостью создан и апробирован в практике напыления программный комплекс для численного моделирования формирования слоистой структуры плазменных покрытий, напыляемых из металлических порошков. В основу алгоритма положены экспериментально апробированные теоретические решения, характеризующие диаметр и толщину сплэтов при соударе-
СПИСОК ЛИТЕРАТУРЫ
1. Кудинов В.В., Пекшев П.Ю., Белащенко В.Е. и др. Нанесение покрытий плазмой. - М.: Наука, 1990. - 408 с.
2. Knotek O., Elsing R. Monte Carlo Simulation of the Lamellar Structure of Thermally Sprayed Coatings // Surf. Coat. Technol. -1987. - №32. - P. 261-271.
3. Ghafouri-Azar R., Mostaghimi J., Chandra S., Charmchi M. A Stochastic Model to Simulate the Formation of a Thermal Spray Coating // Journal of Thermal Spray Technology. - 2003. - V. 12. -№ 1. - P. 53-69.
нии капель металлических расплавов с подложкой. Укладка сплэтов на поверхность с произвольной шероховатостью производится с использованием сплайн-сглаживания, учитывающего гидродинамические особенности растекания капель. В качестве иллюстрации приводятся результаты расчетов формирования покрытия из порошка никеля, напыляемого на медную подложку.
Работа выполнена при частичной поддержке Президиума СО РАН, Междисциплинарный интеграционный проект № 1 на 2009-2011 гг., а также гранта РФФИ № 09-01-00433-а.
4. Steinke T., B?ker M. Monte Carlo Simulation of Thermal Sprayed Coatings // Proc. ITSC'06, Seattle, USA, May 15-18, 2006. -P. 329-334.
5. Солоненко О.П., Алхимов А.П., Марусин В.В. и др. Высокоэнергетические процессы обработки материалов. - Новосибирск: Наука, 2000. - 425 с.
6. Шикин Е.В., Плис А.И. Кривые и поверхности на экране компьютера. Руководство по сплайнам для пользователя. - М.: ДИАЛОГ-МИФИ, 1996. - 240 с.
Поступила 02.09.2010г.
