К развитию приближенной теории прямой и обратной задач для неоднородной прямоугольной области

В литературе по механике материалов опубликовано большое количество статей, в которых представлены результаты теоретических и экспериментальных исследований функционально-градиентных (ФГ) материалов и конструктивных элементов из них. В последнее время для таких структур все активнее развиваются теории связанных полей [1] . Однако эти исследования чаще всего ограничиваются прямым конечно-элементным моделированием с последующим анализом прочностных и термомеханических характеристик, электрического потенциала. Вместе с тем область обратных задач (ОЗ) идентификации их реальных переменных свойств, несмотря на острую востребованность, изучена недостаточно [2] .

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

В качестве основных и немногочисленных методик изучения ОЗ для материалов сложной структуры отметим следующие: наноиндентирование, фотоакустику, построение изображений с помощью диэлектрических наночастиц, нано-оптические методы, компьютерную томографию (на основе построения 3D изображений), методы рентгеновской и нейтронной дифракции, различные разрушающие методы. Большинство неразрушающих подходов является дорогостоящим, требующим значительных временных затрат, и позволяет определять искомые переменные свойства в материалах сложной структуры лишь в приповерхностном слое.

Зачастую при индентировании принимаются гипотезы об однородности и изотропности образцов, поэтому напряженно-деформированное состояние не зависит от глубины вдавливания. Авторы [6] провели эксперименты для оценки механических свойств различных фаз неоднородных материалов. Их метод заключается в использовании большой сетки наноиндентирования на репрезентативной поверхности в сочетании со статистическим анализом данных. В рамках такой методики авторы [7, 8] оценили свойства нескольких

Статья опубликована в открытом доступе по лицензии CC BY 4.0

неоднородных материалов различной микроструктуры. Обнаружено, что небольшие фазовые фракции могут быть скрыты другими фазами, что затрудняет их точное обнаружение. Данные сетки наноиндентирования также применялись для построения карт механических свойств на поверхности образца [9, 10] .

В [11] авторами предложен метод решения ОЗ распределения тепла в ФГ материалах, основанный на функциях гомогенизации. Для краевой задачи вводится семейство функций усреднения и применяется метод их суперпозиции. В рамках этого метода ОЗ определения источника тепла решаются напрямую, путем расчета линейной матричной системы. Эта схема не включает в себя создание сетки, численное интегрирование, итерационные процессы, регуляризацию и построение фундаментальных решений.

При производстве ФГ структур, в силу специфики технологии изготовления, появляется ряд проблем, часто связанных с несоответствием проектируемых и достигнутых реальных свойств. В [12] проведен обзор технологий создания и обработки ФГ материалов с помощью порошковой металлургии, искрового плазменного спекания, центробежного литья, аддитивного производства и электрофоретического процесса. В статье подробно описываются этапы изготовления различных ФГ материалов, а также приводится анализ получаемых на выходе реальных свойств путем применения различных дорогостоящих способов, включая рентгеновский анализ для оценки микроструктуры материала в приповерхностном слое, а также оптическую микроскопию; зачастую обнаруживается присутствие «побочных» материальных фаз, несоответствие объемных долей спроектированным, образующихся вследствие особенностей технологии изготовления ФГ материалов. Помимо этого, в [13] отмечается образование трещин из-за появления побочных материальных фаз. Значения твердости получаемых ФГ материалов в различных производственных процессах, показаны в таблице в работе [12] .

Одним из способов минимизации критических напряжений, улучшения свойств покрытий и предотвращения растрескивания за счет снижения различия в коэффициенте теплового расширения покрытий и подложек является использование слоев с изменяющимся составом. Для производства ФГ термобарьерных покрытий между керамическим верхним покрытием и металлическим связующим покрытием добавляется несколько смешанных металлокерамических промежуточных слоев [14] . Как известно, такие структуры могут снизить концентрацию напряжений на границах раздела и решить проблему короткого срока их эксплуатации. При этом для вязкоупругих материалов с переменными свойствами решение задач в основном производится численно. Например, в случае трехслойной композитной пластины из вязкоупругого ФГ материала в статье [15] рассмотрен специальный метод расчета колебаний на основе метода Фурье–Ритца. Получены решения в виде собственных частот и собственных форм для различных типов нагружения и граничных условий.

Отдельного внимания заслуживают работы, посвященные двумерным задачам для ФГ материалов. Авторы [16] использовали балочную теорию Тимошенко для изучения колебаний наностержня, изготовленного из ФГ материала со свойствами, изменяющимися в двух направлениях. В [17] рассмотрены колебания 2D ФГ наностержня с помощью балочной теории Бернулли–Эйлера и нелокальной теории упругости. В [18] для исследования распространения трещин в ФГ двумерных структурах осуществляется моделирование фазового поля методом конечных элементов (МКЭ). В соответствии с правилом смеси Фойгта, свойства ФГ материалов в конструкции непрерывно меняются в двух направлениях.

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

В работах [19 –21] приведены некоторые теоретические результаты решения ОЗ определения коэффициентов Ламе λ и µ по неполным данным о полях перемещений на поверхности упругого тела. Подвергнуты обсуждению вопросы существования и единственности решения ОЗ, доказаны теоремы и леммы, представлены формулы обращения. При наличии полных данных Дирихле–Неймана (при измерениях перемещений на всей границе) произведено восстановление значений коэффициентов на границе, а также их нормальных производных. Изучены постановки ОЗ восстановления параметров упругости λ , µ по частичным данным Коши (по измерениям на произвольном открытом подмножестве границы) при рассмотрении ограничений: коэффициент µ постоянен либо близок к постоянному. В работах [22 –26] представлены результаты решения различных ОЗ идентификации переменных материальных и физических характеристик ФГ материалов. При этом считалось, что дополнительная информация задана на том же участке тела, где происходит зондирующее нагружение. Такая постановка ОЗ наиболее традиционна и хорошо изучена.

Актуальным направлением является рассмотрение ОЗ в обобщенной постановке, если области зондирования и съема информации о полевых характеристиках тела, вообще говоря, различны. Пример такой реконструкции в задаче о продольных колебаниях неоднородного стержня, где амплитудно-частотная характеристика задана во внутренней точке стержня, а нагружение реализовано на торце, представлен в работе [27] . С целью развития этого нового подхода к ОЗ в настоящей работе исследуется статическое деформирование двумерной области с переменными механическими характеристиками под действием некоторой механической нагрузки. Прямая задача (ПЗ) состоит в определении поля перемещений области по известному закону неоднородности материальных параметров, граничным условиям и геометрии области. В ОЗ по информации об измеренных смещениях некоторой ненагруженной части границы области необходимо восстановить закон неоднородности упругой характеристики.

• 2.    Общая трехмерная постановка задачи для неоднородного упругого тела

  

  Рис. 1. Изображение трехмерного тела с обозначением частей границы c известными на них условиями

  
  

Рассмотрим задачу деформирования упругого линейного изотропного тела объемом V , ограниченного поверхностью S = S u U S 0 U S f U S _a , под действием статической механической нагрузки P , приложенной к части поверхности S σ . Будем считать, что тело жестко защемлено на части поверхности S u , а часть поверхности S 0 ∪ S f свободна от напряжений. При этом на части границы S f дополнительно заданы перемещения f = u | S^ (Рис. 1) .

Постановка краевой задачи имеет вид [23, 28] :

V-о = 0, а = AE trE+2^е, е = |(Vu+VuT),

ulsu =0> ulsf = f, n • а0 = P, n • °lso =0, где σ — тензор напряжений, ε — тензор линейных деформаций, λ , µ — коэффициенты Ламе, E — единичный тензор, u — вектор перемещений, n — вектор единичной нормали к поверхности тела. В случае неоднородного материала коэффициенты Ламе λ, µ являются функциями координат.

Рассмотрим также слабую постановку задачи, которая понадобится при вычислительной реализации. Для этого введем пробные функции v , удовлетворяющие главным граничным условиям задачи v | S =0 . Уравнение слабой постановки задачи принимает вид [20] :

! (av

V

• u V • v +2 ^ E ^u 0 E v

)dV-!

P • v dS = 0 .

S σ

Здесь E v = ( V v + V v T )/ 2 ; ⊙ — операция полного тензорного умножения [29] . Слабая постановка дает возможность искать «слабое решение» как приближение в конечномерном пространстве. При известных материальных параметрах на основе (2) с помощью МКЭ нетрудно получить приближенное решение задачи.

Будем считать, что в ПЗ заданы: геометрические и материальные параметры тела; условие закрепления u | s =0 ; внешняя нагрузка P , приложенная к части поверхности S _a . При этом дополнительное граничное условие на части S f границы в (1) игнорируем. Требуется определить компоненты поля перемещений u в каждой точке тела x ∈ V . Класс таких задач в механике деформируемого тела является традиционным и чаще всего подразумевает использование стандартных техник решения, включающих полуаналитические и численные методы, такие как методы пристрелки, МКЭ и другие.

В ОЗ необходимо идентифицировать законы неоднородности материальных параметров A ( x ) , ^ ( x ) по дополнительной информации об измеренных смещениях f = u l _S/ на части границы S f (Рис. 1) . Нагрузку P на части S _σ полагаем зондирующей. В общем случае части границы S _σ и S f не совпадают; в настоящей работе рассматривается именно такое обобщение. Постановки частных ОЗ о колебаниях неоднородных тел, когда съем перемещений происходит непосредственно в месте приложения нагрузки ( S _σ ≡ S f ), исследованы во многих работах, в том числе в [22 –26] . Отметим, что при задании дополнительной информации на некоторой части границы ОЗ становятся нелинейными некорректными и требуют специальных подходов к их решению (например, итерационной регуляризации, различных проекционных методов, направленных на алгебраизацию этих задач) [22] .

• 3.    Деформирование неоднородной плоской области

На основе общей модели (1) рассмотрим задачу о статическом деформировании двумерной области Ω с переменными материальными характеристиками. Граница области по аналогии с предыдущим разделом разбита на части: дС = l u U 1 0 U l f U l _a . Считаем, что область испытывает плоское напряженное состояние. Сформулируем постановку краевой задачи:

{ & ij,j —⁰ ,   & ij — ^{A d} ij u k,k ⁺ d ^{( u} i,j ^{+ u} j,i ) ,

^ ^^{A u} i,i v j,j ⁺2 № ⁽ u i,j ⁺ u j,i ⁾⁽ v i,j ⁺ v j,i )^ ^dV   j P i v i ^dS — ⁰‘

Ω                                                        l σ

Подобные задачи для неоднородных тел чаще всего решаются численно, например, с помощью МКЭ. Однако, с точки зрения анализа решений и применения к ОЗ, больший интерес представляют приближенные аналитические решения и создание на их основе упрощенных моделей. Опишем далее подход к построению приближенного аналитического решения ПЗ (3) .

Для областей канонической формы (например, прямоугольных) и определенных видов силового нагружения часто возможен ввод функции напряжений Эри Φ в форме, отвечающей граничным условиям задачи. В общем случае удобно использовать ее полиномиальное представление, как это сделано в классических трудах Менаже, Тимошенко и других [28, 30] :

NN ф = XXAmnXmxn,                            (5)

m=2n=2

где коэффициенты A mn подбираются исходя из граничных условий задачи. Соответствующие компоненты поля напряжений определяются согласно формулам: ^ 11 = Ф, 22 , ^ 22 = Ф, 11 , а ₁₂ = — Ф, 12 . Далее представим коэффициенты Ламе через упругий модуль E и коэффициент Пуассона ν , воспользовавшись формулами:

E         Eν

^М =2(1+ V ) ^{, А} =(1+ v )(1 — 2 v ) .

На основании обобщенного закона Гука выразим компоненты деформаций через напряжения:

( i,j,k =1 , 2) ,

^£ ij = E 1[(1+ V ) ^a ij -V^ ij V kk ]

£ 11 = ^( ^ 11 —V^ 22 ) ,  £ 22 = —( ^ 22 - VOn ) ,   £ 12 = ——^ 12 .

EEE

Производные перемещений, а затем и сами компоненты перемещений найдем через компоненты тензора деформаций [28] :

x1

u i,j = ^{u 0} ,j +    ^{( £} ij,1 + ^£ i1,j - ^£ j1,i ^{) l} x ₂ = 0 ^dx 1 +    ^{( £} ij,2 + ^£ i2,j - ^£ j2,i ^)dx 2    ^{( i} = ^j) ,

x1

U i = u ⁰ + У U i,1 | _x 2 = ₀ dx 1 + У U i,2 dx 2   ( i,j = 1 , 2) .

При этом константы u _i ⁰ , u _i ⁰ _,j определим исходя из граничных условий задачи.

4. ПЗ для прямоугольной области

Рис. 2. Схема к задаче для 2D неоднородной области с консольным закреплением

Рассмотрим в декартовой системе координат O x 1 x 2 задачу об изгибе вытянутой прямоугольной области с одной защемленной гранью (Рис. 2) . К противоположной грани приложена касательная изгибающая нагрузка P . Будем считать, что область неоднородна в направлении продольной оси x ₁ : упругий модуль является переменным и зависящим от координаты E = E ( x 1 ) >  0 , V x 1 G [0 ,l ] , в то время как коэффициент Пуассона сохраняет постоянство. Для удобства дальнейших построений введем функцию податливости: s ( x 1 ) = E ^{- 1} ( x 1 ) .

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

• 4.1.    Приближенная процедура решения 2D прямой задачи

Уравнения равновесия выглядят так же, как и в системе (3) . Запишем граничные условия задачи:

' u1 | x i =l = ^0,

< ^11 \ x i =0 = ^0,

u2|xi = l =0, h/2

У Mxi = 0 ^dx 2 = P,

- h/ 2

. ^a22 \ x i = ± h/2 = ⁰ ,    ^a 12 \ x i = ± h/2 = 0

Применим методику, описанную в разделе 3 и определим поле напряжений с использованием функции Эри в полиномиальном представлении (5) . Для удовлетворения граничным условиям (9) ограничимся в (5) двумя слагаемыми: Ф = b ₂ x 1 x ₂ + d 4 x 1 x 2 . В результате приходим к выражениям для напряжений: 7 11 = Ф ,₂₂ = 6 d 4 x 1 x ₂ , а ₂₂ = Ф ,11 = 0 , (7 12 = — Ф , ₁₂ = - (b ₂ + 3 d 4 x 2) . Из граничных условий (9) находим коэффициенты b ₂ , d ₄ и получаем следующие формулы для вычисления напряжений:

P

7 11 = JX 1 X 2 ,

7 22 = 0 ,

P h ²

^{7 12} =J i

в которые введен момент инерции J . Далее определяем компоненты тензора деформаций:

. P .           . Pv ( у _e    (1+ v )P

£ 11 = JS ( X 1 ) X 1 X 2 , £ 22 = JS ( X 1 ) X 1 X 2 , £ 12 = J

Подставляя (11) в формулы (8) и дважды производя интегрирование, выражаем компоненты перемещений:

x i

PP u1 = u^+U0 2x2 +x2 s(z)zdz +(vs,1x1 -(v+2)s)x2,

’ J J           ⁶ J

( x i                 \       xi z

P             Pν 2

s ( z ) dz — s (0) x 1 I ——      s ( y ) ydydz — ~^ysx 1 x ₂ .

Константы u 1 , u 2 , u 12 , u 2 1 в (12) могут быть найдены из условий закрепления, причем не единственным образом (аналогичная ситуация для однородной области обсуждается в [28] ). Например, при условиях 2 ( U 1, 2 + u 2,1 ) = £ 12 1 x 1 =0 ,x ₂ =0 , U 1 ( l, 0) = 0 , U 2 ( l, 0) = 0 , U 1,2 ( l, 0) = 0 получаем:

l z                                i

I

u ⁰ ₁ = 0 ,

u 0 ₂ =

— УУ s ( n ) ydydz — l У s ( z ) zdz — -h ² (1+ v )У s ( z ) dz ,

I

u

■ 1,2 = ^— P У ^{s ( z ) zdz},

u

Г i

.₂₁ = — js(z)zdz + -s(0) h ² (1 + v ) .

С учетом этих выражений поле перемещений представим в виде:

P

U 1 = jX 2

x i

/ s ( z ) zdz + — ( vx 1 s ' ( x 1 ) — ( v +2) s ( x 1 )) x 2 , 6

l

l

P

U2 = J I   S(z) I X1z xi

z ² — -h ² (1+ v ) ]dz — -vx 1 x 2 s ( x 1 ) I .

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

• 4.2.    Балочная 1D модель Тимошенко

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

U 1 = 0 ( X 1 ) X 2 , U 2 = w ( X 1 ) , U 3 = 0 , (15)

где θ — угол поворота нормали к оси балки, w — прогиб. Постановка краевой задачи в этом случае принимает

вид [30] :

( ( Ejey — nF ( е + w ‘ )=0 , 1 ( nF ( е + w' )' = 0 , w ( l ) = 0 , i ( l ) = 0 ,                                                            (16) I Eje ‘ \ x i =o =0 , nF ( e + w ‘ ) \ x i =0 = P.

Здесь µ — модуль сдвига; F — площадь поперечного сечения балки; штрих обозначает дифференцирование по координате.

Рассмотрим пробные функции W (l) = 0 ,G (l) = 0 и запишем уравнение слабой постановки для последующего

КЭ решения задачи:

l J [ Eje ’ O ’ + nF ( e + w ‘ )( O + W ‘ )] dx + PW (0) = 0 .                        (17) 0

Осуществив замену n = E/ (2(1 + v )) , преобразуем уравнение (17) к виду:

l

IE ( x ) je'O' +   F^ ( е + w')( O + W ‘ ) dx + PW (0) = 0 .

Для сравнения приведем выражения для компонент перемещений (14) из приближенной 2D модели, сгруппировав их по степеням x ₁ , x ₂ :

u i

P

u ¹ , ^{2 +} J

x 1

I sl.z t zdAx ² ₊

|

e(x i )

}

P

—^(vs, 1 x 1 ^{- ( v} + ^2)s)x 3 , O J

|---------------V--------------- '

Q(x i )

U 2 = u ⁰ + u 2, ₁ X 1 +

Ph ² (1 + v ) 4 J

x 1

s ( z ) dz - s (0) x 1

W (X 1 )

x 1 z

- J// ^{s (} n ) ndnd'

' z -PJ^sX 1 x 2 ,

}    R(x i )

где константы u 2 , u 1 ₂ , u 2 ₁ определяются согласно (13) . В отличие от структуры поля перемещений в рамках модели Тимошенко (15) , здесь в выражениях для u ₁ и u ₂ , присутствуют дополнительные члены при степенях более высокого порядка: Q ( x ₁ ) x 2 и R ( x ₁ ) x 2 .

• 4.3.    Балочная 1D модель Бернулли–Эйлера

Рассматривая в условиях (15) в качестве угла поворота функцию вида: е = — w, ₁ , с помощью вариационного принципа получаем постановку задачи об изгибе балки согласно упрощенной модели Бернулли–Эйлера:

( EJw ” )" = 0 ,

w ( l ) = 0 , w ‘ ( l ) = 0 ,

EJw"\ _x1 = ₀ =0 ,   ( EJw ")'^_^ = P.

Для формулировки слабой постановки введем дополнительную функцию y ( x ₁ ) = EJw " , а также пробные функции, удовлетворяющие главным граничным условиям: W ( l ) =0 , Y (0) = 0 . Уравнение слабой постановки принимает вид:

l

! w ’ Y' + y'W' + -JY

dx — PW (0) = 0 .                           (20)

• 4.4.    2D КЭ модель

Для численной оценки предложенной в разделе 4.1 теоретической модели также использован МКЭ в 2D постановке. Для КЭ расчетов применен пакет FreeFEM++ на основе слабой постановки. Представим слабую постановку задачи (2), выразив коэффициенты Ламе через E и ν с помощью формул (6). В условиях плоского напряженного состояния получаем уравнение:

[ f V               1

j E \ -V 2 U ^Uiij ■ 1(1- V )

Ω

^{( u} i,j + ^u ji ^{)( v} i,j + ^v j,i )

dV -

l

σ

P i v i dS = 0 ,

которое справедливо для любых пробных функций v i | l^ = 0   ( i,j = 1 , 2 ).

• 5.    Вычислительные эксперименты по решению ПЗ. Сравнение моделей

  Проведены вычислительные эксперименты по решению ПЗ для прямоугольных областей с различной неоднородностью, характеризуемой податливостью, описываемой разными законами. Полученные на основе предложенной упрощенной 2D модели результаты сравниваются с результатами решений аналогичных задач, осуществленных по 1D балочным теориям и посредством 2D МКЭ-реализации. Расчеты выполнены при следующих параметрах задачи: l = 1 м, E 0 = 2 • 10 ¹¹ Па, s Q = 1 /E 0 , v = 0 . 3 , P = 10 6 Н. Размер балки h принимал разные значения (см. конкретную задачу).

• 5.1.    Описание экспериментов и рассмотренных расчетных моделей

В разделе 5 на всех рисунках расчетные данные приведены для продольной переменной x = l — x 1 Е [0 ,l ] , то есть у балки защемлена левая грань, а поперечная нагрузка приложена к правой грани. Линии на рисунках отвечают расчетам по предложенной 2D модели (модели 1), точки — по модели 4 (2D КЭ модели). В легендах к рисункам для обозначения выбранной модели использован ее номер (1 или 4), а в скобках указана расчетная координата x 2 . Ниже перемещений графики демонстрируют погрешность их вычислений по модели 1, отнесенных к результатам по модели 4, при этом δ 1 и δ 2 — соответственно, относительные погрешности горизонтальной и вертикальной составляющих вектора перемещений.

Итак, в вычислительных экспериментах задействованы 4 модели:

• –    модель 1 — 2D аналитическая приближенная теория, полученная с использованием функции напряжений Φ . Определялось поле перемещений в виде (14) . Расчеты проводились в пакете Maple;

• –    модель 2 — 1D балочная модель Тимошенко. В пакете FreeFEM++ выполнялась одномерная КЭ-реализация на основе слабой постановки (17) ;

• –    модель 3 — 1D балочная модель Бернулли–Эйлера. Как и для модели 2, использовалась КЭ-реализация соответствующей слабой постановки (18) ;

• –    модель 4 — стандартная численная 2D КЭ модель, которая базировалась на постановкее (3) . Численное решение получалось из слабой постановки (21) .

• 5.2.    Сравнительный анализ решений для однородной области

Для 1D КЭ-реализации по моделям 2 и 3 осуществлялось разбиение области на 100 квадратичных элементов в продольном направлении. Для 2D КЭ-реализации строилась равномерная сетка вдоль оси x 1 с числом квадратичных элементов не менее 100. В качестве законов для функции податливости при получении нижеследующих результатов использованы постоянный закон s = s Q , моделирующий однородный материал, и квадратичный закон s ( x i ) = s Q (1 + 2( x i / l )²) , описывающий в материале плавный функционально-градиентный переход.

На рисунке 3а , б представлены компоненты вектора перемещений u 1 ,u ₂ в срезах х ₂ = 0 , x ₂ = ±h/ 2 (см. Рис. 2) , найденные с помощью предложенной 2D приближенной теории (модель 1) и 2D МКЭ-реализации (модель 4). Рассмотрено соотношение сторон h/l = 0 . 3 . Расчетные данные воспроизведены в 20 точках вдоль оси x 1 Е [0 ,1 ] .

Рис. 3. Перемещения u i ( а ) и U 2 ( б ) из решения ПЗ для однородной области в срезах Х 2 =0 , Х 2 = ± h/ 2 , рассчитанные по моделям 1 и 4; ( в ) и ( г ) – их погрешности вычисления в рамках модели 1, отнесенные к решению МКЭ

( в )

Рис. 3. Продолжение

Максимальная относительная погрешность результатов по модели 1 по сравнению с КЭ моделью 4 составляет менее 1.4% для компоненты u 1 и менее 2.3% для компоненты u 2 . И хотя размеры области не вписываются в рамки балочной теории, тем не менее наблюдается достаточно хорошее согласование результатов.

На рисунке 4 для сравнения представлены вертикальные перемещения и ₂ в двух срезах х ₂ = 0 , x ₂ = — h/ 2 , полученные с помощью всех рассматриваемых моделей, для узкой (1:20) и широкой (1:2) прямоугольных областей. Отметим, что для узкой области (Рис. 4а ) максимальная относительная погрешность решения по модели 1, по сравнению с КЭ моделью 4 (Рис. 4в ), составляет менее 0.1%. 1D модели дают близкие погрешности, не превышающие 3%; расхождение результатов наблюдается в окрестности среднего продольного сечения области, при этом на концах области ( х = 0 , x = 1 ) погрешность не достигает и 0.1%. Для широкой области максимальная погрешность решений по моделям 1 и 2, по сравнению с моделью 4, составляет около 6%, в то время как погрешность данных, рассчитанных по модели 3 (в рамках теории Бернулли–Эйлера), равняется почти 16%.

( в )

Рис. 4. Перемещения из решения ПЗ на основе моделей 1–4 для прямоугольной области при разном соотношении сторон: 1:20 ( а ) и 1:2 ( б ); ( в ) и ( г ) – погрешности решений ( а ) и ( б ) относительно решения МКЭ

• 5.3.    Сравнительный анализ решений ПЗ для неоднородной области

Отметим, что выводы, сделанные для однородных областей исходя из сравнительного анализа расчетных данных, установленных по разным моделям, сохраняются для неоднородных областей при различных монотонных и немонотонных законах материальной неоднородности.

Ниже представлены результаты для неоднородной прямоугольной области при использовании квадратичного закона для функции податливости: s ( x 1 ) = s 0 (1 + 2( х 1 / 1 )2) , x 1 6 [0 ,1 ] . В узкой области, как и в однородном случае, погрешность решения по модели 1, по сравнению с погрешностью МКЭ решения, составляет менее 0.1%. Результаты для более широкой области содержат рисунки 5 и 6.

Результаты проведенных вычислительных экспериментов демонстрируют достаточно хорошее согласование предложенной 2D приближенной теории с КЭ моделью при решении 2D задачи для прямоугольных областей с различными соотношениями сторон, особенно, если область узкая, то есть ее форма соответствует балочным

Рис. 5. Перемещения u i ( а ) и U 2 ( б ) в ПЗ для неоднородной широкой области (соотношение сторон 1:3) в срезах Х 2 = 0 , x ₂ = ± h/ 2 ; ( в ) и ( г ) - их погрешности по сравнению с МКЭ

Рис. 6. Перемещения u 2 ( а ) из решения ПЗ для неоднородной широкой области (соотношение сторон 1:2) в срезах х ₂ = 0 , x ₂ = — h/ 2 ; относительная погрешность решения по модели 1 по сравнению с решением МКЭ ( б )

теориям. Следует также отметить хорошую точность балочной модели Тимошенко, по сравнению с МКЭ, даже для широких областей. При этом модель Бернулли–Эйлера показывает самую низкую точность.

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

• 6.    ОЗ реконструкции закона неоднородности материальной характеристики прямоугольной области

6.1. 2D приближенная теория (модель 1)

Далее при решении ОЗ будем использовать две из приведенных выше моделей — предложенную 2D приближенную (модель 1) и 1D балочную модель Тимошенко (модель 2). Модель Бернулли–Эйлера (модель 3), как наименее точную, исключим из рассмотрения.

Рис. 7. Схема области к ОЗ и граничные условия

ОЗ состоит в следующем: необходимо восстановить закон неоднородности функции податливости s ( x ₁ ) по заданным смещениям на нижней грани области (Рис. 7) . В качестве дополнительной информации введем функции:

f 1 ( x i )= U 1

f 2 ( X 1 ) = U 2 ( x—-hj , X 1 6 [0 ,1 ] . (22)

Рассмотрим два подхода, основанные на упрощенных теоретических моделях 1 и 2.

Сформулируем уравнение для ОЗ исходя из представлений предложенной приближенной 2D теории. Выразим дополнительные функции (22) через заданные компоненты смещений (14) :

f 1 ( x i )= U 1

f 2 ( x i )= U 2

x 1

Ph ¹          h ²

J I s(z Zzd z + — ( vx i s ( x i ) - ( v +2) s ( x i )) ,

l

l

x 1 z

-

z ²

-

¹ h ²

(1+ v ) 1 dz

-

- vh²x^s ( x^ ) . 8

Дифференцируя один раз первое уравнение из (23) и дважды второе, получаем:

48J h2vx1s (x1) — 2h2s (x1) + 24x1s(x1) = — -j-j^f1(x1),

h ² vx ₁ s' ‘ ( x ₁ ) — 2 h ² s' ( x ₁ ) + 8 x ₁ s ( x ₁ ) = — ^pf ^ ‘ ( x ₁ ) .

Как видно, левые части уравнений (24) содержат общий член h ² vx ₁ s' ‘ ( x ₁ ) — 2h ² s' ( x ₁ ) . Вычитая из первого уравнения второе, приходим к явному выражению искомой функции податливости:

₍ . J 1 s (x-| ) =--

^{( 1 )}    2 PX 1

(^f 2 ' ^(x 1 ) ^— h ^f 1 ^(x1^) .

Основная вычислительная трудность использования формулы (25) состоит в нахождении производных численно заданных функций (в решаемой ОЗ дополнительная информация известна в наборе точек (см. Рис. 7) ). Численное дифференцирование подобного рода функций представляет собой, вообще говоря, некорректную проблему, требующую регуляризации. В данной ОЗ можно прибегнуть, например, к сплайн-интерполяции.

6.2. 1D модель Тимошенко (модель 2)

На основе условий балочного деформирования (15) запишем функции перемещений на нижней грани области x ₂ = — h/ 2 : u ₁ = —№ /2, u ₂ = w . С их учетом получаем явные представления для функций f ₁ ( x ₁ ) , f ₂ ( x ₁ ) в решаемой ОЗ: f ₁ ( x ₁ ) = — h^ ( x ₁ ) / 2 , f ₂ ( x ₁ ) = w ( x ₁ ) . Отсюда связь между этими функциями и решением соответствующей ПЗ для балки имеет вид:

в ( х 1 ) = — 2 f 1 ( x 1 ) /h, w ( x 1 )= f 2 ( x 1 ) .

Интегрируя второе уравнение системы (16) ирассматривая его вместе с граничным условием pF ( в + w¹ ) |_xi =о = P , получаем:

pF ( в + w ¹ ) = P, V X 1 6 [0 ,1 ] .                                         (27)

Подстановка выражения (27) в первое уравнение системы (16) приводит к тривиальной краевой задаче, описывемой обыкновенными дифференциальными уравнениями 1-го порядка:

Г ( EJ0¹ ) = P,

[ ^Eje l x i =0 = ^0.

Ее решение позволяет выразить закон неоднородности для упругого модуля через функции f ₁ ( x ₁ ) , f ₂ ( x ₁ ) : E ( x ₁ ) = Px ₁ / ( Jв’ ) . С учетом соотношений (26) имеем:

E ( x 1 ) = —hPx 1 / (2 Jf 1 ) .

Формула (29) дает возможность вычислять искомую характеристику всюду, кроме малой окрестности нагруженного конца балки.

• 7.    Вычислительные эксперименты по решению ОЗ

  • 7.1.    2D приближенная теория (модель 1)

На рисунках 8,9 представлены результаты реконструкции различных законов неоднородности s ( x 1 ) с помощью формулы (25) . Принято, что l = 1 , h = 0 . 05 , остальные параметры задачи такие же, как в разделе 5. При численном дифференцировании функций, входящих в (25) и имеющих вид набора значений, для регуляризации использована сплайн-интерполяция методом, встроенным в пакет Maple (тип сплайна «not-a-knot spline», порядок N ). В

Рис. 8. Примеры реконструкции слабой неоднородности, описываемой разными законами: точный закон s 0 , N =2 ( а ) и

N =3 ( б ); квадратичные законы s 01 ( в ) и s 02 ( г ) при N =5 ; число точек восстановления M = 10

Рис. 9. Примеры реконструкции существенной неоднородности при разных законах представления функции податливости, а также различных порядках сплайна и числах точек восстановления: s 01 ( а ), N = 4 , M = 10 ; s 02 ( б ), N =5 , M =10 ; s 03 ( в ), N =5 , M =20 ; s 04 ( г ), N =7 , M =15

Рис. 9. Продолжение

восстановлении участвует M = 10 узлов. Сплошная линия соответствует точному закону s ₀ ( x _i ) , точки — восстановленному s ( x _i ) (в обозначениях на рисунке нижний индекс у x опущен).

Рисунок 8 содержит примеры реконструкции сплайнами разного порядка N слабой неоднородности. Функция податливости описывается точным s ₀ ( x _i ) = s Q (Рис. 8а , б ) и квадратичными s _0i ( x _i ) = s Q (1+0 . 3( x _i / l )²^ (Рис. 8в ) и s ₀₂ ( x 1 ) = s Q ^1 — 0 . 3( x 1 / l ) ²^ (Рис. 8г ) законами.

На рисунке 9 показаны результаты восстановления функции податливости, характеризующей существенную неоднородность, при различных порядках сплайн-интерполяции и разном числе точек восстановления. Использованы представляющие ее законы: S oi ( x i ) = s 0 (1+2( x i / l )) (Рис. 9а ); S 02 ( x i ) = s Q (1+2( x i / l )²^ (Рис. 9б ); s ₀₃ ( x _i ) = s Q ( 1+2 e ^{8 x- 8} l ) (Рис. 9в ); s ₀₄ ( x _i ) = s 0 (1 + 1 . 4sin( nx _i / l )) (Рис. 9 г ).

Из анализа результатов вычислительных экспериментов следует, что для получения приемлемой реконструкции существенно неоднородных немонотонных законов необходимо задавать достаточно высокий порядок сплайн-интерполяции. Случай на рисунке 8а приведен исключительно с целью демонстрации качества реконструкции сплайном низкого порядка. Следует отметить, что, не считая концов ( x _i = 0 и x _i = l ) при существенной неоднородности s ₀ ( x _i ) , погрешность реконструкции в большинстве рассмотренных примеров составляет менее 0.1%.

• 7.2.    1D модель Тимошенко (модель 2)

Рисунок 10 демонстрирует два реконструированных закона неоднородности — линейный s ₀ ( x _i ) = s Q (1+2( x _i / l )) (Рис. 10а ) и квадратичный s ₀ ( x _i ) = s Q (1+2( x _i / l ) ²) (Рис. 10б ); приведены графики соответствующих функций упругого модуля E ( x _i ) = s ^{- 1} ( x _i ) . Прочие параметры задачи описаны выше.

( a )

( б )

Рис. 10. Примеры реконструкции упругого модуля E ( x 1 ) при линейном ( а ) и квадратичном ( б ) законах неоднородности;

( в ) и ( г ) – соответствующие погрешности реконструкции E относительно точного решения E 0

Рис. 10. Продолжение

X, м

X, м

ПЗ решена МКЭ на основе слабой постановки (18) . Рассмотрено 100 элементов вдоль оси x ₁ . При численном дифференцировании функции смещения f 1 ( x 1 ) (ее значения известны в 20 точках) в формуле (19) использована аппроксимация производных на КЭ сетке, осуществленная с помощью алгоритма, реализованного в пакете FreeFEM++. В построенных графиках первый узел, соответствующий свободному концу прямоугольной области ( x 1 =0 ), исключен. Сплошная линия изображает точный закон, точки — восстановленный. Также показана погрешность реконструкции относительно точного решения.

• 8.    Заключение

Предложен подход к построению аналитических представлений двумерного поля перемещений в ПЗ, базирующийся на введении функции напряжений. Исследована задача статического нагружения прямоугольной области с материальной характеристикой, являющейся функцией координат. Проведен сравнительный анализ решений, полученных с помощью предложенной приближенной теории, одномерных балочных моделей Бернулли– Эйлера и Тимошенко, 2D КЭ модели. Сравнение продемонстрировало хорошую точность аналитических представлений для различных соотношений сторон прямоугольных областей и законов неоднородности упругой характеристики. Опробован подход к решению ОЗ реконструкции одномерного закона неоднородности упругого модуля на основе разработанной приближенной теории и балочной модели Тимошенко. Получены явные формулы для определения искомых характеристик через функции смещений, заданные на одной из граней прямоугольной области в виде набора значений.

Сформулированный подход может быть применен для нахождения приближенных модельных решений прямых задач и анализа решений обратных задач при 2D-законах материальной неоднородности.

Исследование выполнено при поддержке гранта Российского научного фонда (проект № 22-11-00265-П, в Южном федеральном университете.