CC BY
58
2/2006
я • т
МОДЕЛИРОВАНИЕ РАБОТЫ ШВОВ И
ТРЕЩИН В РАСЧЕТАХ НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ БЕТОННЫХ ПЛОТИН
В.В. Толстиков
расчетах напряженно-деформированного состояния бетонных сооружений и их оснований существенную роль играют прочностные и деформа-тивные свойства различных нарушений сплошности (швов, контактов сооружений с массивом основания и крупных трещин в массиве основания) по которым могут развиваться подвижки. На кафедре гидротехнических сооружений накоплен достаточный опыт решения различных задач в рамках метода конечных элементов с учетом нарушений сплошности. В основном он связан с применением специальных контактных элементов, основное преимущество которых перед другими способами решения таких задач состоит в том, что их использование не нарушает основной алгоритм МКЭ и не ведет к увеличению временных затрат на получение решения по сравнению с затратами на решение упругопластических задач. Несмотря на широкое применение контактных элементов, методика и алгоритмы моделирования нелинейных эффектов контактного взаимодействия в публикациях представлены весьма скудно, в том числе и в публикациях работ выполненных на кафедре. В данной работе подробно рассмотрен механизм контактного взаимодействия, используемый в исследовательских программных комплексах разработанных на кафедре.
Деформационные свойства нарушений сплошности (и в частности трещин) в настоящее время принято рассматривать с двух точек зрения:
1) перемещения в направлении, нормальном к нарушению сплошности (нормальная деформация);
2) перемещения в направлении, параллельном нарушению сплошности (сдвиговая деформация).
Поведение трещины под действием нормальных напряжений нелинейно, как показано на рис. 1а.
Освобождение от первоначального сжатия вызывает взаимное разделение точек контактирующих поверхностей трещины, которая раскрывается и достигает состояния отсутствия нормальных напряжений ап после того, как достигнута прочность на растяжение Яр. При сжатии контактирующие поверхности будут смыкаться за счет упругого (и неупругого) смятия шероховатостей или деформирования заполнителя. По мере увеличения нормальной нагрузки величина смыкания асимптотически стремится к пределу, максимально возможному закрытию трещины Ктах, которое должно быть меньше «толщины» трещины. Зависимость оп=/(у) на участке смыкания трещины в общем случае нелинейна, но для расчетов, где не исследуются детально вопросы деформации трещины, можно принять линейную аппроксимацию, как принято во многих исследованиях. Тогда модель поведения трещины при нормальной деформации можно построить используя три параметра: прочность на растяжение (Кр), максимально возможное закрытие трещины (^тах), удельный модуль деформации (£п) в направлении
Рис. 1. Поведение трещины под действием нормальных и касательных напряжений
нормальном к трещине, который далее будем называть нормальной жесткостью. Необходимо также учитывать и начальное сжатие в трещине (о0).
Поведение трещины под действием касательных напряжений показано на рис. 1б. Характерным является приблизительно линейная податливость контакта под действием касательных напряжений т5 не достигающих предельного сопротивления сдвигу (т^,). Угол наклона графика на участке упругого деформирования определяет единичную жесткость трещины при сдвиге (кД Следует отметить, что для нарушений сплошности кп и к3 являются независимыми параметрами, в отличие от континуальной среды.
При разработке модели поведения трещины при сдвиге использовался подход, предложенный в [7] к учету рельефа поверхностей трещин на их сопротивление сдвигу и нормальные к плоскости трещины перемещения, которые могут существенно изменять напряженно-деформированное состояние в окрестности трещины. В соответствии с этим подходом трещины подразделяются на трещины первоначально сомкнутые и трещины, которые уже частично претерпели сдвиг. Первоначально сомкнутые трещины всегда имеют тенденцию к дилатансии (расширению), независимо от направления сдвига. Поведение же трещин уже подвергавшихся сдвигу зависит от направления сдвигающего усилия. Указанный подход, дающий возможность учета в рамках МКЭ, как дилатантного, так и контрак-тантного поведения трещины был использован при разработке алгоритма контактного взаимодействия.
Используя зависимости оп=Ау) и т8=/(м) связь между напряжениями и относительными перемещениями для гладкого контакта на участке упругого деформирования выражается следующим образом:
"к, 0"
т„ - т0
оп -о0
0 кп
(1)
Учет угла дилатансии или контрактансии г (фактически поворот на угол г главной плоскости трещины) дает:
т„ -т0
оп -о0
Кп
кпп
(2)
где к88 = к8 соб2/ + кп б1п21 Кп = к-кп) г эт / (3) кпп = кп соб2/ + к8 эт2/ .
Угол дилатансии определяется по зависимости аппроксимирующей экспериментальные данные и предложенной в работе Б.Е.Ладани [8]:
I — (¡о - ¡т)
От — О
1 = 1„
От _
при О >От
+1
к'.
т при О <От
2 - 4,
(4)
где ¡о - угол дилатансии при о = 0;
От
- напряжение, при котором угол дилатансии становится постоянным. Предел упругости достигнут, когда или равно предельной величине тр или когда вертикальные деформации превосходят значения максимально возможного закрытия или открытия трещины. Таким образом, дополнительные параметры, необходимые для моделирования шероховатых трещин, в этом случае: тип трещины; величина от; угол дилатансии /о; максимальное вертикальное смещение, допустимое при дилатансии Уар (высота выступа шероховатости).
При разработке программы МКЭ, учитывающей нарушения сплошности материала, был взят за основу контактный элемент предложенный Р. Гудманом, Р. Тэйлором, Т. Брекке [6]. Элемент учитывает физическую сущность взаимодействия на контактах, согласно которой значение имеют не абсолютные, а относительные смещения границ блоков и усилия на соприкасающихся границах, в области упругих деформаций, остаются непрерывными и зависят лишь от разности смещений. Формирование матрицы жесткости и вектора сил модифицированного контактного элемента подробно рассмотрены в [1, 5].
Численная процедура учета нелинейных эффектов контактного взаимодействия состоит в следующем. На первом шаге перемещения и напряжения вычисляются в точках интегрирования по характеристикам упругой части зависимостей о„=/(у) и т5=/(и). При этом, если растягивающее напряжение оп(1) (рис. 2, правая часть графика зависимости о„=/(у)) превышает предел прочности на растяжение Яр, контакт нарушается и образуется раскрытая трещина, которая уже не несет никакой нагрузки.
Это приводит к перераспределению напряжений в окружающих элементах. Такое перераспределение достигается путем вычисления узловых сил, ликвидирующих напряжения в трещине, которые прикладываются на следующей итерации и являются, по существу, новым
случаем нагружения.
Рис. 2 Схема итерационного процесса для зависимости оп=/(у).
Очередное распределение напряжений дает новое значение нормального напряжения оп(2) , которое равно:
Оо + А°И(1> (5)
Оп(2) = кп у(2) + О0 + Аоп(1>
аналогично касательное напряжение т^) равно:
т^(2) = К и(2) + т0 + Ат^(1)- (6)
Невязки напряжений вычисляются как разница между упругими (точки оп(1), 2', 3', 4' и т.д. рис. 2) и истинными (в данном случае нулевыми) напряжениями при достигнутых перемещениях:
Аоп(2) = - кп у(2) - о0
Ат^(2) = - к и(2) - т0 (7) Для достижения истинного напряжения требуется несколько итераций. Зависимости (5), (6) и (7) можно записать в общем виде для «/ -ой» итерации:
an(i) = kn v(i) + aC Aan(i) = - kn v(i) -Ts(i) = ks u(i) + TC + Als(i-1> Axs(i) = - ks u(i)
+ Aa
n(i-1);
aC
(8)
Tc. (9)
Сходным образом учитываются нелинейные свойства среды при сдвиге (рис. 3). Невязка касательных напряжений вычисляется по следующим зависимостям:
- при «правостороннем» сдвиге (ts > C ):
ATs(i) = Tr- ks u(i) - Tc; (1C)
- при «левостороннем» сдвиге (ts < C ):
ATs(i) = -Tr - ks u(i) - TC (11)
где Tr = Cr - an(i) tgtyr - остаточная прочность при сдвиге.
Если в процессе деформирования перемещения превышают заданную величину Vmax (рис. 2, левая часть графика зависимости an=f (v)), то на каждой итерации при достигнутом уровне напряжений an(i) корректируются относительные перемещения V(i). Невязка напряжений при этом равна:
Aan(i) = an(i) - kn Vmax - aC (12)
При моделировании дилатантных свойств трещины, находящейся в поле сжимающих напряжений an(;), сначала определяется угол дилатансии i по зависимости (4) с учетом типа трещины и направления сдвига, а затем учитывается поворот на угол i главной плоскости трещины. Связь между компонентами напряжений и относительными смещениями берегов трещины определяется соотношениями (2), (3). Отсюда, при смещении в касательном к трещине направлении на величину u(i), касательные напряжения удовлетворяющие соотношениям (2) равны:
Рис. 3 Схема итерационного процесса для зависимости т s=fU).
Ts(i ) =
k -
U( i) + ^ (с k. V
an ( )-aC
(13)
а величина относительных перемещений в нормальном к трещине направлении (обусловленная перемещением по плоскости наклоненной под углом 1) равна:
kn 1 /
V(i) =--— U(i ) +-IG„ (i) - Oo
(i) 7, 0 ) к V n (i )
k.
(14)
Невязки напряжений, которые необходимо приложить для коррекции полученных на «i-ой» итерации напряжений и перемещений при этом равны:
Aü«(o = ü«(i) - kn v(i)-
Axs(i) = Ts(i) ksu(i)- T0- (15)
Если в ходе решения напряжения xs(i) превысят предельные значения, то невязки напряжений определяются зависимостями (10), (11). Нормальные деформации ограничены, по прежнему, величиной Vmax и дополнительно величиной максимально возможного расширения трещины Vap, обусловленного высотой зубцов шероховатости стенок трещины. В случае если V(i)>Vap, невязка напряжений для коррекции относительных перемещений равна:
A°n(i) = °n(i) - kn Vap - °0 (16)
При уровне сжимающих напряжений on(i) превышающем критическую величину om в зависимости (4), угол дилатансии i принимается равным нулю, а прочность на сдвиг равной прочности на сдвиг скального массива:
= С0 — (1 + ^Е—ЩО*, (17)
т p
n C,
'0
где Со - прочность на сжатие;
п - отношение прочности на сжатие к прочности на растяжение;
т = (1 + п )°'5.
Напряжения (Аоп(г) + о°) и (Ат^) + т°) являются начальными напряжениями для следующей «¡+1» итерации. Вектор узловых усилий {Я} элемента, подсчитанный по значениям этих напряжений, добавляется к вектору сил системы и производится следующее упругое решение с прежней матрицей жесткости, но с новым набором узловых сил. Добавление сил, обусловленных начальными напряжениями, увеличит упругие напряжения в элементе на следующей итерации, однако на величину меньшую, чем начальные напряжения, по которым были рассчитаны узловые силы, поскольку в ансамбле элементов добавленные силы распределяются также и на другие элементы расчетной области. В физическом смысле это означает итерационный поиск таких дополнительных нагрузок, которые сообщают линейно-деформируемому телу перемещения, равные перемещениям нелинейно-деформируемого тела при заданной нагрузке. Итерации продолжаются до стабилизации решения. Достоинством метода является постоянство матрицы жесткости системы, что позволяет лишь однократно формировать и триангулировать ее.
При решении задач в блочных средах, при большом числе элементов «вышедших в пластику», для ускорения сходимости оказывается целесообразным осуществлять коррекцию матрицы жесткости на одной или нескольких промежуточных итерациях, а затем вести итерационный процесс методом начальных напряжений в порядке изложенном выше. Удобно производить коррекцию жесткости при переходе на следующий шаг нагружения или этап при решении задач с поэтапностью, и наделить модель элемента «памятью».
Коррекция жесткости в элементе осуществляется следующим образом. В случае открытой трещины, если нормальная жесткость элемента, соответствующая
достигнутому напряженно-деформированному состоянию на «i» итерации
—- Оо-) -а *
к„!) = „ (■)—- больше жесткости к„ - минимально допустимой по условиям обусловленности системы уравнений МКЭ, первоначальная жесткость снижается до величины к„1+1\ которая принимается равной:
к„!+1) = -0,/^, если k„i+1) >к* или к„+1) = к„ , если к„+1) <к„ , (18)
аналогичным образом корректируется и касательная жесткость открытых трещин
к(г+1) = -т0/и(i), если к(!+1) > к* или к(г+1) = к*, если к?+1) <к* . (19)
s Of s s s s s s
В случае сомкнутых трещин перемещения берегов, которых превысили заданную величину максимально возможного закрытия Vmax, новое значение нормальной жесткости к„1+1 равно:
к„!+1) = < к* (20) max
и не должно превосходить максимального значения к„ , по условиям обусловленности системы уравнений МКЭ.
Касательная жесткость при нарушении сдвиговой прочности (|> тp ) в элементе может быть скорректирована при переходе к следующей итерации следующим образом:
- при "правостороннем сдвиге" (т^} > 0)
т(■) -т
к(!+1) --TR_1° (21)
s ~ (i) u '
- при "левостороннем сдвиге" (т_(г) < 0)
-т('■) -т
к^ц = %r то . (22)
s u(i)
В зависимостях (18)^(22) ао и то - значения нормальных и касательных напря-
«_» 7*7 **
жений с предыдущего шага нагружения; к„ и к„ - пределы варьирования жесткости по условиям обусловленности системы уравнений МКЭ. Рекомендуется принимать к„*>0,001к„ , к„**<1000к„; где к„- первоначальная жесткость в элементе.
При новых значениях к„1+1\ к}1+1 в элементе формируется его матрица жесткости и дальнейшее решение, после сборки и повторной триангуляции глобальной матрицы жесткости, осуществляется с откорректированной матрицей жесткости элемента.
Модель контактного элемента наделена "памятью". Так при переходе на следующий этап нагружения, в элементах, моделирующих открытые трещины, понижается жесткость, а прочность на растяжение Rp и сдвиг Tp принимаются равными нулю. В элементах, моделирующих трещины претерпевшие на предыдущих этапах сдвиговые нарушения, Tp принимается равной остаточной прочности Tr и Rp принимается равной нулю. В случае разгрузки и закрытия первоначально открытых трещин в них восстанавливается касательная и нормальная жесткости и сдвиговая прочность на уровне Tr.
Рассмотренная выше методика и алгоритм расчета сооружений, с учетом контактных взаимодействий по имеющимся нарушениям сплошности, были использованы при исследовании статической работы ряда высоких бетонных плотин
[3, 4, 5]. На основании анализа полученных результатов выполнялось научное обоснование проектных решений на различных стадиях проектирования сооружений, которое заключалось в обосновании выбора типа и конструкции сооружения при вариантном проектировании, выявлении запасов надежности выбранного варианта сооружения, определении резервов его прочности и устойчивости, обосновании технологических процессов, обусловленных схемой возведения сооружения.
Рассмотрим в качестве примера использования методики результаты расчета напряженно-деформированного состояния и устойчивости станционной секции бетонной плотины Богучанской ГЭС.
В плотине высотой 87 метров контактными элементами воспроизводились 11 горизонтальных строительных швов, два вертикальных межстолбчатых шва и контактный шов между сооружением и основанием. Имеющиеся в теле сооружения водовод, отверстия и проемы характеризовались приведенными значениями модуля деформации и объемного веса (рис. 4).
Размеры массива основания принимались равными 380x180 метра. Верхняя более податливая часть массива основания мощностью от 13 до 26 метров имела модуль деформации 6000 МПа, нижняя - 15000 МПа.
Расчетная сетка содержала 693 элемента. С учетом внеузловых степеней свободы в элементах сплошной среды и граничных условий общее число степеней свободы при решении задачи составило 2466.
Строительные швы в теле плотины имели следующие значения деформатив-ных и прочностных характеристик:
- горизонтальные швы к8 = 100 МПа/см; кп=250 МПа/см; С=1,0 МПа; 1§ф=1; Яр=0,1 МПа; Ут£К=0,5 мм;
- вертикальные швы к3 =10 МПа/см; кп=25 МПа/см; С=0,2 МПа; 1§ф=0,85, Яр=0; Утах=1,0 мм с учетом возможного ослабления за счет некачественной цементации.
Контактный шов к3 =50 МПа/см; кп=130 МПа/см; С=0,15 МПа; 1§ф=0,75 Яр=0,1 МПа; Утах=1,0 мм.
Модуль упругости бетона составлял 20000 МПа; объемный вес бетона -24 кН/м3. Прочность на растяжение была снижена (за счет наличия швов не моделируемых явно) и принималась равной 0,38 МПа.
Расчеты выполнялись для совместного воздействия основных статических нагрузок эксплуатационного периода (вес, гидростатическое давление, фильтрационное и взвешивающее давление) и температурного воздействия зимнего периода. Расчетные значения температур в теле сооружения задавались на основе решения температурной задачи, выполненного проектной организацией.
Рис. 4 Фрагмент расчетной сетки МКЭ (копия экрана программы визуализации расчетной схемы)
Среднегодовая температура воздуха принималась равной -3,20С. Максимальное значение расчетной температуры для зимнего периода составляло -330С; со стороны верхнего бьефа на отметке НПУ она равнялась нулю, а на дне водохранилища +3,20С. При расчете температурного поля учитывался фактор проникновения отрицательных температур в полости водовода (при штрабленом профиле в период возведения), которые не были закрыты. Коэффициент линейного расширения бетона а(^) принимался равным 1-10-5 1/град.
Рис. 5 Фрагмент деформированной схемы сооружения в приконтактной области со стороны напорной грани (копия экрана программы визуализации)
Моделирование возможного разуплотнения скальных пород в зоне растягивающих напряжений (область у верховой грани плотины) и трещинообразова-ния в бетоне (вне зоны раскрытия швов) на напряженно-деформированное состояние плотины реализовывалось с помощью упруго-пластической модели работы материала и в элементах сплошной среды [2]. Фактически реализовывалась модель среды не воспринимающей растягивающих напряжений превосходящих прочность на растяжение.
Рассмотрим кратко результаты выполненных расчетов. В упругой стадии
работы системы "плотина-основание" (1 итерация) в плотине фиксируются растягивающие напряжения, вызванные в основном действием отрицательных температур. На низовой грани плотины растягивающие напряжения у достигают величин 13,3 кг/см2 (главные О! до 20,5 кг/см2). На напорной грани максимальное растяжение оу составляет 8,8 кг/см2 в зоне порога водоприемника и 1,12 кг/см2 в бетоне выше контактного шва. В контакте растяжение составляет 13,6 кг/см2. Максимальное сжатие в бетоне плотины достигает величины о2 = -28,7 кг/см2 в начале наклонного участка подошвы плотины.
В основании с верховой стороны плотины формируется зона двухосного растяжения с максимальной величиной растяжения 01 до 15,6 кг/см2. Полученная картина напряженного состояния показывает, что и в сооружении и в основании (при заданных прочностных параметрах) возможны раскрытия швов и контакта и трещинообразование, как в бетоне, так и в массиве основания.
Для достижения сходимости решения при одновременном моделировании контактного взаимодействия в швах и трещинообразования в бетоне и скале потребовалось 25 итераций (на 5 и 10 итерациях выполнялась коррекция жесткости в контактных элементах).
Максимальное горизонтальное перемещение гребня плотины составило 5,57 см (5,1 см на 1 итерации). На рис. 5 приведен фрагмент деформированной сетки (в увеличенном масштабе перемещений) в приконтактной области плотины со стороны напорной грани на котором отчетливо видно произошедшее раскрытие контактного шва. Глубина раскрытия составила 5,0 метров при величине раскрытия до 3,1 миллиметра. На рисунках 6 и 7 в виде цветных заливок приведены главные напряжения в плотине и примыкающей области основания. Максимальное растяжение в плотине не превышает 3,4 кг/см2
в пригребневой зоне плотины со стороны низовой грани. Максимальное сжатие в бетоне плотины несколько увеличилось и достигло величины Ст2 = -29,8 кг/см2. Вертикальные напряжения вдоль напорной грани вблизи основания сжимающие. В области водоприемника (где фиксировалось максимальное растяжение в упругой задаче) напряжения практически нулевые, так как произошло раскрытие строительных швов в этой зоне и трещинообразование в бетоне расположенном между швами. Раскрытие строительных швов и трещинообразование в бетоне произошло и со стороны низовой грани плотины, где растягивающие напряжения, вызванные действием температуры, превысили прочность на растяжение в бетоне и строительных швах. На рис. 8 показано состояние строительных швов, контактного шва и области тре-щинообразования в бетоне плотины и скальном основании. Глубина раскрытия швов и зон трещинообразования со стороны низовой грани составляет от 1,5 до 2,5 метров при величинах раскрытия швов до 0,05 миллиметра (если не моделировать трещинообразование в бетоне раскрытие швов составляет 1-2 мм). На напорной грани глубина раскрытия швов составляет от 0,8 до 5 метров при величине раскрытия до 0,4 миллиметра. В вертикальных ослабленных швах имеются локальные зоны сдвигового нарушения.
В основании область разуплотнения распространяется под подошву сооружения на 5 метров. На наклонной части подошвы контактный шов также раскрыт но величина раскрытия не превышает 0,4 мм.
На рис. 8 приведены информационные панели, показывающие деформации и напряжения в контактном элементе (в данном случае элемент в контакте со стороны напорной грани) и позволяющие подсчитать коэффициент устойчивости по любому шву. На рисунке
видно, что при заданных сдвиговых параметрах в контактном шве коэффициент устойчивости составил 1,31, что удовлетворяет требованиям норм.
Таким образом, по результатам расчета, устойчивость станционной секции от воздействия основных статических нагрузок и зимней температуры обеспечивается с коэффициентом устойчивости выше нормативного.
Рис. 6. Напряжения Ст! (кг/см2) в плотине и примыкающей области основания (копия экрана программы визуализации)
Рис. 7. Напряжения Ст2 (кг/см2) в плотине и примыкающей области основания (копия экрана программы визуализации)
Рис. 8 Характер работы строительных швов, контактного шва и положение зон трещинообразования
Выводы
1. Выполненные расчеты напряженно-деформированного состояния бетонной гравитационной плотины Богучан-ской ГЭС и ряда других показали, что использование разработанной методики, учета в расчетах МКЭ нелинейных эффектов контактного взаимодействия (с использованием контактных элементов), является весьма эффективным способом оценки влияния на статическую работу сооружения различных строительных и конструктивных швов, ослабляющих тело сооружения, прочностных и деформационных характеристик швов и контакта сооружения с основанием, а также отдельных трещин или систем трещин в скальном основании.
2. Получаемые оценки предельного состояния строительных и деформационных швов при их различных прочностных и деформативных характеристиках, количественные оценки глубины и величины раскрытия различных швов, в том числе и контактного при статических и температурных воздействиях, размеры зон разуплотнения скального основания и их влияние на состояние
контактного шва позволяют более объективно и детально оценить напряженно-деформированное состояние плотины и основания, а также реализовать в проектных разработках гидротехнических сооружений инженерные идеи с надежным расчетным обоснованием.
Литература
1. Орехов В.Г., Зерцалов М.Г., Толстяков В.В. и др. Исследование схемы разрушения системы «бетонная плотина - скальное основание». - Известия ВНИИГ им. Б.Е. Веденеева, 1987, т. 204, с.71-76.
2. Зерцалов М.Г., Толстиков В.В. Учет упруго-пластической работы бетонных плотин и скальных оснований в расчетах с использованием МКЭ. - Гидротехническое строительство, №8, 1988, с.33-36.
3. Орехов В.Г., Захаров В.Ф., Зерцалов М.Г., Толстиков В.В., Шимельмиц Г.И. Напряженно-деформированное состояние плотины Бурейской ГЭС с учетом очередности ее возведения. - Материалы конференций и совещаний по гидротехнике «ПТТС-88». - Л.: Энергоатомиздат, 1989, с.46-49.
4. Орехов В.Г., Захаров В.Ф., Толстиков В.В., Шимельмиц Г.И. Напряженно-деформированное состояние бетонной плотины гидроузла Капанда с учетом работы ослабленных строительных швов. - М.: Энергетическое строительство, №1, 1993, с.26-29.
5. Орехов В.Г., Зерцалов М.Г. Механика разрушения инженерных сооружений и горных массивов. - М: АСВ, 1999.
6. Goodman R.E., Taylor R.L., Brekke T.L. A model for mechanics of jointed rock. - Journal of the soil mechanics and foundation division. Proceeding of the ASCE, 1968, SM3, №5, p.637-659.
7. Goodman R.E., Dubois J. Duplication of dilatancy in analysis of jointed rocks. -Journal of the soil mechanics and foundation division, ASCE, 1972, SM4, p.399-422.
8. Ladany B., Archambault G. Simulation of Shear Behavior of a jointed Rock Mass. Proc. of Symp. On Rock Mech., ASME, 1970.
