Численное моделирование механического поведения модельных хрупких пористых материалов на мезоуровне | Вестн Том. гос. ун-та. Математика и механика. 2013. № 5(25).

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

Для описания механического отклика мезообъемов хрупких пористых материалов применен эволюционный подход, а сам материал рассматривается как нелинейная динамическая система. Определяющие нелинейные соотношения описывают накопление повреждений и их влияние на деградацию прочностных и упругих свойств каркаса. В рамках иерархического моделирования пористость мезомасштаба учтена явно, а на более мелких масштабах - интегрально как мера поврежденности. Этот тип мелкомасштабных повреждений учтен в определяющих уравнениях как дилатансия, в том числе сдвиговой природы. Рассмотрено два типа морфологии модельных пористых структур: перекрывающиеся сферические поры и перекрывающиеся сферические тела. Полученные результаты численного моделирования свидетельствуют о слабом влиянии морфологии пор на особенности развития поврежденности в каркасе. А усредненная диаграмма нагружения оказывается чувствительной не только к значению пористости, но и к морфологии пористой структуры. Аналогичную зависимость имеет и предел прочности.

Numerical modelling of mechanical behaviour of model brittle porous materials.pdf Хрупкие пористые материалы широко встречаются как среди природных, так и среди искусственно создаваемых материалов. Большинству природных геоматериалов присуща исходная пористость или трещиноватость. Костные ткани обладают сложной пористой структурой. При изготовлении керамических материалов в результате особенностей технологических процессов они приобретают ту или иную степень пористости. Оксидные керамики находят широкое применение в разных технических устройствах, например в виде фильтров, теплоизоляционных материалов, а в биомедицинской сфере - как протезы (заменители) костной ткани. Во многих случаях пористые керамики, костные ткани и геосреды представляют собой иерархически организованные многомасштабные системы. В полях внешних приложенных сил эволюция такой многомасштабной системы приводит к накоплению повреждений разных масштабов и формированию новых структур, т.е. усложнению самой системы. Одной из фундаментальных проблем при моделировании таких сред является построение определяющих уравнений, описывающих все аспекты механического поведения, включая деформационный отклик и особенности разрушения этих систем. Такие определяющие уравнения должны учитывать как характерные черты локального накопления повреждений, так и эволюцию порового пространства в различных условиях нагружения. Моделированию деформации пористых материалов посвящена обширная литература, например [1-8]. Ее можно разделить на несколько направлений: а) моделирование самой пористой структуры; б) определение эффективных упругих модулей пористых материалов; в) непосредственно численное моделирование процессов деформирования и разрушения пористых материалов на разных масштабных уровнях. Эффективным подходом к изучению подобных систем является развиваемая авторами эволюционная методология [9-12]. С математической точки зрения описывающая нестационарное поведение твердых деформируемых тел полная система дифференциальных уравнений механики сплошной среды вместе с определяющими соотношениями представляет собой систему нелинейных динамических уравнений. Как показано в работах [9-12], решения этих уравнений механики деформируемого твердого тела демонстрируют все основные особенности эволюции нелинейных динамических систем, включая особенности медленной динамики, бифуркации и смену сценариев эволюции. Показано также, что разрушение на заключительной неустойчивой стадии развивается как катастрофа в режиме с обострением. С физической точки зрения, нагружаемый материал рассматривается как нелинейная динамическая многомасштабная система, в которой явно задаются нелинейные положительные и отрицательные обратные связи, регулирующие взаимодействия между формирующимся в материале напряженно-деформированным состоянием и его откликом на нагружение (накопление повреждений на разных масштабах, деградацию прочностных свойств). Подобный подход позволяет эффективно описать эволюцию свойств пористых материалов на разных уровнях, включая локализованное накопление повреждений и неупругих деформаций микро- и мезомасштаба, деградацию прочности, формирование трещин разных масштабов и макроскопическое разрушение (образование магистральных трещин). Цель настоящей работы состоит в обсуждении подходов к моделированию пористой структуры реальных материалов (как геометрической задачи) и механического поведения (деформации, накопления повреждений и разрушения) мезообъ-емов пористых материалов как нелинейных систем. Геометрические модели пористой среды Следуя работам [1-3], рассмотрим две разные морфологии модельных пористых структур: перекрывающиеся сферические поры (ПСП) и перекрывающиеся сферические тела (ПСТ). В первом случае модельный образец пористой среды состоит из сплошного тела, включающего в себя сферические пустоты разного радиуса в случайно выбранных точках заданного объема, и иногда называется моделью «швейцарского сыра». Согласно работе [2], поры становятся проницаемыми при пористости p ~ 0,3, в то время как каркас остается связанным до пористости p s 0,97. Этой морфологии соответствует в реальных материалах система изолированных пор при низкой пористости. Подобная пористая структура в керамике может быть обусловлена процессами коалесценции и сфероидизации пор на последней стадии спекания. Во втором случае геометрическая модель строится заполнением объема сферическими частицами сплошного материала, которые располагаются в случайно выбранных точках объема с разными радиусами и могут перекрываться. Эта модель является отражением процесса получения керамического материала спеканием идеальных сферических частиц (порошка). Поровое пространство, как указано в [2], является проницаемым, если пористость превышает p s 0,03, а каркас остается связанным до p ~ 0,7. Однако отметим, что создать связанные модельные структуры оказывается затруднительно уже при пористости более 60 %. Для дальнейшего исследования методами численного моделирования были созданы модельные объекты указанных двух типов морфологии и для двух значений пористости 15 и 31,5 %. Изображения модельных структур в дискретном представлении расчетной сеткой из кубических ячеек (вокселов) представлены на рис. 1. Рис. 1. Изображения модельных объемов с разной морфологией пористой структуры и пористостью: а - 15 % ПСП, б - 15 % ПСТ, в - 31,5 % ПСП, г - 31,5 % ПСТ Насколько эти модельные структуры соответствуют реальности можно представить, если сравнить их с фотографиями поверхностей образцов пористой циркониевой керамики, представленными на рис. 2. Керамика получена спеканием наноструктурного порошка диоксида циркония [13]. Видно, что в реальности встречаются поры разной морфологии, причем они могут присутствовать в одном и том же материале. Это обусловлено разными факторами, например изменением режимов технологического процесса, разной формой отдельных частиц порошка и т.п. Тем больший интерес представляет изучение идеальных структур в численном эксперименте. Рис. 2. Примеры пористых структур в керамике из диоксида циркония Математическая постановка задачи в рамках эволюционного подхода В соответствии с эволюционной концепцией описания процессов деформирования и последующего разрушения материалов [9-10] полная система уравнений при лагранжевом описании движения сплошной среды включает в себя фундаментальные законы сохранения массы (1), импульса (2) и энергии (3), геометрические соотношения (4) и две группы определяющих соотношений (5) и (6) - (9). pV = PqVq ; (1) (2) Р1; = j; РЕ = a.ёij; (3) 2sj = vi,j + vj,i , 2сСij = vi,j - vj,i ; (4) a. =-Р5у + sj , Р = -K(е-eP), sij = 2^[(sy -sP)-1 (е-eP )5ij ]-+ j ik, e = sii; (5) SP = X , если f (a.) > 0; (6) J daj J f (a.) = -aP + j J2 - Y ; (7) g (a. ) = J2 -ЛР(2Y +aP) + const; (8) Y = Y0(1 - D). (9) Здесь p0, p - начальное и текущее значение плотности материала; V0, V - начальное и текущее значение объёма некоторой малой области материала; vi - компоненты вектора скорости; с. - компоненты тензора напряжений; s. - компоненты тензора скорости деформации; Е - внутренняя энергия единицы начального объёма; точка над символом означает материальную производную по времени, а запятая в нижнем индексе - частную производную по соответствующей пространственной координате. В первой группе определяющих уравнений (5) содержится разложение тензора напряжений на шаровую (давление Р) и девиаторную si. части, а также связь этих составляющих с параметрами деформированного состояния. При этом используется коротационная производная по времени Яуманна (a j =a j +aik со. -aj cb ik), учитывающая поворот элемента среды как целого. 5. означает символ Кронеккера, J2 = -2 s^s. - второй инвариант девиатора тензора напряжений, K и ц - модули всестороннего сжатия и сдвига соответственно, a -коэффициент внутреннего трения, Л - коэффициент дилатансии. Y обозначает сцепление (сдвиговую прочность среды при нулевом давлении), которое уменьшается от начального значения Y0 по мере накопления поврежденности D. Пластический множитель X в (6) определяется из условия удовлетворения напряжений условию текучести (7). Эволюционные определяющие уравнения первой группы (5) задают связь между напряжениями и деформациями в релаксационной (скоростной) форме: скорости приращения напряжений a. пропорциональны скоростям приращения упругих деформаций sE . Поскольку скорость полных деформаций s. складывается из скоростей упругих sE и неупругих составляющих sP , то скоростная форма записи определяющих уравнений описывает рост напряжений, вызываемый увеличением полных деформаций, и релаксацию напряжений за счет развития неупругих деформаций. При sP > s. в (5) a. < 0, идёт релаксация, реализуется отрицательная обратная связь, стабилизирующая деформационный процесс в состоянии динамического равновесия [10], определяемого текущим значением прочностных параметров. При sP < s^ напряжения растут ст^ > 0, в этом случае непосредственно реализуется положительная обратная связь, приводящая к катастрофе -макроскопическому разрушению, в том числе в режиме с обострением, если учтена деградация прочности [12]. Определяющую роль в переходе эволюции прочности среды в сверхбыстрый катастрофический режим будут играть особенности накопления средой локализованных неупругих деформаций и соответствующая деградация физико-механических параметров. Задачей эволюционных уравнений второй группы (6) - (9) является определение скоростей неупругих деформаций в уравнениях первой группы (5). В общем случае это кинетические уравнения, задающие скорости неупругой деформации и обеспечивающие релаксацию упругих напряжений в (5) [10]. В настоящей работе компоненты тензора скоростей неупругой деформации определены в соответствии с теорией пластического течения (6). В этом случае имеет место мгновенная релаксация напряжений на каждом временном шаге при численном решении уравнений. Предельная поверхность напряжений /(о^) записана в форме Мизеса -Шлейхера (7) и является обобщением критерия текучести Кулона - Мора. Она позволяет учесть зависимость сдвиговой прочности от давления, что в свою очередь определяет разную прочность материала в условиях сжатия и растяжения. За основу взята модель Друккера - Прагера - Николаевского с неассоциированным законом течения, позволяющая описывать процесс дилатансии как независимый от внутреннего трения. В этом случае пластический потенциал g(Oj) не совпадает с функцией текучести и записывается в виде (8) [14, 15]. Разрушение материла в развиваемом подходе описано как процесс обвальной деградации прочности материала до нуля при формировании трещин на сверхбыстрой катастрофической стадии эволюции напряженно-деформированного состояния, следующей за стадией предразрушения, на которой постепенно копятся повреждения. Следует отметить, что среда остаётся консолидированной при 7^0. Следовательно, остаются справедливыми все уравнения неупругого деформирования (1) - (8), нет также необходимости вводить в модель прочностные параметры, определяющие «предельное» состояние материала. По развиваемым в настоящей работе представлениям предельное состояние должно формироваться в нагружаемой среде по мере накопления в ней неупругих деформаций - повреждений. Необходимо только задать величину исходной прочности материала 70. Функцию деградации среды (поврежденность) D = D(t, цо, s) < 1 представим в виде зависимости от накапливаемой средой неупругой деформации s^ и вида напряжённого состояния (коэффициента Лоде - Надаи): D =f [H(Цст )(sтеК -sC)2 + (1 - H(Ц0 ))(sтеК -sT)2]dt tQ s?[ H (Цс )tC + (1 - H (Цс ))t*T] , s.=Sq. (1.01 + Ца )2, Ца= -1, (10) Здесь H(x) - функция Хевисайда, sтек = sHm + аб , sHm, - интенсивность тензора полных деформаций, sj, sT - начальные степени деформации на упругой стадии, по достижении которых в материале начинается накопление повреждений в областях сжатия и растяжения соответственно. Причём sT 0 в областях сжатия-сдвига. Скорости накопления повреждений ]/tj для локальных областей, где цс < 0 также существенно выше, чем 1/t*C в областях сжатия-сдвига (Цс > 0). Этот процесс дополнительно регулируется параметром е* = е*(Цс) в (10). Таким образом, отклик среды на вид напряжённого состояния (её текущая прочность) формируется в среде в процессе её нагружения. Следовательно, прочностные параметры будут деградировать существенно быстрее в тех областях среды, где коэффициент Лоде - Надаи пониженный, Цо < 0, т.е. преимущественный характер напряжённо-деформированного состояния определяется растяжениями-сдвигами. Этот отклик зависит также от истории нагружения. Если какая-то частица среды прежде испытывала, например, растяжение-сдвиг, а потом вид напряжённо-деформированного состояния изменился на сжатие-сдвиг, то процесс дальнейшего разрушения будет развиваться уже при пониженных прочностных параметрах, хотя и более медленно, в соответствии с кинетикой (10). е0* - параметр модели, t* имеет смысл характерного времени процесса. S1, S2, S3 - главные значения девиатора тензора напряжений. Любое твёрдое тело под нагрузкой рано или поздно разрушается, причём макроскопическое разрушение обычно представляется как мгновенное, в то время как устойчивая стадия неупругого деформирования (стадия предразрушения) может быть весьма продолжительной [12]. Однако для квазихрупких сред в стеснённых условиях деформирования за счёт дилатансионных процессов сверхбыстрая за-критическая стадия также может быть существенно затянута во времени [11, 12]. В модели эти процессы описаны зависимостью для меры накопления повреждений D и законом деградации. Приведённые ниже результаты моделирования процессов неупругого деформирования и разрушения хрупких и квазихрупких сред получены решением выписанной системы уравнений в трехмерной постановке методом конечных разностей по схеме второго порядка точности, подробно описанной в работе [16]. Результаты расчётов и их обсуждение Было проведено численное моделирование одноосного сжатия созданных модельных структур пористых материалов. Граничные условия, определяющие на-гружение образцов, были выбраны следующие. На гранях, перпендикулярных оси х, задавались на начальном периоде плавно растущие и затем постоянные значения скорости, направленные в сторону образца. Другие компоненты вектора скорости не ограничивались, что позволяло этим граням расширяться при сжатии. На остальных гранях, параллельных оси х, поддерживались условия свободной поверхности (отсутствия напряжений). Следует отметить, что такое плавное увеличение прилагаемой нагрузки, чтобы звуковая волна несколько раз прошла по образцу, позволяет в динамических задачах получить квазистатические решения [16]. Поскольку модель материала не учитывает скоростную чувствительность, то для экономии вычислительных затрат выбиралась такая скорость, чтобы при условии плавного наращивания скорости достигнуть заданной деформации за минимальное время. Расчеты были выполнены со следующими физико-механическими свойствами модельной сплошной среды: модуль сдвига G = 83 ГПа, модуль объёмной упругости K = 200 ГПа (модуль Юнга Е = 220 ГПа), плотность p = 6 г/см3, Y0 = 2200 МПа, a = 0,4, Л = 0,1, что соответствует частично стабилизированному диоксиду циркония в тетрагональной фазе. Значения для параметров выбранной кинетики накопления повреждений были подобраны следующие: е0* = 0,1, t*T = 10-6 , t*C = 10-4, eC = 0,005 , eT = 0,00125 . Размер расчетной области во всех случаях составил 9x9x9 мкм3. Расчётная сетка - равномерная 150x150x150 узлов. При этом в местах расположения пор задавался эффективный упругий материал, масса и упругие свойства которого были на несколько порядков меньше, чем соответствующие характеристики модельной сплошной среды. На рис. 3 показаны усредненные диаграммы нагружения для всех четырех образцов. Видно, что для образцов с меньшей пористостью диаграммы идут выше, что согласуется с экспериментальными данными. Для одной и той же пористости ниже расположены диаграммы для ПСТ-морфологии. Эту же информацию для чисто упругого поведения сообщают авторы статей [1, 2]. Боле того, согласно [1] зависимость усредненного модуля Юнга E от пористости p для разных морфоло-гий подчиняется степенной зависимости E = Ed (1 - p)m (11) где Ed - модуль Юнга плотного материала скелета, а показатель m = 2 для ПСП-структур и m = 4 для структур ПСТ. Этот показатель степени m авторы статьи называют показателем или индексом морфологии пористости. Полученные в наших расчетах значения этого показателя равны m =2,9 для структур ПСП и m = 4,0 - для ПСТ. Отличие для случая сферических пор может быть вызвано тем, что объемы не были достаточно представительными, а также тем, что в модели поры распределены квазиравномерно, в то время как в реальности наблюдается ярко выраженная кластеризация (см. рис. 2). Хотя другие авторы для тех же пористых структур приводят иные зависимости, например [2] E = f1 Г, E = f1 _Р_ Г Ed ^ 0,652) Ed ^ 0,818) для структур ПСТ и ПСП соответственно. 4 1000 800 600 200 ^■I1I■I■I1г^ 0;0 0,2 0,4 0,6 0,8 1,0 Деформация, % Рис. 3. Диаграммы нагружения для разных морфологий и значений пористости: 1 - 31,5 % ПСТ, 2 - 31,5 % ПСП, 3 - 15 % ПСТ, 4 - 15 % ПСП Экспериментальные исследования на пористых материалах из диоксида циркония с различной пористостью в диапазоне 10-75 % и разным соотношением размеров пор и размеров зерна керамики свидетельствуют о том, что зависимость модуля Юнга керамических образцов от пористости лучше аппроксимируется экспоненциальной функцией E = Ed exp(-b-p) [17]. Сравнение степенной и экспоненциальной аналитических зависимостей, а также полученные в наших расчетах данные изображены на рис. 4. Видно, что отличие экспоненциальной зависимости, полученной в результате аппроксимации экспериментальных данных, и степенной зависимости для ПСТ-морфологии незначительное в выбранном диапазоне пористости. Отличие экспоненциальной аппроксимации и степенной зависимости с показателем m = 4 для ПСТ-морфологии можно объяснить тем, что для большой и малой пористости морфология полученных образцов меньше соответствовала структуре перекрывающихся сферических тел. Полученные в наших расчетах данные также хорошо совпадают со степенной зависимостью для ПСТ-морфологии из [1]. E/Ed 0,9 0,8 0,7 0,6 0,5 0,4 0,3 0,2 0,1 о я V о Л с ч Л с 0 0,1 0,7 N \ \ \ V V \ к \ ч ч ♦ ' ч ч ч ч \3 чЧ. ч Ч 1 ч ч ч S 2 0,2 0,3 0,4 0,5 0,6 Пористость, усл. ед. 1400 200 0,1 0,4 1200 1000 800 600 400 А \\ А V N \ \\2 ....... ....... ....... V ч .4 \\ 0,15 0,2 0,25 0,3 0,35 Пористость, усл. ед. Рис. 4. Зависимость приведенного эффективного модуля упругости для пористых образцов: 1 - экспоненциальная аппроксимация для пористой керамики Zr02(Y203) [17]; 2 - степенная аппроксимация для ПСТ (m = 4); 5 - степенная аппроксимация для ПСП (m =2); + - данные расчетов для ПСТ; ♦ - данные расчетов для ПСП Рис. 5. Влияние пористости на предел прочности при сжатии пористых образцов: 1 -аппроксимация для керамики Zr02(Y203) с размером пор, соизмеримым со средним размером зерна (а = 1900exp(-5p)); 2 - то же с размером пор, значительно превышающим средний размер зерна (а = 2400 х х exp(-6p)); + - данные расчетов для ПСТ; ♦ - данные расчетов для ПСП Представляет интерес также оценка зависимости прочностных характеристик от пористости. Экспериментальные данные для керамик из диоксида циркония в работе [17] также были аппроксимированы экспоненциальной функцией с разными параметрами в зависимости от соизмеримости размеров пор и размеров зерна. На рис. 5 показаны эти аппроксимирующие зависимости и нанесены точки, соответствующие данным, полученным в наших расчетах. Видно, что расчетные данные для ПСП-морфологии хорошо ложатся на экспериментальные кривые. Отметим, что о количественном совпадении говорить не стоит, поскольку прочностные характеристики и параметры модели накопления повреждений не подобрались из условия строгого соответствия определенным значениям прочности для беспористых образцов. На рис. 6 представлены распределения поврежденности в разных образцах на момент, соответствующий последней точке расчетной диаграммы, т.е. непосредственно перед макроразрушением. Анализируя характер накопления повреждений в каркасе на мезоуровне, стоит отметить, что однозначных выводов о влиянии общей пористости и морфологии пор на форму и количество образующихся трещин выявить не удалось. Трещины преимущественно начинают образовываться около пор в местах высокой концентрации напряжений. Можно отметить, что в более тонком каркасе и большей кривизне пор ПСТ-структуры концентрация Рис. 6. Распределение поврежденности на поверхностных гранях исследуемых объемов для разных морфологий пор и пористости: а - 15 % ПСП, б - 15 % ПСТ, в - 31,5 % ПСП, г -31,5 % ПСТ напряжений выше. В связи с этим трещины начинают образовываться и прорастать через весь образец раньше. Заключение Проведенные расчеты показали, что для пористых керамических материалов, полученных спеканием наноструктурных порошков, хорошо проявили себя модельные структуры, созданные объединением перекрывающихся сферических тел. Такие геометрические модели лучше соответствуют реальным образцам пористой циркониевой керамики как по морфологии пор, так и по получаемой в моделировании механического поведения зависимости усредненного модуля Юнга от общей пористости. Хорошее качественное и количественное совпадение расчетов с экспериментальными данными свидетельствует о том, что учет двух масштабов пористости: явный учет крупных пор на мезоуровне и интегральный учет пор и накапливаемых повреждений на микроуровне вполне достаточен в рамках иерархического моделирования. Некоторое несовпадение расчетов с данными эксперимента для случая сферических пор следует объяснить тем, что в модельных пористых структурах поры расположены квазиравномерно, в то время как в реальности наблюдается их кластеризация. Подобное моделирование с различной степенью кластеризации пор будет целью дальнейших исследований. Применение эволюционного подхода, в рамках которого на стадии макроскопического предразрушения формируются области локализованных повреждений, позволило правильно описать образование областей повышенной поврежденности (трещин) на мезоуровне и влияние деградации в локальных областях материала на общий характер диаграммы нагружения, включая ее ниспадающую ветвь. Далее макроскопическое разрушение развивается как катастрофа в режиме с обострением. В рассматриваемой модели отклик среды на вид напряженно-деформированного состояния формируется в среде в ходе эволюции ее локальной прочности. Нет необходимости специально задавать различные прочности для условий сжатия и растяжения. Расчеты выполнены на высокопроизводительном кластере СКИФ-Cyberia Томского государственного университета.

Ключевые слова

хрупкие пористые материалы, напряженно-деформированное состояние, поврежденность, разрушение, влияние структуры, морфология пор, brittle porous materials, stress-strain state, damage, fracture, structure effect, pore morphology

Авторы

ФИООрганизацияДополнительноE-mail
Евтушенко Евгений ПавловичИнститут физики прочности и материаловедения Сибирского отделения Российской академии наукпрограммист лаборатории механики структурно-неоднородных средeugene@ispms.tsc.ru
Кульков Сергей НиколаевичТомский государственный университет; Институт физики прочности и материаловедения Сибирского отделения Российской академии наукдоктор физико-математических наук, профессор, заведующий кафедрой прочности и проектирования физико-технического факультетаkulkov@ispms.tsc.ru
Буякова Светлана ПетровнаТомский государственный университет; Институт физики прочности и материаловедения Сибирского отделения Российской академии наукдоктор технических наук, профессор кафедры прочности и проектирования физико-технического факультетаsbuyakova@ispms.tsc.ru
Макаров Павел ВасильевичТомский государственный университет; Институт физики прочности и материаловедения Сибирского отделения Российской академии наукдоктор физико-математических наук, профессор кафедры прочности и проектирования физико-технического факультетаpvm@ispms.tsc.ru
Еремин Михаил ОлеговичТомский государственный университетаспирант кафедры прочности и проектирования физико-технического факультетаbacardi@sibmail.com
Смолин Игорь ЮрьевичИнститут физики прочности и материаловедения Сибирского отделения Российской академии наук; Томский государственный университетдоктор физико-математических наук, доцент, ведущий научный сотрудник лаборатории механики структурно-неоднородных сред; профессор кафедры прочности и проектирования физико-технического факультетаsmolin@ispms.tsc.ru
Всего: 6

Ссылки

Wilkins M.L. Computer Simulation of Dynamic Phenomena. Berlin, Heidelberg: Springer-Verlag, 1999. 246 p.
Буякова С.П., Кульков С.Н., Масловский В.И. Структура, фазовый состав и механическое поведение керамики на основе диоксида циркония // Вестник Томского государственного университета. Бюллетень оперативной научной информации. 2003. № 13. С. 28-34.
Гарагаш И.А., Николаевский В.Н. Неассоциированные законы течения и локализации пластической деформации // Успехи механики. 1989. Т. 12. № 1. С. 131-183.
Макаров П.В., Смолин И.Ю., Стефанов Ю.П. и др. Нелинейная механика геоматериалов и геосред. Новосибирск: Академич. изд-во «Гео», 2007. 235 с.
Буякова С.П., Кульков С.Н. Фазовый состав и особенности формирования структуры в нанокристаллическом Zr02 // Российские нанотехнологии. 2007. Т. 2. № 1-2. С. 119- 132.
Макаров П.В., Еремин М.О. Модель разрушения хрупких и квазихрупких материалов и геосред // Физич. мезомех. 2013. Т. 16. № 1. С. 5-26.
Костандов Ю.А., Макаров П.В., Еремин М.О. и др. О разрушении хрупких тел с трещиной при сжатии // Прикладная механика. 2013. Т. 49. № 1. С. 124-132.
Макаров П.В. Математическая теория эволюции нагружаемых твердых тел и сред // Физич. мезомех. 2008. Т. 11. № 3. С. 19-35.
Макаров П.В. Эволюционная природа деструкции твердых тел и сред // Физич. мезо-мех. 2007. Т. 10. № 3. С. 23-38.
Смолин А.Ю., Коноваленко Иг.С., Кульков С.Н., Псахье С.Г. Моделирование разрушения хрупких пористых сред с различной внутренней структурой // Изв. вузов. Физика. 2006. Т. 49. № 3. С. 70-71.
Коноваленко Иг.С., Смолин А.Ю., Коростелев С.Ю., Псахье С.Г. О зависимости макроскопических упругих свойств пористых сред от параметров стохастического пространственного распределения пор // ЖТФ. 2009. Т. 79. Вып. 5. С. 155 - 158.
Скрипняк В.А., Скрипняк Е.Г., Козулин А.А. и др. Влияние поровой структуры хрупкой керамики на разрушение при динамическом нагружении // Известия Томского политехнического университета. 2009. Т. 315. № 2. С. 113 - 117.
Скрипняк Е.Г., Скрипняк В.А., Кульков С.С. и др. Моделирование механического поведения керамических композитов с трансформационно-упрочненной матрицей при динамических воздействиях // Вестник Томского государственного университета. Математика и механика.
Кондауров В.И. Механика и термодинамика насыщенной пористой среды. М.: МФТИ, 2007. 310 с.
Roberts A. Garboczi E. Elastic properties of model porous ceramics // J. Am. Ceram. Soc. 2000. V. 83. № 12. P. 3041 - 3048.
Torquato S. Random Heterogeneous media: microstructure and improved bounds on elastic properties // Appl. Mech. Rev. 1991. V. 44. P. 37 - 76.
Bruno G., Efremov A.M., Levandovskyi A.N., Clausen B. Connecting the macro- and microstrain responses in technical porous ceramics: modeling and experimental validations // J. Mater. Sci. 2011. V. 46. P. 161 - 173.
 Численное моделирование механического поведения модельных хрупких пористых материалов на мезоуровне | Вестн Том. гос. ун-та. Математика и механика. 2013. №  5(25).

Численное моделирование механического поведения модельных хрупких пористых материалов на мезоуровне | Вестн Том. гос. ун-та. Математика и механика. 2013. № 5(25).