О наноотколе после воздействия ультракороткого лазерного импульса



бет2/3
Дата11.06.2016
өлшемі0.85 Mb.
#127706
1   2   3

4. Поведение около порога абляции
На рис. 5 приведена схема звукового распада и отражения от вакуумной границы начального профиля термомеханического напряжения . Схема построена в линейном акустическом приближении, когда гидродинамическая задача сводится к волновому уравнению , const. Общее решение задачи Коши , (начальные данные) этого уравнения дается решением Даламбера , , описывающим распад начального профиля на бегущие вправо и влево звуковые волны.
При распаде начального профиля у границы с вакуумом (рис. 5) в начальные данные по давлению , необходимо включить отраженную волну , . При этом скорость границы с вакуумом дается выражением . Сравним зависимости скорости границы от времени в термомеханическом и ударноволновом [23,24] случае. В случае с ударной волной на границу в момент приходят профили давления и скорости . Тогда скорость границы равна удвоенной [23,24] скорости за фронтом ударной волны . Как видим, в линейном приближении термомеханическое и ударноволновое отношения отличаются в два раза.
На рис. 5 показан пример с начальным профилем давления треугольной формы, распределенный при на отрезке . Мишень находится слева ( ) от границы с вакуумом . Начальный профиль разбит на сумму из двух половин из разбегающихся волн Даламбера. Они обозначены цифрами 2 и 3. Волна 3 распространяется налево вглубь мишени, а волна 2 движется направо. Волна 2 отражается от границы с вакуумом. Отраженная волна отмечена цифрой 1 на рис. 5. В момент она находится в области . На отрезке времени , к которому относится момент , идет процесс отражения. Мгновенный профиль в мишени есть сумма волн 1, 2 и 3, показанных по отдельности на верхнем из пары рисунков . Из-за отраженной волны 1 в мишени возникает область отрицательных давлений. После окончания отражения ( ) в мишень уходит пара из положительного и отрицательного импульсов, показанных на нижнем из рисунков на рис. 5. В ударноволновом случае после отражения в мишени распространяется только импульс отрицательного давления.
Интересно сопоставить линейное решение с численным. Мгновенные профили давления показаны на рис. 6-8 при следующих значениях параметров ( , ): 0.86, 0.3 (рис. 6), 0.86, 0.6 (рис. 7) и 1.08, 0.7 (рис. 8). Давление приводится в МД единицах давления . Рассмотрим сначала допороговые случаи (рис. 6 и 7). Линейное решение, состоящее из секций 1-2, j и 2-3, построено по экспоненциальному начальному профилю давления. Секции 1-2 и 2-3 составлены из волн (1 и 2) и (2 и 3), ср. с рис. 5. В решении Даламбера секции 1-2 и 2-3 разделяет скачок j. На самом деле из-за зависимости скорости звука от давления этот скачок размывается в подъем J. Численное решение представлено флуктуирующей кривой. В случае автомодельной волны разрежения [25] подъем J соответствует участку с веером характеристик. Течение становится автомодельным, если начальные данные однородны: const, const, 0. Давление, вычисленное по формулам линейной акустики, обращается в нуль в точке, в которой находилась граница с вакуумом в момент . Поэтому в области флуктуирующая кривая на рис. 6-8 находится левее линейного решения 1-2. Это вызвано смещением вещества.
Небольшая разница между линейным решением и численным, как сказано, заключается, во-первых, в размывании скачка j, и, во-вторых, в смещении границы с вакуумом. Более важным является отклонение от линейного поведения , const, связанное с пластическими эффектами возле предела прочности . Линейное приближение пригодно при малых возмущениях давления , . Поскольку объемный модуль намного превышает , то главное ограничение обусловлено пределом . Из-за пластичности вблизи порога абляции ( ) описывающее область линейное решение 1-2 (рис. 6-8) сильно отличается от численного. Имеется даже любопытная немонотонность в виде выброса давления в области . Эта немонотонность распространяется вглубь вещества вместе со звуковой волной, ср. рис. 6 и 7.
Аналогичная движущаяся с волной немонотонность наблюдается и в запороговом случае (рис. 8). Она отмечена цифрой 1 на рис. 9. Немонотонность 1 существует на стадии до образования разрыва вещества. Интересно, что точка разрыва не связана с местом нахождения немонотонности. Если немонотонность движется с волной, то разрыв «зацеплен» за вещество. Формирование разрыва является продолжительным процессом. Начало разрыва отмечено стрелкой 2 на рис. 9 на профиле плотности . Видно, что в момент в точке будущего разрыва имеется излом на профиле давления («слабое место»). Позже давление в этой точке приближается к нулю (профиль ).
На рис. 10 представлены линейная и МД зависимости для границы мишени. На стадии торможения границы (стадия 2, рис. 3, 4) эти зависимости примерно одинаковы, хотя из-за пластических эффектов профили тормозящего давления в глубине мишени отличаются сильно (см. рис. 6 и 7). Дело в том, что торможение границы и, следовательно, зависимость на стадии 2 определяются градиентом давления в непосредственной окрестности границы. Акустический и МД профили примерно одинаковы при небольших по абсолютной величине значениях отрицательного давления (рис. 6-8). Поскольку на границе выполняется условие (давление пара в области между кривыми 1 и 2-E-A мало, рис. 1), то как раз в приграничной области конденсированного вещества значения являются небольшими.
В линейной акустике в момент (рис. 5) отражение звука заканчивается. После этого в линейном приближении граница останавливается, сместившись за время торможения на расстояние . Подъем и спад МД функции на рис. 10 при , видимо, обусловлен приходом сигнала, отраженного от тыльного промежутка «atm-atm » расчетного слоя, см. рис. 7 и 8. На этом промежутке расположен «термостат». Он тормозит вещество промежутка так, чтобы смещение тыльной поверхности «atm » было мало и давление на границе «atm » равнялось нулю. Динамически действие такого термостата эквивалентно действию гравитационного «поля» с примерно однородным на промежутке зависящим от времени ускорением свободного падения, направленным налево на рис. 7 и 8. Действие эффективного «веса» увеличивает градиент на промежутке.
Начальная толщина слоя вещества составляет 1420 . Толщина промежутка равна 210 . Промежуток «atm-atm » размазывает и ослабляет отраженный сигнал по сравнению с расчетом с граничным условием в виде жесткой стенки в точке «atm ». Влияние термостата не сказывается на зависимости скорости внешней границы конденсированной части мишени от времени , если порог абляции превышен (кривые b и c на рис. 3). Дело в том, что в наших расчетах отрыв откольного слоя происходит до прихода отраженного сигнала в область отрыва. После отрыва слоя отраженный сигнал не может повлиять на динамику полета откольного слоя.
5. Переход от откольного к безоткольному истечению
На рис. 1 (b) пятью стрелками (i-v) указаны позиции характерных случаев, совокупность которых формирует наблюдаемую в экспериментах морфологию лазерного выброса (п. 3).

Описанию случаев (i) и (ii) около порога абляции посвящен п. 4. Теперь необходимо представить оставшиеся три случая, относящиеся к утоньшению купола (iii, iv) и течению в области отверстия купола (v).


Рассмотрим область между кривыми 1 и A-E-2 на рис. 1 (b), (c). Для случая (iv) с = 2.6 и (v) = 4.3 этой области соответствует отрезок между стрелками 1 и A-E на рис. 11 и отрезок между стрелками 1 и 2 на рис. 12 («хвостик»). В МД расчетах эта область заполнена переохлажденным паром. Из-за относительно небольшого времени МД моделирования и ограниченности числа Кнудсена по толщине этой области в МД расчетах не наблюдается конденсации пара. Количество пара в анализируемой области и ее толщина возрастают с ростом отношения , ср. рис. 11 и 12.
«Хвостик», похожий на распределение плотности в области переохлажденного пара на рис. 11 и 12, имеется и на профилях плотности, вычисленных из уравнений гидродинамики с равновесным уравнением состояния [25]. Равновесное уравнение состояния получается из неравновесного с метастабильной областью с помощью правила Максвелла. В расчетах [25] «хвостик» соответствует участку равновесной адиабаты, находящемуся в двухфазной области и соединяющему точку пересечения бинодали и адиабаты и точку . На равновесной двухфазной адиабате относительная доля пара по объему двухфазной смеси монотонно возрастает по мере смещения от точки пересечения бинодаль-адиабата к точке .
Рассмотрим двухфазную область. На рис. 1 (b) она ограничена кривой 2-E-cr и дном кратера cr. На рис. 1 (c) это область под кривой 2 и нижней границей купола. См. также рис. 11 и 12. Вещество, попадающее в двухфазную область, не сразу переходит в двухфазное состояние. Сначала вещество растягивается, а его плотность уменьшается без потери сплошности. Снижение плотности начинается под действием положительного давления и продолжается за счет сил инерции при . Смена знака давления происходит при прохождении равновесного значения плотности. При некотором растяжении вещества ниже равновесной плотности в среде зарождаются первые очень маленькие пузырьки. Период зарождения пузырьков короткий. Появление пузырьков уменьшает степень растяжения жидкости вне пузырьков, повышает плотность этой жидкости и соответственно уменьшает абсолютную величину отрицательного давления в жидкости, растянутой ниже равновесной плотности. Поэтому новых пузырьков более не образуется. Далее по мере растяжения происходит расширение пузырьков. При этом объемная относительная доля смеси, занятая паром в пузырьках, увеличивается, .
При образуется пена, в которой толщина стенок мала по сравнению с размером пузырей. Увеличение объема пузырей происходит, во-первых, из-за растяжения, и, во-вторых, вследствие разрыва стенок и слияния пузырей. Второй процесс уменьшает число пузырей на единицу массы смеси. Кроме того, число пузырьков уменьшается вследствие коллапса самых мелких пузырьков при уменьшении по абсолютной величине напряжений, вызванных растяжением. При малых значениях отношения доминирует первый процесс, при больших – второй. С приближением отношения к единице вместо пены формируется система из капель жидкости. Капли возникают в местах схождения нескольких стенок пузырей пены.
На рис. 13 показан процесс формирования и эволюции двухфазной системы, включающий растяжение сплошной среды, зарождение и увеличение радиуса и объема пузырей до образования пены. С течением времени система трансформируется от состояния пузыри в жидкости (перколяция по жидкости) к состоянию капли в паре (перколяция по пару). На правую из четырех фигур рис. 13 попал момент разрыва стенки между пузырями. Прорывается стенка между вторым сверху справа пузырем и соседним к нему средним пузырем.
МД расчеты выполнены для атомов, взаимодействующих по закону Леннарда-Джонса (ЛД). Использование многопроцессорного алгоритма [26] и простота ЛД потенциала позволили провести расчеты с большим количеством атомов (10-100 миллионов). В результате была исследована сложная картина пены с очень большим количеством взаимодействующих пузырей. С пеной связан важный эффект длительного торможения откольного образования. Эффект обусловлен слабым, но конечным сопротивлением пены растяжению из-за поверхностного натяжения стенок пузырей.
В МД расчетах (рис. 13-15) моделировалось воздействие ультракороткого лазерного импульса на полубесконечную мишень. Поэтому начальный размер кристалла по оси разлета был велик. Начальный размер слоя составлял 1420 на 36 на 151 . На рис. 13, 14 и верхнем из рис. 15 показана плотность вещества в плоскости . До начала нагрева атомы были упакованы гранецентрированную кубическую (ГЦК) равновесную решетку с при . Для укрупнения масштаба на рис. 13-15 показана только та часть слоя по , которая охвачена движением.
На рис. 13 начальная температура поверхности и глубина прогрева , определенные в п. 2, равны и 590 . Слева направо показаны моменты времени (50, 0.12); (150, 0.36); (450, 1.07) и (1175, 2.79) в МД единицах времени (МДЕ) и в единицах . На рис. 14 последовательно слева направо значения четверки параметров (МДЕ), равны (1.25, 770, 2100, 3.81); (2, 590, 1000, 2.37); (3, 500, 570, 1.6); (5, 414, 180, 0.61). На рис. 15 показан расчет при значениях параметров равных 2, 590, 1000, 2.37. Толщины откольного образования равны в случаях с =1.25, 2 и 3 значениям 250, 80 и 30 соответственно.
На рис. 13 и 14 ось основного движения направлена по вертикали рисунка. Такая ориентация выбрана с целью соответствия направлению лазерного луча на рис. 1. Рис. 13 показывает, что происходит в поверхностном слое мишени с течением времени («развертка по времени»). Назначение рис. 14 – проиллюстрировать, как устроена структура факела в поперечном лазерному лучу направлении из-за изменения флуенса внутри лазерного пучка («развертка» по оси , рис. 1).
Как говорилось выше, над откольным образованием располагается область однофазного пара, см. рис. 1, 11 и 12. Снизу эта область ограничена верхней из трех черточек на правом из рис. 13. Черточка на правом кадре на рис. 14 отделяет область однофазного пара от двухфазной области, в которой переход в двухфазное состояние происходит в результате механического разрыва вещества. На трех левых кадрах на рис. 14 область однофазного пара находится над верхней границей A-E откольного слоя, см. рис. 1 b, 11. Правый из кадров на рис. 14 размещен по оси так, чтобы в поле этого кадра попала большая часть области однофазного пара между метками 1 и 2 на рис. 12.
На рис. 15 повторен второй слева кадр с рис. 14 и приведены соответствующие ему профили средней плотности и давления. Всплески плотности связаны с присутствием крупных пузырей. На профиле плотности выделяется летящее впереди откольное образование из сплошного вещества. На профиле давления ясно видна область растяжения с отрицательным давлением, связанная с неразрушенным веществом вблизи дна будущего кратера. Она находится справа, ср. кадры на рис. 15. Обратим внимание на зону пены на профиле давления. Здесь имеется небольшое отрицательное давление. Это свидетельствует обо все еще продолжающемся сопротивлении пены растяжению, хотя средняя плотность двухфазной смеси составляет меньше половины от исходной плотности вещества. Это связано с действием поверхностного натяжения в стенках пузырей, см. верхний из рисунков 15.

6а. Кольца Ньютона и микроинтерферометрия
Сравним развитые выше теоретические представления с экспериментом.

Источником излучения служила хром-форстеритовая фемтосекундная лазерная система [27,28], генерирующая импульсы длительностью ≈ 100 фс на длине волны λ1 = 1240 нм. Лазерный импульс разбивался на два импульса (схема разделения не показана): мощный нагревающий (1) и слабый диагностический (2). Основная длина волны λ1 использовалась для нагрева образца. Нагревающий импульс p-поляризации фокусировался на мишень под углом падения 45є в эллиптическое пятно с диаметром порядка 50 мкм по малой оси. Его энергия плавно варьировалась с помощью поляризационного ослабителя и контролировалась калиброванным фотоприемником. Пространственное распределение плотности энергии в пятне на мишени соответствовало гауссову. Параметры распределения – плотность энергии в центре пятна Fc и радиус пучка на уровне e-1 определялись с помощью методики, описанной в [ ] .

Диагностический импульс задерживался относительно греющего по времени прибытия на мишень на промежуток времени . Для диагностики применялась вторая гармоника лазерного излучения с длиной волны λ2 =620 нм. Преобразование частоты осуществлялось в кристалле LBO.

В качестве экспериментальных образцов использовались золотые пленки толщиной 200 нм, напыленные на стеклянные подложки. После каждого воздействия нагревающего импульса образец смещался на новое место. Перемещение образца осуществлялось с помощью трехкоординатного позиционера с шаговыми двигателями, управляемого с компьютера. Наведение на фокус по одной из координат производилось с точностью до 1 мкм.

Исследование динамики образования откольного купола проводилось с помощью двух методик: микроинтерферометрической и наблюдением за интерференцией при разлете поверхностного слоя (кольца Ньютона) .

В первом случае регистрация интерферограмм возбужденной поверхности образца осуществлялась с помощью микроинтерферометра Линника (рис. 16). На входе интерферометра диагностический импульс разделялся на два луча – объектный (3) и опорный (5) с помощью делителя пучка. Луч (3) через микрообъектив осуществлял подсветку мишени в области нагрева. Этот же объектив использовался для переноса изображения поверхности мишени в плоскость ПЗС матрицы примерно с десятикратным увеличением. Опорный луч (5-6-7), отраженный от опорного зеркала, интерферирует в плоскости матрицы с объектным лучом (3-4), образуя систему интерференционных полос (рис. 17). В фазе и амплитуде объектного импульса содержится информация о состоянии мишени в области нагрева на момент времени после воздействия греющего импульса (1). Из смещения полос извлекаются сведения о форме поверхности, отражающей объектный луч. Размер пятна диагностического импульса примерно на порядок больше нагревающего. Поэтому на интерферограммах видна как сильно нагретая область с выпучиванием поверхности, так и периферия пятна, на которой возмущение поверхности отсутствует.

Компьютерная обработка интерферограмм проводилась с помощью алгоритма двумерного Фурье-преобразования, обеспечивающего точность измерения смещения поверхности до нескольких нанометров. Пространственное разрешение в плоскости мишени соответствовало примерно 2 мкм. Временное разрешение интерферометрической методики составляло ~10-13 с и определялось длительностью лазерных импульсов. В одном опыте можно использовать только один диагностический импульс вместе с греющим импульсом, поскольку время экспозиции ПЗС приемника (~10-5 с) намного больше, чем требуемые для рассматриваемых целей задержки .

Сравним эксперименты с кольцами Ньютона (рис. 18) и микроинтерферометрию [1,2,27-29]. При опытах с кольцами опорный импульс 5-7 отсутствует (см. рис. 16). В этом случае диагностический луч 4 состоит из двух лучей. Первый из них - это луч , отраженный от купола, а второй луч складывается в результате отражения от дна кратера. Эти лучи интерферируют. Поэтому на изображении, регистрируемом ПЗС камерой, в области нагрева видны чередующиеся кольца с центром на оси греющего пучка. Купол и дно кратера показаны на рис. 1, 21.


Для проявления интерференции лучей в виде колец необходимо, чтобы расстояние от купола до дна было достаточно большим. А именно, это расстояние должно превышать половину длины волны диагностического импульса. При использовании второй гармоники лазера на хром-форстерите имеем λ2 /2 = 310 нм. Возможности методов с кольцами Ньютона и микроинтерферометрией разные. Кольца появляются на далеких временах (порядка нс) при значительных удалениях купола (порядка мкм). Тогда как микроинтерферометрия позволяет начать слежение за движением поверхности облученной мишени с интереснейшей ранней стадии, на которой удаления будущего купола еще малы. Речь идет о временах 10 пс и расстояниях порядка нескольких десятков нм. На этих временах только начинается формирование двухфазной области.
7а. Анализ экспериментальных данных
Были проанализированы результаты ряда экспериментов, сочетающих микроинтерферометрию и измерения с помощью колец Ньютона. Для изучения развития процесса во времени были проведены серии опытов с одинаковыми значениями параметров и и разными задержками .

Диапазон вариации задержек состоит из шести характерных поддиапазонов: (i) ультракороткие задержки 0.1, 0.2 пс и т.д., (ii) короткие задержки 1, 2 пс, (iii) ранняя стадия 10, 20 пс, (iv) переходная стадия 50, 100 пс, (v) средняя стадия 300, 500 пс и (vi) далекая стадия 1, 2 нс. На первых стадиях (i)-(iii) происходит электрон-ионная релаксация, которая здесь в данной работе не рассматривается.


Из интерферометрических измерений следует (рис.19), что глубина кратера слабо зависит от флуенса в диапазоне 1.5ч5 и составляет примерно

=110 10 нм. (3)

Из расчетов с учетом (3) получаем, что глубина прогрева медленно убывает от значения 300 нм на пороге до значений 150 нм при 4-5 Дж/см . Следовательно, продолжительность переходной стадии (iv) порядка звукового времени (скорость звука =3.3 км/с для золота). На этой стадии формируется откольный купол и двухфазное облако. Переходной стадии соответствует «звуковой зубец» (рис. 3), в течение которого продолжается торможение откольного слоя. При этом скорость слоя снижается от значения до значения (рис. 2). Между переходной (iv) и средней (v) стадиями постепенно прекращается пополнение факела новыми порциями вещества и устанавливается режим разлета по инерции.


Измерение зависимости диаметров кратера от энергии лазерного импуьса (рис. 19а) в случае гауссова пространственного распределения нагревающего импульса позволяет с высокой точностью измерить порог абляции по падающему флуенсу [ПЖЭТФ-2006]. Его величина для золота в наших экспериментах составила:

≈ 1.35 Дж/см (4)

Очень важным параметром для расчетов является коэффициент поглощения мишени в рассматриваемом диапазоне флуенсов. В данной работе с помощью системы калиброванных фотоприемников было проведено измерение интегрального коэффициента отражения R мишени в диапазоне от 1 до 7 Дж/см . Соответствующая величина коэффициента поглощения оказалось равной:



0.3. (5)

С учетом (4) и (5) поглощенный флуенс на абляционном пороге составляет ≈ 0.4 Дж/см . Падающий (4) и поглощенный флуенсы приводятся в расчете на единицу поверхности мишени.


В инерционном режиме разлета скорость перестает изменяться во времени (рис. 3). При достаточно больших задержках начинает выполняться асимптотическая формула , связывающая смещение откольного слоя и его скорость . И МД моделирование, и оба экспериментальных метода (микроинтерферометрия и кольца Ньютона) указывают на рост скорости с ростом флуенса. Соответственно вещество откольного слоя образует вздымающийся в сторону больших флуенсов купол, как показано на рис. 1 (b). Об этом однозначно свидетельствуют искривление интерференционных полос и чередование колец Ньютона.
По данным эксперимента (рис. 17) находится форма купола. На рис. 20 построены профили h(y) для двух опытов с разными значениями центрального флуенса при = 1.5 нс. Там же приведены профили (пунктир) соответствующих кратеров абляции. С помощью формулы (1) находится зависимость h(Finc) приведенная на рис. 21 для опыта с =7.2 Дж/см2.
Важно, что с ростом флуенса толщина купола убывает. Следовательно, существует значение такое, что выполняется =0 (порог испарения). Нас интересуют данные по характеристикам факела при возможно больших превышениях порога (4), чтобы в одном факеле присутствовали и купол, и безоткольное течение, как на рис. 1. При больших превышениях течение разлета в окрестности порога абляции около края кратера сильно отличается от течения в центральной области, см., например, МД иллюстрацию на рис. 14 с примерами по параметру от 1.08 до 4.3. Диапазон МД вариантов на рис. 14 накрывает порог испарения. Основная цель данной работы заключается в доказательстве существования порога в эксперименте на золоте.
Эксперимент однозначно свидетельствует об убывании толщины купола с ростом флуенса. Уже этот факт означает, что опытным путем доказано существование порога испарения. Интересно оценить величину , и изучить купол вблизи этого порога и безоткольное течение за этим порогом.
Полосы на интерференционном изображении на рис. 17a начинают искривляться на границе кратера и следуют в соответствии с поверхностью куполообразного выпучивания. На рис. 17b с большим центральным флуенсом плавно искривляющиеся полосы пропадают по достижении границы внутренней «дыры», которой (см. рис.21) соответствует флуенс:

≈ 5 Дж/см . (6)
Исчезновение плавно продолжающихся полос означает, что толщина купола сильно уменьшилась по сравнению с большим значением у границы кратера =110 нм (3). Тонкий купол пропускает свет, отраженный от структур снизу. Отражение от него самого уменьшается с уменьшением его толщины и при малых толщинах становится незаметным на фоне света, отраженного от дна кратера. Конечно, это не означает, что купол исчез. Однако величина (6) дает первое ограничение снизу для истинного значения порога .
Кольцам Ньютона соответствуют осцилляции коэффициента отражения. По мере того, как купол становится тоньше, амплитуда осцилляций убывает. Кольцам Ньютона в отличие от микроинтерферометрии «не мешает» опорная волна 5-7 (рис. 16). Поэтому ослабляющиеся осцилляции прослеживаются на более тонком куполе, где интерферометрические полосы от купола уже не видны. На рис. 19 осцилляции исчезают при значении

≈6 Дж/см . (7)

Таким образом, приходим к другой более высокой, чем (6), оценке порога испарения. Порог (7) превышает порог абляции (4) примерно в четыре раза. Это указывает на высокую прочность жидкой фазы золота, поскольку именно прочность обеспечивает существование откольного слоя на краю облака разлета.


Из общих соображений следует, что значение флуенса , при котором достигается термодинамическая критическая точка, выше, чем порог испарения . МД расчеты дают значение для ЛД системы, где - температура на изохоре начальной плотности. Из расчетов ЛД адиабат на фазовой плоскости получается, что адиабата с начальной температурой = 11 проходит через критическую точку при пересечении с бинодалью – кривой равновесия жидкость-пар. Таким образом, для ЛД системы отношение значительно (равно 3). Значение этого отношения служит некото

рым указанием для оценки отношения в случае золота. Из этих соображений и оценки (7) порога испарения снизу следует, что величина является очень значительной.





Достарыңызбен бөлісу:
1   2   3




©dereksiz.org 2024
әкімшілігінің қараңыз

    Басты бет