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

УДК 681.3.06:624.131

ИНТЕГРАЛЬНЫЙ МЕТОД ПОСТРОЕНИЯ МАТЕМАТИЧЕСКОЙ МОДЕЛИ И АЛГОРИТМА ИССЛЕДОВАНИЯ ВЯЗКОУПРУГИХ ДЕФОРМАЦИЙ ГРУНТОВЫХ ОСНОВАНИЙ

Докт. техн. наук, проф. БЫХОВЦЕВ В. Е., канд. техн. наук СЕСЬКОВ В. Е., канд. техн. наук, доц. КУРОЧКА С. К.

Гомельский государственный университет имени Ф. Скорины, РУП «Институт БелНИИС», Гомельский государственный технический университет имени П. О. Сухого

В механике сплошной среды рассматриваются твердые тела, жидкости и газы. Грунты по своей структуре и физико-механическим свойствам не являются простыми телами, которые можно было бы отнести к одной из указанных групп. В большинстве случаев грунты включают в себя твердую, жидкую и газообразную фазы. Твердая фаза грунта (скелет) представляет собой пористую среду. Поры между твердыми частицами в грунте заполнены жидкостью (вода и водные растворы) и газом (воздух). Если поры в грунте заполнены целиком жидкостью (водонасыщенный грунт) или только газом (сухой грунт), то такие грунты называют двухфазными. Грунт называют трехфазным, если поры в нем содержат одновременно и воду и воздух. Двух- и трехфазные грунты воспринимают внешнюю нагрузку иначе, чем сплошные твердые тела. Для грунтовых оснований характерно изменение во времени напряженно-деформированного состояния при постоянной внешней нагрузке. Распределение нагрузок в толще грунтов и их сопротивление внешним силам обусловлено силами сцепления и трения частиц грунта, а также сжатием поровой воды и ее отжатием из пор. Сами частицы при этом практически не сжимаются, происходит только их сдвиг, вследствие чего уменьшаются поры в определенном объеме [ 1-4].

В настоящее время существует несколько подходов к определению деформаций грунтов. В зависимости от конкретных условий для определения деформаций грунтов применяют [5-7]:

• теорию ползучести;

• теорию ползучести с одновременным учетом фильтрационной консолидации;

• теорию ползучести с одновременным учетом фильтрационной консолидации и сжимаемости поровой воды.

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

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

Механико-математическая модель вязко-упругих деформаций грунтовых оснований.

Принято, что скелет грунта изменяется в соответствии с законом линейной наследственной ползучести Больцмана - Вольтерра [1, 5, 8]. В этом случае уравнение состояния в интегральной форме имеет следующий вид:

s(t ) = — En

j(t ) + J K (t, т)о(т)^ t

K(t -т) =

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

Если деформация B(t) известна, то (1) должно решаться относительно напряжения c(t). B этом случае интегральное уравнение второго порядка Вольтерра имеет следующий вид:

c(t ) = Eo

B(t ) -J R(t, t)b(t)J t

где R(t, т) - ядро релаксации, является резольвентой ядра ползучести.

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

Для внесения ясности в представление о природе ядра ползучести K(t, т) необходимо напряжение в (1) принять постоянным, т. е. с(т) = с0 = const и продифференцировать обе части уравнения по t

-1 к (t)=-! dB.

E

с0 dt

Тогда функция K(t) будет соответствовать скорости ползучести при единичной нагрузке, увеличенной в Е0 раз.

Особые формы ядра ползучести имеют следующий вид:

• для реологической модели Кельвина

K(t -т) = 5exp [-SjCt -т)];

• для логарифмического закона вторичной компрессии

( t -в)Y

при 0 < у < 1;

• для комбинации из приведенных выше видов

= Sехр[-51 (t- т)] при о < у < j,

(t-s)Y Р Г

где S - коэффициент ядра ползучести; S1 - то же затухания ползучести. Эти коэффициенты определяются опытным путем в течение нескольких дней.

Интегральная форма уравнения деформирования в настоящей работе принята определяющей. В дальнейшем для краткости будем пользоваться операторной формой записи уравнений (1) и (2):

s^-1(1 + K)а, а = Е(1 -R)s;

E

где K, R - операторы ползучести и релаксации;

t t

Ка = j K (t, т)а(т)^т; Rs = j R(t, l)s(l)d T.

о о

Построение математической модели напряженно-деформированного состояния вяз-коупругого грунтового основания. Физические уравнения для вязкоупругих тел представим в следующем виде:

(t ) = 2G (t ) (1 -Г )(sx (t )-s^ (t ) ) +

+ 3K(t ) (1 -Го Кр (t ); , (t) = 2G (t ) (1 -Гс )(s„ (t)-s^ (t))+ (3) + 3K(t ) (1 -Г0 )sCр (t );

с

xy

(t ) = G (t ) (1 -Гс )Yxy (t ),

где с х

мации; вср (V) =-—; О(V) - модуль

сдвига; К(V) - то же объемной деформации; Го, Гс- операторы объемной и сдвиговой релаксаций.

В скалярной форме выражение (3) будет иметь вид:

5

а

а

:(0 = 20(0) (Вх(0-5Ср (0) )-2^(0, Ъ) (ех-5Ср(^))^ ^ + 3К(^)еср(^)-3| ГД Ъ;

0 о

t 0

,(() = 20(0) ^(0)-5Ср (t) )-2^(0, Ъ) (еу(^)-Еср^У Ъ +3^0^(0)-3| ГД ¡^(¡^ Ъ

0 0 t

т^ (0 = 0(1)1 у (t)-{ Г с (t, ^ (№ Ъ

где Гс(0, Ъ) и Го(t, Ъ) - ядра сдвиговой и объемной релаксаций.

Предположим, что при вязкоупругом деформировании объемные деформации будут упругими, тогда К(0)(1 -Го) = К(t). Выразим модуль сдвига О (0) и модуль объемной деформации К(0) через модуль Юнга Е(0) и коэффициент Пуассона ц(0), тогда выражение (3) с учетом принятой гипотезы об упругости объемных деформаций можно записать в виде:

аХ0 =

а (о =

Е

Е

+ -

(1 -Гс )(Б((00-8Ср (0) ) еср (0);

+

1 -ц

Е

+

+

1 -ц

еср (0);

т

(1 -Гс )У(у (00

Рассмотрим первое уравнение

вим бСр (0) =

а(0) =

2

Е

2

(1-Гс )(б((0}-б((0})-

Е

+--—(б((0)+б((о) =

2

1 -ц

-(б((0+иб((0)-

Г (6 (0-6 (0).

ОЛ , ^ЛЛ ^ "(V//'

2(1 + ц( 0}}

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

Е

2(1 + ц( 0}} а (0} Е(0}

ау(0} = "-—7

( 1 -ц(0 )2

- -^Т^ГТ^ ^ ((0} + 6 у(0}}; 2(1 + ц( 0}}

т (0} Е(0} 1 -ц(0} ^ (0} Е(0} т(((0} = -—тгг—^— У (((0}--

1 -ц(0 }2 2

2(1 + Ц( 0}}

Гс 1 (( (0.

В скалярной форме выражение (7) будет иметь вид:

а (0) =

Е

(6( (0 ) + Ц(0 )6у (0 ) )-

Е

2

Е

-{ Г^, Ъ) (6( (Ъ)-6у (Ъ)) d Ъ

или

(0) = 1-402 (6( (0) + Ц(0 )6( (0) )-

( 0 0 ^ (8) | Гс (0, Ъ)6 ((^ Гс (0, Ъ)6 у (^ Ъ

Рассмотрим интеграл вида

I = | Гс (0, Ъ.

Временной отрезок [0, 0] разделим на п интервалов, тогда

0 И-1 0!+1

i = |гс (0, = x | гс (0, фода.

0 /=0 0,

X

X

По теореме о среднем получим

Я-1 Ь,+1

I =

Хв(п) | Гс (ь, %) с

,=0 Ь

11

где П^; ^+1 ].

В качестве точки п приближенно возьмем одну из границ отрезка, например ^. Тогда

я-1

I «Хв(ь,) I Гс (V,

причем

I = Нш Х

,=0

КЬ) | Гс (V, %

Ядро релаксации примем в виде [1, 8] Г(Ь, %) = 5О (Ь )ехр(-51 (Ь-%)),

где 5 и 51 - параметры релаксации.

Подставляя (9) в предыдущее выражение и интегрируя, получим:

я-1

^О) 15 ехр(-5^-%))С %

1=0 Ь

5

я-1

= - Х(в(/,.)О(/1.)(ехр(-51 (Ь - ^)) - ехр(-51 (Ь - ,))).

51 1=0

В случае другого ядра релаксации (9) интегрирование в (10) может быть проведено численно.

Используя (10), выражение (8) можно переписать в виде

С (V) = Т^^т (вх (V) + (Ь)) -- Х °(<1) х

( я-1

1 -|(02у ......^ М 1=0

Л

х

Рассмотрим последнее слагаемое (11). Для начального времени t0 = 0 последнее слагаемое будет равно 0.

Для времени t1:

с

-О(tl)(в,(tl)(1 -exp(-5l(tl -to)))-

-в.

(Ь1 )(1 - ехр(-5^ - ^)))).

Для времени Ьш:

( т-1

Х 0(Ь] ) (в, (Ь] ХехрН^т - Ь]+1))-

}=0

- ехр(-5! Ь - V]))) - в у (V] )(ехр(-5! (^ - Ь}+1)) -

- ехрН^ - Ь]))))

Таким образом получим

Е

с, (Ь+1) =

с

-Х[О

1 -|(Ь )2

(в, (Ь +) + |(Ь)ву (Ь +))-

ч ]=0

- ехр(-51(Ь+1 - Ь]))) (в,(Ь)-в у }

Аналогично рассуждая, можно найти:

Су

с _

-1тХ[о (Ь] )(ехр(-51(ь,+1 - Ь]+!)) -51 ]=0

- ехр(-5х (Ь,+! - Ь]))) (-в, (Ь]) + в у (Ь]))];

_5

5

т и ) = Е(Ь) 1 -|(Ь) у (Ь ) — х

хХ[О (Ь] )(ехр(-51(Ьг+1 - Ь]+1)) -

]=0

- ехр(-51 +1 - )))у ху (Ь;)].

Геометрические уравнения (уравнения Ко-ши) в случае плоской деформации будут иметь вид [2, 8-10]:

= ди = _

х дх' у ду '

У ху

ди д£ ду дх:

,=0

я-1

I

X

где |6}Г ={6( 6у у((} - вектор деформаций;

и, 3 - соответственно перемещения вдоль осей X, 7.

Закон Гука [5, 9, 12] для линейно упругого тела будет иметь вид

1 -Ц2

1 ц 0 ц 1 0

0 0 ^

Используя (12), (14) можно переписать в виде

|а}г+1=[Е 1+1{6},+1-^^ОХ

Н+1 ={а((0-+1) с( (^+1) т(((0-+:)}Г ; И+ ={6(+!) 6(^+,) у((^+!)}Г ;

где

[ Е] =Е0+1 1 ] 1+1 1 -ц(0 )2

1 ц(0) 0 ц(0) 1 0

0 0 Ь40)

[ О ] = (ехр(-51(0г+1 - +>}} -

- ехр(-51(0г+1 - 01}}}

1 -1 0 -1 1 0 0 0 1

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

времени деформации {6}0 могут быть найдены

из решения упругой задачи.

Плотность энергии деформации будет выражаться площадью диаграммы деформирования материала [9, 10, 12]

6

и

,(0г +1} = | а ^61.

При конечно-элементной дискретизации области существования исследуемой системы энергия деформации для одного конечного элемента будет равна

и+1 = Л и °(0г+1)^,

где S - площадь конечного элемента.

Тогда согласно принятой гипотезе в соответствии с принципом возможных перемещений для момента времени 0г+1 будем иметь

{¿С {*} = Я {6}Т+! {а}г+1 ^,

г+1

вектор узловых перемещений в

где

(1 + 1)-й момент времени, символ «~» означает вариацию; {Я} - вектор узловых усилий.

Рассмотрим треугольный конечный элемент [12-15]. Каждый узел такого конечного элемента в любой момент времени будет иметь две степени свободы

{§}] ={и(0г ) 3(0, )}.

Для аппроксимации перемещений внутри конечного элемента воспользуемся линейными полиномами:

'"(и ) = а1 (0- ) +а2 (0- )( + а3 (0г )(; . З^ ) = «4 (0г ) +а5 (0г )( + «6 (^ )(,

где а1 (0) - а6 (0) - параметры линеаризации, постоянные для элемента в фиксированный момент времени.

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

Ыг=[ ^н,

где Ы ={и(0) и,(0) ик(0) 3г(0) 3,(0) 3^(0)}

вектор узловых перемещений;

"1 х1 У1 0 0 0

1 х, у 0 0 0

[ Л] = 1 0 хк 0 Ук 0 0 1 0 х 0 у1

0 0 0 1 х, у,

0 0 0 1 хк ук

матрица; (а}, = (а1 а 2 а

координатная

а

.} -

вектор параметров; х, у, х,, у,, хк, ук - координаты узлов.

Из (17) могут быть найдены неизвестные параметры

н =[ лг ы, •

Деформации определяются из уравнений Коши (2), (8), (13):

н=

В х & ) В у & )

У ху

ди

дх дЩ)

дУ

ди

ду

дх

Подставляя в уравнения (20) равенства (17) и производя дифференцирование, получим

н=

а 2

Выражение (21) можно переписать в виде

Н =[ВМ,

где

[ В ] =

0 1 0 0 0 0 0 0 0 0 0 1 0 0 10 10

Подставляя (19) в (22), получим выражения, связывающие три составляющих деформаций с перемещениями:

Н =[В][л]- Ыг .

Подставляя далее (23) в (15), для произвольного момента времени будем иметь

И,+1 =[Е 1.1 [ В][ Л]- (Я у,л }

'+1 5!

0ь)[в],, [в][л]-1 (яузл},).

]=0к '

Подставим

(я „л }Гм 1*1 = 1 / (я „ },+1 [ л-1 ]г [ В ]

0 0

[ Е ]„1 [ В ][ Л]-1 (Я „ },+1 -

Щ) [ В] „ [ В][ л]-1

1 ,=0 у '

dxdy,

где а, Ь - размеры конечного элемента вдоль осей X и У соответственно.

Учитывая, что матрицы [л] и [в] не зависят от координат, их можно вынести за знак интеграла, тогда, проинтегрировав, будем иметь

(*} = *[Л]' [В]Т [Е],+1 [В]х

*[Л](яузл },+ - * [Евязк ],

где * - площадь конечного элемента;

( \

[ Е ] = -У

[ вязк ] Г» / , 51 ,=0

0(/,)(ехр(-51(/,.+1- ,)) --ехр(-51 (/,.+1 - /,))) [Л ]Т [ В]

1 -1 0 -1 1 0 0 0 1

[ "МЯл},

Последнее выражение можно переписать в виде

К Ь =[*](Яузл },+1,

где [Я] = *[Л1]Т[В]Т[Е]г+1[В][Л1] - матрица жесткости конечного элемента; (Лф} ^ - вектор

узловых усилий в момент времени +1, причем

{Яф Ь ={Я} + [ Евязк ].

и моментов я

{яузл },

вектор известных гранич-

Матрица жесткости для всей деформируемой системы будет равна

N

КГл =Х КГ,

г=1

где N - количество конечных элементов, на которые разбита рассматриваемая область; К^ -

элемент матрицы [Я^л], характеризующий вклад ,-го единичного перемещения в г-й ком-

понент узловых сил системы;

КГ -

элемент

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

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

Алгоритм компьютерного моделирования вязкоупругих деформаций можно представить в виде ряда следующих процедур.

1. С помощью графического интерфейса, реализующего методику визуального объектно-ориентированного моделирования [12, 16, 17], осуществляется задание геометрической конфигурации рассматриваемой системы дисперсных твердых тел, задание физических свойств (модуля упругости, коэффициента Пуассона, параметров ядра релаксации) элементов рассматриваемой системы, граничных условий и действующих внешних сил.

2. Согласно изложенному алгоритму формируются матрица жесткости [Я^1] и вектор узловых усилий {Я}.

3. Учитываются граничные условия, для этого осуществляется корректировка матрицы [Я^л] и вектора {Я}. Результатом корректировки матрицы [ЯГл] является понижение ее размерности за счет исключения граничных значений перемещений, при этом не нарушается ленточная и симметричная структуры матрицы. В итоге будут получены: новая матрица жесткости

|К/^, новый вектор узловых усилий {Я0}, вектор неизвестных узловых перемещений

{я"}

ных перемещений и моментов я и вектор

номеров граничных узлов.

4. Осуществляется предобусловливание Хо-

лецкого для матрицы | КГ ^.

5. Задаем номер шага г = 0 и временной интервал 0 е[00; tN ], на котором исследуется

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

6. Осуществляется решение системы линейных алгебраических уравнений (26) с предобу-

словленной матрицей |К^ ^ методом сопряженных градиентов.

7. Для каждого конечного элемента по (25) вычисляется матрица [ Евязк ].

8. Для каждого конечного элемента по (27) вычисляется новый вектор {Я0} . При этом

' ' г+1

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

перемещения {я}.

9. Сохраняется найденный вектор узловых перемещений и моментов {яу°зл } .

10. Проверяется, если

0г+1 е [00; tN ],

то осуществляется переход на шаг 6.

11. Осуществляются визуализация и сохранение результатов моделирования.

Моделирование напряженно-деформированного состояния исследуемой системы. На плитный фундамент шириной 20x40 м, работающий в условиях плоской деформации, и толщиной 0,4 м действует равномерно распределенная нагрузка интенсивностью q =114 кПа. Модуль упругости плиты Е = 27 • 103 МПа, коэффициент Пуассона ц = 0,20. Грунтовое основание представляет собой слой мергеля мощностью 6,5 м. Общая толщина сжимаемой зоны - 9,5 м. Ядро релаксации имеет вид (9), где 5 = 0,05 1/сут; = 0,02 1/сут; Е = 9 МПа; ц = 0,45 [18, 19].

При решении задачи методом конечных элементов рассматриваемая область дискрети-зировалась 7200 конечными элементами в фор-

ме равностороннего А размером стороны 0,3 м. Временной шаг выбирался равным 0,5 сут. Время нахождения решения на отрезке от 0 до 1000 сут. составило менее 30 с. Стабилизационная осадка 5,55 см приходится на точку 395 сут.

В Ы В О Д Ы

1. Согласно результатам проведенного моделирования предлагаемая математическая модель, алгоритм и программное обеспечение могут быть использованы для исследования вяз-коупругого деформирования неоднородных систем механики грунтов.

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

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

Л И Т Е Р А Т У Р А

1. Цытович, Н. А. Механика грунтов / Н. А. Цыто-вич. - М: Стройиздат, 1982. - 264 с.

2. Вялов, С. С. Реологические основы механики грунтов / С. С. Вялов. - М.: Высш. шк., 1978. - 446 с.

3. Шукле, Л. Реологические проблемы механики грунтов / Л. Шукле. - М.: Стройиздат, 1976. - 486 с.

4. Рейнер, М. Реология / М. Рейнер. - М.: Наука, 1965. - 280 с.

5. Александров, А. В. Основы теории упругости и пластичности / А. В. Александров, В. Д. Потапов. - М.: Высш. шк., 1990. - 400 с.

6. Малинин, Н. Н. Прикладная теория пластичности и ползучести / Н. Н. Малинин. - М.: Машиностроение, 1976. - 400 с.

7. Ржаницын, А. Р. Теория ползучести / А. Р. Ржани-цын. - М.: Стройиздат, 1968. - 416 с.

8. Фадеев, А. Б. Метод конечных элементов в геомеханике / А. Б. Фадеев. - М.: Недра, 1987. - 224 с.

9. Старовойтов, Э. И. Основы теории упругости, пластичности и вязкоупругости: учеб. для строит. спец. вузов / Э. И. Старовойтов. - Гомель: БелГУТ, 2001. -344 с.

10. Победря, Б. Е. Численные методы в теории упругости и пластичности / Б. Е. Победря. - М.: Изд-во Московского университета, 1981. - 344 с.

11. Зарецкий, Ю. К. Теория консолидации грунтов / Ю. К. Зарецкий. - М.: Наука, 1967. - 320 с.

12. Быховцев, В. Е. Компьютерное моделирование систем нелинейной механики грунтов / В. Е. Быховцев, А. В. Быховцев, В. В. Бондарева. - Гомель: УО «ГГУ им. Ф. Скорины», 2002. - 215 с.

13. Зенкевич, О. Метод конечных элементов в технике / О. Зенкевич. - М.: Мир, 1975. - 540 с.

14. Журавков, М. А. Математическое моделирование деформационных процессов в твердых деформируемых средах / М. А. Журавков. - Минск: БГУ, 2002. - 456 с.

15. Амусин, Б. З. Метод конечных элементов при решении задач горной геомеханики / Б. З. Амусин, А. Б. Фадеев. - М.: Недра, 1975. - 144 с.

16. Курочка, К. С. Формализация процесса компьютерного моделирования линейной системы: здание-фундамент-основание / К. С. Курочка // Новые математические методы и компьютерные технологии в проектировании, производстве и научных исследованиях: материалы IV респ. науч. конф. студ. и асп. - Гомель: ГГУ, 2001. -С. 29-31.

17. Быховцев, В. Е. Визуальное объектно-ориентированное моделирование зданий с фундаментами на грунтовых основаниях / В. Е. Быховцев, А. В. Быхов-цев, К. С. Курочка // Тр. междунар. науч.-техн. конф., Минск, 10-12 окт. 2001 г.: Т. 2. - Минск: БНТУ, 2001. -С. 5-16.

18. Сеськов, В. Е. Биогенные грунты Беларуси и использование их в качестве оснований под здания и сооружения / В. Е. Сеськов. - Минск: БелНИИНТИ, 1989. -48 с.

19. Сеськов, В. Е. Физико-механические характеристики погребенных биогенных грунтов / В. Е. Сеськов // Современные архитектурно-конструктивные системы зданий и сооружений, новые строительные материалы и технологии: сб. науч. тр. - Минск: Стринко, 2000. -С. 119-128.

Поступила 19.11.2007