CC BY
102
Исследование особенностей разрушения хрупких керамических покрытий на основе метода подвижных клеточных автоматов
С.Г. Псахье, Д.Д. Моисеенко, А.Ю. Смолин, Е.В. Шилько, А.И. Дмитриев
Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия
В данной работе на основе метода подвижных клеточных автоматов проведен анализ особенностей развития процесса генерации и развития повреждений в хрупких покрытиях. Применение метода подвижных клеточных автоматов позволяет проводить исследования, начиная от генерации первых повреждений и заканчивая формированием макротрещин и потерей целостности структуры. Результаты моделирования подтверждают экспериментальные данные о том, что граница раздела в материале с покрытием, являясь макроконцентратором напряжений, служит причиной квазипериодического растрескивания. Было показано, что при выходе трещины на поверхность покрытия происходит формирование волны сдвига, инициирующей образование вторичных трещин.
1. Введение
Задача исследования механизмов деформации и разрушения гетерогенных сред и структур является достаточно общей и касается самых разнообразных приложений физики и механики деформируемого твердого тела. В частности, изучение таких механизмов является важной задачей для исследования особенностей разрушения композиций типа «керамическое покрытие -эластичная подложка».
В соответствии с основными концепциями физической мезомеханики, поведение гетерогенных систем в процессе нагружения представляет собой сложный нелинейный процесс, самосогласованно развивающийся на разных масштабных уровнях и включающий в себя движение мезообъемов, фрагментацию материала и, наконец, разрушение как последнюю стадию процесса [1-3]. Как показано в работах [1, 2], в результате нагружения гетерогенного материала каждый локальный объем деформируется по схеме «сдвиг + поворот» и упругие напряжения внутри образца распределяются весьма сложным образом. Вследствие сильной неоднородности упругих модулей на границах раздела (например в переходном слое между покрытием и подложкой) в материале формируются градиенты напряжений, причем уровень этих напряжений в локальных областях может значительно превышать уровень средних приложенных напряжений [1-4].
Способ упрочнения материала посредством нанесения керамических покрытий является перспективной технологией, которая позволяет повышать не только износостойкость обрабатываемых поверхностей, но и в ряде случаев усталостную прочность изделий. Как показано в работах [3,4], для описания прочности ком-
позиций данного класса важным фактором является учет формирования блочной структуры, приводящей к множественным локальным нарушениям сплошности материала и в результате к глобальному перераспределению упругой энергии моделируемого объекта. Поскольку поведение таких композиций под нагрузкой всегда сопровождается генерацией и развитием повреждений, то возможности континуальных методов для компьютерного моделирования отклика подобных систем являются ограниченными [5]. Данная проблема может быть решена в рамках дискретного подхода, например на основе созданного в последние годы метода подвижных клеточных автоматов (метода MCA — от английского «movable cellular automata») [6-8].
В рамках данного метода моделируемый материал представляется как ансамбль элементов, взаимодействующих между собой по определенным правилам. В отличие от классического подхода взаимодействие подвижных клеточных автоматов приводит не только к изменению их состояния, но и к пространственному перемещению, при этом окружение элементов может меняться. Благодаря мобильности отдельных элементов, такой подход позволяет моделировать различные процессы, имеющие место в реальном гетерогенном материале, включая эффекты проникания, массопе-ренос, эффекты фрагментации, формирование повреждений, развитие трещин и т.д. При этом эффекты потери сдвиговой устойчивости могут быть учтены явным образом.
Путем задания условий на границах можно имитировать различные режимы механического нагружения (сжатие, растяжение, сдвиговая деформация и др.), имеющие место в реальных условиях эксплуатации.
© Псахье С.Г., Моисеенко Д.Д., Смолин А.Ю., Шилько Е.В., Дмитриев А.И., 1998
Следует отметить, что метод МСА позволяет моделировать поведение образцов и структур различного масштаба в целом, что, в свою очередь, делает этот метод полезным для инженерных расчетов.
2. Постановка компьютерного эксперимента и особенности модели
На основе метода подвижных клеточных автоматов была проведена серия численных экспериментов по исследованию процесса разрушения керамического покрытия. Поскольку в данном исследовании наибольший интерес представляло изучение разрушения хрупкого покрытия, акцент был сделан на описании процессов, протекающих преимущественно в упругой области деформаций.
Моделируемые образцы представляли собой двухслойные структуры, состоящие из сравнительно мягкого слоя, имитирующего подложку, и твердого хрупкого покрытия. В рамках данной работы не ставилось задачи исследовать процессы локализации пластической деформации и разрушения мягкой подложки, а поэтому предел упругости материала подложки был принят заведомо большим чем предел его прочности. При этом параметры автоматов задавались таким образом, чтобы подложка не имела возможности разрушаться. Также проводилось моделирование трехслойных структур с дополнительной прослойкой в зоне контакта покрытия и подложки.
Упругие модули и прочностные характеристики моделируемых материалов выбирались следующим образом. Модуль Юнга материала, имитирующего подложку, был равен 71 ГПа и соответствовал модулю Юнга алюминия. Модуль Юнга покрытия был равен 200 ГПа, что соответствует модулю Юнга спеченной из нанокрис-таллического порошка ZrO2 керамики, механический отклик которой достаточно хорошо описывается в рамках данного метода [8]. Все остальные механические параметры керамики также задавались в соответствии с реальным материалом [9]. Размеры образца были равны 4.5х0.5 см. Толщина подложки при этом составляла 0.17 см, а толщина покрытия - 0.33 см. Нагружение структуры имитировалось посредством задания на торцевых площадках подложки постоянной скорости 10 см/сек (рис. 1).
В целях имитации поведения материала, близкого к реальному, в параметры функции отклика и механические характеристики автоматов был заложен стохастический разброс величиной в 3% по линейному закону распределения.
Необходимо отметить, что вследствие различия упругих модулей зона контакта покрытия и подложки является мощным макроконцентратором напряжений, который оказывает основное влияние на характер всего дальнейшего процесса деформации и разрушения. Условия нагружения (т.е. навязанная деформация подложки), реализующиеся в данном численном эксперименте, накладывают отпечаток на состояние границы раздела таким образом, что вдоль всей зоны контакта покрытия и подложки образуется резкий перепад сдвиговой компоненты напряжений по обе стороны этой границы. Эти эффекты впервые были экспериментально обнаружены и исследованы в работах [3, 4, 10, 11].
В работе [4] было показано, что в зоне границы раздела формируются периодические мезоконцентраторы напряжений, инициирующие осцилляции полей напряжений и сдвиговой деформации вдоль всей границы раздела. Поэтому, прежде чем приступить к анализу результатов теста на растяжение данной структуры, необходимо привести некоторые необходимые сведения из теоретических основ метода подвижных клеточных автоматов, а именно уравнения движения и вид функции отклика автомата при сдвиговой деформации.
В рамках метода подвижных клеточных автоматов изменение состояния пары автоматов определяется их относительным перемещением. Отдельные пары не являются независимыми, так как на перемещения автоматов в паре оказывает влияние их окружение. Совокупность отношений пар подвижных клеточных автоматов образует бистабильную активную среду [12]. Пусть задан потенциал центрального межавтоматного взаимодействия ф(М). Тогда выражение для парной силы взаимодействия автоматов р имеет вид:
= Р1 ■-
(1)
где 1, і — номера автоматов в паре; Ь — расхождение центров автоматов.
При этом необходимым условием к Ф(Ь1) является требование наличия стационарного состояния пары в точке Ы = 0.
Сила центрального многочастичного взаимодействия в паре автоматов і и і (с учетом наличия остальных соседей) определяется по формуле:
Pi1 =
m + m1 1 (( 1
т ІІ- + -71Р* +
m ту
, (2)
Рис. 1. Общий вид структуры образца
где у(агу к) и у(а..,р — некоторые коэффициенты, определяемые углами а., ки а. ^ между направлениями у и Д соединяющими центры элементов пар — и]-1; г и
у — номера автоматов в паре; к и ' — номера соседей каждого автомата пары.
Очевидно, что при моделировании процессов на мезо- и макромасштабных уровнях, наряду с трансляционным движением и изменением объема отдельных клеточных автоматов, необходимо также учесть их повороты. В данном случае необходимо вводить дополнительный параметр, характеризующий состояние отношения пары элементов г и у — угол относительного разворота автоматов 9. Следует отметить, что параметры ¥ и 9. являются полностью независимыми.
Эволюция параметра 9.у описывается дифференциальным уравнением:
(3)
где дг — расстояние от центра элемента г до точки контакта автомата г с автоматом у; Т — парная сила тангенциального взаимодействия элементов; Зг — момент инерции автомата. Необходимым условием для Т является наличие устойчивого стационарного состояния пары в точке 9. =0. Ньютоновские уравнения движения для поворотной степени свободы элементов среды в терминах углов относительных разворотов пар запишутся в виде:
й29г (дгу дуг) у ,Лдгк
но и возможность деформирования по схеме «сдвиг + поворот», наиболее часто реализующейся в задачах, связанных с локализацией деформации и интенсивным разрушением [11].
3. Формирование волны сдвига и генерация вторичных трещин
При вышеописанном моделировании нагружения двухслойной структуры каждый элемент объема под воздействием сдвиговых напряжений деформировался таким образом, что наряду с изменением объема претерпевал разворот. Начиная с момента установления динамического режима нагружения и вплоть до возникновения первых повреждений, все поле автоматов моделируемой структуры разделилось на две равные части с противоположными по знаку углами разворотов: автоматы левой половины образца претерпевали развороты по часовой стрелке (знак «+» на рис. 2, б-г), а автоматы правой половины — против часовой стрелки (знак «-»). Этот факт объясняется тем, что в процессе нагружения подложка, имея более низкие модули упругости, деформируется сильнее чем твердое покрытие. В результате этого на границе раздела возникает резкий перепад продольных деформаций (относительно оси приложенной нагрузки X), т.е. величина де /ду вдоль зоны контакта становится гораздо больше, чем в остальных участках образца, и элементы объема вынуждены испытывать сильный сдвиг и разворот; а в силу самоорганизации элементного поля развороты приобрели автоматы всего образца.
(4)
где 5(гу, гк) и 5(гу,]') — некоторые коэффициенты переноса параметра 9. Отметим, что для случая мгновенного взаимодействия эти коэффициенты равны единице.
Выражение для расчета тангенциального многочастичного взаимодействия элементов г и у в паре с учетом остальных соседей выглядит следующим образом:
+ + + + + + + + + + + +/-) + ++++++++++++ + + + + + + + + + + + + <,-( +
Тгу ■
Г/дг] + Т]/д]г 1
((+у
\
ту +
+ + + + + ++\------г+ + + + + + + }-■
+ + + + + + 3----\ + + + + ++,+( -
+ + + + + ++<'----------\ ++++++?-■
гк у' \
+ У т гк + V ту'
(5)
\+ + + + + + + + + + + + \++ + + + + + + + + + +/ ) ++ + + + + + + + + +
Как видно из этого соотношения, тангенциальное взаимодействие автоматов пары зависит от их смещений относительно общей точки контакта. От этих же смещений зависят взаимные развороты элементов (см. (4)). Таким образом, в рамках данного подхода при моделировании нагруженного состояния в каждый элемент объема естественным образом закладывается не только возможность деформирования по всем осям,
Рис. 2. Распространение волны сдвига при раскрытии трещины: а — структура связей автоматов при выходе на поверхность центральной трещины; б-г —распространение волны сдвига в последовательные моменты времени; д — образование вторичных трещин при прохождении отраженной волны
б
Рис. 3. Последняя стадия разрушения керамического покрытия: а — структура связей образца; б—знаки самоорганизованных разворотов элементов в блоках
Так, в данном численном эксперименте зарождение трещин наблюдалось у границы раздела, причем их пространственное распределение носило квазиперио-дический характер. На первом этапе разрушения наблюдалось образование трещин в центре и у концов образца (рис. 2, а). Очевидно, что в процессе дальнейшего растрескивания керамики пространственный период образования трещин уменьшался.
При выходе трещины на поверхность керамического слоя сдвиговые компоненты напряжений вблизи трещины начали быстро релаксировать, вследствие чего произошло зарождение волны сдвига, распространяющейся от трещины к торцевым граням структуры. При этом направления разворотов автоматов (знак сдвиговой составляющей деформации у ) изменились на противоположные по сравнению с направлениями до генерации повреждений. Три последовательные стадии распространения волны изображены на рис. 2, б-г. При движении отраженной упругой волны сдвига от торцевых граней образца происходит генерация «вторичных» трещин на границе раздела, как это показано на рис. 2, д. Следует отметить, что эти «вторичные» трещины образуются на расстоянии, приблизительно равном толщине покрытия.
4. Формирование блочных структур и квазипериодическое растрескивание покрытия
Теоретические и экспериментальные исследования, выполненные Паниным В.Е. и Гриняевым Ю.В. [3, 4, 11], четко показали, что граница раздела, являясь мощным концентратором напряжений, инициирует формирование осциллирующих концентраторов напряжений и связанное с ними образование блочных структур, движущихся по схеме «сдвиг + поворот».
В описываемом численном эксперименте при дальнейшем нагружении образца наблюдалось квазипе-риодическое растрескивание покрытия по всей длине с периодом, приблизительно равным толщине покрытия. При этом можно видеть, что поле автоматов образца начало делиться на блоки с разворотами, чередующимися по знаку. Причем линии перемены знака совпадали с зарождающимися у поверхности раздела трещинами. Здесь необходимо отметить, что не всем трещинам можно поставить в соответствие границы
н
г иА
2
< ,3 \
Л
ц ■ \
О 0.002 0.006 0.01 1, с
Рис. 4. Диаграммы нагружения различных структур: 1 — керамическое покрытие на алюминиевой подложке; 2 — структура с корундовым дополнительным слоем; 3 — структура с прослойкой из органического волокна
блоков с разными по знаку ориентациями разворотов элементов. Это объясняется тем, что в данной сложно самоорганизованной системе задействованы различные критерии потери сплошности, а достижение некоторого максимума сдвиговой деформации — только один из них. Этот результат схематически изображен на рис. 3.
Очевидно, что в результате действия внешней нагрузки происходит падение способности сопротивляться внешней нагрузке (рис. 4, кривая 1). Вследствие того, что керамический слой и подложка имеют существенно различные упругие модули, при растяжении на поверхности раздела формируется сильный скачок напряжений. При выходе трещины на поверхность эти напряжения релаксируют, что приводит к колебаниям подложки, изгибающим образец в направлениях, перпендикулярных оси приложенной нагрузки. Поля смещений автоматов изображены на рис. 5 (масштаб смещений увеличен).
Таким образом, в данном численном эксперименте показано, что процесс зарождения осцилляции деформаций имеет место уже в упругой области, т.е. на стадии формирования первых нарушений сплошности материала вдоль границы раздела.
5. Влияние параметров границы раздела
Поскольку состояние границы раздела оказывает существенное влияние как на характер распределения деформации, так и на разрушение всей композиции, в данной работе представляло интерес исследовать влияние свойств зоны контакта на прочностные характе-
Рис. 5. Поля смещений в образце при выходе трещины на поверхность
ристики данной структуры в целом. В двух следующих численных экспериментах проводилось моделирование структур с дополнительной прослойкой в зоне контакта покрытия и подложки. Рассматривались два различных случая: 1) промежуточный слой обладает более высокими упругими и прочностными свойствами; 2) слой — с низкими модулями и относительно малой прочностью. Толщина покрытия при этом увеличилась до 0.42 см. Тип и скорость нагружения были идентичны предыдущему случаю.
В первом случае в качестве материала переходного слоя был выбран корунд (А1203). В результате введения такой прослойки прочность данной структуры существенно повысилась (рис. 4, кривая 2). Как и в предыдущем численном эксперименте, в образце с переходным слоем наблюдалось квазипериодическое растрескивание керамического и корундового слоев. Однако пространственный период возникновения трещин существенно увеличился. Более того, скорость образования и распространения трещин в образце была гораздо меньше. В данном численном эксперименте высокие упругие модули и высокая прочность корунда привели к уменьшению градиента деформации вдоль границы «корунд -керамика», что затормозило развитие и образование трещин.
Следует обратить внимание на то, что сильная неоднородность свойств материалов покрытия и подложки приводит к различным специфическим эффектам на границах раздела. Поскольку локальное нарушение сплошности в зоне контакта является дополнительным концентратором напряжений, его релаксация приводит к распространению потоков дефектов, формирующих существенно неоднородное состояние вблизи этого концентратора. Так, в момент прохождения трещины через поверхность раздела «корунд - керамика» наблюдается распространение области локализации деформации, связанное с сильным динамическим перераспределением упругих напряжений вблизи данной локальной области (рис. 6).
Наряду с введением прочного корундового слоя, было проведено моделирование нагружения аналогичной структуры с переходным слоем из органического волокна, модуль Юнга которого в два раза меньше модуля Юнга подложки. Несмотря на довольно низкую
■ ■V. *' 1г»<*м
■* Р 1 II Г р * I'
■ и г ( I ; * * - -*
ГА"
1 Г1 г / г г .* \ / Гп
* н * н м „ н
0
Рис. 6. Эффект распространения локализации деформации при раскрытии трещины на границе раздела. Местоположение трещины указано стрелкой: а — локальный разрыв связей на границе «корунд -керамика»; б — поле скоростей в области локализации деформации
Рис. 7. Структура образца с изначально «гофрированным» профилем зоны контакта
прочность данной структуры, ее предельная деформация по сравнению со структурой без переходного слоя в расчетах оказалась на 50% выше. При этом период квазипериодического растрескивания покрытия был больше, чем в структуре без прослойки. Заметим, что также был проведен численный эксперимент по одноосному растяжению структуры без прослойки с увеличенной толщиной покрытия, равной суммарной толщине покрытия и дополнительной прослойки в двух предыдущих случаях. При сравнении пространственных периодов растрескивания во всех случаях, когда толщина покрытия была увеличена, можно сделать вывод о том, что характер развития повреждений зависит не только от параметров переходного слоя, но и от его толщины. Этот результат хорошо согласуется с экспериментальными данными [13], свидетельствующими о том, что период растрескивания прямо пропорционален толщине покрытия.
Необходимо отметить, что благодаря резкому скачку механических характеристик, зона контакта покрытия и подложки является сильным концентратором напряжений, который оказывает основное влияние на характер всего дальнейшего процесса деформации и разрушения. Поэтому большое значение при анализе вышеописанных эффектов имеет не только состояние границы раздела в смысле тех или иных ее механических характеристик, но и геометрический фактор этой границы, т.е. форма ее профиля.
Было проведено моделирование нагружения структуры с изначально «гофрированным» профилем зоны контакта. Следует обратить внимание на то, что такие «гофрированные» структуры являются не просто предметом чисто теоретического интереса, но также весьма часто имеют место (как правило вследствие различных коэффициентов термического расширения различных материалов) в тех или иных технологиях нанесения покрытий или получения слоистых материалов (например при росте пленок эпитаксиальных структур для различных приборов микроэлектроники). Структура образца показана на рис. 7.
При нагружении структуры также наблюдалось квазипериодическое растрескивание, однако динамика развития повреждений имела иной характер. Если в случае структуры с плоским профилем зоны контакта на первых стадиях разрушения происходило образование трещины в центре образца, то при моделировании данной структуры трещины начинали зарождаться и быстро расти вблизи краев образца на равноудаленных расстояниях от центра (рис. 8, а).
б
Рис. 8. Особенности разрушения структуры с «гофрированным» профилем границы раздела: а — структура связей образца; б — углы разворотов элементов
При исследовании образования волны сдвига и разворотов элементов было также обнаружено существенное отличие от поведения структуры с плоским профилем зоны контакта. В данном случае гофрированный профиль выступал в качестве сдерживающего фактора распространения волны сдвига: при выходе трещин на поверхность покрытия происходило зарождение волны сдвига, которая сразу же затухала.
При детальном анализе этого эффекта было обнаружено, что одновременно с ростом трещины в месте «углубления» подложки происходит зарождение «сопряженной» трещины по другую сторону «выпуклости» подложки. Очевидно, что эта трещина порождает вторую волну сдвига, которая гасит первую. При этом в области между трещинами происходит образование блока с гораздо большими углами разворотов элементов, однако знаки разворотов совпадают с близлежащими участками образца. На рис. 8, б показано поле углов разворотов элементов. Углы графически изображаются в виде окружностей, диаметр которых характеризует абсолютное значение угла. Места, в которых проходят трещины, отображаются жирными черными линиями.
Этот результат согласуется с экспериментальными данными и ранее проведенными численными экспериментами, описанными в [14, 15].
6. Заключение
Проведенные численные расчеты и анализ результатов экспериментов [3, 4, 10, 11, 13] по исследованию особенностей разрушения керамических покрытий на эластичной подложке показал, что в рамках метода подвижных клеточных можно корректно моделировать генерацию и развитие повреждений на границах раздела хрупкого покрытия и подложки. При этом важной особенностью метода является способность автоматов (элементов структуры) перемещаться по схеме «сдвиг + поворот», т.е. реализовывать не только трансляционный, но и ротационный отклик. Расчеты показали, что процесс формирования блочных структур, описанный в [3, 4], в керамическом покрытии приводит к инициированию повреждений в упругой области деформаций.
Проведенные расчеты показали, что при выходе трещины на поверхность покрытия возникает упругая
волна сдвига, причем дальнейшее распространение этой волны вдоль зоны контакта и ее отражение от боковых граней образца может приводить к динамической генерации вторичных трещин.
В заключение хотелось бы поблагодарить академика Панина В.Е. за постановку задачи, весьма полезное участие в обсуждении данной проблемы и дискуссиях относительно полученных результатов, а также за предоставленные экспериментальные данные.
Работа выполнена в рамках Федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития науки и техники гражданского назначения» по проекту 07.08.00600.М.
Литература
1. Физическая мезомеханика и компьютерное конструирование материалов: В 2-х т./ Под ред. В.Е. Панина. - Новосибирск: Наука, 1995. -297 с. и 320 с.
2. Physical Mesomechanics of Heterogeneous Media and Computer-Aided Design of Materials/ Ed. by V.E. Panin. - Cambridge International Science Publishing, 1998. - 339 p.
3. Панин В.Е. Основы физической мезомеханики// Физическая мезомеханика. - 1998.- Т. 1.- № 1.- С. 5-20.
4. Панин В.Е., Плешаное В.С., Гриняев Ю.В. и др. Формирование периодических мезополосовых структур при растяжении поликристаллов с протяженными границами раздела// ПМТФ. - 1998. - Т. 39. -№4.- С. 141-148.
5. ГолоенееИ.Ф., ГолоенееаЕ.И, КонееА.А., ФоминВ.М. Физическая мезомеханика и молекулярно-динамическое моделирование // Физическая мезомеханика. - 1998. - Т. 1. - № 2. - С. 21-33.
6. Псахъе С.Г, Хори Я., Коростелев С.Ю. и др. Метод подвижных клеточных автоматов как инструмент для моделирования в рамках физической мезомеханики// Изв. вузов. Физика. - 1995.- № 11.- С. 58-69.
7. Псахъе С.Г., Коростелев С.Ю., Смолин А.Ю. и др. Метод подвижных клеточных автоматов как инструмент физической мезомеханики материалов// Физическая мезомеханика. - 1998. - Т. 1.- № 1.- С. 95108.
8. Псахъе С.Г., Моисеенко Д.Д., Дмитриев А.И. и др. О возможности компьютерного конструирования материалов с высокопористой и каркасной структурой на основе метода подвижных клеточных автоматов//Письма в ЖТФ. - 1998. - Т. 24. - № 4. - С. 71-76.
9. Физические величины: Справочник. - М.: Энергоатомиздат, 1991. -1232 c.
10. Панин В.Е. Методология физической мезомеханики как основа построения моделей в компьютерном конструировании материалов// Изв. вузов. Физика. - 1995. - № 11. - С. 6-26.
11. Панин В.Е. Физическая мезомеханика материалов// Изв. вузов. Физика. - 1998. - № 9. - С. 8-37.
12. ШилъкоЕ.В. Изучение отклика твердого тела на мезоуровне на основе развития подхода клеточных автоматов с явным учетом эффектов массопереноса: Дис. ... канд. физ.-мат. наук. - Томск, 1997. - 155 с.
13. Панин В.Е., Слосман А.И., Колесова Н.А.Закономерности пластической деформации и разрушения на мезоуровне поверхностно-упрочненных образцов при статическом нагружении // ФММ. -1996. - Т. 82. - Вып. 2. - С. 129-136.
14. Ковалъ А.В., Панин С.В., Трусова Г.В. Влияние поверхностного слоя, упрочненного методом борирования, на характер пластической деформации на мезоуровне стали 15Н3МА// Конф. молодых ученых «Физическая мезомеханика материалов»: Тезисы. - Томск, 1998. -С.17-18.
15. Моисеенко Д.Д., Татаринцев Е.М., Псахъе С.Г. О возможностях компьютерного конструирования материалов с покрытиями на основе метода подвижных клеточных автоматов// Конф. молодых ученых «Физическая мезомеханика материалов»: Тезисы. - Томск, 1998. -С. 40^1.
