Об определении линейных и нелинейных свойств неоднородных материалов Текст научной статьи по специальности «Математика»

УДК 539.3

А.С. Кравчук

Московская государственная академия приборостроения и информатики

ОБ ОПРЕДЕЛЕНИИ ЛИНЕЙНЫХ И НЕЛИН1ЕЙНЫХ СВОЙСТВ НЕОДНОРОДНЫХ МАТЕРИАЛОВ

Abstract

The following identification problems arc considered: electrical conduct] vny

identification, linear and non-linear elasticity parameters, viscosity, thermal parameters. The approximate solution methods are developed, based on the Newton method and variational approach.

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

Источником неоднородности конструкционных материалов являются, прежде всего, технологические процессы их получения и обработки [2]. Кроме того, неоднородность, которую необходимо учитывать при построении солее точных теорий поведения классических конструкционных материалов типа сталей, обусловлена их кристаллическим (поликристаллическим) строением. В современных композиционных материалах неоднородность определяется изначально их гетерогенностью. Осооую роль неоднородность структуры играет при построении теорий поведения биологических тканей, поскольку те или иные отклонения от «стандартной» неоднородности являются основой для построения различных методов диагностики. Знание параметров неоднородности должно быть заложено как исходная информация в мезомеханику [3].

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

Модельная задача

Рассмотрим в качестве модельной краевую задачу для обобщенного уравнения

Лапласа

V ■ (о(х)Уф(х)) = 0 ,

(1)

которое описывает распределение электрического потенциала ф = ф(д") в области . материал в которой имеет электрическую проводимость <7, С помощью уравнения (1) можно построить математическую модель одного из современных методов диагностики - метода импедансной или потенциальной томографии [4]. Электрические схемы измерений в данном методе описаны в [4, 7].

Схема этих измерений такова. К поверхности диагностируемого тела (технической конструкции или биологического объекта) прикладывается некоторое распределение электрического потенциала, и измеряется соответствующее распределение электрического тока через поверхность. Естественно, что при технической реализации может быть использовано только конечное количество электродов, а при построении математической модели необходимо рассматривать сложную физическую задачу о контактном взаимодействии двух электропроводных тел. Для упрощения построения теории предполагают, что известно распределение электрического потенциала на всей поверхности X = д£1 диагностируемого объекта:

а также известно соответствующее распределение электрического тока 7у£на той же поверхности:

Соответствие между' граничными токами и потенциалами устанавливается посредством либо решения задачи Дирихле (1)-(2), либо решения задачи Неймана для уравнения (1) с граничным условием

где V - вектор единичной внешней нормали к поверхности X (точка означает скалярное произведение векторов), V - опера гор Гамильтона.

Поскольку основной задачей импедансной томографии является нахождение пространственного распределения проводимости с(х), а известной из эксперимента информацией являются граничные значения тока и потенциала, то математическая модель данного вида диагностики должна включать уравнение (1) и граничные условия (2) и (4), в которых неизвестными будут две функции ст(л'), ф(У). Следовательно, задача является, нелинейной. Более того, единственного измерения- оказывается недостаточным для однозначного ее решения. Данное утверждение сначала было установлено для конечномерного аналога поставленной задачи на основании известной из теоретической электротехники теоремы Телегена (по этому поводу см. работу [8]). Применительно к рассматриваемой здесь континуальной постановке задачи высказанное утверждение означает, что необходимо (теоретически) иметь бесконечное число граничных пар фг, /уТ для того, чтобы задача могла быть решена единственным образом. К счастью (и это впервые отмечено в работе [6]), пространство граничных потенциалов имеет счетный базис, что позволяет поставить задачу о доказательстве

ф(А’) 1е ~ Фх >

ф(х) -» ЛЕ(х), хеЪ .

О)

(4)

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

Приведем три метода решения поставленной задачи.

]. Лагранжев метод последовательных приближений

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

Пометим верхним индексом "V" истинное распределение проводимости и соответствующие ему распределения потенциала внутри и на границе и распределение граничных токов:

Приведем описание основных шагов итерационного процесса в варианте метода Бубнова - Галеркина. Процесс начинается с задания некоторого нулевого приближение: проводимости:

Кроме того, предполагается, что для всех функций базиса ф,-' из эксперимента

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

функцией % разность истинного распределения и нулевого приближения граничных токов (невязку):

Заменяя знак приближенного равенства в (1 П точным, получим два вариационных уравнения, второе из которых связывает две нейзвестные функции (10) в об.тас ш 12:

Пусть имеется набор соответствующих друг другу граничных пар фт'. /.V}

а = а^°-(х), іей

і (>)

найдены соответствующие граничные токи :

а) для каждой базисной функции решаем задачу Дирихле с граничным

условием:

І

б) выбираем произвольную функцию Xе #4^) и вычисляем взвешенную с

Введем разности:

(ИІ)

и линеаризуем правые части в выражениях (9) по разностям (10):

в) считая формально функцию 8у известной, выразим все функции 5о' ' через 8у. Из второго из уравнений (11) видно, что эта операция эквивалентна решению задач Дирихле для уравнения

V • [<т(0)У(5ф(0)] = - V ■ [5оУ(фП'0))] (12)

с нулевым граничным условием для 8ф('). Записывая решения этих задач в виде Бф*'1 =1{1>(&о), где 1со - некоторые линейные операторы, и подставляя полученные решения в первое из вариационных уравнений (11), придем к системе

. Эф(,> 0> Э(11,)(5а))

Г[бо——- + 0 ----------]%с1Ъ = . (13)

^ ду ду

Для решения этой системы предварительно производится ее дискретизация, тип которой определяется выбором базиса в пространстве функций %. Обозначим вектор неизвестных, полученный в результате дискретизации, через 8о, вектор правой части -через §1. Тогда система вариационных уравнений (13) приведется к следующему-набору систем линейных алгебраических уравнений:

которые в целом представляют собой переопределенную систему для неизвестных 5с;

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

о11) = а(0) + Р8о<0\ (14)

где р - экспериментально подбираемый числовой параметр - длина шага. Процесс повторяется с новым значением проводимости 0.

2. Двойственный метод последовательных приближений

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

3. Вариационный метод

Данный метод основан на утверждении о том, что задача Дирихле (1)-(2) эквивалентна задаче минимизации функционала

= || ^оУф|2<Я2 (15)

при ограничениях:

V і = 0, (16)

ФІї=Фї, (17)

3v

(18)

где j - вектор гоков внутри оолаети

j = cV(p

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

Рассмотрим алгоритм решения данной задачи, получивший название метода переменных направлений, в англоязычной литературе - "alternating direction implicit method', сокращенного/:

а) задается распределение проводимости в нулевом приближении cr=oJ и строится решение 2N задач минимизации функционала (15) - N задач по аргументу / с ограничениями (16), (18)-(19) при г=1,2,...,Лт и N задач минимизации по аргументу ф с ограничениями (17) при г—1

Нетрудно видеть, что первая серия задач - это задачи Неймана для уравнения (1), вторая серия - это задачи Дирихле для того же уравнения;

б) просуммируем функционалы вида (15) по всем значениям индекса /. Считая в получившемся функционале проводимость неизвестной, произведем минимизацию по данной переменной. После этого произведем коррекцию проводимости по формуле (14) и перейдем к этапу а) с уточненным значением проводимости о.

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

• проводимость материала заменяется тензором упругих модулей ’а - ча(.т). индекс «4» определяет ранг данного тензора;

• место потенциала ф занимает поле перемещений и(х):

• закон Ома (19) заменяется законом Гука, при этом плотность тока заменяется тензором напряжений ст, а градиент потенциала - тензором деформаций є:

• граничные условия типа Дирихле и Неймана будут иметь вид

где и£ - заданные на границе перемещения, Рх - заданные на границе усилия.

После этого алгоритмы типа I) и 2) предыдущего п.п. будут содержать в точности те же этапы, что и алгоритмы дня решения модельной задачи. Для построения вариационного метода решения следует использовать следующий аналог функционала (15):

Описание подробностей работы с функционалом (20) приведено в работе [9] и поэтому на них не останавливаемся и переходим к более интересным результатам, связанным с идентификацией свойств материалов с нелинейным поведением.

Идентификация параметров неоднородных материалов в линейно упругой области

» >~Uv. a(i/;-V v l\.

СО і

Идентификация параметров неоднородных материалов в нелинейно-упругой области

Начнем с вариационной формулировки основной рассматриваемой здесь задачи

- об определении физико-механических параметров материала неразрушающими методами. Прежде всего, заметим, что аналогом функционалов (15) и (20) будет следующий функционал:

JNL (а, е) = ^ || а - Е(е, х) ? сК1, (21)

2 а

где 5 - тензорная функция тензорного аргумента, определяющая нелинейное

поведение неоднородного материала.

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

Упростим задачу идентификации, предположив, что структура искомой функции 5 задается конечным числом параметров А{, Л2Лу, являющихся функциями только пространственных задач, Тогда функционал (21) при фиксированных распределениях напряжений и деформаций будет зависеть только от этих параметров и мы приходам к классической задаче нелинейного программирования: найти минимум функционала:

]Ж(А, Л2,...,4)Л£ ||с(0 -Е(ес0; А, А2,...,Ах)\2 с!П (22)

■ и

по переменным А/ при соблюдении указанных выше требований на поля напряжений и деформаций.

Рассмотрим конкретизацию возникшей задачи ла примере деформационной теории пластичности без разгрузок. В соответствии со сложившимися традициями решения прямых задач деформационной теории пластичности - методом упругих решений А.А. Ильюшина и методом переменных параметров упругости И.А. Биргера -можно предложить два метода решения задач идентификации. Для этого запишем определяющее уравнение в следующих двух формах:

Ъц = +2Щ ~ 2\ий(е(е))е° =

, 2 - (2-^ (К --цю((е))еХе,и )Ь0 + 2ц(1 - со(е(е)))е(/ а а.(£и )5? + 2Д£г/,

где ю(е) - функция А.А. Ильюшина интенсивности е(н), X, р. - параметры Ламе материала в упругой области.

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

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

Задавая теперь определенную параметризацию функции А. А. Ильюшина ох с) при помощи параметров А/, получаем прежнюю формулировку задачи идентификации с параметрами X, р., А,■.

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

Идентификация параметров вязкости и теплопроводности

Основные соотношения

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

1) для малых амплитуд ультразвуковых воздействии материал можно считать физически линейным;

2) при воздействии ультразвука температура материала повышается, и это повышение полностью описывается вязкой диссипацией;

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

Принятые гипотезы приводят (в рамках аксиоматики сплошной среды) к следующей разрешающей системе уравнений - уравнениям движения:

й , г, ^ . сісд,(Т) сЮ(Т). , г)"г;

-г-Н —Н--------------Ф„(Г-Т)———\с(х) -р—- (24)

ох/ і ох Эт Эг

и уравнению теплопроводности:

Э А Л ■ д Лг , 30(Х) , . <>.. . Т), ,

т—-(-—-0., ) =—■{ \т(‘ - Т)—— + <ру. (I - х) -с/т} <■ рг . (25)

ОХ/ Т() ді - Эт с)х

где С 1;у {т) - тензор функций релаксации, 0,; - тензор функций теплового расширения, / - временная переменная, к- - тензор коэффициентов теплопроводности, 0 = Г~Т0 -перепад температур, т(/) - функция теплоемкости, г=г(х,/) - плотность заданных источников тепла; остальные обозначения прежние. Выражение в фигурных скобках в левой части уравнения (24) представляет собой компоненты тензора напряжений о;;.

Уравнения (24)-(25) представляют собой математическую модель термомеханических процессов в области Г2 с границей X. На границе необходимо задать соответствующие физической постановке граничные условия, при исследовании переходных динамических процессов необходимо задать начальные условия.

Постановка и методы решения задач идентификации

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

а) материал изотропный,

б) нагружение гармоническое,

в) процессы в материале установившиеся.

Принятые предположения позволяют перейти от системы уравнений (24)-(25) к более простой системе:

^ г \ / \ ^ г- \ г / Э0(Т) с , , , Э ~и

— { Г [С?1 (г - т)—Г---+ с2(г-т)— -----5;; - ср (/ - х) 6 • ]с1 т} - р . (2о)

Эх,- дх сЗг ' Эх Зг

^ 1И* ~ ^ + Р>'. (27)

ёхI Т0 О! с)Т г)т

где е£ - девиатор деформаций. Дальнейшие упрощения, основанные на данных эксперимента, заключаются в следующем. Полагаем

С2 (0 = сош1 -ЗК, ф(?) = 0. (28)

Вводится для удобства функция релаксации;

ед = 2ц. (29)

о/

Предположение о том, что ф(/) = 0, отделяет систему уравнений движения 01 уравнения теплопроводности (27). Отметим, что ситуация здесь противоположна гой. которая обычно исследуется в теории термоупругости, когда в несвязанных задачах сначала решается уравнение теплопроводности, а затем - уравнения равновесия или движения, в которых температура дает вклад в виде фиктивных массовых сил; таким-образом, в классической теории разогрев материала является следствием его сжатия, подобно разогреву газа при его сжатии. Из опыта известно, что в твердых телах сложной структуры типа биокомпозитов эффект разогрева за счет сжатия пренебрежимо мал по сравнению с разогревом, вызванным переходом работы вязких напряжений в теплоту (естественно, все это справедливо при большой вязкости и невысоких давлениях).

Гипотеза о том, что вся работа вязких напряжений переход»т в теплоту, приводит к следующему выражению для величины pr [] 1]:

1 V V 3 „ Л Эе? (т) Эе" 01) , ,

рг=2 i] ^ ц-аГ"”"1 -

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

э 1

— {[2ц- f G{ (t - т) ехр(-;сот)*/т]е£ (х) -f К(Еа )8ц} = -p<o~и, . (31)

Эх,

Выражение

Г

J G, (/ - т) ехр(-гсот)с?~ = G.(ito) (.32,)

приводится к одностороннему преобразованию Фурье функции -G', (/); подставляя выражение (32) в систему уравнений (31), получаем следующую систему.

{[2|д(х) - G,*(х, /co)Je,v К(x)(Ett)8й.} + p«)V = 0 . (33)

ОХ/

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

1) каждую из серий задач, возникающих в методах последовательных приближений, необходимо реша й, в комплексных неременных или же реша л, каждый раз две отдельные задачи в вещественных переменных;

2) для восстановления ядра сдвиговой релаксации G;(f ) необходимо построить

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

3) уравнение (33) имеет структуру уравнения Гельмгольца, что при малых вязкостях может привести к проблемах! при построении решений вблизи резонансных частот.

Идентификация теплофизических характеристик

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

0(х,() = 0(х)ехр(-2/н>0 . (3-г >

Примем дополнительное предположение о том, что теплопроводность не зависит от температуры, т.е.

m(t) = pcp=Cp, (35)

где С - удельная теплоемкость, рассчитываемая на единицу объема. Подстановка

выражений (34)-(35) в уравнение теплопроводности приводит к следующему

результату';

А(±е.)=_2/а)С;)0-у^ J jj^G, (2/-т-л)хехр(1(в(2г-х-т1))«яГт£/т1. (36) Примем предположение (основанное на данных экспериментов) о том, что

<?! (2f - т - п) = £ Н\к) (* - ^Н2к) “ Л) • (37)

к

Внося выражение (37) в уравнение (36) и повторяя выкладки, проведенные при выводе уравнения (33), получаем уравнение

e,.) + 2icoC. +со2 тн\к)' , (38)

dr, Т0 2 ■ к

где, как и ранее, символом «*» помечено преобразование Фурье соответствующих функций.

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

Работа выполнена при финансовой поддержке гранта РФФИ №99-01-00021 и гранта,Минобразования РФ №ТОО-8,0-2574.

Библиографический список

1. Поздеев А.А.. Няшин КЗ.И,, Трусов П.В. Большие упруго-пластические деформации. - М.: Наука, 1989. - 240 с.

2. Поздеев А.А., Няшин Ю.И., Трусов П.В. Остаточные напряжения. Теория и приложения. - М.: Наука, 1982. - Г11с.

3. Физическая мезомеханика и компьютерное конструирование материалов / Под ред. акад. В.В. Панина. - Новосибирск; Наука, 1995. Т. 1 - 297с., Т. 2 - 320с.

4. Физика визуализации изображений в медицине / Под ред. С. Уэбба; Пер. с англ. -М.: Мир, 1991. Т. 1 - 407с,, Т. 2 - 406с.

5. Натерер Ф. Математические аспекты компьютерной томографии / Пер, с англ. - М.: Мир, 1990. - 267с.

6. Кравчук А.С. Алгоритмы томографии в теории упругости // Прикладная математика и механика. - 1999. - Т. 63. - Вып. 3. - С. 491-494.

7. Webster J.G. Electrical Impedance Tomographie. - Adam Hilger, Bristol, 1990. - 240 p.

8. Dines K.A., Lytle R.J. Analysis of electrical conductivity imaging // Geophysics. - 1981. -Vol. 46. - №7.-P. 1025-1036.

9. Кравчук А.С. Применение вариационного метода решения задач устойчивости в томографических методах контроля качества // Устойчивость, пластичность, ползучесть при сложном нагружении / Тверской гос. техн. ун-т. - Выи.2. - Твсоь. 200. - С.44-48.

10, Nishimura М. A numerical method of crack determination by the boundary integral equation method'/' Некорректно поставленные задачи в естественных науках: -Тр Между нар. конфер. Москва, 19-25 августа 1991, - М.: ТВГ1, 1992,- С.553-562.

11. Ильюшин А.А., Победря Б.Е. Основы математической теории термовязкоупругости. -М.: Наука, 1970.-280с.

Получено 20.03.2001