Бураго Н.Г. Вычислительная механика (1185926), страница 34
Текст из файла (страница 34)
главу проинтерполяцию). остальные переменные аппроксимируем кусочнопостоянными функциями со значениями в центрах элементов.Интегралы в вариационных уравнениях (1) вычислим с помощьюквадратурных формул прямоугольников. Точки численногоинтегрирования для интегралов по области V от виртуальнойработы массовых сил и сил инерции будем брать в узлах сетки, дляинтеграла виртуальной работы напряжений – в центрах элементов.Для интегралов по границе точки численного интегрированиявозьмем в центрах граничных элементов.18.3.
Явные лагранжевы схемыДля понимания работы явных лагранжевых схем достаточнорассмотреть конечно-элементное представление классической явнойсхемы Уилкинса (1967). В соответствии с этой схемой каждый шагпо времени рассчитывается в три этапа: на первом этапе материалсчитается упругим, определяются новые скорости и координатыузлов, на втором этапе учитывается влияние контактныхвзаимодействий на скорости и координаты узлов, на третьем этапе всоответствии с уравнениями теории пластичности определяютсянапряжения, пластические деформации и внутренняя энергия.Соотношения первого этапа имеют вид:M i ( v in +1 − γ in v in − (1 − γ in ) vɶ in ) = Fin ∆tn ( i ∈ Ω \ ΩV )vn +1i=vn +1*i ,( i ∈ ΩV )197(6)(7)Глава 18.
Методы для задач упругопластичностиxin +1 = xin + v in +1∆tn ,(i∈Ω )(8)гдеNC M CM i = ∑∑Vk0 [ ρ ]0k M C−1 Hɶ (i − C (k , l )) , ( i ∈ Ω )(9)vɶ in = 0.5(max v nj + min v nj )(10)k =1 l =1j∈Ωi(i∈Ω )j∈Ωiγ nj = min(1, a1 + a2 | vɶ nj − v nj |)(11)NC M CFin = ∑∑ [V ]nk [σ ]nk ⋅∇ nkl Hɶ (i − C (k , l )) +k =1 l =1NB M B+ ∑∑ M B−1[P]kn Skn Hɶ (i − B(k , l ))(12)k =1 l =1Здесь (9) – формула для определения "узловых масс" M i , (10) –формула для осреднения скорости узла i по шаблону, (11) – формуладля коэффициента гибридности, (12) – формула для определения"узловых сил".Поясним формулу для коэффициента гибридности.
Модульразности | vɶ nj − v nj | пропорционален второй производной отскоростей по пространственным переменным в окрестности узла i ,коэффициенты a1 и a2 принимались равными a1 = vS / c иa2 = (0.2vS ) −1 , где vS = max | v in | − min | vin | - диапазон измененияi∈Ωi∈Ωмодуля скорости в области решения, c = dp / d ρ + 4 / 3µ / ρ скорость звука.В соответствии с формулой (11) сглаживание производится вузлах с сильным отклонением вычисленного значения скорости отсреднего по шаблону значения.Поясним значения новых величин в формулах (6)-(12): [ S ]nk площадь граничного элемента k , определяемого узлами B (k , l ) ,[V ]nk - объем k -й ячейки в актуальной конфигурации. Векторы ∇ nklпредставляютдискретныйоператорпространственногодифференцирования в актуальной конфигурации (см.
главу прочисленное дифференцирование). M i - узловые массы, вычисляемыеодин раз перед началом расчета, [V ]0k - начальный объем элемента(ячейки), M C - число узлов в элементе, N C - число элементов,функция Hɶ равна единице для нулевого значения аргумента инулю в противном случае, она указывает адрес рассылки вкладов от198Глава 18. Методы для задач упругопластичностимассы элемента в приузловые массы, Ωi - номера узлов соседейузла i (шаблон узла i).На втором этапе сначала определяется зона контакта. Дляэтого исследуется положение каждого граничного узла поотношению к каждому граничной ячейке. Если по результатампервого этапа расчета граничный узел проник сквозь граничнуюячейку внутрь деформируемого тела, то это явилось следствиемтого, что в расчете не были приняты во внимание контактные силы.Эти контактные силы определяются так, чтобы при учете ихвлияния на скорости и положения данного узла и ячейки устранитьобнаруженное счетное проникание.Определение зоны контакта подразумеваетопределениемножества контактных пар k ∈ Ωcont .
Имеются в виду пары"граничный узел (C(k,0)) – граничная плоская ячейка(C(k,1),C(k,2),C(k,3))", где C(k,i) – номера узлов в контактных парах.Алгоритм поиска имеет вид:цикл по граничным узламцикл по плоским граничным ячейкам- если номер узла совпадает с номером узла в ячейке,то контакта нет;- если узел находится от центра ячейки далее чемхарактерный размер ячейки по одной из осей координат, то контактанет;- если объем пространственной ячейки образованныйграничным узлом и граничной ячейкой больше нуля, то контактанет;- если нормаль, опущенная из узла на плоскостьграничной ячейки, не пересекает ячейку, то контакта нет;- граничный узел и граничная плоская ячейкаобразуют контактный конечный элемент или контактную пару.конец цикла по граничным плоским ячейкамконец цикла по граничным узламРасчет контакта.
Нелинейная система уравнений дляопределения скоростей, координат и нормальных контактныхусилий имеет видn +1 ( k )v Cn +(1k ,i ) = v 0nC+1( k ,i ) + FkCLi / M C ( k ,i ) ∆tn ( k ∈ Ω cont ; i = 0,1, 2, 3 )xCn +(1k ,i ) = xCn ( k ,i ) + v Cn +(1k ,i ) ∆tn( k ∈ Ω cont ; i = 0,1, 2, 3 )[(xCn +(1k ,2) − xCn +(1k ,1) ) × (xCn +(1k ,3) − xCn +(1k ,1) )] ⋅ (xCn +(1k ,0) − xCn +(1k ,1) ) = 0 ( k ∈ Ωcont )где контактное усилие представлено нормальными и касательными кгранице составляющими199Глава 18. Методы для задач упругопластичностиn +1FkC= FkNn +1n kn +1 + FkTn +11 τ kn1+1 + FkTn +21τ kn +21Усилия контактного трения FkTn +11 и FkTn +21 определяются скачкомкасательной скорости и контактным давлением по закону трения, аконтактное давление FkNn +1 определяется из условия непроникания(из равенства нулю объема контактного элемента).В алгоритме решения совершается обход контактнойграницы и условие непроникания линеаризуется по неизвестнымнормальным контактным усилиям с учетом малости контактныхпоправок координат узлов по сравнению со значениями самихкоординат.
Условие непроникания разрешается по очереди длякаждого контактного элемента. Таким образом, используютсяитерацииГаусса-ЗейделявместесметодомНьютона.Соответствующие найденным контактным усилиям поправкивносятся тут же в величины координат и скоростей контактныхузлов.
Итерирование контактных усилий оканчивается приустранении счетного проникания. В явных схемах для этогодостаточно двух обходов контактной границы.На третьем этапе в каждой ячейке k области решенияопределяютсяновыескоростидеформаций,напряжения,пластическая работа и внутренняя энергия:n +1k[ L]MC= ∑ ∇ kln vCn +(1k ,l ) ,[e]nk +1 = 0.5([ L]nk +1 + [ LT ]nk +1 )l =1[e m ]nk +1 = [e]nk +1 : I / 3 ,[e ']nk +1 = [e]nk +1 − [e m ]kn +1 I[σɶ ']nk +1 = [σ ']kn + 2 µ[e ']kn +1 ∆tn , [σ ']nk +1 = [σɶ ']nk +1α kn +1α kn+1 = σ s / [σɶ ']nk +1 :[σɶ ']nk +1 ,[a p ]kn +1 = [a p ]nk + σ s (1 − α kn +1 )α kn +1[U ]nk +1 = [U ]nk + [σ ]nk +1 :[L]kn +1 ∆tn /[ ρ ]kn( k ∈ ΩC )(13)Описанный расчет напряжений, пластических деформаций ипластической работы предложен Уилкинсом и известен какалгоритм “посадки напряжений на круг текучести”.Метод Уилкинса устойчив при выполнении критерияКуранта200Глава 18.
Методы для задач упругопластичности1∆tn = min nk∈ΩC c max(| ∇ n ⋅ e |,| ∇ n ⋅ e |) kl1kl2 k l(14)и условия точности, выражающего требование достаточной малостиприращений деформации(∆tn ≤ ∆ε max max ([e]nk :[e]kn )где1/ 2k∈ΩC∆ε max = 0.1ε Y)−1(15)- максимально допустимое приращениедеформации на шаге по времени, ε Y - деформация, отвечающаяпределу текучести.В оригинальной схеме Уилкинса для устранениямелкомасштабных осцилляций решения используется явнаяискусственная вязкость в виде вязкостного давления. Эта вязкостьтребует подбора коэффициентов в зависимости от решаемой задачи.Поэтому здесь она заменена не требующим регулировкинелинейным лаксовым сглаживанием.В оригинальной схеме Уилкинса дискретизированныеуравнения движения получены исходя из интегральнойформулировки уравнений сохранения импульса и внутреннейэнергии интегро-интерполяционным методом.
Результирующиеуравнения практически совпадают с выписанными здесь.В качестве варианта рассмотренной схемы приведем схемуквазивторого порядка точности. Для гиперболического уравненияпереноса она описана в книге Поттера (1975). Она отличается отсхемы Уилкинса расчетом первого этапа, который в этом случаеимеет вид:M i ( v in +1 − v in ) = ((1 + γ )Fin − γ Fin −1 )∆tn( i ∈ Ω \ ΩV )v in +1 = v*ni+1 ,( i ∈ ΩV )xn +1in +1i= x + v ∆tn ,ni(i∈Ω )(16)где 0.0 < γ ≤ 1.0 - малое положительное число. Аналогом этойсхемы интегрирования по времени для случая задач Коши дляобыкновенных дифференциальных уравнений является схемаАдамса-Башфорта (трехслойная схема).
Условия устойчивости дляэтой схемы, как показывает анализ Фурье для модельной задачи суравнением переноса, являются вдвое более ограничительными,нежели условия Куранта (13). Достоинством схемы являетсяпростота управления вязкостью, осуществляемая единственным201Глава 18. Методы для задач упругопластичностипараметром γ (не надо "химичить" с искусственной вязкостью, т.е.не надо изменять вид вязких членов и значения коэффициентовискусственной вязкости в зависимости от условий задачи).Дополнительно требуется держать в памяти ЭВМ вектор узловыхсил со старого слоя Fin −1 . На первом шаге по времени (n=0)полагают Fi−1 = Fi0 . При γ = 0.5 схема имеет второй порядокточности.Другим примером видоизменения явной лагранжевой схемыУилкинса является полностью консервативная явная схема,предложенная Самарским и Поповым (1980) применительно кзадачам лагранжевой газовой динамики и выведенная для расчетаупругопластических сред Тишкиным (1984) из вариационногопринципа наименьшего действия Гамильтона.
В нашихобозначениях она может быть записана в следующем виде:M i ( v in +1 − v in ) = Fin ∆tn ( i ∈ Ω \ ΩV )v in +1 = v*ni+1 ,( i ∈ ΩV )xin +1 = xin + v in +1/ 2 ∆tn , ( i ∈ Ω )k ∈ ΩC :[σɶ ']nk +1 = [σ ']nk + 2 µ[e]nk +1/ 2 ∆tn , [σ ']nk +1 = [σɶ ']kn +1α kn +1α kn+1 = σ s / [σɶ ']kn+1 :[σɶ ']kn +1 ,(17)[a p ]kn +1 = [a p ]nk + σ s (1 − α kn +1 )α kn +1[U ]nk +1 = [U ]nk + [σ]nk +1/ 2 :[L]nk +1/ 2 ∆tn /[ ρ ]nkгде дробный верхний индекс подразумевает следующее действие[ f ]n +1/ 2 = ([ f ]n +1 + [ f ]n ) / 2(18)Для устойчивости в схеме используется малая искусственнаявязкость в виде тензора вязких напряжений:[σ v ]nk = −I[ pv ]nk + [σ 'v ]nk , [ pv ]kn = − K [∇ ⋅ v ]kn ∆tn / 2[σ 'v ]nk = 2 µ [e ']nk ∆tn / 2(19)Устойчивость схемы при наличии вязких членов легкоустанавливаетсяизанализапервогодифференциальногоприближения, который и подсказывает вид выражений (19).202Глава 18. Методы для задач упругопластичностиОтметим, что при больших скоростях удара в схемы (15),(16)на ударных волнах для устранения осцилляций приходится вводитьлаксово сглаживание (см.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
