Abaqus-之线性粘弹性 –2 求解方程
介绍
在本系列的第 1 部分中，推导出了如何使用以下积分方程计算任何应变历史的应力响应：
正如在第 1 部分中展示的那样，这个方程可以在像 julia 这样的高级数学语言中使用像 quadgk 这样的函数来解决。这种方法的问题在于非常慢。为了说明有多慢，设N是计算σ（Δt）所需的浮点运算（FLOP）数，然后计算应力σ（M⋅Δt）需要MN浮点运算。更重要的是，需要 O（N⋅M2） 浮点运算来计算由 M 个点组成的整个应力-应变曲线。因此，如果应力-应变曲线有 1，000 个数据点，那么计算整个曲线所需的时间大约是第一个数据点的 100 万倍。这很糟糕！我们需要一个更好的方法。在本文中，我将展示如何将此计算速度提高 100 倍到 1000 倍。
单项Prony系列
在计算应力响应之前，我们需要定义如何量化应力松弛响应。在这里，我将使用以下 1 项 Prony 级数方程：
在此方程中，g=0 对应于线弹性响应，g=1 对应于完全应力松弛。
分部积分
从概念上讲，使用以下众所周知的方程按分部积分进行积分是一个好主意：
这得到了以下应力方程：
这个方程很酷的地方在于，我们可以看到应力由瞬时部分（σi）和粘弹性部分（σv）组成。但是，这个等式并没有帮助我们解决运行时问题。为此，我们需要一种不同的方法。
将积分分成两部分
现在是时候应用一些数学“魔法”了。首先重写σ（t+Δt）的等式（4）。请注意，我将积分分为两部分：
如果你考虑了一段时间，你应该能够说服自己，第一个积分可以简化很多：
为了计算最终积分，我们需要假设应变在时间间隔 [t，t+Δt] 内随时间线性变化：
有了这个假设，就可以计算方程（7）中的积分。经过 2 页代数推导，能够得到以下应力的最终方程：
计算应力
下面的 Julia 代码计算了给定应变历史的应力，图 2 显示了运行代码的结果。此实现的运行速度明显快于原始公式 （1）。
预测的应力
using QuadGK, Plotsfunction strain(t)t1 = maxStrain / strainRateif t <= t1return strainRate * telsereturn maxStrain - (t - t1) * strainRateendendfunction calcStress_direct(time, strain)N = length(time)stress = zeros(N)stressV = zeros(N)for i in 2:Ndt = time[i] - time[i-1]de = strain[i] - strain[i-1]stressE = E0 * strain[i]stressV1 = stressV[i-1] * exp(-dt/tau0)stressV2 = E0*g * (strain[i-1] - de/dt * time[i]) * (1 - exp(-dt/tau0))stressV3 = E0*g * de/dt * (time[i]+dt-tau0 + (tau0-time[i])*exp(-dt/tau0))stressV[i] = stressV1 + stressV2 + stressV3stress[i] = stressE - stressV[i]endreturn stressendE0 = 10.0tau0 = 4.0g = 1.0maxStrain = 0.1strainRate = 0.01N = 100time = range(0, 2*maxStrain/strainRate, length=N)ee = zeros(N)for i in 1:Nee[i] = strain(time[i])endss3 = calcStress_direct(time, ee)plot!(ee, ss3, label="direct", linewidth=2)plot!(size=(1000,800), framestyle=:box)plot!(tickfontsize=14,guidefontsize=16,xlabel="Strain", ylabel="Stress")
