Transient dynamics of 3D inelastic heterogeneous media analysis by the boundary integral equation and the discrete domains methods

Full Text

Введение. Изучение переходных процессов в трехмерных нелинейно деформируемых средах с локальными особенностями в геометрии, структуре или режимах нагружения является актуальной и достаточно сложной проблемой. Под действием ударных волн происходят изменения геометрии и объема, образование и локализация повреждений в зонах концентрации с последующим катастрофическим разрушением, страгивание и распространение существующих дефектов или образование множественных очагов разрушения для хрупких материалов и др. [1]. В большинстве случаев реальных нагрузок изменения в объеме (плотности) сред оказываются при этом незначительными и возникновение указанных предельных состояний происходит вследствие нелинейных деформаций сдвига. Очевидно, что анализ таких состояний возможен только с использованием численных методов, иногда в сочетании с экспериментальными методами, если ударным воздействиям подвергаются особо ответственные объекты техники. Наиболее широко для решения подобных задач применяется метод конечных элементов (МКЭ), его достоинства неоспоримы и хорошо известны. Однако связанные с ним вычислительные затраты оказываются чрезмерными, если требуется детальное описание геометрии и представление решения с большими градиентами изменения полевых функций, например для составных конструкций или их систем при наличии зон концентрации в виде локальных особенностей в геометрии или дефектов, трещин и т. п. В силу малости этих зон соответствующие краевые задачи часто переходят в разряд внешних для неограниченных областей с граничными условиями «на бесконечности», требующими для их учета в МКЭ неординарных подходов. Более эффективным в этих случаях оказывается применение метода граничных интегральных уравнений (МГИУ), чаще именуемого в зарубежной литературе методом граничных элементов (МГЭ). В основе его лежат интегральные представления соответствующих краевых задач, что позволяет свести их решения к многообразиям меньшей размерности и определять с более высокой точностью, чем в МКЭ, искомые полевые функций и их градиенты внутри расчетных областей [2-5]. Граничные условия на «бесконечности», или условия Зоммерфельда, выполняются при этом автоматически. Основные трудности в применении МГИУ для рассматриваемых задач связаны, прежде всего, с наличием фундаментальных решений, а также необходимостью многократного по времени вычисления граничных интегралов и обращения несимметричных плотно упакованных матриц для формирования и решения систем разрешающих алгебраических уравнений большого порядка. Поэтому очень важна разработка эффективных схем вычисления таких интегралов, особенно почти сингулярных, которые возникают при анализе переменных физического поля вблизи границы или изучении тонких пограничных структур. С этой целью могут быть использованы численные схемы, разработанные для решения не зависящих от времени краевых задач эллиптического типа [2, 5] и др. На гладких явно параметризованных границах обычно применяются численные схемы с квадратурами высокого порядка. Все эти схемы становятся более эффективными при их использовании совместно с так называемыми быстрыми методами формирования и реше138 Изучение переходных процессов в нелинейно деформируемых средах. . . ния больших систем уравнений [3, 6]. Когда эти методы применимы, число операций O(N 2 ), обычно необходимых на формирование, и O(N 3 ) на решение систем уравнений удается снизить до O(N ), где N - порядок системы. Тем самым расширяются возможности МГИУ для моделирования переходных процессов в реальных объектах техники. К настоящему времени сформировались два основных направления в применении интегральных представлений для решения задач волновой динамики нелинейно деформируемых сред. В первом используются формулировки с зависящими от времени фундаментальными решениями Стокса и их аппроксимация по пространству и времени, во втором - с фундаментальными решения Кельвина для задач статики совместно с шаговыми методами по времени. Детальный обзор состояния и перспектив развития этих направлений представлен в недавних работах [3, 7-9]. Строгое математическое обоснование соответствующих формулировок метода, берущее начало с работ Купрадзе [10], изложено в ряде монографий [5, 7]. Однако примеры успешного применения метода ограничиваются в основном двухмерными приближениями. Для ядер, зависящих от времени, формулировка метода оказывается достаточно сложной и имеет проблемы с устойчивостью. Использование не зависящих от времени статических фундаментальных решений совместно с шаговыми методами аппроксимации по времени делает формулировку менее сложной, но трудоемкой с вычислительной точки зрения, поскольку требует многократного вычисления объемных интегралов с инерционными силами и нелинейными функциями физического поля внутри расчетной области. Однако это компенсируется простотой ее реализации на ЭВМ [11]. В предлагаемой статье представлены теория и основные элементы компьютерной реализации МГИУ, определяющие эффективность его применения для моделирования нелинейной волновой динамики неоднородных деформируемых сред, подвергаемых ударным воздействиям. Формулировка метода основана на успешно использованных в [12], не зависящих от времени фундаментальных решениях в интегральном представлении соответствующих краевых задач и шаговых методах интегрирования во времени. В решении используются методы дискретных областей со схемами высокого порядка аппроксимации границ и коллокационного приближения. Такой подход органически вписывается в современные вычислительные технологии для многопроцессорных систем [3]. В подтверждение его эффективности приведены решения модельных задач с результатами, полученными МКЭ, и задач, имеющих практическое значение. 1. Постановка задачи и основные соотношения МГИУ. Рассматривается изотропная вязко-упругопластическая среда объема V , занимающая в R3 произвольную, в общем случае многосвязную область D, ограниченную поверхностью Липшица n+1 Sr, ∂D = r где S n+1 - поверхность, внешняя по отношению ко всем остальным. Область D (тело, конструкция) может быть неоднородной, составленной из различных материалов, и включать в себя локальные геометрические осо139 П е т у ш к о в А. В. бенности, например типа углов, отверстий или поверхностей разрыва типа трещин Ωc (ˆ xα ), где x ˆα , α = 1, 2, - локальная система координат, связанная с поверхностью разрыва, рис. 1. Здесь и далее используются обозначения и соглашения, принятые в тензорном исчислении. Пусть среда деформируется под действием быстро меняющихся во времени массовых Xj (xi , t), xi ∈ D, поверхностных pj (xi , t), xi ∈ ∂Dσ сил и перемещений uj (xi , t), заданных на части поверхности c xi ∈ ∂Du , а также известного температурного поля T (xi , t), xi ∈ D. Следовательно, ∂D = ∂Du ∪∂Dσ , при этом ∂Du ∩ ∂Dσ = 0, а t ∈ Dt = (0, τ ), где τ - продолжительность истории динамического нагружения. Состояние среды в момент времени t = 0 может быть как недеформированным, так и включать в себя начальные или остаточные напряжения и деформации, в том числе технологического происхождения. Поведение среды при указанных воздействиях характеризуется полями перемещений uj (xi , t), деформаций εjk (xi , t) и напряжений σjk (xi , t) или их скоростями, соответственно, u˙ j , ε˙jk , σ˙ jk . Для его описания ограничимся так называемым квазиакустическим приближением. В этом случае актуальное давление p на фронте ударной волны полагается функцией только объема, которую можно принять в следующем виде: p = K(ρ - ρ0 )/ρ, где K = (2µ+3λ)/3 - так называемый объемный модуль, а µ и λ - параметры изотропной среды (параметры Ламе). Ударная адиабата совпадает с изоэнтропой расширения, а необратимость процессов деформирования обусловлена в основном вязкостью и пластическими деформациями сдвига. Полагая деформации малыми, тензор скорости деформаций Коши ε˙jk представим в виде суммы упругой ε˙ejk и начальной ε˙0jk составляющих: ε˙jk = 1 (u˙ j ,k + u˙ k,j ) = ε˙ejk + ε˙0jk , 2 (1) причем ε˙0jk , в свою очередь, может включать в себя как пластические деформации ε˙pjk , так и любые другие, например технологические ε˙∗jk . Отнесем к числу перечисленных и температурные деформации ˙ ε˙θjk = αδjk θ, Рис. 1. Схема к постановке задачи [Figure 1. The problem statement scheme] 140 Изучение переходных процессов в нелинейно деформируемых средах. . . где α - коэффициент температурного расширения, δij - символ Кронекера, θ˙ - разность температур. Соответствующее поле скоростей напряжений определяется обобщенным законом Гука σ˙ jk = Λjklm ε˙lm - ε˙0lm , (2) где тензор упругости для изотропной среды Λjklm = λδij δkl + µ(δik δjl + δil δjk ), а параметры λ, µ в общем случае могут зависеть от температуры. Представим их в форме, удобной для последующего использования: λ= Eν , (1 + ν)(1 - 2ν) µ= E , 2(1 + ν) где E и ν - обычные (технические) модуль Юнга и коэффициент Пуассона. Для вычисления скоростей пластических деформаций, входящих в (1), используем модифицированные соотношения теории течения упругопластических сред с анизотропным упрочнением [13, 14]. В этом случае полная деформация εjk с учетом классических условия течения Мизеса и условий нагружения и разгрузки представляется в виде εjk = εejk + εθjk + εpjk + ∆εpjk , где приращения пластической деформации ∆εpjk определяются из соотношения ∆εpi ∆εpjk = Ejk . (3) Ei Здесь Ejk - девиатор тензора деформаций: Ejk = ejk - 1 1 ρjk = (Sjk - ρjk ) + ∆εpjk , 2µ 2µ Ei - интенсивность так называемой смещенной деформации: Ei = 2 Ejk Ejk 3 1/2 , ∆εpi - интенсивность приращения пластических деформаций: ∆εpi = 3µEi - σin-1 3µ + H n-1 -1 , причем 1 ejk = εjk - εll δjk , 3 ρ˙ jk = H ε˙pjk , 1 Sjk = σjk - pδjk , a p = - σll - давление, 3 H= dσi dεpi σi = ∗ εi =εi 3 Sjk Sjk 2 1/2 . 141 П е т у ш к о в А. В. Интенсивность приращения пластических деформаций ∆εpi определяется на каждом n-ном шаге нагружения деформируемой среды диаграммой деформирования σi ∼ εpi для соответствующей температуры θ и/или скорости деформаций ε˙i для сред, чувствительных к скорости деформирования. При этом полагаются справедливыми следующие соотношения: g(x ˙ i , t) = ∂g/∂t = dg/dt и ∆g = g∆t, ˙ (4) где g(xi , t) - любая переменная физического поля. Следуя законам сохранения, уравнение волновой динамики нелинейно деформируемой среды запишем в обобщенном виде с учетом сил вязкого демпфирования: c21 - c22 ui,ij + c22 uj ,ii + Xj = u ¨j + γ u˙ j , (xi , t) ∈ D × Dt . (5) Здесь запятая перед индексом обозначает производную по соответствующей координате, u ¨j - вектор ускорения, γ = c/ρ - коэффициент вязкого демпфирования. Свойства материала деформируемой среды представлены волновыми скоростями c21 = (λ + 2µ) /ρ и c22 = µ/ρ - соответственно для волн расширения/сжатия и сдвига, а начальные и нелинейные составляющие деформирования отнесены к массовым силам Xj (xi , t): Xj (xi , t) = Xj - 222 ε0kk,j - 2 1 - 222 ε0kj ,i . Для однозначного решения уравнения (4) необходимо располагать граничными и начальными условиями. Граничные условия, записанные в форме Неймана и Дирихле для усилий и перемещений, соответственно имеют вид pj = σjk nk = pbj (xi , t) на ∂Dσ × Dt , uj = ubj (xi , t) на ∂Du × Dt , (6) где nj - компонента вектора внешней нормали к поверхности ∂D в точке xi . Выражения для тензора напряжений σjk и вектора усилий pj можно записать в виде, удобном для последующего использования: σjk = ρ(c21 - 2c22 )δjk ul,l + ρc22 (uj ,k + uk,j ), pj = σjk nk = ρ(c21 - 2c22 )δjk nk ui,i + ρc22 (uj ,k nk + uk,j nk ). (7) (8) Неоднородные, отличные от нуля, начальные условия имеют вид uj (xi , 0) = u0j (xi ), u˙ j (xi , 0) = vj0 (xi ), xi ∈ D. (9) Для описания медленно меняющихся во времени или квазистатических процессов деформирования, когда u˙ j (xk , t) cα , α = 1, 2, соответствующая нелинейная краевая задача следует из (4)-(9), если положить ρ¨ uj = 0. В этом случае время t играет роль только индекса, характеризующего изменения в истории нагружения. При наличии в области D внутренних известной формы границ Γ, обусловленных соединением разнородных материалов или составных тел, условия (5) дополняются следующим соотношением: u˙ 1j = u˙ 2j , 142 p˙1j = p˙2j на Γ × Dt , (10) Изучение переходных процессов в нелинейно деформируемых средах. . . где цифрами 1, 2 обозначены тела (материалы), находящиеся по обе стороны от границы Γ, которая в данном случае не является обычной физической. Если соединение допускает относительное смещение тел, то внутренняя граница Γ между ними рассматривается в виде первоначально совпадающих поверхностей Γ1 , Γ2 тел в зоне контакта, рис. 1. В отличие от статики, в задачах динамического деформирования вместо усилий pj слева и справа от границы Γ используются их импульсы. Отрыву (расслоению) общих точек поверхностей Γ1 , Γ2 соответствуют условия равенства в них нулю усилий: p˙1j = p˙2j = 0, а скольжению с трением - следующие условия на (Γ1,2 × Dt ): u˙ 1n = u˙ 2n , p˙1n + p˙2n = 0, p˙1τ + p˙2τ = 0, -κp˙1,2 n p1τ , p˙2τ κp˙n1,2 , (11) где n и τ отмечены соответственно нормальное и касательные направления к поверхностям Γ1 , Γ2 в точке контакта, κ - коэффициент трения. В этом случае краевая задача (4)-(10) оказывается нелинейной и при упругом (линейном) поведении деформируемой среды. 2. Интегральные представления краевой задачи. Используя метод взвешенных невязок или обобщенную теорему взаимности Бетти (иначе - тождество Грина), краевой задаче (4)-(10) можно поставить в соответствие пространственно-временное интегральное уравнение с зависящими от времени фундаментальными решениями для смещений и сил. В основе его лежит представление решения uj (xi , t) задачи в виде зависящего от времени поверхностного потенциала с заданной плотностью и свойства непрерывности подобных потенциалов [3, 5, 10]: t Cij (x)uj (x, t) = - Pij (x, t; y, τ ) ∗ uj (y, τ )dsy dt+ 0 ∂D t Uij (x, t; y, τ ) ∗ pj (y, τ )dsy dt+ + 0 ∂D t Uij (x, t; y, τ ) ∗ Xj (y, t) - u ¨j (y, τ ) - γ u˙ j (y, τ ) dvy dt+ + 0 V t Σjki (x, t; y, t) ∗ ε0jk (y, τ )(x, t; y, τ )dvy dt. (12) + 0 V Напряжения σij (x, t) для точек внутри области D определяются соответствующим интегральным представлением, однако они могут быть получены с использованием (10) и уравнения состояния для деформируемой среды (2), что оказывается более эффективным с вычислительной точки зрения. Здесь y ∈ ∂D - точка поверхности; x ∈ R3 \∂D - точка (источник) внутри области D; Vp , V - соответственно области пластического деформирования и/или ненулевой разности температур и отличными от нуля объемными силами; Cij (x) = δij , если x ∈ D, и Cij (x) = 0.5δij , если x ∈ ∂D, при условии ее гладкости; знак (∗) означает риманову свертку по времени, которая в этом случае 143 П е т у ш к о в А. В. обладает всеми групповыми свойствами и определяется соотношениями t f (xi , t - τ )g(xi , t)dτ f (xi , t) ∗ g(xi , t) = ∀(xi , t) ∈ D × Dt , 0 f (xi , t) ∗ g(xi , t) = 0 ∀(xi , t) ∈ D × Dt . Динамические фундаментальные решения Pij , Uij и Σjki , входящие в (12), являются сингулярными при |x - y| → 0, обладают свойствами симметрии, взаимности, причинности и не зависят от начала отсчета по времени. Кроме того, они автоматически удовлетворяют условию излучения Зоммерфельда для бесконечных областей [5]. В явном виде выражения для них достаточно громоздкие и здесь в целях экономии места не приводятся, их можно найти, например, в [2]. Решения интегральных уравнений для перемещений uj (xi , t) и напряжений σij (x, t) ищутся в пространстве L2 (D × Dt ) интегрируемых с квадратом функций и их обобщенных производных и строятся численно МГИУ с использованием коллокационного приближения и аппроксимаций по пространству и времени в D × Dt . Однако этот обычный в МГИУ подход, основанный на зависящих от времени фундаментальных решениях, имеет известные проблемы с устойчивостью и оказывается весьма трудоемким и малоэффективным с вычислительной точки зрения [3, 9, 11]. Причина заключается в достаточно сложных выражениях для ядер, необходимости многократного вычисления интегралов свертки Римана (12) с обеспечением при этом причинно-следственной связи и устойчивости в течение всего времени интегрирования, а также объемных интегралов для изменяющихся во времени нелинейных деформаций, инерционных и вязких сил. Более простым и эффективным оказывается подход, основанный на применении фундаментальных решений трехмерной задачи статики и шаговых методов интегрирования по времени [11]. Примеры его успешного использования для решения трехмерной задачи упругой динамики во временной и частотной областях были показаны в [12]. Здесь мы обобщим этот подход, включив учет нелинейного поведения деформируемых сред, как в [14, 15]. Интегральное представление для краевой задачи (5)-(10) в этом случае формально имеет вид, аналогичный (12): Cij (x)uj (x, t) = - Pij (x, y)uj (y, t)dsy + ∂D Uij (x, y)pj (y, t)dsy + ∂D Uij (x, y)ρ u ¨j (y, t) + γ u˙ j (y, t) - Xj (y, t) dv- + V Σjki (x, y)ε0jk (y, t)dv. (13) - Vp Однако входящие в него ядра - тензоры Uij , Pij , Σjki - теперь соответствуют фундаментальному решению Кельвина-Сомильяны и определяются достаточно простыми выражениями: 144 Изучение переходных процессов в нелинейно деформируемых средах. . . Uij = [16πµ(1 - ν)r]-1 [(3 - 4ν)δij + r,i r,j ] , -1 ∂r Pij = - 8π(1 - ν)r2 [(1 - 2ν)δij + 3r,i r,j ] - ∂n - (1 - 2ν)(r,i nj - r,j ni ) , Σjki = - 8π(1 - ν)r2 r = x-y и r,i = -1 (14) (1 - 2ν)(δij r,k + δik r,j - δjk r,i ) + 3r,i r,j r,k , ∂r = r-1 (yi - xi ). ∂xi Существование и единственность решений уравнения (13) с условиями (6) для задач упругости, в нашем случае на каждом шаге по времени, установлено в классе непрерывных по Гельдеру функций C 1,α (D) ⊂ W21 (D), где W21 (D) - пространство Соболева с параметром 0 < α 1 [7, 16]. Для нерегулярных областей, в том числе с поверхностями разрыва типа трещины, такие исследования выполнены, например, в [16] на основе теории псевдодифференциальных операторов и метода Винера-Хопфа. Следует отметить, что уравнение (13) с выражениями для ядер (14) применимо только для ограниченных сред, поскольку условия излучения Зоммерфельда, как известно, в этом случае не выполняются. Для рассматриваемых неоднородных сред уравнение (13) записывается отдельно для каждой из подобластей, входящих в область D и имеющих общие границы. В целом для области соответствующие уравнения объединяются с помощью условий (10), (11) на этих границах. Выражение для напряжений σjk (xi , t) может быть получено дифференцированием (13) с учетом (2)-(4) и для точек внутри области D имеет вид uijk (x, y)pk (y, t)ds - σij (x, t) = ∂D pkij (x, y)uk (y, t)ds+ ∂D uijk (x, y)ρ u ¨k (y, t) + γ u˙ k (y, t) - Xk (y, t) dV + + V Σijkl (x, y)ε0kl (y, t)dV - + Vp 4µ(1 + ν) 7 - 5ν ˙ (15) 2µεpij - δij αθ, 15(1 - ν) 3(1 - ν) где ukij = -Σkij = λδij ukl,l , pkij = λδij pkl,l . Тензор Σijkl определяется соотношением Σijkl = 2µ 3(1 - 2ν)(δij r,k r,l + δkl r,i r,j )+ + 3ν(δli r,j r,k + δjk r,l r,i + δik r,l r,j + δjl r,i r,k ) - 15r,i r,j r,k r,l + + (1 - 2ν)(δik δlj + δjk δli ) + (1 - 4ν)δij δkl (8π(1 - ν)r3 )-1 . Для вычисления напряжений на границе ∂D области D выражение (15) неприменимо из-за сингулярного характера входящих в него интегралов. Однако они могут быть получены в точках xi ∈ ∂D через граничные усилия и производные от граничных перемещений по касательным направлениям [14]. Для 145 П е т у ш к о в А. В. этого необходима параметризация геометрии границы ∂D и аппроксимация решения на ее поверхности. Аппроксимация решения также необходима и внутри области D для определения последних двух интегралов по внутренним объемам V , Vp в выражениях (13), (15) и параметров механики разрушения в окрестности трещины Ωc (ˆ xα ). Полагая рассматриваемые поверхности области D кусочно-гладкими, численное решение уравнения (13) выполним с использованием M граничных элементов ∆k так, чтобы M h ∂D ≈ ∂D = ∆k и ∆k ∩ ∆l = ∅, k = l. k=1 Параметрическое представление геометрии и искомых функций u˙ i (xj , t) и p˙i (xj , t) на каждом из них выберем подобно МКЭ в следующем виде: xhi :=Nr (ζγ )xri ; u˙ hi :=Nr (ζγ )u˙ ri , p˙hi := Nr (ζγ )p˙ri , ζγ ∈ ∆ k ⊂ R 2 . (16) Здесь Nr (ζγ ), r = 1, 2, . . . , 8, - квадратичная функция формы локальных координат ζγ ∈ (-1, 1), принадлежащая к лагранжеву семейству конечных элементов Shk,m (∂Dh ) в смысле [17]. Аналогичное параметрическое представление используется и для вычисления объемных интегралов внутри области D. В этом случае функция формы Nr в (16) включает в себя третью локальную координату η ∈ (-1, 1), нормальную к двум первым ζγ , и принадлежит тому же семейству: L Nr (ζγ , η) ∈ Shk,m (D), D ≈ Dh = ∆k , r = 1, 2, . . . , 8, . . . 20, k=1 для объемных элементов с 8 и/или 20 узлами. Следует отметить, что размер области с пластическими деформациями заранее не известен и устанавливается в процессе решения. Заменяя далее интегралы по границе и объему соответствующими суммами по граничным и объемным элементам, вместо (13), (15) получаем их численные аналоги: ¨ (t) = BX(t) + GP (t) + Dε0 (t), HU (t) + C U˙ (t) + M U ¨ (t), σ(t) = G P (t) - H U (t) + B X(t) + D ε0 (t) + Qε0 (t) - C U˙ (t) - M U (17) где M и C - соответственно матрицы масс и вязкого демпфирования, остальные матрицы и векторы являются стандартными в МГИУ. Вычисление входящих в эти выражения интегралов основано на численных квадратурах Гаусса. Объем и точность вычислений при этом зависят от порядка формулы и минимального расстояния между точкой наблюдения и узлом рассматриваемого элемента ∆k . Если точка наблюдения совпадает с одним из узлов граничного или объемного элементов, т. е. r → 0, интегралы в (13), (15) становятся сингулярными. Для их вычисления, как и в задачах 146 Изучение переходных процессов в нелинейно деформируемых средах. . . статики, используются переход соответственно к полярным и сферическим координатам и квадратуры Гаусса по радиальным и угловым координатам. Сильно сингулярные интегралы с особенностью O(r-3 ) рассматриваются в смысле главного значения Коши. Их вычисление по объемным элементам строится путем перехода к сферическим координатам и выделения локальной зоны вблизи сингулярной точки, по которой проводится расчет аналитически по радиусу и численно с использованием квадратуры Гаусса по угловым координатам [14]. Прямое численное интегрирование уравнений (17) во времени может быть осуществлено с помощью стандартных разностных схем; среди них выберем неявную двухслойную разностную β-схему Ньюмарка, обычно применяемую в МКЭ для решения задач нелинейной динамики. Предпочтение такому выбору обычно отдают из-за наличия у схемы контролируемых характеристик устойчивости и возможности простого подключения различных итерационных схем для учета нелинейностей. Согласно этой схеме на Dt вводится равномерная разностная сетка с шагом ∆t таким образом, что Dt = tn : tn = n∆t; n ∈ Z, 0 n∆t τ . На каждом интервале времени [tn , tn+1 ] ускорение, входящее в (17), обычно полагается постоянным и равным среднему значению (осредненному ускорению): ¨ (t) = (1 - γ)U ¨n + γ U ¨n+1 , U где γ - свободный параметр, принимаемый здесь равным 1/2. Ускорения и скорости на временном слое (n + 1) могут быть вычислены через перемещения U (t) следующим образом: ¨n+1 = 1 Un+1 - 1 Un + U˙ n ∆t + 1 - β ∆t2 U ¨n , U β∆t2 β∆t2 2 (18) 1¨ 1 1 1 2¨ ˙ Un+1 + U˙ n + U ∆t - U + U ∆t + - β ∆t U , U˙ n+1 = n n n n 2β∆t 2 2β∆t 2 где β - второй свободный параметр, управляющий точностью и устойчивостью схемы Ньюмарка. Для безусловной устойчивости схемы необходимо, вообще говоря, выполнение условия 2β γ 1/2. Выражения (18) включают в себя текущие (актуальные) значения вектора перемещений. Тем самым учитывается нелинейный, зависящий от решения характер изменения инерционных и вязких сил в процессе динамического деформирования среды. Подставим эти выражения вместе с граничными условиями в уравнения (17) и сведем в левую часть все входящие в них неизвестные. Полученную разрешающую систему нелинейных уравнений МГИУ можно переписать в следующем виде: [A]X˙ n+1 = F˙n+1 + [D] ε˙0 σ˙ = -[A ]X˙ n+1 + F˙n+1 + [D ] ε˙0 n+1 n+1 , + [G]ε˙0n+1 . (19) (20) 147 П е т у ш к о в А. В. Таким образом, на каждом шаге по времени ∆t требуется решение обычной нелинейной задачи квазистатики, детали которой подробно обсуждаются в [14, 15]. Входящие в (19), (20) [A] и [A ] - матрицы коэффициентов при векторе граничных неизвестных, X˙ n+1 - вектор, содержащий неизвестные усилия и перемещения на границе, векторы F˙n+1 и F˙n+1 содержат члены, определяемые отличными от нуля заданными значениями граничных перемещений и усилий, а также действием объемных сил, включая инерционные и вязкие. Матрицы [D] и [D ] включают в себя объемные интегралы с начальными и нелинейными деформациями в (13), (15), квазидиагональная матрица [G] определяется наличием свободных от интеграла членов в (15). Выражения для перечисленных матриц и векторов стандартные в МГИУ и здесь не приводятся. Их можно найти, например, в [2]. Для неоднородных (составных) сред уравнения (19), (20) записываются отдельно для каждой из сред и на общей поверхности между ними Γ = Γ1 ∪Γ2 вводится два слоя узлов и элементов. В случае идеального контакта на первом из них в качестве граничных условий задаются неизвестные перемещения и нулевые усилия, на втором - неизвестные усилия и нулевые перемещения. В случае проскальзывания с трением на первом слое задаются неизвестные компоненты перемещений - нормальные и касательные, а на втором - неизвестные нормальные усилия и касательные перемещения. При этом условия (10), (11) учитываются при формировании матриц [A] и [A ] . Моделирование нелинейного деформирования сред осуществляется здесь, как и в [13, 14], с использованием на каждом шаге по времени уравнений состояния (2), (3) и классического метода последовательной линеаризации в форме начальных деформаций с применением известной схемы «упругий предиктор - пластический корректор». Приращения пластической деформации на каждой итерации определяются из решения (19), (20) с использованием соотношений (3) и экспериментально полученных диаграмм деформирования для материала каждой из рассматриваемых сред. Большая размерность, плотная упакованность и несимметричная структура матриц и векторов в уравнениях (19), (20) делают вычислительные операции над ними чрезвычайно трудоемкими и представляют серьезное препятствие для применения МГИУ в математическом моделировании переходных процессов в реальных крупномасштабных объектах. В связи с этим для ускорения формирования и решения систем уравнений МГИУ в последние годы разработаны так называемые быстрые методы, такие как FMM, АСА и другие [6], а также используемый ниже метод дискретных областей [14, 18]. В основе этого метода лежит представление исходной расчетной области D в виде совокупности подобластей с искусственными границами между ними и граничными условиями типа (10). В результате матрицы уравнений становятся блочно-симметричными с ленточной структурой, что позволяет эффективно использовать прямые алгоритмы гауссовского типа для их обращения [19]. 3. Вычисление напряжений на границе тела, метод дискретных областей. Интегральное выражение (14) неприменимо для вычисления напряжений на границе нелинейно деформируемого тела из-за сингулярного характера входящих в него интегралов с особенностью O(r-2 , r-3 ). Вместо него 148 Изучение переходных процессов в нелинейно деформируемых средах. . . воспользуемся вычисленными на границе усилиями и касательными производными от граничных перемещений и соотношениями (7), (8). Перемещения uj (xi , t) и их производные по локальным координатам ζγ с учетом аппроксимации (16), а также граничные усилия pj (xi , t) определяются из решения уравнений (19), (20) для каждого шага по времени. Через них с использованием (7) могут быть вычислены напряжения на границе ∂D рассматриваемой области D. Такой прием оказывается весьма эффективным с вычислительной точки зрения, поскольку обладает высокой точностью и не требует дополнительного вычисления производных от интегралов по области пластического деформирования, нет также необходимости вычисления матрицы D в выражении (20) в процессе решения нелинейной задачи. Его легко обобщить на нелинейно деформируемые среды, расчетные области которых состоят из элементов с сильно различающимися свойствами и/или размерами или включают в себя локальные особенности, например, типа произвольно ориентированных трещин. В этом случае область D, занимаемая деформируемой средой, рассматривается в виде совокупности N дискретных подобластей Ωi N D= Ωi и Ωl ∩ Ωm = ∅, l = m, i=1 образованных произвольно ориентированными фиктивными границами (см. рис. 1). На каждой из них задаются условия типа (10) или (11) и при необходимости учитываются, например, взаимодействия берегов трещины [20]. Кусочно-неоднородные среды с различными свойствами материалов и реальными разделяющими их границами учитываются при этом естественным образом. При таком подходе появляется возможность объединения предлагаемого метода с другими, например, МКЭ для областей с сильной нелинейностью, или МГИУ в его прямой формулировке (12) с областью, включающей «бесконечность», для учета граничных условий Зоммерфельда [11], или даже с экспериментальным методом определения переменных физического поля на поверхности. Изложенные в статье методы и разработанные на их основе алгоритмы реализованы в компьютерной программе, как развитие программы [11] разработанной ранее применительно к изучению квазистатических процессов нелинейного деформирования трехмерных сред - конструкций и их элементов. В качестве пре- и постпроцессора в пакет включена программа GID [21] с широкими возможностями автоматического построения расчетных схем МГИУ и визуализации результатов. 4. Численные результаты. Для обоснования вычислительной эффективности предложенного подхода и компьютерной программы [10], развитой применительно к моделированию нелинейной волновой динамики, вначале приведем результаты решения МГИУ классической задачи о распространении упругой и упругопластической волны вдоль неоднородного стержня при его мгновенном нагружении в момент времени t = 0. Задача рассматривается в объемной постановке (рис. 2), где показаны геометрия, условия закрепления и функция нагружения стержня. Длина стерж149 П е т у ш к о в А. В. ня L = 10 см, поперечное сечение - квадрат с длиной стороны l = 0.2L. Стержень состоит из двух подобластей Ω1 и Ω2 с общей границей B, каждая со своим материалом, диаграммы деформирования которых приведены на рис. 2, b. Механические свойства материалов: модули Юнга E1 = 2.15 · 105 МПа, E2 = 2.05 · 105 МПа; коэффициенты Пуассона ν1 = 0.3, ν2 = 0.33, начальные плотности материалов ρ1 = ρ2 = 7.85 · 103 кг/м3 . Функция нагружения - равномерно распределенное на торце давление p(t) = p0 H(t), где H(t) - функция Хевисайда, значение p0 = 220 МПа, рис. 2, c. В случае упругого поведения материала задача имеет аналитическое решение. Для нелинейного деформирования воспользуемся сравнительными результатами, полученными МКЭ с использованием программы [22]. Для аппроксимации решения на поверхности стержня выбраны квадратичные восьмиузловые граничные элементы, при этом объемные интегралы вычислялись с использованием аналогичных двадцатиузловых элементов. Изучались вопросы точности и вычислительной устойчивости МГИУ при последовательном увеличении числа граничных элементов для каждой подобласти (материала) от 88 до 1426 и выборе соответствующих значений шагов по времени с использованием известного критерия КФЛ - c1 ∆t/h < 1/3, где h - размер равномерной сетки, и β, выбираемого из диапазона β = 0.25 ÷ 0.10, в соотношениях (18). Результаты, представленные ниже, получены с использованием в аппроксимации стержня 352 квадратичных элементов, β = 1/6 и шагом по времени 1.5 · 10-4 с. Рис. 2. Ударное нагружение неоднородного стержня: a - геометрия, b - диаграммы нелинейного деформирования материалов, c - функция нагружения [Figure 2. Impact loading of the inhomogeneous rod: a-the geometry, b-the diagrams of nonlinear deformation of the materials, c-the loading function] 150 Изучение переходных процессов в нелинейно деформируемых средах. . . Распространение упругопластической волны напряжения вдоль неоднородного стержня, полученное с использованием уравнения состояния (3) и диаграмм деформирования составляющих его материалов (рис. 2, b), показано на рис. 3. Здесь же приведены результаты упругого решения в сечениях A и C стержня и для сравнения результаты, полученные МКЭ. Для демонстрации объемности возникающих в стержне напряженных состояний в точке C приведены также другие значимые компоненты тензора напряжений. Из сравнения полученных результатов следует хорошее качество получаемых МГИУ решений. Деформированная геометрия стержня и распределение в нем упругопластического напряжения σ33 для момента времени t = 6.0 · 10-3 c от начала ударного нагружения показаны на рис. 4, a. Видно, что основная концентрация напряжения имеет место в зоне соединения разнородных материалов (сечение B на рис. 2, a). Здесь, по-видимому, и происходит накопление повреждений и разрушение стержня. Практический интерес представляет и характер взаимодействия распространяющих волн напряжения в стержне с границей этого соединения (рис. 4, b), где также приведено аналогичное решение, полученное МКЭ. Сравнение результатов и в этом случае свидетельствует о хорошей точности и разрешающей способности МГИУ в описании больших градиентов изменения полевых функций. С использованием результатов вычислительного эксперимента, выполненного на объемном стержне, МГИУ решена более сложная задача о кинетике напряженно-деформированных и возникновении предельных состояний в зонах резкого локального изменения геометрии в трехмерных конструкциях при ударных воздействиях. В объемной постановке выполнено моделирование переходных нелинейных процессов деформирования в биметаллическом образце с острым надрезом под действием внезапно приложенной к одному из торцов ударной нагрузки (рис. 5, a). Надрез расположен вблизи зоны соединения разнородных материалов, которые сильно отличаются друг от друга своими деформационными свойствами (рис. 5, b). Условия нагружения представлены на рис. 5, c. Уровни нагружения и механические свойства материалов аналогичны приведенным в предыдущем примере. На этом же рисунке показаны сечения образца, представляющие интерес для анализа указанных состояний, и распределение зон пластических деформаций в центральном сечении ABCD образца для момента времени t = 8.83 · 10-4 c от начала ударного нагружения. В аппроксимации расчетной области полосы использовано 784 объемных и 398 граничных элементов. На рис. 6 представлена информация о максимальном раскрытии берегов надреза (по существу трещины) в точке N образца (см. рис. 5, d) для упругого и упругопластического поведения материалов. На рис. 7, a приведено распределение упругопластического напряжения в самом образце с сечением вблизи вершины надреза для момента времени t = 8.83 · 10-4 c и изменение этого напряжения вдоль траектории S, построенное для различных моментов времени от начала нагружения (рис. 7, b). Распределение по форме характеризует собой известный режим «с обострением», возникающий в зоне разнородного соединения материалов с концентратором напряжений, и результаты его очень важны для оценки прочности 151 П е т у ш к о в А. В. Рис. 3. Распространение волн напряжений в точках A и C стержня: 1 - упругое решение МГИУ для σ33 , 1 - пластическое решение МГИУ для σ33 ; 2 - упругое решение МКЭ для σ33 , 2 - пластическое решение МКЭ для σ33 ; 3 - решение МГИУ для σ11 , 3 - решение МГИУ для σ22 (онлайн в цвете) [Figure 3. Stress waves distribution in the points A and C of the rod: the line 1 is an elastic BIEM solution for the stress σ33 , the line 1 is a plastic BIEM solution for the stress σ33 ; the line 2 is an elastic FEM solution for the stress σ33 , the line 2 is a plastic FEM solution for the stress σ33 ; the line 3 is a BIEM solutions for the stress σ11 , the line 3 is a BIEM solutions for the stress σ22 (color online)] Рис. 4. Распределение упругопластических напряжений σ33 (МПа) в стержне: a - геометрия и напряжения в момент времени t = 6 · 10-3 c от начала нагружения, b - изменение напряжения во времени в точке B по МГИУ (линия 1) и МКЭ (линия 1) (онлайн в цвете) [Figure 4. Distribution of the elastoplastic stresses σ33 (MPa) in the rod: a-the geometry and the stresses at time t = 6.0 · 10-3 s from the beginning of loading; b-stress evolution in time at point B by BIEM (line 1), and by FEM (line 2) (color online)] 152 Изучение переходных процессов в нелинейно деформируемых средах. . . Рис. 5. Образец с надрезом при мгновенном нагружении: a - геометрия и условия закрепления, b - функция нагружения, c - диаграммы деформирования материалов, d - зоны пластических деформаций в центральном сечении для момента времени t = 8.83 · 10-4 c (онлайн в цвете) [Figure 5. A bi-material plate with edge notch/crack along interface at instantaneous loading: a-the geometry and the fixing conditions, b-the loading function, c-the diagrams of deformation of the materials, d-the plastic deformation zones in the central section at time t = 8.83 · 10-4 s (color online)] Рис. 6. Максимальное раскрытие берегов надреза uN 3 во времени: пластическое (линия 1) и упругое (линия 2) решения по МГИУ [Figure 6. Maximum opening of the bi-material sample edge notch uN 3 in time: the plastic (line 1) and elastic (line 2) solutions by BIEM] 153 П е т у ш к о в А. В. Рис. 7. Максимальное напряжение σ33 в вершине надреза образца: a - распределение для момента времени t = 8.83 · 10-4 c, b - кинетика распределения напряжений для различных моментов времени t вдоль траектории S, c - кинетика распределения напряжений для различных моментов времени t вдоль сечения S1 ; 1 - t = 3.65 · 10-4 c, 2 - t = 7.28 · 10-4 c, 3 - t = 8.83 · 10-4 c, 4 - t = 1.46 · 10-3 c, 5 - t = 2.18 · 10-3 c (онлайн в цвете) [Figure 7. The maximum stress σ33 in the top of the edge notch of the bi-material specimen: a-distribution at time t = 8.83 · 10-4 s, b-stress distribution kinetics for various times t along the trajectory S, c-stress distribution kinetics for various times t along the cross-dection S1 ; line 1 for t = 3.65 · 10-4 s, line 2 for t = 7.28 · 10-4 s, line 3 for t = 8.83 · 10-4 s, line 4 for t = 1.46 · 10-3 s, line 5 for t = 2.18 · 10-3 s (color online)] 154 Изучение переходных процессов в нелинейно деформируемых средах. . . самого соединения. Не меньший интерес представляет распределение напряжений и в другом сечении образца на границе Γ соединения разнородных материалов. Изменения напряжений во времени вдоль выбранной траектории S1 (рис. 7, c) также свидетельствуют о возможности разрушения образца в этом сечении. Заключение. Для изучения переходных процессов в объемных нелинейно деформируемых средах при динамических, в том числе ударных, воздействиях разработаны методы моделирования, основанные на интегральных представлениях краевых задач, численных схемах высокого порядка аппроксимации границ и коллокационного приближения. Представлены соответствующие формулировки МГИУ с фундаментальными решениями статической упругости и разностными методами по времени, позволяющие учитывать неоднородности и локальные особенности сред. Также обсуждаются основные аспекты численной реализации выполненных разработок. Изучены вопросы точности, вычислительной устойчивости и эффективности метода и компьютерной программы. Приведены примеры моделирования линейных и нелинейных процессов распространения волн деформаций и напряжений в трехмерных средах и дано их сравнение с численными решениями, полученными МКЭ. Они подтверждают вычислительную эффективность разработанных алгоритмов, а также общность и полезность для практических целей предлагаемого подхода. Конкурирующие интересы. У меня нет конкурирующих интересов. Авторская ответственность. Я несу полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи мною одобрена. Финансирование. Исследование не имело финансирования.

About the authors

Vladimir A Petushkov

A. A. Blagonravov Mechanical Engineering Institute RAS

Email: pva_imash@bk.ru
Dr. Phys. & Math. Sci.; Professor; Lab. of Mathematical Simulation 4, M. Khariton’evskii per., Moscow, 101990, Russian Federation

References

Supplementary files