Оптимизация конструкции проводящего слоя в зоне изгиба гибких органических светоизлучающих диодов
С целью решения проблемы обрыва в проводящем слое, вызванного складкой в зоне изгиба гибких органических светоизлучающих диодов FOLED, и повышения эффективности оптимизирована схема проводящего слоя в этой зоне. На полимерной подложке были приготовлены образцы FOLED, в которых зона изгиба сформирована путем наложения полимера, контакта, металлической проволоки и других слоев. Проанализировано эквивалентное рабочее состояние процесса изгиба и проведено трехмерное конечно-элементное моделирование области изгиба FOLED. По результатам имитационного анализа области изгиба предлагается оптимальная толщина слоя полимера. При этом должны быть строго учтены пороговые условия повреждения материала в проводящем слое: необходимо, чтобы кумулятивная пластическая деформация в положениях А - А и В - В была меньше, чем в слое без повреждений. Показано, что оптимальная схема применения проводящего слоя заключается в использовании оптимального распределения механических напряжений для уменьшения угла скрутки, приводящего к разрушению зоны изгиба FOLED.
Optimized design of wiring scheme for conductive layer in bending area of FOLED.pdf Введение Органические светоизлучающие диоды (OLED) представляют одну из наиболее перспективных технологий для отображения информации благодаря своей высокой яркости, широкому углу обзора, низкому энергопотреблению и широкому диапазону рабочих температур [1, 2]. Одно из самых больших преимуществ OLED заключается в том, что они могут быть реализованы, например, в виде гибкого дисплея, что показывает преимущества органического полупроводника [3]. Гибкие органические светоизлучающие диоды (FOLED), представляющие собой тип OLED, состоят из анодов, катодов и органических функциональных слоев между ними. По крайней мере, одна сторона электродов прозрачна для получения светоизлучающей поверхности [4, 5]. При этом подложка FOLED гибкая, она легче, тоньше и более ударопрочна, чем OLED-дисплей на стеклянной подложке [6]. Кроме того, ожидается, что FOLED будет производиться в виде рулонов, что значительно снизит их стоимость. Разработка OLED началась в 1982 г. Была использована проводящая подложка из стекла, содержащая оксиды индия и олова (ITO) [7, 8], на которую напылялась пленка фталоцианина с медью толщиной 0.1 мкм (CuPc) с хорошей термостойкостью, низкой рабочей функцией и высокой эффективностью инжекции на аноде ITO для введения в проводящий слой. Для предотвращения появления точечных пробоев, приводящих к повреждению изоляции, на слой CuPc был нанесен слой полимерной пленки, содержащей флуоресцентный пигмент трифенилбутадиен. При использовании катода из Ag при постоянном напряжении 30 В был получен синий люминесцентный материал (EL) с яркостью 170 кд/м2. С этого времени начались исследования органической пленки. Позже стала применяться ультратонкая пленочная технология и органический дырочный проводящий слой TPD с лучшим эффектом передачи [9]. В разработанном люминесцентном слое OLED использовались катоды из сплавов Alq3 (гидроксихинолинат алюминия) и Mg и Ag с низкой рабочей функцией, что значительно повысило эффективность инжекции носителя [9]. С помощью этих мер светимость EL OLED при напряжении питания 10 В постоянного тока может достигать 1000 кд/м2. В 1988 г. К. Адачи с соавторами предложили многослойный сэндвич-дисплей, который сделал исследования OLED актуальными. Использовали полианилин (PANI) для скручивания проводящей пленки на гибкой прозрачной подложке и полиэтилентерефталат (PET) путем нанесения раствора покрытия в качестве прозрачного электрода [10]. После более чем 20-летнего применения FOLED становится все более востребованным. По сравнению с традиционным OLED, FOLED обладает уникальной гибкостью благодаря своей подложке. Для изготовления FOLED могут использоваться гибкие пластиковые материалы (такие, как PI или PET), металлические фольги (фольга из нержавеющей стали) или ультратонкое стекло для достижения большой кривизны или даже скручивания. Однако в области изгиба проводящий слой легко ломается, образуются нарушенные области. Эффективность электропроводимости в местах изгиба различна, поэтому выбор лучшей схемы и уменьшение нарушений является одной из важных проблем повышения эффективности электропроводности [11]. В данной работе построена трехмерная конечно-элементная модель для симуляции процесса изгиба, а затем вводится двумерная плоская эквивалентная напряжению модель схемы электропроводки. Эффективность подхода проверяется в результате анализа характеристик напряжение - деформация и их сравнения с результатами испытаний. На этой основе предлагается идея создания схем электропроводки для оптимальной конструкции. Материалы и методы Приготовление полимерной подложки В качестве материала подложки для получения слоистых компонентов используется полиэтилентерефталат с проводящей пленкой оксида индия-олова (ITO). Сопротивление блока ITO составляет 60 Ом. Устройство имеет структуру ITO/NPB. Эффективная световая площадь составляет 3×3 мм2. Фотоэлектрические характеристики проверяли полупроводниковым тестером Keithley 4200. В табл. 1 приведены параметры яркости и эффективности OLED-компонентов, приготовленных на полимерной и стеклянной подложках. Как следует из табл. 1, оптические и электрические свойства полимерной подложки, поданной с той же структурой, аналогичны таковым у стеклянной подложки OLED. Таблица 1 Сравнение производительности устройств со стеклянной и полимерной подложками Подложка Рабочее напряжение, В Токовая чувствительность, кд/с Плотность тока, А/м2 Яркость, кд/м2 Стеклянная 7.8 8.44 1431 12621 Полимерная 7.8 10.32 1581 14510 3D-конечно-элементное моделирование области изгиба. Конструктивные особенности многослойной укладки Типичная конструкция со штабельной укладкой по толщине для зоны изгиба получена на основе патентной технологии компании «LG» и показана на рис. 1. Соответствующая информация о толщине слоя и свойствах материалов приведена в табл. 2. Рис. 1. Состояние укладки многослойной стековой конструкции в зоне изгиба компании «LG» Таблица 2 Информация о толщине слоя и свойствах материала в зоне изгиба Параметр Смола PLN PVX Металл PI PI Толщина, мкм 76.43 1.72 0.18 0.68 1.45 19.25 Модуль упругости, ГПа 0.48 6.92 78.25 70.5 6.95 6.95 Схема алмазной поперечной проводки адаптирована для металлической проводки (рис. 2). Рис. 2. Алмазная схема маршрутизации металла в зоне изгиба компании «LG» Эквивалентное рабочее условие процесса изгиба Разные формы изгиба имеют различную специфику процесса [12], но механическое напряженное состояние этой области в основном одно и то же, как показано на рис. 3. Типичный вид процесса изгиба приведен на рис. 3, а. В процессе 1 из начально прямого состояния 1 образец постепенно изгибается к конечному кольцевому. Аналогично, в процессе 2 образец сворачивается из начального состояния 2 в конечное кольцевое состояние. Конечное состояние обоих процессов одно и то же. В анализе в основном фокусируется внимание на характеристиках напряжений средней части испытуемого образца, а эквивалентный процесс изгиба (рис. 3, б) может быть использован для эффективного уменьшения размера модели в соответствии с симметричными характеристиками напряжений при половинном изгибе [13]. Рис. 3. Эквивалентный процесс изгиба участка Разработка трехмерной модели конечных элементов В соответствии со структурными данными и характеристиками напряжений области изгиба для завершения моделирования в HyperMesh 2017 была использована 1/4 модели одиночной металлической проволоки (рис. 4). Рис. 4. Модель конечных элементов 1/4 единичной металлической проволоки в зоне изгиба Рис. 5. Детальное отображение конечно-элементной сетки в области изгиба С целью повышения эффективности и точности анализа для уменьшения сетки моделирования используются шестигранные элементы (рис. 5). Согласно табл. 1 и рис. 5, толщина слоя PVX составляет 0.2 мкм. При проведении анализа блок находился в состоянии изгиба. Поэтому для получения более точных результатов необходимо использовать больше точек расчета в направлении толщины. Как правило, по толщине блок разделяют на три-четыре части, но это приводит к снижению плотности ячеек и эффективности расчета [14]. Поэтому в данной работе для моделирования был взят элемент толстой оболочки (T_ Shell), так что использовался только однослойный элемент и развернуто больше точек интегрирования, чтобы рассмотреть эффективность и точность расчета. Выбрано квазистатическое неявное решение. Использовался программный комплекс ANSYS 19.0 с учетом геометрии, приведенной на рис. 3. Конец А закреплен, а нормальное (направление Z) и вертикальное (направление Y) смещения нижнего узла ограничены. Конец В движется методом смещения при нагружении для придания жесткости торцу [15]. Соединение MPC пилотного узла с поверх¬ностью осуществляется с помощью точечного контакта поверхности, и к пилотному узлу прикладывается смещающая нагрузка при повороте на 90 (рис. 6). Рис. 6. Граничные условия изгиба Результаты Анализ результатов моделирования области изгиба Как несложно видеть из анализа эквивалентного процесса изгиба (см. рис. 3), в этом процессе материалы в каждом слое будут иметь значительную деформацию. Для пластических материалов кумулятивная эквивалентная пластическая деформация является важным показателем для измерения степени повреждения материалов и кумулятивной величины приращения эквивалентной пластической деформации: , (1) где - кумулятивная эквивалентная пластическая деформация; - ее приращение. Тогда , (2) где - приращение пластической деформации. Поэтому . (3) Рис. 7. Кумулятивная эквивалентная пластическая деформация каждого слоя в конечном состоянии Кумулятивная эквивалентная плас¬тическая деформация каждого слоя материала определяется так, как показано на рис. 7. Согласно рис. 7, проводящий слой имеет наибольшую кумулятивную эквивалентную пластическую деформацию. Кривая роста для двух типичных положения A-A и B-B показана на рис. 8. Положение А-А более подвержено повреждениям, что согласуется с результатами испытаний, показанными на рис. 9. Рис. 8. Кривая роста кумулятивной эквивалентной пластической деформации А-А и В-В Рис. 9. Микрофотография разрушения при A-A-позиции в испытаниях на изгиб На стадии изгиба (рис. 10) позиция А-А проводящего слоя находится в опасной зоне из-за склонности к разрушению. Изменение кривизны позиции В-В невелико, что безопасно. Рис. 10. Изменение состояния А-А и В-В при изгибе На стадии прессования, хотя деформация в позиции А-А проводящего слоя происходит при уменьшении общей деформации, эквивалентная пластическая деформация будет увеличиваться [16-18]. Скорость будет уменьшаться, а кривизна при деформации в позиции B-B будет постепенно расти в процессе превращения полукруга в полуэллипс. Из-за большого диапазона изменения постепенно позиция B-B войдет в опасное состояние, что приведет, скорее всего, к разрыву в проводящем слое. Оптимальная толщина слоя полимера При оптимизации толщины слоя полимера принимается во внимание, что обычно она используется в качестве одной из переменных для регулировки положения нейтрального слоя, снижения уровня напряжений токопроводящего слоя, снижения скорости повреждения области изгиба [19]. Переменная толщины слоя смолы задается как T. Кумулятивная эквивалентная пластическая деформация положений А-А и В-В в проводящем слое в зоне изгиба принимается в качестве целевого значения для завершения оптимизации (условия оптимизации в табл. 3). Таблица 3 Условия оптимизации толщины слоя смолы Переменная Интервал изменения, мкм Цель оптимизации T (полимер) 80-140 min{plAA, plBB} Кривые отклика T (полимера), plAA и plBB, показаны на рис. 11. С увеличением T (полимера) plAA увеличивается, а plBB уменьшается. Основная причина заключается в следующем. С увеличением толщины слоя полимера T жес¬ткость образца вблизи ограничительного конца, т.е. положения А-А и его левой части, будет увеличиваться, что приведет к изменению формы в позиции А-А. Сравним состояния при трех типичных толщинах (T = 100, 120, 140 мкм) после изгиба (рис. 12). Угол в положении А-А является наибольшим, значит и деформация является наибольшей. Следовательно, проводящий слой в этом месте склонен к разрушению. В то время как при T = 100 мкм угол наименьший, а значит, и состояние в этом месте является наиболее безопасным. В позиции B-B кривизна при Т = 140 мкм является самой маленькой и деформация самая маленькая, т.е. в это время в этом месте электропроводящий слой склонен к разрушению. Искривление при Т = 100 мкм является наибольшим и провод в токопроводящем слое может разрушиться. Эти характеристики также согласуются с кривыми отклика позиций A-A и B-B. Рис. 11. Кумулятивный эквивалентный отклик на плас¬тическую деформацию в позициях А-А и В-В Тенденция кумулятивной пластической деформации в положениях А-А и В-В противоположна. Хотя кумулятивная пластическая деформация в положении А-А меньше, чем в положении В-В, диапазон изменений А-А больше. При оптимизации толщины слоя для предотвращения разрушения проводящего слоя следует принимать во внимание условия разрушения материала [20]. Кумулятивная пластическая деформация в положениях А-А и В-В должна одновременно быть меньше достаточной для разрушения материала. Рис. 12. Сравнение трех состояний изгиба при трех толщинах Оптимизированная конструкция токопроводящего слоя. Стратегия оптимизации конструкции Для повышения эффективности оптимизации и упрощения модели используется небольшой участок конечно-элементной модели в области изгиба и применяется такая же продольная деформация, как на рис. 13. Рис. 13. Упрощенная модель оптимальной конструкции Распределение напряжений и пиковое значение упрощенной модели согласуются с данными исходной модели, поэтому упрощенная модель может быть использована для сравнения различных схем. Кроме того, поскольку на оценку эквивалентного номинального напряжения будут влиять сингулярности напряжений, постольку используется метод линеаризации для извлечения суммы мембранного напряжения и напряжения изгиба для оценки эквивалентного напряжения опасной точки в расчетной схеме, устранения влияния полосы пиковых напряжений и снижения чувствительности сетки модели в процессе оптимизационного проектирования [21]. Если взять в качестве примера исходную схему, то влияние сетки на номинальное значение эквивалентной эффективной силы, изгиб и напряжение мембраны после линеаризации в опасной точке показано на рис. 14. Рис. 14. Сравнение различных методов оценки эквивалентного напряжения Эквивалентное номинальное напряжение понижается со 165 до 95 МПа, а линеаризованное эквивалентное напряжение - со 105 до 90 МПа с увеличением размера ячейки сетки (рис. 14). Причем амплитуда уменьшения меньше эквивалентного номинального механического напряжения. В качестве критерия текучести используются номинальное и линеаризованное эквивалентные напряжения. Контур напряжений используется для установления их распределения в упрощенной расчетной модели, которая может четко описать изменение результата, так что можно быстро определить наиболее опасную область [22]. Поэтому при уменьшении эквивалентных номинального и линеаризованного напряжений возможна реализация упрощенной схемы. Опасность зоны снижается, поэтому за счет оптимизации дизайна безопасность FOLED в зоне изгиба повышается. Стратегии сравнения схем подключения могут быть стандартизированы следующим образом: 1. Приняты четыре конструктивных элемента для схем проводки в горизонтальном сравнении, а конец нагрузки смещения увеличен для уменьшения конечного эффекта. 2. Используется одна четверть симметричной плоской модели напряжения для нагружения одной и той же продольной деформации. 3. Используется метод эквивалентной линеаризации напряжений для оценки вероятности разрушения, а в качестве цели оценки берется сумма мембраны и напряжения изгиба. Некоторые идеи и результаты оптимизации монтажных схем В качестве примера взят один проводник, площадь которого имеет прямоугольную форму. В этой работе приведены и классифицированы несколько схем (рис. 15). Конструктивные особенности различных схем сводятся к следующему: Класс А: полностью заполненная конструкция, материалы равномерно распределены при высоком общем уровне напряжения мембраны. Класс B: конструкция поперечной и наклонной линий распределяет продольное и поперечное напряжения, эффективно снижает общий уровень напряжения мембраны, но концентрации напряжений легко возникают в локальном месте. Конструкция класса С по сравнению с конструкцией класса А снижает количество материалов, но незначительно изменяет основное механическое напряжение, что приводит к невозможности снижения уровня напряжений при одновременном увеличении концентрации напряжений в локальном месте. Рис. 15. Три типа схемы конструкции По сравнению с тремя видами расчетных схем, с точки зрения оптимизации распределения напряжений для уменьшения трещиноватости проводка должна выполняться в соответствии со следующим порядком: поперечный тип прокладки > полный тип заполнения > внутренний тип отверстия. Обсуждение В данной работе в качестве подложки для получения слоистых компонентов используется полиэтилентерефталат с проводящей пленкой из оксида индия-олова. На основе FOLED-компо¬нентов выполнено трехмерное конечно-элементное моделирование области изгиба. В качестве объекта исследования была выбрана конечно-элементная модель. Проведены анализ оптимизационного проектирования схемы проводящего слоя, моделирование области изгиба и оптимизация толщины слоя полимера. FOLED-материал каждого слоя будет сильно деформирован и при этом будет генерироваться соответствующая кумулятивная эквивалентная пластическая деформация, которая для проводящего слоя является наибольшей [23]. Проанализированы типичные положения А-А и В-В проводящего слоя в зоне изгиба. На стадии изгиба позиция А-А более опасна, потому что проводящий слой в этом положении с большей вероятностью оборвется. Кривизна в положении B-B относительно невелика, поэтому проводящий слой в этом положении находится в относительной безопасности. На стадии сжатия эквивалентная пластическая деформация позиции А-А медленно повышается, а деформационная кривизна положения В-В постепенно увеличивается и переходит в опасное состояние. В это время позиции А-А и В-В становятся непрочными, растет вероятность перелома. В аспекте оптимизации толщины слоя полимера с увеличением его толщины пластическая деформация в положении А-А растет, а в положения В-В снижается. Когда толщина составляет T = 140 мкм, деформация в положении A-A является самой большой, а проводящий слой с наибольшей вероятностью сломается. При толщине Т = 100 мкм деформация в позиции А-А наименьшая, а состояние является безопасным. Когда толщина составляет T= 140 мкм, тогда кривизна в положении B-B является наименьшей, т.е. деформация наименьшая, а состояние является наиболее безопасным. При толщине T = 100 мкм кривизна в позиции В-В является самой большой, а деформация - самой серьезной. В этом положении проводящий слой легко нарушить. Оптимизация схемы проводящего слоя: при приложении продольной деформации к небольшому участку в области изгиба обнаружено, что распределение напряжений и их пик согласуются с исходной плоской моделью. Поэтому необходимо учитывать только поперечную деформацию. При ее приложении к малому сечению конечно-элементной модели области изгиба эквивалентные номинальное и линеаризованное напряжения уменьшаются с увеличением размера ячейки сетки, причем эквивалентное напряжение понижается со 165 до 95 МПа, а линеаризованное - со 105 до 90 МПа. При понижении эквивалентных номинального и линеаризованного напряжений уменьшается возможность обрыва в проводящем слое зоны изгиба, поэтому предлагаемый метод эффективен при проектировании проводящего слоя в зоне изгиба [24]. Необходимо также дополнительно стандартизировать стратегии горизонтального сравнения различных схем маршрутизации, которые заключаются в следующем. Для различных схем маршрутизации в горизонтальном сравнении выбраны четыре конструктивных элемента, увеличен конец нагрузки смещения, уменьшен конечный эффект, принята модель симметричного плоского напряжения 1/4, предполагается та же продольная деформация. Применены метод эквивалентной линеаризации напряжений для точек риска разрушения и метод добавления мембраны. Целью оценки является суммарное напряжение при изгибе. Разработаны три вида схем проводки, преимущества и недостатки которых анализируются. При полностью заполняющем дизайне материалы расположены равномерно, а общий уровень напряжения мембраны высок. При поперечной прокладке наклонная конструкция разлагает продольное и поперечное распределение напряжений, эффективно снижает общий уровень напряжения мембраны, но возникают концентрации напряжений в локальных положениях. Во внутреннем отверстии это уменьшает использование материалов по сравнению с полностью заполненным типом, но эффективно не изменяет основных напряжений, что приводит к невозможности снижения уровня напряжений при одновременном увеличении концентрации напряжений в локальной точке. Поэтому с точки зрения оптимизации распределения напряжений, для уменьшения трещиноватости, электропроводящий слой следует прокладывать в следующем порядке: тип поперечной прокладки > тип полного заполнения > тип внутреннего отверстия. Выводы Выполнен имитационный анализ и исследованы механические свойства проводящего слоя зоны изгиба в его процессе. Основные результаты заключаются в следующем: 1. Полиэтилентерефталат (PET) с проводящей пленкой оксида индия-олова (ITO) в качестве материала подложки используется для получения слоистых компонентов. Вводятся и анализируются типичные характеристики многослойной стековой структуры FOLED с полимером, PLN, PVX, металлическим проводником, PI в качестве основного слоя стека, а также эквивалентное рабочее состояние в процессе изгиба для завершения трехмерного конечно-элементного моделирования области изгиба FOLED. 2. Моделирование процесса изгиба позволило проанализировать кумулятивный процесс пластической деформации проводящего слоя, который не может быть проконтролирован в процессе испытания, что расширило возможности последующего исследования механизма его разрушения. 3. Изучаются и обсуждаются проектные спецификации различных схем электропроводки с целью их оптимизации и сравнения различных характеристик распределения механических напряжений. 4. В соответствии с конструктивными характеристиками различных схем сравниваются их преимущества и недостатки. 5. Исследована оптимальная конструкция толщины слоя эпоксидной смолы, получены кривые отклика различных положений проводящего слоя с толщиной слоя смолы и проанализированы влияющие факторы.
Ключевые слова
FOLED,
область изгиба,
проводящий слой,
конечно-элементное моделирование,
схема проводаАвторы
Liang Ma | Jiangsu University | Master Degree, Engineer Jiangsu University | maliang72@163.com |
Jinan Gu | Jiangsu University | Doctoral Degree, Professor, Doctoral Supervisor Jiangsu University | gujina2013@163.com |
Всего: 2
Ссылки
Huan L., Xin J.L., and Xing B.C. // IEEE. Trans. Electron. Devices. - 2018. - V. 48. - P. 1-5.
Zhamaletdinov A.A., Velikhov E.P., and Shevtsov A.N. // Dokl. Earth. Sci. - 2017. - V. 474. - P. 641-645.
Kim T., Lee K.H., and Lee J.Y. // J. Mater. Chem. C. - 2018. - V. 72. - P. 88-91.
Peng W., Zhang J., and Lei D. // Chem. Mater. - 2017. - V. 29. - P. 112-118.
Ma L.L., Yang L., and Zong H.C. // Chemelectrochem. - 2017. - V. 4. - P. 1443-1449.
Cao R.R., Liu S., and Liu Q. // IEEE. Electron. Dev. Lett. - 2018. - V. 38. - P. 1371-1374.
Gruber M., Müller A., and Jovanov V. // IEEE J. Photovoltaics. - 2017. - V. 39. - P. 1-8.
Liu H.Y., Liu G.J., and Huang R.C. // IEEE. Sensors. J. - 2017. - V. 17. - No.16. - P. 5087-5092.
Jin X., Sendek A.D., and Cubuk E.D. // Acs. Nano. - 2017. - V. 11. - P. 19-70.
Yan R., Ma X.J., and Liang Y.Z. // J. China Acad. Electron Inform. Technol. - 2017. - V. 12. - P. 7-13.
Pan L., Wang H.C., and Chu M.Y. // J. Pow. Supply. - 2019. - V. 17. - P. 140-147.
Xia K., Fu X.L., and Li Z.R. // Chin. J. Pow. Source. - 2019. - V. 43. - P. 116-119.
Qu J.J. // Autom. Instrum. - 2017. - V. 15. - P. 226-227.
Sa R.N. // J. Jilin Univ. (Sci). - 2017. - V. 55. - P. 180-184.
Li Y.Y., Liang K.F., and Yuan Z.Y. // Comput. Simul. - 2017. - V. 34. - P. 55-62.
Brown T.S., Du S., Eruslu H., and Sayas F. // Appl. Math. Nonlinear Sci. - 2018. - V. 3. - P. 55-96.
Fernández-Pousa C.R. // Appl. Math. Nonlinear Sci. - 2018. - V. 3. - P. 23-32.
Xiong Z., et al. // Multimed Tools Appl. - 2019. - V. 78. - No. 22. - P. 31035-31055.
Khellat F. and Khormizi M.B. // Appl. Math. Nonlinear Sci. - 2018. - V. 3. - P. 15-22.
Lakshminarayana G., Vajravelu K., Sucharitha G., and Sreenadh S. // Appl. Math. Nonlinear Sci. - 2018. - V. 3. - P. 41-54.
Bozduman H.C. and Afacan E. // Appl. Math. Nonlinear Sci. - 2020. - V. 5 (1). - P. 479-484.
Gao N., Hou H., and Wu J.H. // Int. J. Mod. Phys. B. - 2018. - V. 32. - No. 20. - P. 1850005.
Gao N., Hou H., Zhang Y., et al. // Mod. Phys. Lett. B. - 2018. - V. 32. - No. 4. - P. 1850040.
Ruiz-Fernandez J.P., Benlloch J., Lopez M.A., et al. // Appl. Math. Nonlinear Sci. - 2019. - V. 4. - No. 1. - P. 21-34.