弹性材料曲线问题

倪举鹏

前面的帖子发得不严谨，可以删了。大家看看这个问题

倪举鹏

我突然发现了，图中的F直接这样将头尾拉到一起，形成的形状就是后面8字曲线的上半段或下半段，大家再也不用产生歧义了，计算吧

倪举鹏

还有个问题  F在这么拉到头尾相连的过程中 ， F是一直在减小吗？如果不是，拉到什么形状的时候F最大

yinhow

yinhow

倪举鹏 发表于 2014-10-22 20:28
我突然发现了，图中的F直接这样将头尾拉到一起，形成的形状就是后面8字曲线的上半段或下半段，大家再也不用 ...

弹性势能极小，变分法求的EL方程

yinhow

yinhow 发表于 2014-10-22 22:07
弹性势能极小，变分法求的EL方程

你先要知道弹性势能的表达式，300多年前的Euler也想过。

倪举鹏

eq := x*(1+(D(y))(x)^2)^(3/2) = .18*(diff((D(y))(x), x));
print(outputredirected...); # input placeholder
(3/2)
/ 2\
x \1 + D(y)(x) / = 0.18 @@(D, 2)(y)(x)
s := dsolve({eq, y(0) = 0, (D(y))(0) = -.86}, y(x), numeric);   很接近你的曲线  感觉这个水滴形状是固定的，不会随着弹性系数变化而变化，尖凸地方切线夹角大约81度

葡萄糖

8字曲线为贝努利双扭线！

倪举鹏 发表于 2014-10-22 19:51
前面的帖子发得不严谨，可以删了。大家看看这个问题 ...

为何不严谨？何必删去？

倪举鹏 发表于 2014-10-22 20:28
我突然发现了，图中的F直接这样将头尾拉到一起，形成的形状就是后面8字曲线的上半段或下半段，大家再也不用 ...

http://bbs.emath.ac.cn/thread-5892-1-1.html
http://www.eecs.berkeley.edu/Pubs/TechRpts/2008/EECS-2008-103.pdf

kastin

让楼主去看材料力学和弹性力学的教材他不愿意，单靠网上自学也就能学个半调子。

没学过就不知道，不知道就会想当然，想当然就会出现让人啼笑皆非的错误。

任何均匀线弹性材料（这里强调“均匀”，是为了让弹性模量是各项同性的标量，否则它就是一个二阶张量；“线弹性”也是为了利用广义胡克定律），产生弯曲是因为弯矩（注意，应力都是内力）的存在，而材料之所以有不同程度抵抗弯曲变形的能力，是因为它们有弹性大小和几何尺寸。这个能力我们在力学中用“刚度”来称呼，其大小取决于材料的截面惯性矩 `I` 和材料弹性模量 `E`，大小等于二者之积。

具体来说，设有线弹性材料，受到使其产生纯弯曲（没有扭转与剪切变形）的外力（包括力偶矩）作用，内部产生内力（弯矩）。根据受力分析可得任意位置截面弯矩函数 `M(x)`，记该截面处的挠度为 `w(x)`（也就是材料弯曲产生偏离水平切线的位移），该段弧微元的曲率半径为 `\rho`。根据广义胡克定律以及一些几何性质可推出$$\frac{1}{\rho}=\frac{M}{EI}$$上面这个式子说明，对于给定线弹性材料，其弯曲的曲率半与该截面处的弯矩成反比。

根据数学中曲率公式$$\kappa=\frac{1}{\rho}=\frac{|w''|}{(1+w'^2)^{3/2}}$$可导出变形与受力的关系$$\frac{|w''|}{(1+w'^2)^{3/2}}=\frac{M}{EI}$$若是小变形，那么上式左端分母部分近似有 `w'^2<<1`，于是弯曲变形的解法是：通过受力分析给出沿着材料长度方向的弯矩，根据上面二阶微分方程，求出`w(x)`即。若不是小变形，则需要数值求解微分方程。

首先，楼主不给出刚度，根本没法求。并且前面楼主说让材料无穷细，那么截面惯性矩 `I\to 0`，即刚度 `EI \to 0`，此时材料跟水流有啥区别？那就不可能有弹性了，如同一根软绳，可以摆成随意形状，谈何弯曲？

其次，能量法也需要用到弹性模量，截面惯性矩。能量法有好两种，最小势能原理和最小余能原理。这是材料力学和弹性力学中经常用的解法，其本质与通常非能量解法并不有什么区别，只是侧重角度不同。但注意这两种方法的结论只是适用于线弹性的保守系统。因为 1） 保守系统中，势能对位移的偏导数就等于应力，因为应力都是保守。2）线弹性则满足力正比于变形（位移）。这两点在能量法的证明过程中要用到。在线弹性系统中势能就是应变能。

最小势能原理是是势能驻值原理在特殊情况。势能驻值原理说的是，任何系统真实位移状态总是使其势能取驻值（即一阶变分为零）。注意，是驻值（极值）而不是最值。不过这一原理用在线弹性问题中，恰好变为取最小值，因此被称作最小势能原理。

同样可以证明出最小余能原理——在线弹性系统中，对于稳定的平衡状态，真实应力使总应变余能取最小值。因此，在线弹性保守系统中，最小势能原理和最小余能原理是等价的表述。
这里"余能"是材料力学和弹性力学中一个特殊的概念。它并不是真实存在的一种能量，而是一种数学计算出来具有能量量纲的一个量，因为它与真实存在的应变能存在"互余"的数学关系。
概括来说，应变余能密度与应变能密度的表达式有一种对偶关系——应变能和应变余能的微分表达式中，应力与应变的位置是互换的，前者是应变作为微分量，后者是应力作为微分量。应变（余）能密度作体积分便得到应变（余）能，从而应变余能加上应变能就是应力力乘以应变的体积分（前面不需要加上1/2），这就是“互余”的意思。

最为广义、反应本质的则是虚功原理（它适用任何情况，不一定非得是线弹性）：如果弹性体是处于静力平衡状态的，外力在虚位移上所做的虚功之和，等于真实应力（内力）各分量在对应的虚应变上所做的虚功（即虚应变能）之和（即所有可能的位形中，真实位形满足内外力虚功之和为零）。
注：虚位移（或虚应变）指的是弹性体在满足约束条件及连续性条件下的无限小可能位移，它与受力无关，只是一种符合条件可能性（相当于位移函数的变分）。

这个自然法则确定了平衡后的真实位形。势能极值原理只是虚功原理的特殊情形。

无论是用哪种能量法，必须知道系统的属性，比如剪切模量，弹性模量，截面惯性矩（给出几何形状就能计算），不给就无法给出应变（余）能密度，无从计算能量表达式，就无法利用能量法求解。这也就是为什么yinhow说必须知道弹性势能的表达式了。

kastin

弹性力学倒是需要较强的数学功底，毕竟需要解出应力函数（平衡方程是一组偏微分方程组），而且很多解都是级数解或者带半猜测性质的，不过做起来挺有意思的，而且复变函数中的很多知识也能用到求解应力函数的中。但材料力学就简单多了，它是弹性力学的简化版，做了很多假设，所以用的都是很简单微积分知识。

考虑一个长为 `l` 的等截面均匀水平悬臂梁（或细杆），在自由端施加竖直向下的集中力 `F`，经变形后达到静平衡。

以悬点为原点，水平向右为 `x` 轴，竖直向下为 `y` 轴建立坐标系。设距离原点水平距离为 `x` 处竖直位移为 `y(x)`。在沿竖直方向截断梁（或杆），显然截面处内力只有弯矩。对右边一段受力分析可知，截面处的弯矩与外力平衡，所以弯矩 `M(x)=F(l-x)`.
所以$$\kappa=\frac{y''}{(1+y'^2)^{3/2}}=\frac{F(l-x)}{EI}\tag{1}$$下面考虑边值条件，在固定端，位移为零：`y(0)=0`；固定端不会发生转动，仍然保持原来水平角度：`y'(0)=0`.

事实上，上述公式左端（也就是9楼给出的公式）有很强的几何意义。因为曲率的定义是，曲线切线的转角 `\theta` 相对于弧长 `s` 的变化率，并注意到弧微分关系 `ds=\sqrt{1+y'^2}\,dx` 和几何关系`\theta=\arctan y'`，因此就有$$\kappa=\frac{d\theta}{ds}=\frac{dy'}{(1+y'^2)\sqrt{1+y'^2}dx}=\frac{y''}{(1+y'^2)^{3/2}}\tag{2}$$

1. 若变形较小，最大转角 `y'\to 0`，故(1)可近似为 `EI y''=F(l-x)`，积分可得 `\D y=-\frac{F}{6EI}x^3+\frac{Fl}{2}x^2+c_1x+c_2`，代入边界条件得到挠曲线（变形后的曲线）方程$$y=-\frac{F}{6EI}x^3+\frac{Fl}{2}x^2$$这是一个三次曲线。

2. 这次我们换个稍微精确一点的近似方法来处理（思想来自于8楼资料）。注意到弧微分关系可变形为$$\frac{dy}{dx}=\pm \frac{dy/ds}{\sqrt{1-(dy/ds)^2}}\tag{3}$$如果我们能从等式(1)和(2)中得到 `dy/ds` 的近似表达，那么代入方程(3)后，相当于将方程(1)降为了一阶。

为方便起见，记 `\D a=\frac{F}{2EI}`。因为有小角度近似 `\arctan y'\sim y'\,(y'\to 0)` 代入(2)得到$$\frac{d\theta}{ds} \sim \frac{dy'}{ds}=\frac{d(dy/dx)}{ds}\stackrel{\color{green}{(a)}}{=}\frac{d^2y}{dxds}\stackrel{\color{green}{(b)}}{=}\frac{d(dy/ds)}{dx}=2a(l-x)\tag{4}$$因为 `s=0` 时 `s=x`，故`\D \frac{dy}{ds}|_{s=0}=\frac{dy}{dx}|_{s=0}\*\frac{ds}{dx}|_{x=0}=0\* 1=0`因此$$\frac{dy}{ds}=\frac{dy}{ds}|_{s=0}+\int_0^x 2a(l-u)\dif u=a(2lx-x^2)$$于是，将结果代入(3)便得到稍精确降阶微分方程$$\frac{dy}{dx}=\pm \frac{2alx-ax^2}{\sqrt{1-(2alx-ax^2)^2}}\tag{*}$$使用Mathematica求解(物理解满足转角 `y'>0`，前面取正号 ). 上述方程都没有解析解，通过Mathematica数值求解 `a=3,l=1` 得到结果

1. s1 = NDSolveValue[{y''[x] == 3 (1 + y'[x]^2)^(3/2) (1 - x), y[0] == 0,
   
2.      y'[0] == 0}, y, {x, 0, 0.5}];
   
3. s = NDSolveValue[{y'[x] == 3 (2 x - x^2)/Sqrt[1 - 3^2 (2 x - x^2)^2],
   
4.     y[0] == 0}, y, {x, 0, 0.5}];
   
5. Plot[{s1[x], s[x]}, {x, 0, 0.4}, PlotLegends -> {"精确解", "近似解"}]

复制代码

a(l-x)情形的解

但是可以肯定的是，它们的解与椭圆积分有关，比如上面的积分，如果方程中的 `(l-x)` 换成 `x` 那么结果恰好就是第一类不完全椭圆积分 `F(x|m)` 和第二类不完全椭圆积分 `E(x|m)` 在参数 `m=-1` 时的复合函数 $$y=\frac{1}{\sqrt{a}}E(\arcsin(\sqrt{a}x)|-1)-\frac{1}{\sqrt{a}}F(\arcsin(\sqrt{a}x)|-1)\tag{5}$$Mathematica代码

1. sol1 = DSolve[{y'[x] == a x^2/Sqrt[1 - a^2 x^4], y[0] == 0}, y[x],
   
2.    x] // Simplify[#, a > 0] &

复制代码

2. 若各个位置的位移（挠度）尺度跨度较大（非小变形）。
    a. 但在原点附近，位移较小，并且转角接近零，即 `x \to 0，y' \to 0`，上述近似解都成立。
    b. 在全局范围，则不能适用上述近似。直接求解方程(1)，解有两个（一正一负），考虑到物理解必须满足 `y>0`，使用Mathematica计算出结果仍然是(5)：

1. sol = DSolve[{y''[x] == 2 a (1 + y'[x]^2)^(3/2) x, y[0] == 0,
   
2.    y'[0] == 0}, y[x], x] // Simplify[#, a > 0]&

复制代码

令`a=3`，绘图比较

1. Plot[{y[x] /. sol1 /. a -> 3, y[x] /. First[sol] /. a -> 3}, {x, 0,
   
2.   1}, PlotLegends -> {sol1[[1]][[1]] // TraditionalForm,
   
3.    First[sol][[1]] // TraditionalForm}]

复制代码

结果如下

ax情形的解

很显然，(*)与(1)不等价，这在上面的图像中给出了
我猜测，问题出在：要么Mathematica直接求解(1)给出的解只是一个分支，或者是某种情况下退化的解；要么(4)式中的步骤(a)或步骤(b)有问题（因为全微分不等于偏微分，弧长 `s` 也是坐标 `x` 的函数）。


下面回到能量法。弹性势能在弹性系统中指的就是应变能（条件是“完全弹性“，即能量不会损耗），因为势能对应保守系统，保守系统能量不会耗散，而是相互转化，且与路径无关，就意味着可全微分。弹性势能来自于弯矩做功。由于假设纯弯曲，所以只有弯曲应变能，单位弧长上转角为`\D\frac{1}{\rho}=\kappa`单位长度梁段弯矩做功(弯曲应变能密度)为`\D V_{\varepsilon}=\int_0^k M\dif k=\int_0^{\kappa} EIk\dif k=\frac{1}{2}EI\kappa^2=\frac{1}{2}EI\kappa^2=\frac{M^2}{2EI}`总能为整个梁的弧长积分$$V=\frac{1}{2}EI\int \kappa^2 ds $$即找到上述积分的最小值，而`\kappa`与变形曲线的关系就是$$\kappa=\frac{y''}{(1+y'^2)^{3/2}}$$这就相当于一个变分问题了。