Конечно-разностное решение задачи о росте трещины антиплоского сдвига в среде с поврежденностью Текст научной статьи по специальности «Физика»

УДК 539.376

КОНЕЧНО-РАЗНОСТНОЕ РЕШЕНИЕ ЗАДАЧИ О РОСТЕ ТРЕЩИНЫ АНТИПЛОСКОГО СДВИГА В СРЕДЕ С ПОВРЕЖДЕННОСТЬЮ1

© 2004 Л.В.Степанова, Т.Н.Макарова2

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

Введение

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

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

Терминология "в связанной постановке" означает, что параметр поврежденности входит в определяющие соотношения рассматриваемого материала, что позволяет учесть как влияние процесса накопления повреждений

1 Представлена доктором физико-математических наук профессором В.И.Астафьевым.

2Степанова Лариса Валентиновна ([email protected]), Макарова Татьяна Николаевна ([email protected]), кафедра механики сплошных сред Самарского государственного университета, 443011, г. Самара, ул. Акад. Павлова, 1.

на напряженно-деформированное состояние, так и описать влияние эволюции напряженно-деформированного состояния на процесс накопления рассеянных повреждений.

Следует отметить, что вопросам определения напряженно-деформированного состояния в окрестности как стацинарной, так и растущей трещины в связанной постановке задачи в последнее время посвящается значительное количество исследований [1-17], которые показывают следующее:

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

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

Так в [8, 9] установлено, что к берегам растущей полубесконечной трещины нормального отрыва примыкает область полностью поврежденного материала, а перед вершиной трещины существует зона активного накопления повреждений, граница которой определяется в первом случае — в случае усталостного роста трещины — из кинетического уравнения, задающего степенной закон накопления повреждений при действии растягивающих нагрузок, превышающих некоторое пороговое значение, и отсутствие процесса накопления повреждений в противном случае, а во втором — из условия непрерывности интенсивностей напряжений во внутренней области (в области накопления повреждений) и во внешней области (в зоне отсутствия процесса накопления повреждений).

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

Целью настоящего исследования является получение численного решения задачи, сформулированной в [15] для автомодельной переменной, предложенной в [18], методом конечных разностей.

1. Автомодельная переменная в задаче о росте трещины в среде с поврежденностью

Рассмотрим прямолинейную полубесконечную трещину в подвижной системе координат с началом в вершине трещины (рис. 1).

Определяющие соотношения теории установившейся ползучести — уравнения степенного закона Нортона — для среды с поврежденностью, предложенные Л.М. Качановым [19] и Ю.Н. Работновым [20], имеют следующий

У

поверхность трещины

X

Рис. 1. Вершина продвигающейся полубесконечной трещины

вид

/■> / \п— 1

3 „/О/

^ =о В\- (1-26)

2 \у/ у

где Щ/ — компоненты тензора скоростей деформаций, вг/ — компоненты деви-атора напряжений, ое = ^/3^7/^7/72 — интенсивность напряжений, В, п — константы материала, определяемые экспериментально.

Параметр сплошности у = 1 - ш (ш — параметр поврежденности) изменяется в соответствии с кинетическим уравнением [20]

\ т

е

% = ^ , у(, = 0) =1, а«)

где А, т — константы материала, определяемые экспериментально-

Начальные условия для компонент тензора напряжений имеют вид

/ с* \ (п+Г)

Оц(т, ф, I = 0) = I — I Ог/ф, п), (1.28)

где Ог/(ф, п) — функции, определяемые асимптотическим решением Хатчинсона—Райса—Розенгрена (НИИ) [21, 22], т, ф — полярные координаты с полюсом в вершине трещины, С*-инвариантный интеграл теории установившейся ползучести- Для двумерного случая С*-интеграл определяется как контурный интеграл по кривой Г, охватывающей вершину трещины [23]

С* = £ W*dx2 - Ог/пг-^йв , (1.29)

где W* = f Ог/й&г/, 1п —функция, зависящая от п и определяемая как безразмерный С*-интеграл.

Асимптотическое условие в бесконечно удаленной точке (т —> ю), где у = 1, может быть сформулировано как требование асимптотического сближения искомого решения с решением НИИ, поскольку рассматриваемые

определяющие соотношения (1.26) при у = 1 представляют собой традиционные соотношения степенного закона теории установившейся ползучести:

/ C* \

Oij(r га, ф, t) И Bi~r I 5¿Дф' ")• (1.30)

Установлено, что для определяющих соотношений (1.26) с начальными и граничными условиями (1.28), (1.30) существует автомодельная переменная [18]

r

R =---, (1.31)

k(n)(At)(n+1)/m v 7

где k(n) = C*/BIn.

Действительно, перейдем к безразмерным величинам согласно формулам:

r t 5ij

r =—' t = —' 5ij = -77"-—'

L T 4 (k(n)/L)1/(n+r>

где L — некоторая характерная длина, T —характерное время.

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

1 I k(n)\-m/(n+1'

T = A (т2) • (132)

В этом случае безразмерные напряжения 5ij примут форму

1

^Г'Ф'° = (k(n)/L) 1/(»+1)Ч L' tA\—)

r i k(n) X

\m/(n+1)'

(1.33)

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

r/L

R =--ri-,-. (1.34)

(tA)(n+1)'mk(n)/L v 7

В этом случае напряжения и параметр сплошности могут быть представлены в виде

Oij(r, ф, t) = (At)(n+r>/m^ij(R, Ф); (1.35)

V(r, Ф, t) = y(R, Ф). (1.36)

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

Oij(r га, ф, t) = Crsöij(<f, n), (1.37)

где s < 0.

Если принять, что искомое решение при r —> га стремится к решению HRR, то C = C* и s = -l/(n + 1). Если предположить, что в бесконечно удаленной точке решение сближается с решением линейной теории упругости, то C = Km и s = -1/2. В [15] показано, что граничное условие

(1.30) не может быть принято в качестве граничного условия в бесконечно удаленной точке. Асимптотический и численный анализ полей напряжений и сплошности, проведенный в [15], приводит к необходимости постановки граничного условия в более общей форме (1.37), где значения показателя 5 определяются в результате решения задачи на собственные значения (см. таблицу).

Таблица

п = т = 1 5 = -1,5

п = 2, т = 0,7п 5 = -1,2303

п = 3, т = 0,7п 5 = -1,1830

п = 4, т = 0,7п 5 = -1,1648

п = 5, т = 0,7п 5 = -1,1553

п = 6, т = 0,7п 5 = -1,1495

п = 7, т = 0,7п 5 = -1,1455

п = 8, т = 0,7п 5 = -1,1425

п = 9, т = 0,7п 5 = -1,1405

п = 10, т = 0,7п 5 = -1,1390

В самом общем случае автомодельная переменная и компоненты тензора напряжений представлены как

Я = г[гАС]1/5т, (1.38)

вц(г, ф, г) = С£(п)5г-1/то,7(Я, ф), (1.39)

у(г, ф, г) = у(Я, ф). (1.40)

2. Автомодельное решение

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

уравнения равновесия

д — доф7

Ш ^ + -ф = 0. (2-1>

д - дуЯ

дЯ (ЯУФ7) = Цф

(ЯУф7) = , (2.2)

кинетического уравнения

Я= -sm О) , (2.3)

дЯ Ы у 7

где скорости деформаций ползучести имеют вид

(~\п—1~ /— \п—1 _

г , ^ = г , (2.4)

г/ V V

интенсивность напряжений определяется как

2 = л!°Яг + • (2.5)

Решение этой системы уравнений должно удовлетворять граничным условиям:

условию отсутствия поверхностных усилий на верхнем берегу трещины

офг(Я, ф = п) = 0, (2.6)

Ояг (Я, Ф = 0) = 0, (2.7)

асимптотическому условию в бесконечно удаленной точке

О/Я, ф) = Я"о,7(ф, п) (Я —> то, 5 < 0). (2.8)

Решение системы уравнений (2.1)-(2.3), подчиняющееся граничным условиям (2.6)-(2.8), разыскивается во всей плоскости за исключением полностью поврежденной зоны, примыкающей к вершине трещины и внутри которой материал не удовлетворяет сформулированной системе уравнений. Предполагается, что внутри области полностью поврежденного материала все компоненты тензора напряжений и сплошность обращаются в нуль, а на границе введенной области разыскиваемое решение должно удовлетворять условиям непрерывности

V = 0, Оу = 0. (2.9)

В дальнейшем надстрочный символ опускается.

3. Разностные уравнения задачи о трещине антиплоского сдвига

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

АФ

м

ш,м)

(г-1,3-1)

АН

Н

Рис. 2. Рассматриваемая область, представленная равномерной сеткой

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

Перед нами стоит задача отыскания решения системы уравнений: уравнения равновесия (2.1), условия совместности (2.2) и кинетического уравнения (2.3). Решение этой системы уравнений должно удовлетворять граничным условиям: условию отсутствия поверхностных усилий на верхнем берегу трещины (2.6), условию симметрии на ее продолжении (2.7) и асимптотическому условию в бесконечно удаленной точке (2.8).

Решение этой задачи рассматривалось на плоскости с координатными осями Я и ф. Для решения на этой плоскости введена равномерная сетка с шагами ДЯ и Дф соответственно (рис. 2). Для замены дифференциальных операторов была использована неявная левая разностная схема [24, 25]. Неявная разностная схема относится к классу безусловно устойчивых схем. Безусловно устойчивая схема — это такая схема, которая устойчива при любых соотношениях между шагами. Левый разностный оператор выбран в соответствии с заданными граничными условиями. Условию (2.6) в таком случае будут удовлетворять все точки, принадлежащие верхней границе выбранной плоскости. Условию (2.7) удовлетворяют все точки, относящиеся к нижней границе плоскости. Для правой границы плоскости

п

0

= Т^п ^ ф/2 °фг = 7^ с°8 Ф/2' (3.1)

граничные условия известны и имеют следующий вид: для п = 1

Я1/2з1п Ф/2, = КП2

для остальных п — констант материала — распределение этих функций определяется численно (методом Рунге—Кутта—Фельберга) [26] для введенных в таблице.

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

ОЯт.(1, У - 1) - Оят.(1 - 1, У - 1)

Оя =-дя-' (3.2)

у = у(1,] - 1) - у(» - 1,] - 1), (з

у дя у 7

ОЯ1(и У) - ОЯг^, ] - 1) ,„„4

°Яг = ---, (3.4)

дФ

ОфгО, У) - ОфZ(i' ] - 1)

Офг = ---, (3.5)

дФ

. у(1, У) - У(1, У - 1) „ч

V/ =---• (3.6)

дФ

С учетом этих обозначений получим систему алгебраических уравнений следующего вида:

ДЯ / я• \

ОЯг( - 1, У - 1) = я (офг + ОЯz(i' У) --¿Я ОЯг (Ь У - 1)] , (3.7)

ОфгО' - 1, У - 1) = Оф/г; У - 1)+

+_Ое(', У)_

Яг(п - 1)дЯофг(г', У) + Яю1(и У)

X

X

фz^

Офг(г, ])Ое(1, ]) + ОяЖп - 1)-—-+

авг(и ])ОфZ(i' ]Х / _ (п - 1)

+ (п - 1)-ГГТТ-;:-Офг - ОЯг | Ое0, ]) -

(Ш У) (3.8)

Уа У) ^ \ Ое(/, У)

Ое(l' ])офZ(г' ]) „ Ое(Ь ])Оя7.(1, У)

- УЯп-—--ууп-

дЯ

У(1 - 1,У- 1) =

Я:

У(1, У) У(1, У)

я; / Оe(i' ])

;У(;,] - 1) + т&

хч-,, , , (3.9)

ДЯ т ' \ У(1, ])! ] У 7

где О обозначает производную по Я, О обозначает производную по ф, i — номер текущего узла по оси Я; ] — номер текущего узла по оси ф; т = = 0,7п — константа кинетического уравнения [27]; Ояг, Офг — искомые напряжения; Ое = у!аЯг + Офг; у — параметр сплошности.

Начиная с верхнего правого угла рассматриваемой области (где i = N, ] = М), можно, разрешая эти уравнения, осуществить движение по направлениям сверху вниз и справа налево. Вычисления производятся от у = 1 до у = 0, что соответствует достижению границы полностью поврежденного материала.

0.9 0.8 0.7

0.6 0.5 0.4 0.3 0.2 0.1

0 ^

Рис. 3. Геометрия области полностью поврежденного материала для п = 1. Конечно-разностное решение

4. Результаты вычислений. Сравнительный анализ геометрии областей полностью поврежденного материала

Программа, реализующая метод конечных разностей, включающий метод Рунге—Кутта—Фельдберга для определения угловых распределений функций, входящих в граничные условия в бесконечно удаленной точке для различных значений материальных констант, позволяет найти конфигурации областей полностью поврежденного материала. Эти области показаны на рис. 3-6.

Для нахождения границ этих областей, как и в [15], были использованы собственные числа 5, приведенные в таблице.

На рис. 7-10 представлена геометрия области полностью поврежденного материала, полученная в результате аналитического решения задачи методом разложений по собственным функциям [15].

Сравнение графиков, представленных на рис. 3-6 и рис. 7-10 для различных значений т и п, показывает, что конфигурации областей полностью поврежденного материала имеют значительное сходство. Основываясь на этом утверждении, можно сделать вывод о качественном совпадении областей, полученных в результате численного и аналитического исследований. Зададимся вопросом о количественном совпадении графиков. На рис. 6 и 10 представлены области полностью поврежденного материала для матери-

Рис. 4. Геометрия области полностью поврежденного материала для п = 3. Конечно-разностное решение

П=7 ш=0,7П -- , Y

1

-3 -2.5 -2 -1.5 -1 0.5 0 0.5

Рис. 6. Геометрия области полностью поврежденного материала для п = 7. Конечно-разностное решение

2

X

0

3 / / 2 1

Рис. 8. Геометрия области полностью поврежденного материала для п = 3. Асимптотическое решение задачи

-------------

------------- 3 21 V

1 /

Рис. 10. Геометрия области полностью поврежденного материала для п = 7. Асимптотическое решение задачи

альной константы п = 7. При рассмотрении рис. 6, где показана граница области, полученная методом конечных разностей, видно, что протяженность области полностью поврежденного материала в отрицательном направлении оси X составляет 2,5; а на рис. 10, для одночленного разложения параметра сплошности эта же величина имеет значение, равное 2,5. Рассматривая далее графики, полученные для остальных значений материальных констант, видно, что области полностью поврежденного материала имеют не только качественное, но и количественное совпадение.

Выводы

1. Получено численное решение методом конечных разностей задачи о стационарной трещине типа III в связанной постановке (в связке "ползуче-сть-поврежденность") с введением автомодельной переменной. Установлено, что к вершине трещины примыкает область полностью поврежденного материала. Конфигурация данных областей для различных констант позволяет заключить, что найденное решение отвечает стационарной трещине (с введением автомодельной переменной задачи о растущей и стационарной трещине перестают быть различимыми), поскольку область полностью поврежденного материала охватывает вершину трещины и ее берега вблизи вершины (рис. 3-10).

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

Литература

[1] MouY., HanR.P.S. Influence of damage in the vicinity of a macrocrack tip // Engn. Fracture Mechanics. 1996. V. 55. No. 4. P. 617-632.

[2] Ravichandran G., LiuC.T. Crack tip shielding in elastic particulate composites undergoing damage // Engn. Fracture Mechanics. 1998. V. 59. No. 6. P. 713-723.

[3] FengX.Q., YuS.W. Analyses of damage localization at crack tip in a brittle damaged material // Engn. Fracture Mechanics. 1996. V. 53. No. 2. P. 169-177.

[4] JinZ.H., BatraR.C. Crack shielding and material deterioration in damage materials: an antiplane shear fracture problem // Arch. Appl. Mech. 1998. No. 68. P. 247-258.

[5] LeeS.B., LuM., KimJ.Y. An asymptotic analysis of a tensile crack in creeping solids coupled with cumulative damage — Part I. Small damage region around the crack tip // Int. J. Solids Structures. 1997. V. 34. No. 24. P. 3163-3178.

[6] LeeS.B., LuM., KimJ.Y. An asymptotic analysis of a tensile crack in creeping solids coupled with cumulative damage — Part II. Small damage region around the crack tip // Int. J. Solids Structures. 1997. V. 34. No. 10. P. 1183-1197.

[7] WangT., KishimotoK. Higher order fields for damaged nonlinear antiplane shear notch, crack and inclusion problems // Eng. J. Mech. A/Solids. 1999. No. 18. P. 963-986.

[8] Zhao J. and ZhangX., The asymptotic study of fatigue crack growth based on damage mechanics // Engn. Fracture Mechanics. 1995. V. 50. No. 1. P. 131-141.

[9] Zhao J. and ZhangX., On the process zone of a quasi-static growing tensile crack with power-law elastic-plastic damage // Int. J. of Fracture. 1995. V. 108. P. 383-395.

[10] MurakamiS., HiranoT. and LiuY., Asymptotic fields of stress and damage of a mode I creep crack in steady-state growth // Int. J. Solids Structures. 2000. V. 37. P. 6203-6220.

[11] MurakamiS., LiuY. and MizunoM., Computational methods for creep fracture analysis by damage mechanics // Comput. Methods Appl. Mech. Engrg. 2000. 183. P. 15-33.

[12] Астафьев В.И., ГригороваТ.В., Пастухов В.А. Влияние поврежденности материала на напряженно-деформированное состояние в окрестности вершины трещины при ползучести // ФХММ. 1992. Т. 28. №1. С. 5—11.

[13] АстафьевВ.И., ГригороваТ.В. Распределение напряжений и поврежденности у вершины растущей в процессе ползучести трещины // Изв. РАН. МТТ. 1995. №3. С. 160-166.

[14] Степанова Л.В., ФединаМ.Е. Автомодельное решение задачи о трещине антиплоского сдвига в связанной постановке (связка "ползучесть-повре-жденность") // Вестник СамГУ. 2000. №4(18). C. 128-145.

[15] Степанова Л.В., ФединаМ.Е. О геометрии области полностью поврежденного материала у вершины трещины антиплоского сдвига в связанной постановке задачи (связка "ползучесть-поврежденность") // Вестник СамГУ. 2001. №2(20). C. 87-114.

[16] LuM., MaiY.-W., YeL. Crack-tip field for fast fracture of an elas-tic-plastic-viscoplastic material incorporated with quasi-brittle damage. Part 1. Large damage regime // Int. J. Solids Structures. 2001. V. 38. P. 9383-9402.

[17] LuM., MaiY.-W., YeL. Crack-tip field for fast fracture of an elastic-plastic-viscoplastic material incorporated with quasi-brittle damage. Part 2. Small damage regime // Int. J. Solids Structures. 2001. V. 38. P. 9403-9420.

[18] RiedelH. Fracture at High Temperature. Berlin: Springer, 1981. P. 504-519.

[19] КачановЛ.М. О времени разрушения в условиях ползучести // Изв. АН СССР. ОТН. 1958. С. 26-31.

[20] РаботновЮ.Н. О механизме длительного разрушения / Вопросы прочности материалов и конструкций. М.: Изд-во АН СССР, 1959. С. 5-7.

[21] RiceJ. R., Rosengren G.F. Plane strain deformation near a crack tip in a power-law hardening material //J. Mech. Phys. Solids. 1968. V. 16. P. 1-12.

[22] HutchinsonJ.W. Singular behavior at the end of tensile crack in a hardening material //J. Mech. Phys. Solids. 1968. V. 16. P. 13-31.

[23] LandesJ.D., BegleyJ.A. A fracture mechanics approach to creep crack growth // In: Me^anics of Crack Growth. ASTM STP 590. 1976. P. 128-148.

[24] КалиткинН.Н. Численные методы. М.: Наука, 1978. 512 c.

[25] Рябенький В.С. Введение в вычислительную математику. М.: Физмат-лит, 2000. 296 c.

[26] Форсайт Дж., Малькольм М., МоулерК. Машинные методы математических вычислений. М.: Мир, 1980. 280 c.

[27] БойлДж., СпенсДж. Анализ напряжений в конструкциях при ползучести. М.: Мир, 1986. 360 c.

Поступила в редакцию 24/77/2004; в окончательном варианте — 24/77/2004.

FINITE DIFFERENCE METHOD APPROACH TO A CREEP-DAMAGE COUPLED MODE III CRACK

PROBLEM3

© 2004 L.V. Stepanova, T.N. Makarova4

An attempt to describe stress, creep strain rate and continuity (integrity) fields near a Mode III crack tip in the framework of the continuum damage mechanics by the aid of finite difference method is presented. The scalar damage parameter is involved into the power creep constitutive equations. The coupled system of continuum damage mechanics-creep theory equations is considered. The solution is based on introducing a self-similar variable proposed by Riedel for the creep power constitutive law coupled to the power kinetic law of the damage evolution. The stress, strain rate and damage fields near the crack tip are obtained as the functions of the similarity variable assuming the existence of a totally damaged zone. It is presumed that there is the totally damaged zone in the crack-tip region where the damage parameter reaches its critical value and all the stresses are equalled to zero.

All the equations of the considered boundary value problem are then formulated as the finite difference equations resulting in the linear algebraic system of equations. The left difference scheme is employed for the approximation. The shape and the characteristic length of the totally damaged zone are found from the obtained solution.

Paper received 24/77/2004. Paper accepted 24/77/2004.

3Communicated by Dr. Sci. (Phys. & Math.) Prof. V.I. Astafiev.

4 Stepanova Larisa Valentinovna ([email protected]), Makarova Tatjana Nickolaev-na ([email protected]), Dept. of Continuum Mechanics, Samara State University, 443011, Samara, Russia.