Numerical modeling of nonlinear dynamic and static compression of the metal mesh

Abstract


Simulation of elastoplastic compression of a symmetric fragment of a bundle of woven steel grids was carried out. The simulation was carried out under static and dynamic loading modes in ANSYS and ANSYS LS-DYNA computing systems. A porous mesh package is formed by overlaying the layers on top of each other while maintaining the directions of wires. Such a packet has a quasiperiodic structure; therefore, a certain symmetric fragment can be distinguished. The compression was carried out by a pair of absolutely rigid plates moving symmetrically towards each other. In the calculations, a multilinear plasticity model with isotropic hardening was used. The diagram of static deformation of the material was used which was obtained experimentally. The calculations were carried out according to the algorithm of a perfect symmetric contact of bodies without friction and with friction. The dynamic loading of a fragment of the mesh packet was carried out at a constant speed. The characteristics of the pulse and the loading rate correspond to those observed in previous experimental studies. The dynamic strain diagram was assumed to be similar to a static one with an increased yield strength. Calculations showed that in all loading modes there is a high level of internal stresses, and complex inhomogeneous stress-strain state. We investigated two factors that could cause differences in the behavior of the curves of deformation of a porous medium - the finite length of the acting pulse and the differences in the dynamic and static compression diagrams of the initial mesh material. The main influence on the dynamic behavior of a fragment of a porous packet of a steel mesh is exerted by the dynamic properties of the wire material. The strain curves of the porous fragment qualitatively change in accordance with the behavior observed in static and dynamic experiments. The final duration of the loading pulse and the friction between the wires for this type of mesh do not have a significant effect. The numerical dependences of the relative area of the normal and lateral through sections of the symmetric fragment of the grid packet on the compression deformation are obtained.

Full Text

Введение Для снижения импульсных и вибрационных нагрузок в конструкциях применяют пористые элементы в виде гранулированных слоев, проволочных решеток, экранов, сеток, перфорированных перегородок [1-19]. В частности, для защиты силовых корпусов взрывных камер от осколочного повреждения в настоящее время применяются многослойные металлические сетки тканого плетения [20, 21]. Защитный пакет сеток обычно формируется путем наложения слоев друг на друга с сохранением направлений проволок, поэтому многослойный пакет можно считать высокопористым деформируемым элементом конструкций, обладающим ортотропными свойствами [21]. В [22-24] приведены некоторые результаты экспериментальных исследований деформационных свойств многослойных пакетов сеток на сжатие по направлению нормали к слоям сеток. Показано, что при интенсивных нагрузках диаграммы деформирования носят нелинейный и необратимый характер. Экспериментальные исследования [22] показали существенные отличия кривых деформирования пакетов сеток в статическом и динамическом режимах сжатия (рис. 1). Динамические испытания проведены на газовых пушках в системе разрезного стержня Гопкинсона по методикам Кольского [25-27]. Исследования были проведены для сетки с диаметром проволоки 0,5 мм и шагом плетения 2,5 мм, изготовленной из стали 3. Скорость деформации в экспериментах изменялась в пределах 1500-3500 с-1. Деформирование образца проходило в несколько циклов, связанных с распространением волн сжатия и растяжения в мерных стержнях. На рис. 1 показаны диаграммы деформирования образцов при сжатии по нормали к слоям сеток, включающие участки активного нагружения и разгрузки. Напряжение рассчитывалось как сила, деленная на начальную площадь образца в плане (всего пакета), а деформация как изменение толщины пакета к начальной толщине образца. На рис. 1 введены обозначения: 1 - образец из 10 слоев; 2 - образец из 20 слоев; 3 - результаты статических испытаний [22]. Статические испытания на сжатие по нормали проводились на сервогидравлической установке Zwick/Roell Amsler HA 100 для пакетов сеток также из 10 и 20 слоев. Видно, что деформирование при активном нагружении носит нелинейный характер, в то время как разгрузочные ветви близки к прямым линиям. Кривые деформирования для образцов, состоящих из более чем 10 слоев, слабо зависят от количества слоев и формы образца. Полученные при активном нагружении динамические кривые деформирования располагаются значительно выше статической кривой. Одной из целей настоящей статьи является выявление факторов, обусловливающих различие кривых деформирования в разных режимах упругопластического статического и динамического сжатия симметричных фрагментов многослойного пакета плетеных стальных сеток. Исследуются два фактора, которые могли бы вызвать различия в поведении кривых, конечная длина воздействующего импульса при динамическом нагружении и различия динамической и статической диаграмм сжатия исходного материала сетки (сталь 3). В ранее проведенных аналогичных исследованиях при сжатии гранулированных слоев из свинцовых шариков [26] оба фактора оказывали заметное влияние на кривые деформирования в различных режимах нагружения. Другой целью статьи является определение зависимости изменения характеристик пористости и проницаемости деформируемых пакетов сеток в зависимости от степени обжатия для учета в задачах взаимодействия ударных волн с проницаемыми элементами конструкций. Рис. 1. Экспериментальные диаграммы сжатия плетеной металлической сетки при статическом (3) и динамическом (1, 2) режимах нагружения: 1 - 10 слоев; 2 - 20 слоев Fig. 1. Experimental compression diagrams of a woven metal mesh under static (3) and dynamic (1, 2) loading conditions: 1 - 10 layers; 2 - 20 layers Результаты численного моделирования статического и динамического сжатия Пористый пакет сеток формируется, как правило, путем наложения слоев с сохранением направлений проволок. Такой пакет имеет квазипериодическую структуру, и можно выделить некоторый симметричный фрагмент. Для моделирования была выбрана типовая симметричная ячейка одного слоя плетеной сетки (рис. 2). Численное моделирование процесса статического деформирования плетеной металлической сетки при ее сжатии перпендикулярно плоскости слоя проводилось в вычислительной системе ANSYS v. 17.2 методом конечных элементов (лицензия Academic Research, Customer #623640). Расчетная область состоит из четырех трехмерных цилиндрических тел (проволок). В силу симметрии рассматривается половина поперечного сечения проволок. На торцах проволок также выполняются условия симметрии. В начальный момент времени напряжения и деформации отсутствуют. В данной задаче использовался алгоритм расчета идеального симметричного контакта тел без трения и с трением, когда в каждой контактной области используются две контактные пары. Коэффициент трения принимался равным 0,3. Сжатие проводилось парой абсолютно жестких пластин (плоскостей симметрии), движущихся в направлении оси Z навстречу друг другу (кинематическая схема нагружения). Сжатие проводилось до 0,7 начальной толщины фрагмента (далее возникала неустойчивость при моделировании). Рис. 2. Расчетная модель Fig. 2. Numerical model Геометрическая модель проволоки получена вытягиванием полукруглого сечения вдоль оси проволоки. Ось состоит из двух сопряженных дуг окружности. Для построения конечно-элементной модели использовался 20-узловой КЭ второго порядка SOLID185 с сокращенным (2×2×2) интегрированием. Модель каждой из проволок состоит из 6144 КЭ, всего 24576 КЭ. Общее количество неизвестных составило 321660. Задача решалась на последовательности сеток до тех пор, пока различие в вычислении напряжений не превосходило 5 %. В расчетах использовалась диаграмма деформирования на растяжение стали 3 (рис. 3), которая была получена авторами статьи экспериментально при квазистатическом режиме нагружения на экспериментальной установке Zwick Amsler HA 100. Скорость деформации 0,01 1/с. Модуль упругости равен 200 ГПа, коэффициент Пуассона - 0,3, предел текучести - 200 МПа, предел прочности - 482 МПа, максимальное относительное удлинение - 20 %. В квазистатическом расчете описания поведения материала использовалась мультилинейная модель пластичности с изотропным упрочнением. Диаграмма деформирования материала разбивается на 10 линейных участков различного наклона. Первый участок связывается с упругим поведением материала. При превышении в расчете максимальной заданной деформации (20 %) диаграмма автоматически продолжается с линейным модулем упрочнения, равным 1ГПа, продолжение диаграммы с нулевым модулем упрочнения не приводит к заметному изменению результатов. В ANSYS вычислительный процесс решения нелинейных задач квазистатического деформирования конструкций разбивается на шаги, чтобы процесс нагружения передавался наилучшим образом. Это внутренний параметр программ. В рассматриваемой задаче достаточно было разбить процесс нагружения на 200 шагов для получения требуемой точности. Рис. 3. Диаграмма деформирования материала Fig. 3. Material deformation diagram На рис. 4 приведено распределение эквивалентных пластических деформаций, то есть значений второго инварианта тензора пластических деформаций, после снятия нагрузки. Наблюдается развитое пластическое течение во всех проволоках, максимальные значения эквивалентных пластических деформаций в некоторых точках узла плетения достигают 1,6. Имеет место высокая неоднородность деформированного состояния, особенно в окрестности сжатого узла плетения. Сравнение результатов моделирования с учетом и без учета трения на контактных поверхностях показало незначительное влияние сил трения. На рис. 5 приведено распределение эквивалентных напряжений по Мизесу в конце обжатия, когда смещение нагружающих пластин максимально. Эквивалентные напряжения во всех проволоках, за исключением небольших локализованных зон, превышают 440 МПа, что свидетельствует о высоком уровне внутренних напряжений. Диапазон их изменения от 100 до 480 МПа показывает неоднородность распределения. Численная диаграмма деформирования пористого фрагмента пакета плетеных сеток вдоль оси Z в квазистатическом режиме нагружения показана линией 1 на рис. 6. Диаграмма связывает среднее вертикальное напряжение с деформацией фрагмента. Вычисление средних напряжений проводилось по формуле (1) где скобки обозначают среднее значение; Vc.cell - текущий объем деформируемой области вместе с пустотами, равный 2(d-Uz)l2; Uz - текущее смещение горизонтальной плоскости; - напряжения в конечном элементе; - текущий объем конечного элемента; l - шаг плетения (l = 2,5 мм); d - диаметр проволоки (d = 0,5 мм). Характерный перегиб на диаграмме при деформации 50 % связан с тем, что произошло локальное сжатие узлов плетеной металлической сетки и далее происходит сжатие всех проволок по всей их длине. Рис. 4. Распределение эквивалентных пластических деформаций Fig. 4. Distribution of equivalent plastic strains Рис. 5. Распределение эквивалентных напряжений по Мизесу Fig. 5. Distribution of equivalent von Mises stresses Так как пакеты сеток являются газопроницаемыми элементами, большой интерес вызывают характеристики их проницаемости, связанные с изменением площадей проходных сечений вследствие их деформации. На рис. 7 показана геометрия проходного сечения (перпендикулярного оси Z) до начала сжатия (левая часть рисунка) и после его окончания. Относительная площадь проходного сечения изменяется вследствие пластического деформирования проволок от 0,64 до 0,4. Относительная площадь боковых проходных сечений изменяется от 0,41 до 0. Рис. 6. Кривые деформирования статического и динамического сжатия симметричного фрагмента плетеной сетки в вертикальном направлении: 1 - статика; 2 -динамика; 3 - динамическое нагружение с использованием динамической диаграммы стали Fig. 6. Strain curves of static and dynamic compression of a symmetrical fragment of a woven mesh in the vertical direction: 1 - static diagram; 2 - dynamic diagram; 3 - dynamic loading using a dynamic steel diagram Рис. 7. Формы вертикального проходного сечения до начала (слева) и в конце процесса сжатия (справа) Fig. 7. Shapes of the vertical bore before the start (left) and at the end of the compression process (right) Вычисление относительной площади минимального проходного сечения проводилось по формуле (2) где - свободная площадь в плане, которая рассчитывается как , где - полная площадь симметричного фрагмента в плане; - площадь деформированных проволок в плане. На рис. 8 показано изменение относительной площади минимальных проходных сечений от вертикальной средней по объему деформации сжатия одного слоя плетеной металлической сетки (1 - изменение сверху; 2 - изменение сбоку). Резкое изменение относительной площади при деформации 50 %, как и излом на диаграмме деформирования (см. рис. 6), связано с увеличением площади контакта сжимаемых проволок. Численное моделирование динамического деформирования фрагмента плетеной сетки выполнено в программном обеспечении ANSYS LS-DYNA (явная схема интегрирования по времени, лицензия Customer #244793). Нагружение аналогичного фрагмента пакета сетки, рассмотренного выше (см. рис. 2), проводилось с постоянной скоростью 1,25 м/с в течение времени 0,6 мс. Характеристики импульса и скорости нагружения соответствуют экспериментально наблюдаемым и измеряемым воздействиям в исследованиях [22, 27]. При численном моделировании не воспроизводились циклы разгружения в динамическом эксперименте, связанные с пробегом и отражением волн в мерных стержнях, нагружение было непрерывным. При моделировании использовалась модель Piecewise Linear Plasticity (при этом зависимость от скорости деформирования игнорировалась). То есть в динамическом режиме использовалась та же модель деформирования материала, что и при моделировании в статическом режиме. Используемый численный метод второго порядка точности в лагранжевых переменных условно устойчивый, ограничения на шаг по времени определяются условием Куранта - самая быстрая упругая волна не может пробегать расстояние, большее, чем размер разностной ячейки (элемента). Для обеспечения устойчивости и сходимости решения размер временного шага при явной схеме интегрирования ограничен условием Куранта: dt ≤ f×(h/c), где dt - шаг по времени, f - коэффициент запаса к шагу по времени (0,9 по умолчанию), h - характерный размер конечного элемента, с - локальная скорость звука материала конечного элемента. Для данной задачи это время порядка 0,01 мкс. Вычислительный пакет автоматически генерирует устойчивый шаг по времени в зависимости от степени изменения (смятия) расчетных ячеек. Также вычислительный пакет генерирует оптимальное значение коэффициента искусственной вязкости, которая сглаживает нефизические осцилляции, характерные для численных схем второго порядка точности. При выполнении данных условий в расчетах не наблюдалось неустойчивостей и численных осцилляций высокой амплитуды, за исключением стадии расчета, когда формировались очень сильные искажения форм и размеров элементов вследствие больших деформаций. Для оценки практической сходимости численного решения расчеты в вычислительном пакете ANSYS LS-DYNA проводились с использованием конечно-элементных моделей размерами порядка 12000, 25000, 50000 элементов. Результаты расчетов с количеством 25000 и 50000 элементов фактически совпали. Распределение напряжений и деформаций в целом соответствует распределению при статическом нагружении за небольшим исключением областей в окрестности узлов плетения. На рис. 6 показаны зависимости среднего по расчетной области вертикального (по оси сжатия Z) напряжения от средней вертикальной деформации, статическая кривая 1, динамическая - 2. Эти кривые практически совпадают. Таким образом, для одной и той же статической диаграммы сжатия стали 3 (рис. 3) влияние режима нагружения оказывается незначительным. Материал успевает необратимо деформироваться в течение действия нагружающего импульса. Рис. 8. Относительная площадь проходных сечений в зависимости от деформации обжатия: 1 - изменение сверху; 2 - изменение сбоку Fig. 8. The relative area of the bore sections, depending on the compression deformation: 1 -change from above; 2 - change from the side Другим фактором, который может повлиять на изменение характера деформирования фрагмента является отличие динамической диаграммы материала (сталь 3) от статической. В работах [28, 29] при обобщении известных экспериментальных данных по зависимостям характеристик деформирования от скорости нагружения показано, что для многих металлов и их сплавов наибольшую зависимость от изменения скорости деформирования от 10-2 до 104 1/с проявляет предел текучести, который может измениться в несколько раз. В частности, для стали 3 предел текучести при скорости деформации 103 1/с достигает 500 МПа [29]. Другие характеристики ведут себя относительно стабильно. С целью учета этого фактора проведены численные исследования сжатия фрагмента стальных плетеных сеток с диаграммой деформирования, у которой предел текучести 500 МПа. На кривой деформирования стали (см. рис. 3) изменяется только значение предела текучести, угол наклона упругого участка и поведение участков кривой деформирования выше предела текучести не изменяются. На рис. 6 цифрой 3 показана полученная в результате численных расчетов динамическая диаграмма деформирования фрагмента пакета сеток как пористого элемента конструкции, то есть зависимость среднего вертикального напряжения (1) от средней деформации фрагмента. Наблюдается сильная зависимость кривых деформирования от значений предела текучести и, следовательно, от скорости деформирования. С ростом предела текучести материала кривые деформирования фрагмента качественно изменяются в соответствии с поведением, наблюдаемым в статических и динамических экспериментах. Таким образом, основным фактором, который обусловливает отличие динамических диаграмм деформирования от статических для пористых пакетов плетеных сеток, является различие динамической и статической диаграмм материала, из которого изготовлена сетка. Заключение В результате численных исследований статического и динамического деформирования фрагмента плетеной сетки как пористой среды были получены кривые упругопластического сжатия. Исследовались два фактора, которые могли бы вызвать различия в поведении кривых - конечная длина воздействующего импульса и различия динамической и статической диаграммы сжатия исходного материала сетки (сталь 3). Показано, что главное влияние на динамическое поведение фрагмента пористого пакета стальной сетки оказывает динамическая диаграмма деформирования материала проволок. Конечная длительность нагружающего импульса для данного вида сеток существенного влияния не оказывает. Трение при моделировании контактного взаимодействия проволок сетки при сжатии симметричного фрагмента существенного влияния на результаты также не оказывает. Получены численные зависимости относительной площади нормального и бокового проходного сечения симметричного фрагмента пакета сеток от деформации обжатия. Эти зависимости могут быть использованы при моделировании взаимодействия проницаемых пакетов плетеных сеток с интенсивными ударными волнами, вызывающими большие деформации сеток.

About the authors

A V Kochetkov

National Research Nizhniy Novgorod State University N.I. Lobachevsky

I A Modin

National Research Nizhniy Novgorod State University N.I. Lobachevsky

N V Leontev

National Research Nizhniy Novgorod State University N.I. Lobachevsky

I A Turigina

National Research Nizhniy Novgorod State University N.I. Lobachevsky

E Y Poverennov

Afrikantov OKB Mechanical Engineering

References

  1. Hermann W. Constitutive equation for the dynamic compaction of ductile porous materials // J. Appl.Phys. - 1969. - Vol. 40. - No. 6.
  2. Гельфанд Б.Е., Сильников М.В. Фугасные эффекты взрывов. - СПб.: Полигон, 2002. - 272 с.
  3. Гельфанд Б.Е., Губанов А.В., Тимофеев Е.И. Взаимодействие воздушных ударных волн с пористым экраном // Изв. АН СССР. МЖГ. - 1983. - № 4. - С. 79-84.
  4. Передача ударно-волновой нагрузки насыпными средами / Б.Е. Гельфанд, С.П. Медведев, А.Н. Поленов, С.М. Фролов // Прикладная механика и техническая физика. - 1988. - № 2. - С. 115-121.
  5. Shock waves attenuation by granular filters / A. Britan, G. Ben-Dor, O. Igra, H. Shapiro // International Journal of Multiphase Flow. - 2001. - Vol. 27 (4). - P. 617-634.
  6. Britan А., Ben-Dor G. Shock tube study of the dynamical behavior of granular materials // International Journal of Multiphase Flow. - 2006. - Vol. 32 (5). - P. 623-642.
  7. Модин И.А. Упругопластическое деформирование высокопористых элементов конструкций при квазистатическом и импульсном нагружениях. - Н. Новгород: Изд-во Нижегород. гос. ун-та им. Н.И. Лобачевского, 2018. - 87 с.
  8. Simulation the dynamics of a composite cylindrical shell with a gas-permeable layer under the internal impulse loading / A.Yu. Konstantinov, A.V. Kochetkov, S.V. Krylov, I.V. Smirnov // Materials Physics and Mechanics. - 2016. - Vol. 28. - No 1/2. - P. 39-42.
  9. Кругликов Б.С., Кутушев А.Г. Ослабление ударных волн экранирующими решетками // ФГВ. - 1998. - Т. 24, № 1. - С. 115-118.
  10. Stolz A., Ruiz-Ripoll M.L. Experimental and computational characterization of dynamic loading and structural resistance of tunnels in blast scenarios // Fire Technology. - 2015. - P. 24. doi: 10.1007/s10694-015-0496-8
  11. An inverse estimation of high strain rate properties of composite material constituents / S. Chacko, A. Jones, R. Brooks, M.J. Lidgett // 20th International Conference on Composite Materials Copenhagen. 19-24th July 2015.
  12. Splichal J., Pistek A., Hlinka J. Dynamic tests of composite panels of an aircraft wing // Progress in Aerospace Sciences. - 06/2015. doi: 10.1016/j.paerosci.2015.05.005
  13. Cadoni E., Forni D. Strain rate effects on reinforcing steels in tension // EPJ Web of Conferences. 04. 2015. DOI: http://dx.doi.org/10.1051/epjconf/20159401004
  14. Zhu H., Pierron F. Exploration of Saint-Venant’s principle in inertial high strain rate testing of materials // Experimental Mechanics. - 07/2015. doi: 10.1007/s11340-015-0078-1
  15. Hu D., Meng K., Jiang H. Experimental investigation of dynamic properties of AerMet 100 steel // Procedia Engineering - 12/2015; 99:1459-1464. doi: 10.1016/j.proeng.2014.12.685
  16. Альтшулер Л.В., Кругликов Б.С. Затухание сильных ударных волн в двухфазных и гетерогенных средах // ПМТФ. - 1984. - № 5. - С. 24-29.
  17. Ruan H.H., Gao Z.Y., Yu T.X. Crushing of thin-walled spheres and sphere arrays // Int. J. Mech. Sci. - 2006. - No. 48. - Р. 117-133.
  18. The Effect of voids and inclusions on wave propagation in granular materials / M.H. Sadd, A. Shukla, H. Mei, C.Y. Zhu // Micromechanics and Inhomogeneity. - 1989. - P. 367-383.
  19. Shukla A., Damania C. Experimental investigation of wave velocity and dynamic contact stresses in an assembly of disks // Experimental Mechanics. - September 1987. - Vol. 27. - Iss. 3. - Р. 268-281.
  20. Численная модель деформирования противоосколочной сетки при взрывном нагружении / А.И. Абакумов [и др.] // Тр. ВНИИЭФ. Математическое моделирование физических процессов. - 2006. - № 10. - С. 16-30.
  21. Взрывное нагружение деформируемых газопроницаемых осесимметричных элементов конструкций / Е.Г. Глазова, А.Ю. Константинов, А.В. Кочетков, С.В. Крылов // ПМТФ. - 2016. - № 5. - С. 119-126. DOI: 10.15372 / PMTF20160513.
  22. Экспериментальное исследование деформационных свойств пакетов плетеных металлических сеток при динамическом и квазистатическом нагружении / А.М. Брагов, Д.В. Жегалов, А.Ю. Константинов, А.В. Кочетков, И.А. Модин, А.О. Савихин // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2016. - № 3. - С. 252-262. doi: 10.15593/perm.mech/2016.3.17
  23. Исследование деформационных и прочностных свойств металлических плетеных сеток / А.В. Кочетков, Н.В. Леонтьев, И.А. Модин, А.О. Савихин // Вестн. Том. гос. ун-та. Математика и механика. - 2018. - № 52. - С. 53-62. doi: 10.17223/19988621/52/6
  24. Xiao Lijun, Song Weidong. Additively-manufactured functionally graded Ti-6Al-4V lattice structures with high strength under static and dynamic loading // International Journal of Impact Engineering. - 2018. - No. 111. - Р. 255-272.
  25. Bragov A.M., Lomunov A.K. Methodological aspects of studying dynamic material properties using the Kolsky method // Int. J. of Impact Engineering. - 1995. - Vol. 16 (2). - P. 321-330.
  26. Experimental study of deformation properties of a bulk layer from plumbum balls under dynamic and quasistatic loading / A.M. Bragov, A.U. Konstantinov, A.V. Kochetkov, I.A. Modin, A.O. Savikhin // PNRPU Mechanics Bulletin. - 2017. - No. 4. - P. 16-27. doi: 10.15593/perm.mech/2017.4.02
  27. Брагов А.М., Ломунов А.К. Использование метода Кольского для исследования процессов высокоскоростного деформирования материалов различной физической природы. - Н. Новгород: Изд-во Нижегород. гос. ун-та, 2017. - 148 с.
  28. Мержиевский Л.А., Палецкий А.В. Расчет диаграмм динамического деформирования материалов и сплавов // Физическая мезомеханика. - 2001. - Т. 4, № 3. - С. 85-96.
  29. Экспериментальное исследование и математическое моделирование поведения сталей марок Ст.3, 20Х13 и 08Х18Н10Т в широких диапазонах скоростей деформаций и температур / А.М. Брагов [и др.] // Прикладная механика и техническая физика. - 2015. - Т. 56, № 6. - С. 51-58.

Statistics

Views

Abstract - 417

PDF (Russian) - 293

Cited-By


PlumX


Copyright (c) 2019 Kochetkov A.V., Modin I.A., Leontev N.V., Turigina I.A., Poverennov E.Y.

Creative Commons License
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies