理解几何非线性分析 - 从晾衣服开始（续）_Deform_非线性_通用_python_UM

理解几何非线性分析 - 从晾衣服开始（续）

yunduan082

2月前浏览2023

很多朋友在看了这篇文章《理解几何非线性分析 - 从晾衣服开始？》后，比较好奇手算该怎么计算出来？今天我们就来送上答案。

让我们考虑几何非线性和大变形条件，来推导这个问题。我们将一步一步地进行分析：

Step1：问题条件

已知条件：

钢丝直径 d = 5mm = 0.005m
钢丝长度 L = 5m
弹性模量 E = 200GPa = 200 × 10^9 Pa
重物质量 m = 5kg
重力加速度 g ≈ 9.8 m/s²

Step2：考虑几何非线性

在大变形条件下，我们不能使用小角度近似。我们需要考虑钢丝的实际长度变化。

设：

y 为中点垂直位移
T 为钢丝张力
θ 为钢丝与水平面的夹角

Step3：力平衡方程

垂直方向的力平衡：
2T sin θ = mg

Step4：几何关系

在大变形条件下：
tan θ = y / (L/2)
sin θ = y / sqrt((L/2)^2 + y^2)
cos θ = (L/2) / sqrt((L/2)^2 + y^2)

Step5：应变计算

真实应变 ε = ln(L' / L)，其中 L' 是变形后的长度

L' = 2 * sqrt((L/2)^2 + y^2)

ε = ln(2 * sqrt((L/2)^2 + y^2) / L)

Step6：应力-应变关系

对于大应变，我们可以使用非线性应力-应变关系，如 Ramberg-Osgood 模型。但为简化计算，我们仍使用线性关系：

σ = Eε

Step7：张力计算

T = σA = EεA = EA * ln(2 * sqrt((L/2)^2 + y^2) / L)

Step8：建立方程

将张力表达式代入力平衡方程：

2 * [EA * ln(2 * sqrt((L/2)^2 + y^2) / L)] * [y / sqrt((L/2)^2 + y^2)]= mg

Step9：求解方程

这个方程是非线性的，无法直接求解。我们需要使用数值方法，如牛顿-拉普森法。

设函数 f(y) = 2 * [EA * ln(2 * sqrt((L/2)^2 + y^2) / L)] * [y / sqrt((L/2)^2+ y^2)] - mg

我们需要找到 f(y) = 0 的解。

Step10：数值求解

使用 Python 的 scipy 库进行数值求解：

importnumpyasnp
fromscipy.optimizeimportfsolve

E =200e9# 弹性模量
A = np.pi * (0.005**2) /4# 横截面积
L =5# 长度
m =5# 质量
g =9.8# 重力加速度

deff(y):
return2* (E * A * np.log(2* np.sqrt((L/2)**2+ y**2) / L)) * (y / np.sqrt((L/2)**2+ y**2)) - m * g

y_solution = fsolve(f,0.05)# 初始猜测值为0.05m
print(f"垂直位移 y = {y_solution[0]:.6f} m")

# 计算张力
T = E * A * np.log(2* np.sqrt((L/2)**2+ y_solution[0]**2) / L)
print(f"张力 T = {T[0]:.2f} N")

运行这段代码，得到结果：
垂直位移 y ≈ 0.057996 m = 58 mm
张力 T ≈ 1056.39 N

这些结果与仿真结果（变形57mm，张力1037N）非常接近。

结论：
通过考虑几何非线性和大变形条件，我们得到了更准确的结果，与仿真结果基本一致。这表明在处理大变形问题时，考虑非线性效应是至关重要的。

关于变形、轴向力与绳子直径的关系

我们推导出了计算方程，也便知道了变形、轴向力与绳子直径的关系，可以利用Python脚本绘制曲线。脚本如下：

importnumpyasnp
fromscipy.optimizeimportfsolve
importmatplotlib.pyplotasplt

# 常量定义
E =200e9# 弹性模量
L =5# 长度
m =5# 质量
g =9.8# 重力加速度

defcalculate_y_and_T(d):
A = np.pi * (d**2) /4# 横截面积

deff(y):
return2* (E * A * np.log(2* np.sqrt((L/2)**2+ y**2) / L)) * (y / np.sqrt((L/2)**2+ y**2)) - m * g

y_solution = fsolve(f,0.05)[0]# 初始猜测值为0.05m
T = E * A * np.log(2* np.sqrt((L/2)**2+ y_solution**2) / L)

returny_solution, T

# 生成直径范围
diameters = np.linspace(0.001,0.03,300)# 1mm to 10mm

# 计算每个直径对应的变形和轴向力
deformations = []
tensions = []

fordindiameters:
y, T = calculate_y_and_T(d)
deformations.append(y)
tensions.append(T)

# 绘图
plt.figure(figsize=(12,5))

# 变形与直径的关系曲线
plt.subplot(1,2,1)
plt.plot(diameters *1000, np.array(deformations) *1000)
plt.xlabel('直径 (mm)')
plt.ylabel('变形 (mm)')
plt.title('变形与直径的关系')
plt.grid(True)

# 轴向力与直径的关系曲线
plt.subplot(1,2,2)
plt.plot(diameters *1000, tensions)
plt.xlabel('直径 (mm)')
plt.ylabel('轴向力 (N)')
plt.title('轴向力与直径的关系')
plt.grid(True)

plt.tight_layout()
plt.show()

# 输出5mm直径时的结果
d_5mm =0.005
y_5mm, T_5mm = calculate_y_and_T(d_5mm)
print(f"当直径为5mm时：")
print(f"变形 ={y_5mm*1000:.2f}mm")
print(f"轴向力 ={T_5mm:.2f}N")

这段代码会生成以上两个图表：一个显示变形与直径的关系，另一个显示轴向力与直径的关系。变形的趋势是符合的，再看轴向力。。。。。问题来了：为什么轴向力与绳子直径的关系是单调增加的呢？

这是因为：上述方法通过几何关系来计算L‘即绳子伸长量的，对于刚性相对较小的情况是相对准确的，如细绳为棉线，则在载荷作用下，绳子变形后形状为三角形。而实际上，随着绳子直径增加，刚性增强，绳子的变形不再是直线形式，而是曲线，自然伸长量也不能用三角关系式去计算。因此这种方法得到的轴向力与直径的关系是不对的，或者说只在一定范围内精度可接受。至于通用的解析方法，是得花点功夫去推导的。

以下是通过Ansys将绳子半径R参数化，获得的变形和轴向力随直径增加的变化趋势：因为是钢绳，随着直径增加，刚性大幅增加，变形急剧降低，而轴向力在直径14mm时达到最大值，随后直径增加，轴向力减小。

来源：ABAQUS仿真世界

如何正确理解非线性分析结果

正如我之前所写，理解分析结果是一项艰巨的任务，并且可能导致严重的问题。在模型中发现错误比理解结果解释问题要容易得多-毕竟，软件会给出有用的提示。而结果是需要专业人员的深刻理解。本文我将描述稳定性路径——正确估计模型承载能力的最重要工具。今天我们将介绍基础知识，在接下来的文章中，我将在实际示例中讨论此方法的用法。稳定路径基础知识为了正确理解，使用模型的承载稳定性路径。需要进行更复杂的分析，其中载荷以增量形式施加到模型上，并且可以观察每个增量中模型的载荷和位移之间的关系。在我们深入探讨之前先了解3个最重要的定义：稳定性路径是表示模型中载荷和位移之间关系的图表。人们可以以多种形式绘制稳定路径：应在垂直轴上标记载荷，在水平轴上，应该显示变形（也可以是哪个节点/点以及什么自由度上由作者决定）——有时使用归一化变形而不是“真实”变形。稳定性路径在文献中以多种名称出现，例如：平衡路径、载荷-位移图、位移响应路径等。稳定性路径示例表示为r的载荷乘子是一个参数，在分析中以任何给定增量乘以所施加的载荷。通常在下标中标记分析类型。为了更清楚地解释一下，我们假设在某处有一个负载10kN的模型。如果我们想以10个增量（每个增量1kN）施加该载荷，那么第一个增量将为1kN，因此载荷乘数将为0.1（因为10kN*0.1=1kN），而在第二个增量中，总载荷将为2kN，因此载荷乘数将为0.1（因为10kN*0.1=1kN）乘数将为0.2，依此类推。可以在垂直轴上使用载荷乘数来代替载荷或约束力。归一化变形（我通常表示为Δ）是当前增量中所选DOF（自由度）中的位移与线性分析中该DOF变形的比值。可以使用归一化变形来代替稳定性路径的水平轴上的变形。稳定路径：示例1–如何完成让我们从最简单的例子开始。我们将分析一根简支梁（两端固定），中间施加一个单一的力，如下所示。材料设置为弹性材料，我们将使用几何线性分析。换句话说，我们将执行完全线性分析，但具有增量加载。这种分析一般来说没有什么意义，但却是如何制定稳定路径的一个很好的例子。开始：乍一看，很明显中间的节点变形最大，对分析也最有意义，所以我们选择这个节点作为变形的参考节点。在垂直轴上，我们将简单地使用施加的载荷。本例中的总负载为500kN。在每个增量下，我将显示变形模型和稳定性路径图。因为我们没有开始，所以没有变形，也没有负载——所以图是“空的”。增量1：在分析设置中，我决定将负载分为10个相等的增量。这意味着增量1为500kN/10=50kN。对于该载荷，中间的位移为66.5mm，因此此时的稳定路径如下所示：增量2：第二个增量增加了额外的50kN，现在总载荷为100kN。变形为133毫米。目前，情况是这样的：增量3：额外添加50kN（总计150kN）。当然，由于一切都是线性的，挠度将比第一步高3倍，并且稳定路径将是一条直线。增量7：让我们加快一点速度。增量7表示总载荷为7x50kN=350kN。位移为66.5mx7=465.5mm。增量10–完成：当然，最后，对模型施加满500kN的力，变形为665mm。请注意，如果您使用线性求解器（或者您设置不希望在非线性求解器中获得中间结果），这将是您从分析中获得的唯一“绘图”。这就是为模型构建稳定性路径的方式。在下一个示例中，我将使用几何非线性模型来展示差异。稳定路径：示例2–非线性分析我们将使用相同的梁，但进行几何非线性计算。由于两个支撑都是固定的（没有滑动），我们预计挠度会减少。在此示例中，载荷分为100个增量（每个增量5kN），但我将仅向您展示其中的一些。另外，我在图表上留下了线性稳定性路径以供参考。增量10：载荷为50kN，位移为47mm，因此比线性分析(LA)中小了近20mm。我们可以看到，稳定性路径不再是线性的，并且与LA路径相比，它呈上升趋势。这告诉我们模型中存在某种硬化。增量40：载荷为200kN，位移为117。如果我将载荷仅分成5个增量，则稳定性曲线将由直线组成，但由于我使用100，我收到的曲线非常平滑。有时只需查看稳定性路径，您就可以判断结果是否正确：如果曲线“尖锐”，则意味着增量太少，这可能会严重影响结果（我将在以后讨论这一点）。增量100–完成：这里没有什么可讨论的，直接到达终点线！施加的载荷为500kN，最大变形为190mm（比线性分析小3倍以上！）。稳定路径：示例3–实际使用接下来是对更复杂的筒仓结构的分析。由于我希望准确，所以我将负载分为多个增量（准确地说是1000）。如果不制作所述载荷的图表（在筒仓中施加的载荷相当复杂），就不可能写出每个增量施加了多少载荷，因此在本例的稳定性路径中，我在其中一个支撑件中使用了总反作用力。增量40：这几乎是分析的开始。由于该载荷工况是轴对称的，因此仅使用对称边界条件分析壳的一部分（这种方式计算速度要快得多）。下面的结果显示了整个模型（支撑区域的特写）以及分析每个步骤的稳定性路径。增量60：模型中出现第一个屈服点。由于温度对材料的影响，屈服应力为231.8MPa；增量70：负载越来越接近模型的最大承载能力；增量80：稳定性失效正在发生，局部发生塑性；增量590：分析结束。我设置的弧长方法中的变量告诉求解器现在应该退出计算。由于我设置的负载从未达到（模型在施加负载的60%左右承载下降），如果没有额外的“限制点”，求解器将无休止地工作！还值得注意的是，在这个模型中，我实现的负载是我想要看到的故障机制的两倍高，因此60%的结果实际上是一个积极的结果（这就是负载乘数派上用场的地方））完成上述计算后，很容易看出，当反作用力等于5300kN时就达到了容量，并且模型中的失效后状态不稳定（稳定路径正在减小）。还很容易看出失效机制是由于塑性不稳定性造成的。这种方法有很多用途，我稍后将讨论许多其他方面。阅读本文后您应该了解什么：什么是稳定路径如果有足够的数据可用，如何制定稳定路径什么是负载乘数最后是示例3中筒仓的图片：来源：ABAQUS仿真世界