К определению двумерного закона изменения плотности в функционально-градиентной упругой пластине

1. Введение. Развитие новых технологий производства позволяет создавать новые неоднородные материалы, свойства которых изменяются по пространственным координатам. К таким материалам относятся в том числе функционально-градиентные материалы (ФГМ). Благодаря переменным свойствам они находят применение в различных областях науки и техники, таких как автомобилестроение и машиностроение (конструкционные элементы), аэрокосмическая промышленность (защитные покрытия), биомедицина (протезы и импланты) и т.д. [1$–$3]. Основным преимуществом таких материалов по сравнению с классическими композитами в виде сэндвич-панелей является более низкая вероятность возникновения расслоения и трещин [4, 5]. Подобные дефекты могут возникать на этапе производства, а также в результате воздействия эксплуатационных механических или температурных нагрузок.

В настоящее время [6, 7] широко используются многофазные ФГМ, которые описываются переменными свойствами, позволяющими осуществлять переход от одного материала в другой по всему объему. Они обычно создаются в процессе порошковой металлургии, когда требуемый состав порошка укладывается слой за слоем с последующим уплотнением и спеканием. В свою очередь, формирование фаз ФГМ зависит от условий их обработки (скорости охлаждения и процесса термообработки). Поскольку процесс производства является сложным, то следует проверять соответствие проектных свойств реальным.

Одной из важных задач является анализ акустических свойств неоднородных пластин, которые широко используются как конструкционные элементы. В статье [8] исследована задача о свободных колебаниях функционально-градиентной пластины с круговым отверстием. Переменные свойства изменяются по толщине по правилу осреднения Мори$–$Танаки. Решение построено с помощью метода Рэлея$–$Ритца. Получены значения собственных частот, исследовано влияние механических свойств пластины и размера отверстия на компоненты полей перемещений. Равновесие тонкой упругой неоднородной пластины рассмотрено в работе [9]. Определяющие уравнения получены с помощью асимптотического подхода. Проведено сравнение данных, полученных на основе моделей Кирхгофа$–$Лява и Тимошенко$–$Рейсснера. В статье [10] с помощью метода суперпозиции получены аналитические решения для задачи о свободных колебаниях толстой ортотропной пластины. Сформулированная краевая задача сведена к решению бесконечной системы линейных алгебраических уравнений. Рассмотрены пластинки со свободными и жестко закрепленными границами. В статье [11] на основе метода динамической жесткости исследована задача о свободных колебаниях пластины. Свойства материала изменяются по толщине согласно степенному закону. Проведен анализ различных геометрических и материальных параметров на значения собственных частот. В работе [12] получено решение трехмерной задачи теории упругости для однонаправленно армированных композитных пластин в виде рядов. Построены графики смещений и напряжений. Показано, что значительное увеличение жесткости может быть достигнуто благодаря выбору специальной ориентации волокон и значения градиента свойств. В монографии [13] рассмотрены задачи для пластин и оболочек из ФГМ в рамках нелинейной теории упругости.

Помимо изгибных колебаний пластин продолжаются исследования их вынужденных и свободных планарных колебаний. В статье [14] представлен обзор более 230 работ, посвященных разработкам аналитических решений для задач о собственных колебаниях прямоугольных пластин. При этом такие решения могут быть получены только для определенных типов граничных условий [15$–$19]. В статье [20] на основе метода Чебышева$–$Ритца рассмотрены свободные планарные колебания слоистых композитных пластин с криволинейными границами. Проведен анализ влияния структуры, геометрии и вида граничных условий на значения собственных частот. Некоторые результаты сравниваются с данными, полученными с помощью МКЭ. В работе [21] рассматривается задача о планарных колебаниях неоднородной кольцевой области. Свойства материала изменяются по радиальному направлению. Решение построено с помощью усовершенствованного метода Фурье. Представлены значения собственных частот для различных соотношений внешнего и внутреннего радиусов, параметра градиентности и условий закрепления. Точность расчетов продемонстрирована путем сравнения с результатами для однородной области [22]. В работе [23] проведено исследование планарных колебаний ортотропных полых, сплошных круговых пластин и пластин в виде секторов окружностей. В статье [24] с помощью полуаналитического решения исследованы планарные свободные колебания неоднородных пористых пластин. Рассмотрены два вида неоднородности, когда свойства изменяются по толщине или в плоскости пластины. Проведено сравнение некоторых данных с результатами натурного эксперимента. Среди других работ, в которых представлены экспериментальные данные для однородных пластин, можно отметить статьи [25$–$27].

Изучение задач идентификации характеристик упругих функционально-градиентных пластин нашло отражение в современных работах. Статьи [28, 29] посвящены неразрушающей идентификации остаточных напряжений в тонких ФГ пластинах. В первой работе в рамках акустического метода решение обратной задачи сведено к реализации итерационного процесса, представлены численные результаты восстановления. Во второй работе поле остаточных напряжений описывается с помощью функции напряжения Эри. Решение также сведено к итеративной процедуре, основанной на модальном анализе чувствительности. Исследовано влияние зашумления на процедуру реконструкции. Идентификация жесткости пластин осуществлена в работах [30, 31], а в статьях [32$–$34] представлены результаты восстановления нескольких механических характеристик. Например, в работе [32] проводится реконструкция цилиндрической жесткости и плотности неоднородной круглой пластины по дополнительной информации о ее акустическом отклике. Решение построено на основе проекционного метода. В статье [31] предложена процедура идентификации направлений ортотропии и изгибной жесткости. Также авторами осуществлена регуляризация задачи для преодоления высокой чувствительности к зашумлению входных данных.

При исследовании обратных коэффициентных задач в рамках линейной теории упругости в качестве входной информации обычно используются данные, измеряемые внутри области (первая постановка) или на части ее поверхности (вторая постановка) [35]. С практической точки зрения задачи во второй постановке более удобны в реализации при проведении эксперимента, но требуют использования специальных методов, так как являются существенно нелинейными. Для них также характерно более выраженное проявление некорректности. Идентификация искомых двумерных законов изменения из широкого класса функций требует существенных вычислительных затрат [36, 37]. При этом одним из способов снижения таких затрат является сужение класса искомых законов, например до полиномиальных, и задача сводится к отысканию конечного числа параметров [38]. Для задач в первой постановке на практике не для всех исследуемых объектов возможно реализовать съем необходимых данных. Такие задачи являются линейными, поэтому для их решения обычно не требуется использования сложных методов и привлечения специальных численных схем. С другой стороны, их рассмотрение позволяет выявить важные особенности для исследования задач во второй постановке.

В настоящей работе рассматривается задача об установившихся планарных колебаниях функционально-градиентной упругой пластины. Механические свойства описываются модулем Юнга, коэффициентом Пуассона и плотностью, изменяющимися по двум координатам. Влияние переменных свойств на акустические характеристики анализируется на основе ранее апробированного подхода с использованием метода конечных элементов (МКЭ). Сформулирована и решена обратная задача об определении закона изменения плотности по данным о поле смещений, измеренном в конечном наборе точек при фиксированной частоте.

2. Прямая задача. Общая постановка задачи об установившихся колебаниях изотропного тела из неоднородного материала может быть записана в виде [39, 40]:

 $\begin{array}{c}\nabla \cdot s+\rho {\omega }^{2}u=\mathrm{0,}\\ \sigma =\lambda I\text{tr}e+2\mu e,e=0.5(\nabla u+\nabla u^T),\\ {\left.u\right|}_{S_u}\text{}=\mathrm{0,}{\left.n\cdot \sigma \right|}_{S_{\sigma }}\text{}=P,\end{array}$                                                                        (2.1)

где σ – тензор напряжений Коши, ρ – переменная плотность, ω – круговая частота колебаний, u – вектор смещений, μ и λ – параметры Ламе, изменяющиеся по пространственным координатам, ε – тензор деформации, I – единичный тензор, n – единичный вектор внешней нормали, S = S_u ◡ S_σ – поверхность тела, P – амплитуда вектора внешней нагрузки.

На основе этой постановки рассмотрим задачу об установившихся планарных колебаниях прямоугольной пластины малой толщины 2h в рамках обобщенного плоского напряженного состояния. Левая грань жестко закреплена, на правой приложена нагрузка. Ширина пластины $—$ b, длина $—$ l. Срединную плоскость пластины примем за координатную плоскость Ox₁x₂, тогда нулевыми являются следующие компоненты тензора напряжений σ₁₃ = σ₂₃ = σ₃₃ = 0. Компоненты вектора смещений u₁, u₂ являются функциями двух координат u₁ = u₁(x₁, x₂), u₂ = u₂(x₁, x₂).

Пластина изготовлена из ФГМ, переменные свойства которого описываются модулем Юнга E(x₁, x₂), коэффициентом Пуассона v(x₁, x₂) и плотностью ρ(x₁, x₂). Уравнения установившихся колебаний и граничные условия после отделения временного множителя примут вид:

 $\begin{array}{l}{\sigma }_{\mathrm{11,1}}+{\sigma }_{\mathrm{21,2}}+\rho {\omega }^2{u}_1=\mathrm{0,}\\ {\sigma }_{\mathrm{12,1}}+{\sigma }_{\mathrm{22,2}}+\rho {\omega }^2{u}_2=\mathrm{0,}\\ u_1=u_2=\mathrm{0,}{x}_1=\mathrm{0,}{x}_2\in \left[\mathrm{0,}b\right],\\ {\sigma }_{11}=P,{\sigma }_{12}=\mathrm{0,}{x}_1=l,x_2\in \left[\mathrm{0,}b\right],\\ {\sigma }_{21}=\mathrm{0,}{\sigma }_{22}=\mathrm{0,}{x}_1\in \left[\mathrm{0,}l\right],x_2=\mathrm{0,}b,\end{array}$                                                                      (2.2)

где ${\sigma }_{11}=\frac{E}{1-{\nu }^2}\left({\epsilon }_{11}+\nu {\epsilon }_{22}\right),$ ${\sigma }_{22}=\frac{E}{1-{\nu }^2}\left({\epsilon }_{22}+\nu {\epsilon }_{11}\right),$ ${\sigma }_{12}={\sigma }_{21}=\frac{E}{1+\nu }{\epsilon }_{12},$ ${\epsilon }_{11}=u_{\mathrm{1,1}},$ ${\epsilon }_{22}=u_{\mathrm{2,2}},$ ${\epsilon }_{12}={\epsilon }_{21}=0.5(u_{\mathrm{1,2}}+u_{\mathrm{2,1}}).$ Здесь для краткости записи принято обозначение частной производной $\partial F/\partial x_i=F_{,i}.$

Для удобства рассмотрения перейдем к безразмерной постановке задачи:

 $\begin{array}{l}{S}_{\mathrm{11,1}}+S_{\mathrm{21,2}}+{\kappa }^{2}r{w}_1=\mathrm{0,}\\ S_{\mathrm{12,1}}+S_{\mathrm{22,2}}+{\kappa }^{2}r{w}_2=\mathrm{0,}\\ w_1=w_2=\mathrm{0,}{\xi }_1=\mathrm{0,}{\xi }_2\in \left[\mathrm{0,}k\right],\\ S_{11}=p,S_{12}=\mathrm{0,}{\xi }_1=\mathrm{1,}{\xi }_2\in \left[\mathrm{0,}k\right],\\ S_{21}=\mathrm{0,}{S}_{22}=\mathrm{0,}{\xi }_1\in \left[\mathrm{0,1}\right],{\xi }_2=\mathrm{0,}k.\end{array}$                                                                     (2.3)

Определяющие соотношения в безразмерном виде:

 $\begin{array}{c}{S}_{11}=\frac{g}{1-{\nu }^2}({\epsilon }_{11}+\nu {\epsilon }_{22}),S_{22}=\frac{g}{1-{\nu }^2}({\epsilon }_{22}+\nu {\epsilon }_{11}),S_{12}=S_{21}=\frac{g}{1+\nu }{\epsilon }_{12},\\ {\epsilon }_{11}=w_{\mathrm{1,1}},{\epsilon }_{22}=w_{\mathrm{2,2}},{\epsilon }_{12}={\epsilon }_{21}=0.5(w_{\mathrm{1,2}}+w_{\mathrm{2,1}}),\end{array}$              (2.4)

где $x_i=l{\xi }_i$, ${\xi }_1\in [\mathrm{0,1}]$, ${\xi }_2\in [\mathrm{0,}k]$, ${\sigma }_{ij}=E_{\text{*}}{S}_{ij},$ $u_i=l{w}_i$, $E=E_{\text{*}}g\left({\xi }_1,{\xi }_2\right),$ $\nu =\nu \left({\xi }_1,{\xi }_2\right),$ $\rho ={\rho }_{\text{*}}r\left({\xi }_1,{\xi }_2\right),$ E_*, ρ_* $—$ характерные значения соответствующих величин, k² = ρ_*ω²l²/E_* $—$ безразмерная частота колебаний, k = b/l $—$ параметр, характеризующий отношение ширины пластины к ее длине. С учетом линейности задачи значение безразмерной нагрузки выбрано следующим p = P/E_* = 1.

3. Анализ чувствительности. Следует отметить, что для общего случая неоднородности механических характеристик и вида граничных условий решение рассматриваемой задачи может быть построено только численно. Для построения такого решения воспользуемся МКЭ, который является одним из универсальных методов исследования и моделирования поведения объектов конечных размеров. Отметим, что ранее в работах [41, 42] был проведен анализ влияния переменных свойств материала на основные акустические характеристики радиально-неоднородного цилиндра и дана сравнительная оценка точности решений, полученных в КЭ-пакете FlexPDE и с помощью метода, основанного на разложении в ряды. В данной работе в качестве конкретного примера будем рассматривать ФГМ, состоящий из смеси четырех веществ.

Таблица 1. Значения параметров четырех веществ

Переменные свойства такого материала могут быть описаны с помощью двумерных законов смеси [43]:

 $\begin{array}{c}{f}_E\left({\xi }_1,{\xi }_2\right)=E_1{V}_1\left({\xi }_1,{\xi }_2\right)+E_2{V}_2\left({\xi }_1,{\xi }_2\right)+E_3{V}_3\left({\xi }_1,{\xi }_2\right)+E_4{V}_4\left({\xi }_1,{\xi }_2\right),\\ f_{\nu }\left({\xi }_1,{\xi }_2\right)={\nu }_1{V}_1\left({\xi }_1,{\xi }_2\right)+{\nu }_2{V}_2\left({\xi }_1,{\xi }_2\right)+{\nu }_3{V}_3\left({\xi }_1,{\xi }_2\right)+{\nu }_4{V}_4\left({\xi }_1,{\xi }_2\right),\\ f_{\rho }\left({\xi }_1,{\xi }_2\right)={\rho }_1{V}_1\left({\xi }_1,{\xi }_2\right)+{\rho }_2{V}_2\left({\xi }_1,{\xi }_2\right)+{\rho }_3{V}_3\left({\xi }_1,{\xi }_2\right)+{\rho }_4{V}_4\left({\xi }_1,{\xi }_2\right),\end{array}$                       (3.1)

где введены функции

 $\begin{array}{c}{V}_1\left({\xi }_1,{\xi }_2\right)=\left(1-{\xi }_1^n\right)\left(1-{\left(\frac{{\xi }_2}{k}\right)}^m\right),V_2\left({\xi }_1,{\xi }_2\right)=\left(1-{\xi }_1^n\right){\left(\frac{{\xi }_2}{k}\right)}^m\text{},\\ V_3\left({\xi }_1,{\xi }_2\right)={\xi }_1^n\left(1-{\left(\frac{{\xi }_2}{k}\right)}^m\right),V_4\left({\xi }_1,{\xi }_2\right)={\xi }_1^n{\left(\frac{{\xi }_2}{k}\right)}^m\text{}.\end{array}$                                 (3.2)

Положительные параметры n и m характеризуют степень градиентности свойств по координатам ξ₁ и ξ₂ соответственно; E_i, v_i, ρ_i, i = $\overline{1.4}$ $—$ безразмерные значения модуля Юнга, коэффициента Пуассона и плотности соответствующего i-го вещества.

В качестве четырех веществ возьмем следующие: керамики SiC (вещество 1), Al₂O₃ (вещество 2), которые достигают 100% своей концентрации соответственно в точках (0,0) и (0,k) на левой грани области, и металлы Ti6Al4V (вещество 3), Al 1100 (вещество 4), которые достигают 100% своей концентрации соответственно в точках (1,0) и (1,k) на правой границе. Значения модуля Юнга E, коэффициента Пуассона n и плотности r взяты из открытых данных, размещенных на сайте http://www.matweb.com (см. табл. 1). Характерные значения модуля Юнга и плотности выбраны следующими: E_* = 100 ГПа, ρ_* = 1000 кг/м³.

Проведем анализ влияния амплитудных значений модуля Юнга, коэффициента Пуассона и плотности на изменение компонент поля перемещений w_i. Для этого, по аналогии с работой [42], введем в рассмотрение функции:

 $\begin{array}{c}{g}_d\left({\xi }_1,{\xi }_2\right)=\left(1+d\right)\cdot f_E\left({\xi }_1,{\xi }_2\right),\\ {\nu }_d\left({\xi }_1,{\xi }_2\right)=\left(1+d\right)\cdot f_{\nu }\left({\xi }_1,{\xi }_2\right),\\ r_d\left({\xi }_1,{\xi }_2\right)=\left(1+d\right)\cdot f_{\rho }\left({\xi }_1,{\xi }_2\right),\end{array}$                                                                                (3.3)

где параметр d = 0.001, 0.01, 0.1 описывает амплитудные значения. Для исследования влияния конкретной характеристики соответствующий ей закон выбирался из (3.2) при n = m = 1, а два других фиксировались при значении d = 0. Для удобства оценки изменения значений компонент поля перемещений введены функции:

 ${\delta }_{{\xi }_1}=\frac{\left|w_1-w_1^0\right|}{\underset{M}{\text{max}}\left|w_1^0\right|}\cdot 100\text{\%}, {\delta }_{{\xi }_2}=\frac{\left|w_2-w_2^0\right|}{\underset{M}{\text{max}}\left|w_2^0\right|}\cdot 100\text{\%},$                                                 (3.4)

где $M=\{({\xi }_1,{\xi }_2)\in R^2 | {\xi }_1\in [\mathrm{0,1}],{\xi }_2\in [\mathrm{0,}k]\}$ $—$ область рассмотрения, $w_1^0$, $w_2^0$  $—$ компоненты поля перемещений, соответствующего законам (3.3) при d = 0, n = m = 1. Анализ результатов расчетов показал, что наибольшее изменение поля перемещений наблюдалось при изменении амплитудных значений модуля Юнга E, наименьшее $—$ коэффициента Пуассона n.

Рис. 1. График функции δ_ξ₁  для d = 1, n = m = 1, соответствующие амплитудному изменению функции r_d.

Рис. 2. График функции δ_ξ₁  для d = 1, n = m = 1, соответствующие амплитудному изменению функции r_d .

На рис. 1, 2 для примера изображены графики функций  δ_ξ₁, δ_ξ₂ соответствующие функции r_d(ξ₁,ξ₂), d = 0.1. Максимальные значения max δ_ξ₁ = 8.6%, δ_ξ₂ = 6.94%. Для приведенных данных значение безразмерной частоты выбиралось до первого резонанса κ = 1.0. Также были проанализированы аналогичные графики для g_d(ξ₁,ξ₂) (max δ_ξ₁= 17.8%, max δ_ξ₂= 12.6%) и v_d(ξ₁,ξ₂) (max δ_ξ₁= 0.12%, max δ_ξ₂= 7.77%). Для демонстрации изменения значений первой безразмерной резонансной частоты κ в зависимости от значений d для каждой из характеристик в табл. 2 приведены соответствующие данные. При d = 0 и n = m = 1 первый резонанс возникает при κ ≈ 1.39744. Следует отметить, что здесь точность вычислений в пакете FlexPDE задавалась с помощью служебного параметра errlim = 10^$–$6.

Таблица 2. Значения первой безразмерной резонансной частоты κ при различных значениях параметра d и одной из изменяемых функций g_d, ν_d, r_d

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

Для выявления влияния функции изменения плотности выберем законы g и v в виде линейных функций $g\left({\xi }_1,{\xi }_2\right)=f_E\left({\xi }_1,{\xi }_2\right),$ $\nu \left({\xi }_1,{\xi }_2\right)=f_{\nu }\left({\xi }_1,{\xi }_2\right)$ (n = m = 1), а r $–$ в виде кубической функции $r\left({\xi }_1,{\xi }_2\right)=f_{\rho }\left({\xi }_1,{\xi }_2\right)$ (n = m = 3). На рис. 3, 4 показаны графики δ_ξ₁, δ_ξ₂, которые были получены для такого модельного вычислительного эксперимента. Максимальные значения этих функций (max δ_ξ₁ = 4.49%, max δ_ξ₂= 16.1%) показывают, что закон изменения безразмерной функции плотности $r\left({\xi }_1,{\xi }_2\right)$ также может оказывать существенное влияние на поле перемещений.

4. Исследование обратной задачи. Одними из наиболее востребованных методов диагностики свойств материала являются неразрушающие подходы [35]. Рассмотрим коэффициентную обратную задачу (ОЗ) о реконструкции закона изменения плотности. В качестве входной информации используются данные о компонентах поля перемещений, измеренные в конечном наборе точек области M при фиксированной частоте κ [35, 44]. При этом предполагается, что законы изменения модуля Юнга и коэффициента Пуассона являются известными. Например, они могут быть определены из эксперимента на растяжение статической нагрузкой по данным о деформациях, полученных с помощью оптического метода корреляции цифровых изображений (digital image correlation) (см., например, его применение в работах [45$–$47]). Следует отметить, что в некоторых прикладных теориях для ФГМ коэффициент Пуассона полагается постоянным [48, 49].

Рис. 3. График функции δ_ξ₁ для f_ρ при n = m = 3.

Рис. 4. График функции δ_ξ₂ для f_ρ при n = m = 3.

Обратимся к векторному уравнению движения (2.1). Видно, что относительно r оно является линейным и имеет единственное решение для ненулевых полей смещения u≠0:

 ${\rho }_{rec}=-\frac{\left(\nabla \cdot \sigma \right)\cdot u}{{\omega }^2{\left|u\right|}^2}.$                                                                                                                 (4.1)

Запись в таком общем виде позволяет исследовать соответствующие одномерные, двумерные и трехмерные обратные задачи. Для рассматриваемого случая полученное решение может быть записано в виде:

 $r_{rec}=-\frac{\left(S_{\mathrm{11,1}}+S_{\mathrm{21,2}}\right)\cdot w_1+\left(S_{\mathrm{12,1}}+S_{\mathrm{22,2}}\right)\cdot w_2}{{\kappa }^2\cdot \left(w_1^2+w_2^2\right)}$                                                                 (4.2)

и использовано для всех точек области, где $w_1^2+w_2^2\ne 0.$

Докажем, что полученное решение удовлетворяет обоим уравнениям движения. Подставим (4.2) в уравнения (2.3):

 $\begin{array}{c}{S}_{\mathrm{11,1}}+S_{\mathrm{21,2}}-{\kappa }^2\cdot \frac{\left(S_{\mathrm{11,1}}+S_{\mathrm{21,2}}\right)\cdot w_1+\left(S_{\mathrm{12,1}}+S_{\mathrm{22,2}}\right)\cdot w_2}{{\kappa }^2\cdot \left(w_1^2+w_2^2\right)}\cdot w_1=\mathrm{0,}\\ S_{\mathrm{12,1}}+S_{\mathrm{22,2}}-{\kappa }^2\cdot \frac{\left(S_{\mathrm{11,1}}+S_{\mathrm{21,2}}\right)\cdot w_1+\left(S_{\mathrm{12,1}}+S_{\mathrm{22,2}}\right)\cdot w_2}{{\kappa }^2\cdot \left(w_1^2+w_2^2\right)}\cdot w_2=\mathrm{0,}\end{array}$                               (4.3)

После упрощения получаем:

 $\begin{array}{c}{w}_2\cdot \left(\left(S_{\mathrm{11,1}}+S_{\mathrm{21,2}}\right)\cdot w_2-\left(S_{\mathrm{12,1}}+S_{\mathrm{22,2}}\right)\cdot w_1\right)=\mathrm{0,}\\ w_1\cdot \left(\left(S_{\mathrm{12,1}}+S_{\mathrm{22,2}}\right)\cdot w_1-\left(S_{\mathrm{11,1}}+S_{\mathrm{21,2}}\right)\cdot w_2\right)=0.\end{array}$                                                             (4.4)

Покажем, что $\left(S_{\mathrm{11,1}}+S_{\mathrm{21,2}}\right)\cdot w_2-\left(S_{\mathrm{12,1}}+S_{\mathrm{22,2}}\right)\cdot w_1$, присутствующее в обоих выражениях, равно нулю. Для этого достаточно обратиться к исходной системе (2.4), умножить первое уравнение на w₂, второе на w₁ и вычесть одно из другого.

При выборе частоты колебаний до первого резонанса с учетом вида нагружения w₁ будет принимать нулевые значения только на защемленной грани. При этом w₂ помимо левой грани будет обращаться в ноль еще и в других точках рассматриваемой области. Таким образом, при выборе частоты до первого резонанса выражение $w_1^2+w_2^2$, стоящее в знаменателе формулы (4.2), будет обращаться в ноль только в точках, которые принадлежат жестко защемленной грани.

Основная сложность решения (4.2) состоит в численной реализации, а именно $—$ вычислении первых и вторых производных от таблично-заданных функций w₁, w₂. Как известно, задача вычисления производных является в общем случае некорректной [35]. Например, в статье [50] для преодоления этой проблемы решение схожей задачи об определении параметров Ламе неоднородной пластины было построено на основе слабой постановки и процедуры интегрирования по частям. В настоящей работе на основе алгоритма локально взвешенной регрессии LOWESS (Locally Weighted Scatterplot Smoothing) и сплайн-аппроксимации предложен следующий подход, который моделирует численный эксперимент по решению задачи, и состоит из нескольких этапов:

1.  Расчет поля перемещений, соответствующего заданному точному закону изменения искомой функции r_ex в пакете FlexPDE.
2.  Моделирование входной информации путем вычисления значений функций w₁, w₂ в точках (ξ_₁i,ξ₂_j), ξ_₁i = (i – 1)h_ξ₁, h_ξ₁ =1/(n_ξ₁– 1), i  = 1...n_ξ₁, ξ_2i = (i  – 1)h_ξ₂, h_ξ₂ =1/(m_ξ₂ – 1), j  = 1...m_ξ₂.
3.  Построение интерполяционных функций $w_1^s$, $w_2^s$ по соответствующим входным данным с помощью радиальной базисной функции.
4.  Определение значений компонент тензора напряжений согласно закону Гука (2.4) в узловых точках (δ_ξ₁i, δ_ξ₂j).
5. Для “сглаживания” полученных значений компонент тензора напряжений используется алгоритм локально взвешенной регрессии LOWESS (Locally Weighted Scatterplot Smoothing) [51], реализованный в пакете Maple.
6.  Далее по значениям в узловых точках (δ_ξ₁i, δ_ξ₂j) строятся интерполяционные функции $S_{11}^s$, $S_{12}^s$и численно определяются первые производные S  ^{$S_{11.1}^s$, $S_{12.2}^s$}.
7.  С помощью полученных функций согласно решению (4.2) вычисляются значения искомого закона изменения плотности r_rec.

Для проведения вычислительных экспериментов в качестве законов изменения модуля Юнга и коэффициента Пуассона выбраны константы g = 2.407, v = 0.258, которые получены путем осреднения линейных законов (3.1) по площади области M. Для удобства положим h_ξ₁ = h_x2, тогда при заданном параметре k = 0.5 число узловых точек по каждой из координат будет удовлетворять равенству n_ξ₁ = 2m_ξ₂. Здесь при проведении расчетов выбраны n_ξ₁ = 2m_ξ₂= 40. Поскольку наибольшее влияние переменных свойств материала на акустические характеристики проявляется вблизи резонанса, то значение параметра k = 1.2 выбрано вблизи первого резонанса. Выбор такой частоты обусловлен более удобной реализацией на практике, более существенным влиянием переменных свойств на поля смещения и отсутствием дополнительных точек, в которых $w_1^2+w_2^2= 0$ за исключением точек на защемленной грани. Значения параметров алгоритма LOWESS во всех приведенных примерах были заданы следующими: fitorder = 1 (степень полинома, используемого в каждой локальной регрессии), bandwidth = 0.05 (доля точек входных данных, используемых в каждой локальной регрессии), iters = 2 (количество итераций при сглаживании данных одной независимой переменной). Значение параметра bandwidth выбиралось близким к минимальному с учетом величины шага разбиения области (количества узлов сетки).

Рис. 5. Графики сплайнов компонент тензора напряжений: (a) S^s₁₁; (b) S^s₁₂.

Поскольку переменные свойства ФГМ могут быть описаны с помощью различных видов законов [52], то далее представлены примеры вычислительных экспериментов, проведенных для нескольких функций.

Пример 1. Восстанавливаемая функция $—$ квадратичная r_ex(ξ₁,ξ₂) = f_r(ξ₁,ξ₂) (см. закон (5) при n = m = 2). На рис. 5 показаны графики функций S ^s₁₁ и S ^s₁₂, которые имеют особенности у границ области. Поэтому здесь и далее реконструкция осуществляется на подобласти M^* = {(ξ₁,ξ₂) ∊ R²|ξ₁ ∊ [0.1,0.9], ξ₂ ∊ [0.1,0.4]}.

Рис. 6. (a) графики точной квадратичной функции r_ex(ξ₁,ξ₂) = f_ρ(ξ₁,ξ₂) при n = m = 2 (сплошная поверхность) и восстановленной функции r_rec(ξ₁,ξ₂) (сетка); (b) график функции относительной погрешности δ_ρ(ξ₁,ξ₂).

Рис. 7. (a) графики точной экспоненциальной функции r_ex(ξ₁,ξ₂) = 0.5(e^ξ₁ + e^4ξ) (сплошная поверхность) и восстановленной функции r_rec(ξ₁,ξ₂) (сетка); (b) график функции относительной погрешности δ_ρ(ξ₁,ξ₂).

Рис. 8. (a) графики точной тригонометрической функции r_ex(ξ₁, ξ₂) = 3.5 – 0.5cos(1.2π(ξ₁+ξ₂)) (сплошная поверхность) и восстановленной функции r_rec(ξ₁, ξ₂) (сетка); (b) график функции относительной погрешности δ_ρ(ξ₁, ξ₂).

На рис. 6$–$8, демонстрирующих точность реконструкции, сплошной поверхностью обозначается график точного закона r_ex(ξ₁,ξ₂), а в виде черной сетки $—$ восстановленный закон r_rec(ξ₁,ξ₂). Также для удобства анализа распределения относительной погрешности реконструкции приведены графики функции δ_ρ:

 ${\delta }_{\rho }=\frac{\left|r_{ex}-r_{rec}\right|}{\underset{M}{\text{max}}\left|r_{ex}\right|}\cdot 100\text{\%}.$                                                                                                             (4.5)

На рис. 6 представлены графики точной и восстановленной квадратичной функции, погрешность не превосходит 7.1%.

Пример 2. Восстанавливаемая функция $—$ экспоненциальный закон r_ex(ξ₁,ξ₂) = 0.5(e^ξ₁ + e^4ξ₂). На рис. 7 показаны соответствующие графики, max δ_ρ ≤ 7.8%.

Пример 3. Восстанавливаемая функция $—$ тригонометрический закон r_ex (ξ₁,ξ₂) = 3.5 $–$ 0.5cos(1.2π(ξ₁ + 2ξ₂)). Результаты восстановления представлены на рис. 8, относительная погрешность не превосходит 7.7%.

Как известно, для рассматриваемой задачи имеется особенность численного решения в окрестности заделки. Для уменьшения влияния этой особенности при моделировании входных данных на основе решения прямой задачи в представленных вычислительных экспериментах была увеличена плотность расчетной сетки в окрестности заделки. Для этого был использован служебный параметр mesh_density в пакете FlexPDE, позволяющий сгущать сетку по заданному правилу, в частности было задано mesh_density= 50/(ξ₁ + 0.1).

Поскольку анализ чувствительности показал влияние различных законов изменения плотности на значения безразмерной компоненты поля перемещения w₁ менее 5%, то исследование влияния погрешности входных данных проведено для уровня шума 0.1%, изменяющегося по закону $w_i^{noise}({\xi }_{1i},{\xi }_{2j})=w_i^{exact}({\xi }_{1i},{\xi }_{2j})(1+0.001\cdot \Delta ),$ где Δ∊ [$–$1,1] $—$ случайная величина. Погрешность реконструкции для функции из примера 1 составила не более 9.96%, примера 2 $—$ 14.68%, примера 3 $—$ 11.9%. При значении параметра bandwidth = 0.1 погрешность реконструкции на удалении от заделки для функции из примера 1 составила не более 4.65%, примера 2 $—$ 5.53%, примера 3 $—$ 5.42%. Такая нелинейная зависимость погрешности реконструкции от погрешности входных данных характерна для некорректных задач.

Также был проведен анализ влияния количества значений входных данных n_ξ₁ и m_ξ₂ на результаты реконструкции двумерного закона плотности r(ξ₁, ξ₂). Анализ показал, что для успешной идентификации с погрешностью менее 10% требуется снимать смещения не менее чем в 40 узловых точках по координате ξ₁ и 20 узловых точках по ξ₂. При повышении количества узловых значений входных данных до n_ξ₁ = 100 и m_ξ₂ = 50 погрешность реконструкции составляет не более 0.5% на рассматриваемой области M^*.

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

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

Благодарности. Исследование выполнено за счет гранта Российского научного фонда $№$ 18-71-10045, https://rscf.ru/project/18-71-10045/, в Южном федеральном университете.