DETERMINATION OF THE FUNDAMENTAL FREQUENCY OF A SANDWICH PANEL AS CONSTRUCTIONAL PART OF SPACECRAFT BODY


Cite item

Full Text

Abstract

Sandwich panels are widely used constructional parts of modern spacecraft bodies and have various edge-fixing conditions. In designing of such panels there is always a need for an analytical formula to determine with reasonable accuracy the fundamental frequency which is a convenient assessment of sandwich construction efficiency that takes into account the mutual influence of its stiffness and mass. The article observes fundamental frequency determination for rectangular sandwich panel with all edges rigidly clamped. The panel has symmetrical layer package, consists of two identical face-sheets and orthotropic core. The approach uses Galerkin’s method to solve motion equations and the solution of a bending problem of the beam with rigidly clamped edges as a function that approximate panel shape. This reduces the problem of determination of the dimensionless frequency parameter which depends on the panel geometry, inertial and elastic characteristics. The frequency parameter is calculated both including and excluding rotation inertia. In the latter case, an analytical formula for fundamental frequency is the result of this approach. Influence of geometrical, elastic and inertial parameters of the sandwich panel to its fundamental frequency have been analyzed by using the analytical formula. Comparison of calculation data with finite-element package modal solution shows that obtained formula gives high accuracy and low usage of computational resources, which is especially useful in designing process where restrictions are imposed on the first frequency of sandwich panels.

Full Text

Введение. Трехслойные панели широко применяются в составе корпусов современных космических аппаратов и обладают различными способами закрепления краев (рис. 1). Рис. 1. Трехслойные несущие панели Динамический анализ является важной частью расчетов трехслойных панелей. Задачи определения частот и форм колебаний трехслойных пластин привлекали внимание многих исследователей. Впервые уравнения движения трехслойной пластины были получены и решены в статьях [1-4]. В последовавших затем исследованиях были предложены разнообразные динамические модели трехслойных пластин, на основе которых были решены многочисленные задачи колебаний этих конструкций. Разработанные модели и методы определения частот и форм колебаний трехслойных пластин нашли свое отражение в обзорах [5-9]. Результаты динамических расчетов трехслойных пластин представлены в работах [10-16]. Трехслойные пластины как несущие конструкции обладают граничными условиями, соответствующими жесткой заделке. Для такой пластины получить решение динамической задачи весьма сложно. Поэтому большинство исследований, посвященных колебаниям трехслойных пластин, выполнено для граничных условий, соответствующих шарнирному закреплению сторон. Среди задач динамики трехслойных пластин особое место занимает задача определения основной частоты колебаний. Дело в том, что величина частоты зависит от отношения изгибной жесткости пластины к массе единицы поверхности. При этом увеличение изгибной жесткости пластины всегда сопровождается увеличением массы единицы поверхности. Поэтому основная частота колебаний является удобной оценкой эффективности трехслойной пластины, учитывающей взаимное влияние жесткости и массы конструкции. Высокая эффективность трехслойных пластин обусловлена именно тем, что их структура, как никакая иная, позволяет обеспечить значительное увеличение изгибной жесткости при небольшом увеличении массы единицы поверхности. При реальном проектировании защемленных по контуру пластин всегда существует потребность в аналитической формуле, позволяющей с приемлемой точностью определять основную частоту колебаний. Постановка задачи. Для описания движения пластины в статье использованы уравнения сдвиговой теории пластин, основы которой были заложены в работах [17] и [18] и которая получила широкое распространение при анализе слоистых пластин. Рассмотрим трехслойную пластину, состоящую из одинаковых композитных несущих слоев и ортотропного заполнителя. Срединную плоскость пластины отнесем к системе координат x, y, z. Обозначим через a и b размеры пластины в плане, а через t и δ - суммарную толщину несущих слоев и заполнителя (рис. 2). Рис. 2. Трехслойная пластина Для описания колебаний трехслойной пластины воспользуемся системой уравнений, которая включает уравнения движения: (1) физические соотношения: (2) геометрические соотношения: (3) В уравнениях (1)-(3) τ - время; Qx, Qy - перерезывающие силы; Mx, My - изгибающие моменты; χx, χy - кривизны срединной плоскости; χxy - кручение срединной плоскости; ψx, ψy - осредненные по толщине деформации поперечного сдвига; w - прогиб пластины; θx, θy - углы поворота нормали к срединной плоскости; D11, D12, D21, D22, D33 - изгибные жесткости пластины (D12 = D21); Bρ - инерциальный параметр; Dρ - параметр, обусловленный инерцией поворота нормали к срединной плоскости. Получим уравнение движения, содержащие в качестве неизвестных прогиб w и углы поворота θx, θy. Последовательно подставляя (2) и (3) в (1), получим (4) Рассматривая свободные колебания пластины, представим решение уравнений (4) в виде (5) где ω - круговая частота колебаний; . Функции w(x,y), θx(x,y), θy(x,y) определяют форму пластины при колебаниях. Подставляя (5) в (4), получим следующую систему однородных дифференциальных уравнений: (6) здесь и далее w = w(x,y), θx = θx(x,y), θy = θy(x,y). Система уравнений имеет в совокупности шестой порядок по переменным x, y, поэтому на каждом из закрепленных краев пластины надо поставить по три граничных условия, т. е. (7) Определим жесткостные и инерциальные параметры трехслойной пластины. Будем в дальнейшем полагать, что несущие слои являются однородными и ортотропными. Тогда будем иметь (mn = 11, 12, 22, 33), (8) Здесь (9) В уравнениях (8), (9) Amnt, Amnδ (mn = 11, 12, 22, 33) - коэффициенты жесткости несущего слоя и заполнителя; Gxzt, Gyzt, Gxzδ, Gyzδ - трансверсальные модули сдвига несущего слоя и заполнителя; ρt, ρδ - плотности материалов несущего слоя и заполнителя. Для ортотропного материала A11s = , A22s = , A12s = = , A33s = Gxys, где s = t,δ; = Exs/(1 - νxysνyxs), = Eys/(1 - νxysνyxs); Exs, Eys - модули упругости в направлениях x и y соответственно; νxys, νyxs - коэффициенты Пуассона; Gxys - модуль сдвига в плоскости xy. Если материал несущих слоев или заполнителя изотропен, то тогда A11s = A22s = , A12s = , A33s = = Gs = 0,5(1 - νs), Gxzs = Gyzs = Gs, где s = t,δ; = = Es/(1 - ); Es - модуль упругости; Gs - модуль сдвига; νs - коэффициент Пуассона. Выбор функций, аппроксимирующих прогиб и углы поворота. Для решения однородной краевой задачи (6), (7) необходимо выбрать функции, которые аппроксимируют форму пластины, соответствующую основному тону колебаний. В качестве таких функций можно принять решение задачи об изгибе балки с жестко закрепленными концами под действием равномерно распределенной нагрузки (рис. 3). Напряженно-деформированное состояние балки в рамках модели, учитывающей деформации поперечного сдвига, описывается следующими уравнениями: (10) где p - действующая нагрузка; M, Q - изгибающий момент и перерезывающая сила; χ - кривизна продольной оси балки; ψ - деформация сдвига; w - прогиб балки; θ - угол поворота поперечного сечения; D, K - изгибная и сдвиговая жесткости. Рис. 3. Балка под действием равномерной нагрузки Из решения уравнений (10) для угла поворота θ и прогиба w можно получить следующие выражения: (11) где C1, C2, C3, C4 - произвольные постоянные. Для их определения воспользуемся граничными условиями на концах балки при (12) где l - длина балки (рис. 3). Подставляя (11) в (12), найдем (13) Принимая во внимание равенства (13), представим выражения (11) в следующем виде: (14) Здесь (15) где (16) Безразмерный параметр γ характеризует сдвиговую податливость балки. Функции V и U определяют форму балки, которая находится под действием равномерно распределенной нагрузки. Из полученного решения (14), (15) и уравнений (10) следует, что при x = 0,l q = 0, y = -dw/dx = ±pl/(2K). Тогда равнодействующая нагрузки pl уравновешивается краевыми усилиями Q = Ky = ±pl / 2. Решение уравнений движения. Учитывая равенства (15), представим решение уравнений (6) в следующем виде: (17) Здесь F, Tx, Ty - неизвестные числа; Ux, Uy, Vx, Vy - базисные функции, определяемые выражениями (18) где (19) Безразмерные параметры γx и γy характеризуют сдвиговые податливости пластины. Выбранные функции (17) удовлетворяют граничным условиям (7). Подставляя равенства (17) в дифференциальные уравнения (6), получим следующие невязки: (20) Следуя методу Галеркина, потребуем, чтобы невязки R, Rx, Ry были ортогональны к соответствующим базисным функциям, т. е. (21) Подставляя (20) в (21), будем иметь (22) Подставляя Ux, Uy, Vx, Vy из (18) в соотношения (23), после интегрирования найдем (23) где (24) Подставим значения интегралов (23) в уравнения (22) и затем получившиеся соотношения приведем к безразмерному виду, умножив каждое из них на величину . В результате преобразований получим (25) Здесь (26) (27) (28) (29) (30) Величина η, входящая в уравнения (25) и определяемая равенством (30), является безразмерным частотным параметром. Таким образом, задача об определении основной частоты колебаний защемленной по контуру трехслойной пластины сведена к вычислению параметра η, при котором однородная система (25) будет иметь решение, отличное от нуля. Определение частотного параметра η и частоты колебаний ω. Представим систему (25) в виде следующего матричного уравнения: (31) где (32) Традиционный подход к определению собственных значений η предполагает решение характеристического уравнения det(A - ηB) = 0, которое в раскрытом виде представляет собой полином третьей степени. В уравнении (31) матрица A симметричная, а диагональная матрица B очевидно симметричная и положительно определенная. Поэтому все собственные значения обобщенной задачи являются вещественными и, как показывает практика вычислений, положительными. Минимальное из трех собственных значений будет определять требуемый частотный параметр η. Величина η зависит от параметров γx, γy, α, β1, β2, rx, ry, которые содержат всю информацию о жесткостных, геометрических и инерциальных характеристиках трехслойной пластины. Определение частотного параметра требует некоторых вычислительных усилий, которых можно избежать, если учесть, что инерция поворота мало влияет на частоту колебаний основного тона. Получим формулу для η, в которой не будет инерционных составляющих, связанных с поворотом нормали к срединной плоскости трехслойной пластины. Примем Dρ = 0. Тогда rx = ry = 0, и, как следствие, b22 = b33 = 0. Запишем второе и третье уравнения системы (25), учитывая, что a23 = a32, в следующем виде (33) Из (33) найдем (34) Подставляя (34) в первое уравнение системы (25), получим для η следующее выражение: (35) Получим выражение для частоты колебаний трехслойной пластины. Безразмерный частотный параметр η может быть определен из решения обобщенной задачи (31) или из формулы (37). При известном η из формулы (30) найдем частоту колебаний, т. е. (36) где . Используя равенства (8), преобразуем формулу (36) к виду (37) Выражение (37) определяет частоту колебаний основного тона закрепленной по контуру трехслойной пластины. Численный анализ и верификация результатов. Прежде чем приступить к численному анализу, получим выражения, определяющие параметры γx, γy, α, β1, β2, rx, ry, от которых зависит величина η. Подставляя (8) в (19), (28) и (30), будем иметь (38) (39) (40) Если материал заполнителя легкий, т. е. обладает малой жесткостью в направлениях x и y, то можно принять Amnδ = 0 (mn = 11, 12, 22, 33). Тогда из уравнения (9) следует, что ξmn = ξ (mn = 11, 12, 22, 33), и равенства (38), (39) примут следующий вид: (41) (42) Подставляя в формулу (37) ξ11 = ξ22 = ξ, будем иметь (43) Определим частотный параметр η и частоту колебаний основного тона f для квадратной трехслойной пластины с изотропными несущими слоями и легким заполнителем. Коэффициенты жесткости и модули сдвига материала несущих слоев имеют следующий вид: (44) Подставляя (44) в (41), получим (45) Учитывая равенства (45), определим для рассматриваемой пластины сдвиговые параметры γix, γiy (i = 1, 2, 3, 4). Из уравнений (24) и (39) будем иметь (46) Подстановка (44) в (42) при a = b дает (47) Для квадратной пластины инерционные параметры rx и ry (40) равны, т. е. (48) Определим коэффициенты матрицы А и В однородной системы линейных алгебраических уравнений (31). Учитывая равенства (46), (47) и (48), из соотношений (27) будем иметь (49) При известных коэффициентах (49) частотный параметр η определяется как минимальное собственное значение однородной системы (31). Формула для частоты колебаний квадратной трехслойной пластины с изотропными несущими слоями и легким заполнителем может быть получена из (43), если принять a = b и A11t=A22t=. Тогда (50) Для того, чтобы динамический анализ был более конкретным, зададим характеристики материалов и размеры. Пусть несущие слои выполнены из алюминиевого сплава с Et = 70 ГПа, vt = 0,3, Gt = 27 ГПа, ρt = 2700 кг/м3. Рассмотрим два типа заполнителя. Один из них будем называть податливым заполнителем, а другой - жестким заполнителем. Пусть у податливого заполнителя Gδ = 80 МПа, ρδ = 30 кг/м3, а у жесткого заполнителя Gδ = 400 МПа, ρδ = 80 кг/м3. Отметим, что в заполнителях увеличение жесткости, т. е. модуля сдвига Gδ, происходит всегда с одновременным увеличением плотности материала ρδ. Зададим размеры трехслойной пластины. Пусть суммарная толщина несущих слоев t = 0,001 м, толщина заполнителя δ принимает значения 0,02 и 0,1 м, а размер пластины в плане a принимает значения 0,5, 1,0 и 2,0 м. Для вычисления круговой частоты колебаний воспользуемся формулой (50). При этом частотный параметр η будет определен в зависимости от того, учитывается или не учитывается влияние инерции поворота на динамическое поведение трехслойной пластины. В табл. 1 приведены частоты колебаний f = ω/2π в герцах для трехслойных пластин с податливым и жестким заполнителями. Таблица 1 Частота колебаний f (Гц) трехслойной пластины с податливым и жестким заполнителем δ Податливый Жесткий a 0,5 1,0 2,0 0,5 1,0 2,0 0,02 759,2 760,4 247,4 247,5 68,1 68,2 886,1 888,3 240,7 240,9 61,6 61,6 0,10 1571,3 1574,1 656 658,1 222,5 223,1 2023,5 2048,8 664,8 670,6 185,0 185,5 В каждой из ячеек табл. 1 первое значение частоты определено с учетом инерции поворота (Dρ ≠ 0), а второе - без учета инерции поворота (Dρ = 0). Частотный параметр η для первого значения f находится как минимальное собственное число однородной системы (31). Частотный параметр η для второго в ячейке значения f был определен с помощью формулы (50). Анализ данных из табл. 1 позволяет сделать выводы о влиянии размеров трехслойной пластины, упругих характеристик и плотности составляющих материалов на частоту колебаний основного тона. Отметим, что влияние инерционной составляющей Dρ на основную частоту колебаний невелико, даже для толстых пластин (a = 0,5 м, δ = 0,1 м). Из табл. 1 видно, что разница между частотами, вычисленными с учетом и без учета инерции поворота, не превышает 0,3 %. Сравнивая данные в соответствующих ячейках двух таблиц, можно заметить, что увеличение жесткости заполнителя приводит к росту частоты колебаний только для пластин с a = 0,5 м. У пластин с a = 1,0 м, a = 2,0 м увеличение жесткости заполнителя (т. е. увеличение Gδ и ρδ) вызывает уменьшение частоты колебаний. Как показывает численный анализ, для пластин с большим отношением a/(δ + t) такой эффект обусловлен преобладающим влиянием на частоту колебаний плотности заполнителя ρδ по сравнению с влиянием, которое оказывают на частоту колебаний модуль сдвига Gδ. Обнаруженную зависимость f от Gδ и ρδ необходимо учитывать при выборе заполнителя в задачах проектирования трехслойных пластин, когда ограничения накладываются на основную частоту колебаний. Для проверки полученных результатов определим первую частоту колебаний рассматриваемой трехслойной пластины в пакете ANSYS. Воспользуемся при моделировании пластины конечным элементом SHELL181 с опцией, соответствующей расчету трехслойных конструкций. Форма колебаний пластины показана на рис. 4. Рис. 4. Форма колебаний квадратной трехслойной пластины, жестко закрепленной по контуру (МКЭ-расчет) В табл. 2 приведены частоты колебаний трехслойной пластины с податливым и жестким заполнителями, полученные с помощью метода конечных элементов. Таблица 2 Частоты колебаний f (Гц) трехслойной пластины с податливым и жестким заполнителем (МКЭ) δ Податливый Жесткий a 0,5 1,0 2,0 0,5 1,0 2,0 0,02 754,4 245,5 67,9 881,8 240,2 61,6 0,10 1561,9 652,1 221,2 2010,7 661,4 184,5 Сравнение результатов, приведенных в табл. 1 и 2, позволяет сделать вывод о том, что частоты, вычисленные методом конечных элементов, мало отличаются от частот, найденных на основе двух моделей, рассматриваемых в настоящей работе. Максимальная разница между данными в табл. 1 и 2 не превышает 1,9 %. Отметим, что конечно-элементный анализ подтвердил также обнаруженные ранее особенности влияния модуля сдвига заполнителя Gδ и его плотности ρδ на величину частоты колебаний трехслойной пластины. Выполненный выше анализ результатов дает основание утверждать, что определение первой частоты колебаний закрепленной по контуру квадратной трехслойной пластины вполне может быть выполнено с помощью формул (36) и (37). Проектирование. Расчеты по определению основной частоты колебаний трехслойной пластины с защемленными краями с использованием соотношений (36) и (37) не требуют значительных вычислительных затрат. Поэтому полученные формулы могут быть особенно полезны на этапе проектирования трехслойных пластин. Продемонстрируем это на примере задачи по выбору толщин заполнителя δ и несущих слоев t для квадратной трехслойной пластины с заданной частотой колебаний . Пусть несущие слои и заполнитель имеют следующие параметры: Et = 70 ГПа, vt = 0,3, Gt = 27 ГПа, ρt = 2700 кг/м3, Gδ = 400 МПа, ρδ = 80 кг/м3. Размер a = 1,0 м. Определим пространство проектирования в виде неравенств (51) Представим зависимость f(δ,t) набором двумерных контурных линий (рис. 5). Пусть заданная частота = 500 Гц. На рис. 5 контурная линия, соответствующая этой частоте, имеет увеличенную толщину. Рис. 5. Зависимость f(t, δ) в виде набора контурных линий Очевидно, что заданная частота может быть реализована при различных сочетаниях t и δ из пространства проектирования (51). Для каждой пары значений t и δ должно выполняться условие (52) где ω - круговая частота колебаний, вычисляемая по формуле (50). Подставляя (50) в (52), будем иметь (53) где . Для вычисления частотного параметра используются уравнения (30), (35). Задаваясь различными значениями t, решим нелинейное уравнение (53) относительно δ. Для получившейся пары t и δ определим массу пластины, т. е. (54) Результаты вычислений по описанному алгоритму для рассматриваемой пластины приведены на рис. 6 и 7. Рис. 6. Зависимость δ(t) Рис. 7. Зависимость m(t) Из анализа полученных данных следует, что с увеличением толщины несущих слоев t толщина заполнителя δ уменьшается (рис. 6). Масса пластины m в рассматриваемом диапазоне измерения t имеет минимум (рис. 7). Таким образом, mmin = 7,3657кг, topt = 0,0009 м, δopt = 0,0617 м. Рассмотренный выше подход к определению первой частоты колебаний трехслойной пластины с защемленными краями реализован в вычислительных программах. Заключение. Решена задача определения первой частоты колебаний трехслойной пластины с жестко закрепленными краями. Для решения уравнений движения был использован метод Галеркина. В качестве функций, аппроксимирующих форму пластины при колебаниях, было принято решение задачи об изгибе балки с жестко заделанными краями под действием равномерно распределенной нагрузки. Эти функции учитывают сдвиговую податливость пластины и удовлетворяют граничным условиям, предполагающим, что на каждом из краев пластины отсутствуют прогиб и два угла поворота нормали к срединной плоскости. Задача вычисления первой частоты колебаний трехслойной пластины сведена к определению безразмерного частотного параметра, величина которого зависит от геометрических, упругих и инерциальных характеристик трехслойной пластины. Рассмотрены два подхода к определению частотного параметра, которые отличаются друг от друга в зависимости от того, учитывается или не учитывается влияние инерции поворота на динамическое поведение трехслойной пластины. В первом подходе частотный параметр находился как минимальное собственное число соответствующей системы однородных алгебраических уравнений. При втором подходе для определения частотного параметра была получена аналитическая формула. На конкретных примерах был выполнен анализ влияния размеров трехслойной пластины, упругих характеристик и плотности соответствующих материалов на частоту колебаний основного тона. Расчеты проводились как с использованием рассматриваемых в работе подходов, так и с помощью метода конечных элементов. Анализ результатов позволил сделать вывод о том, что определение первой частоты колебаний закрепленных по контуру трехслойных пластин, обладающих разнообразными жесткостными параметрами и размерами, может быть с достаточно высокой точностью выполнено по полученным аналитическим формулам. Эти формулы будут особенно полезны при проектировании трехслойных пластин, когда ограничения накладываются на первую частоту колебаний.
×

About the authors

P. O. Deev

Reshetnev Siberian State Aerospace University

Email: prokhor777@gmail.com
31, Krasnoyarsky Rabochy Av., Krasnoyarsk, 660037, Russian Federation

F. K. Sinkovskiy

JSC Academician M. F. Reshetnev «Information Satellite Systems»

52, Lenin Str., Zheleznogorsk, Krasnoyarsk region, 662972, Russian Federation

References

  1. Yu Y. Y. A new theory of elastic sandwich plates - one dimensional core // Journal of Applied Mechanics. 1959. Vol. 26. P. 415-421.
  2. Yu Y. Y. Simple thickness-shear modes of vibration of infinite sandwich plates // Journal of Applied Mechanics. 1959. Vol. 26. P. 679-681.
  3. Yu Y. Y. Flexural vibrations of elastic sandwich plates // Journal of Aerospace Sciences. 1960. Vol. 27. P. 272-282.
  4. Yu Y. Y. Simplified vibration analysis of elastic sandwich plates // Journal of Aerospace Sciences. 1960. Vol. 27. P. 894-900.
  5. Habip L. M. A survey of modern developments in the analysis of sandwich structures // Applied Mechanics Reviews. 1965. Vol. 18. P. 93-98.
  6. Ha K. H. Finite-element analysis of sandwich plates: an overview // Computers and Structures. 1990. Vol. 37 (4). P. 397-403.
  7. Noor A. K., Burton W. S., Bert C. W. Computational models for sandwich plates and shells // Trans. ASME. Applied Mechanics Reviews. 1996. Vol. 49 (3). P. 155-199.
  8. Librescu L., Hause T. Survey of recent developments in the modeling and behavior of advanced sandwich constructions // Journal of Composite Structures. 2000. Vol. 48 (1-3). P. 1-17.
  9. Vinson J. R. Sandwich structures // Applied Mechanics Review. 2001. Vol. 54 (3). P. 201-214.
  10. Zenkert D. An introduction to sandwich construction. London : Chameleon Press, 1995. 277 p.
  11. Yu Y. Y. Vibrations of elastic plates: linear and nonlinear dynamic, modeling of sandwiches, laminated composites and piezoelectric layers. New York : Springer-Verlag, 1996. 225 p.
  12. Vinson J. R. The behavior of sandwich structures of isotropic and composite materials. Lancaster : Technomic, 1999. 378 p.
  13. Kollar L. P., Springer G. S. Mechanics of composite structures. Cambridge : University Press, 2003. 470 p.
  14. Паймушин В. Н. Точные аналитические решения задачи о плоских формах свободных колебаний прямоугольной пластины со свободными краями // Известия высших учебных заведений. Математика. 2006. № 8. С. 54-62.
  15. Мущанов В. Ф., Демидов А. И. Расчет трехслойной пластинки, лежащей на винклеровском основании, вариационно-разностным методом // Современное промышленное и гражданское строительство. 2010. Т. 6, № 2. С. 77-91.
  16. Коган Е. А., Юрченко А. А. Нелинейные колебания защемленных по контуру трехслойных пластин // Проблемы машиностроения и надежности машин. 2010. № 5. С. 25-34.
  17. Reissner E. The effect of transverse shear deformation on the bending of elastic plates // Trans. ASME. Journal of Applied Mechanics. 1945. Vol. 12 (2). P. 69-77.
  18. Mindlin R. D. Influence of rotary inertia and shear on flexural motions of isotropic elastic plates // Trans. ASME. Journal of Applied Mechanics. 1951. Vol. 18. P. 31-38.

Supplementary files

Supplementary Files
Action
1. JATS XML

Copyright (c) 2017 Deev P.O., Sinkovskiy F.K.

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

This website uses cookies

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

About Cookies