一根绳子，当受到一定的力，弯矩，扭矩时他的形状是什么样的？

  • 预备工作
  • 建立Euler弹性线方程
  • 分析求解
  • 数值求解
  • 黎曼空间中的弹性线
  • 附录Cpp代码

预备工作

分部积分法的妙用

浅谈分部积分公式中通过迭代得到了：

通过这个公式也得到了高阶Lagrange方程，这里简单的总结一下，方便后边使用：
含有高阶导数的Lagrange方程为：

对应守恒量：

几个椭圆函数

一些符号变换而已

建立Euler弹性线方程

不考虑弹性线的伸长时，弹性线的能量可以写为：


这个不考虑伸长有一定的物理意义，并不意味着弹性线的弹性系数很大。比如说一根松弛状态弯曲的橡皮筋，虽然橡皮筋很容易伸缩变形，但是由于绳子并没有绷紧，因此也可以看成不考虑弯曲能的情况。
在不考虑转动时，能量仅仅由曲率进行贡献。
曲率与矢径的关系为：。
利用Lagrange乘数法，得到能量泛函：

这时应用前面的二阶导数的Lagrange方程，得到：

守恒量为：

上面得到的只是矢径需要满足的关系，要想把矢径和曲线的内禀量（什么是内禀量？内禀量是系统本身的性质，不随坐标系选取而改变，这里指的是曲率和挠率。）联系起来，需要把矢径投影到Frenet标架（不一定是Frenet标架，只要是曲线的随体标架都可以）中，并利用Frenet关系进行化简。
（Add：如何理解内禀坐标？ 首先，对于一个空间曲线来说，选取不同位置为基点进行描述，矢径参数方程肯定是不同的，但是这个量与选取基点的方式无关。这也说明了为什么描述曲线自身性质的Frenet标架中不会含有，而含有的各阶导数。）
满足关系式：

投影后得到：（略去大量计算）

要满足必须让三个分量为零：

上式积分得到：

后续只需要通过边界条件求解方程，然后把参数空间转换到欧氏空间中。

分析求解

？？？？？？？？？？？？？

数值求解

弹性线问题中，为了方便起见，总是以内禀参数进行求解的。但是要想进行可视化，必须进行坐标的变换，简而言之就是从内禀参数空间变换到欧式空间。这个变换过程就是求解Frenet方程。
数值求解的关键： 求解Frenet方程，把参数空间转化为欧氏空间。
首先会想到二维情况如下图：
1
猛一看这不是很简单嘛。啪！很快啊！于是通过几何关系很快就写出了参数空间与欧氏空间的变换关系。

离散化：

但是，如果考虑三维空间的情况这个几何关系是什么呢？
也就是说，当从曲线上的一点向前前进了的距离，这个前进的过程中既有弯曲又有扭转，到达的位置怎么表达呢？
要想和二维情况一样通过几何直觉来得到这个关系就很难了。
于是，可以发现，上面的通过所谓的“几何直觉”得到转换关系的过程只不过是参数化的过程，借助中间参数进行转换。

三维可不可以照猫画虎来一遍？
于是想到，同样把弧长进行参数化如下：

与内禀参数的关系怎么写？
通过Frenet方程求解得到：


到这一步，只要求解出和的关系就能把参数空间转化为欧氏空间。但是可以发现上面是一个复杂的非线性方程，离散起来需要花费很大功夫，因此这个方法行不通。
绕了一个大弯，发现绕进了胡同里。
既然参数空间的转化与Frenet方程有关，为何不直接求解Frenet方程？
那么如何把Frenet方程转化到欧氏空间呢？
通过这个关系式。
至此已经有了求解思路，感觉这个思考过程很有意思。
也许前面的弯路并没有白走，二维情况下的几何直觉引导我们想到三维情况。而三维情况下几何直觉不管用了，又得回去归纳前面通过“几何直觉”求解二维问题的流程。好不容易归纳好了，却发现同样的方法三维问题根本没法解。而这个时候发现了三维求解过程中与二维情况的不同：利用了Frenet方程。然后想到可能与Frenet方程有关，最终去求解Frenet方程，得到了参数空间向欧氏空间转换的一般方法。

离散Frenet方程以及得到：

至此，Euler弹性线（仅考虑弯曲能）的求解思路总结如下：

  • 求解弹性线方程
  • 得到内禀参数的演化规律
  • 通过Frenet方程求解标架的演化规律
  • 通过把参数空间映射到欧氏空间

PS: 以上是一个关于曲率的二阶方程，可以通过给定初始曲率以及曲率导数进行迭代（初值问题），也可以通过给定两个端点的导数利用打靶法进行求解（边值问题）。此外，Frenet方程的求解用Euler法会产生很大误差，连求解基本的圆误差都会很大，用预报矫正算法就会好很多。

上面已经给出了求解的方法，下面是几个求解结果的可视化(按照曲率着色)

s
2 1 2 1 50
1.5 1 0.7 0.8 50
0.7 1 2 0.75 50

2
3
4

黎曼空间中的弹性线

？？？？？？？？？？？？？

Cpp代码

主函数(main.cpp)：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#############Elastica类(EulerElasticaRod.h)#############
#include<iostream>
#include<fstream>
#include<math.h>
using namespace std;
struct Para{
	static const int N=2000;
	double kappa[N]{},tau[N]{},coord[3][N];
};
class DER{
public:
	//	Para ks2coord();
	DER(Para &p1,double s1,double* kppa01,double k01,double dk01,double tau0,double ks,double lambda1);//构造函数
	//不需要析构函数！！
	Para SolveER();
	double erfen();
	double f();
	double f1();
private:
	Para p;double s;double* kppa0;double k0;double dk0;double tau0;double ks;double lambda;	
};
DER::DER(Para &p1,double s1,double* kppa01,double k01,double dk01,double tau01,double ks0,double lambda1){
	p=p1;s=s1;kppa0=kppa01;k0=k01;dk0=dk01;ks=ks0;lambda=lambda1;tau0=tau01;
}
//就是不需要析构函数！！！
//定义public函数
Para DER::SolveER(){
	int NN=p.N; 
	double ds=s/NN;
	double c=k0*k0*tau0;
	p.kappa[0]=k0;p.kappa[1]=ds*dk0+k0;p.tau[0]=tau0;
	for (int i =1;i<NN-1;i++) 
	{
		p.kappa[i+1]=2*p.kappa[i]-p.kappa[i-1]+pow(ds,2)*(-0.5*pow(p.kappa[i],3)+0.5*lambda*p.kappa[i]+pow(c,2)/pow(p.kappa[i],3));
		p.tau[i+1]=c/pow(p.kappa[i+1],2);
	}
	return p;	
} 


double DER::erfen(){
	cout<<"ks"<<"         "<<"err"<<endl;
	double ff0=this->f();double ff1=this->f1();
	double eps=1e-6;
	
	while(true){
		dk0=dk0-(ff0-ks)/ff1;
		ff0=DER(p,s,kppa0,k0,dk0,tau0,ks,lambda).f();
		ff1=DER(p,s,kppa0,k0,dk0,tau0,ks,lambda).f1();
		cout<<ff0<<"    "<<abs(ff0-ks)<<endl;
		if (abs(ff0-ks)<eps) break;	
	}
	return dk0;	
}
double DER::f(){
	Para re=SolveER();
	return re.kappa[re.N-1];
}
double DER::f1(){
	double ddk=0.001;
	return (DER(p,s,kppa0,k0,dk0+ddk,tau0,ks,lambda).f()-DER(p,s,kppa0,k0,dk0,tau0,ks,lambda).f())/ddk;
}


#############主函数#############
#include "EulerElasticaRod.h"
Para ks2coord(Para &p, double s) {
	double Frenet[3][p.N][3] {}, ds = s / p.N;
	Frenet[0][0][0] = 1;
	Frenet[0][0][1] = 0;
	Frenet[0][0][2] = 0;
	Frenet[1][0][0] = 0;
	Frenet[1][0][1] = 1;
	Frenet[1][0][2] = 0;
	Frenet[2][0][0] = 0;
	Frenet[2][0][1] = 0;
	Frenet[2][0][2] = 1;
	for (int j = 0; j < 3; j++) {
		p.coord[j][0] = 0;
		for (int i = 1; i < p.N; i++) {
			Frenet[0][i][j] = Frenet[0][i - 1][j] + p.kappa[i - 1] * ds * Frenet[1][i - 1][j];
			Frenet[2][i][j] = Frenet[2][i - 1][j] - p.tau[i - 1] * ds * Frenet[1][i - 1][j];
			Frenet[1][i][j] = Frenet[1][i - 1][j] - p.kappa[i - 1] * ds * Frenet[0][i - 1][j] + p.tau[i - 1] * ds * Frenet[2][i - 1][j];
			//预报校正算法 PS:简单的Euler法会导致比较大的误差，不信试试看！
			Frenet[0][i][j] = Frenet[0][i - 1][j] + p.kappa[i - 1] * ds * (Frenet[1][i - 1][j] + Frenet[1][i][j]) / 2;
			Frenet[2][i][j] = Frenet[2][i - 1][j] - p.tau[i - 1] * ds * (Frenet[1][i - 1][j] + Frenet[1][i][j]) / 2;
			Frenet[1][i][j] = Frenet[1][i - 1][j] - p.kappa[i - 1] * ds * (Frenet[0][i - 1][j] + Frenet[0][i][j]) / 2 + p.tau[i - 1] * ds * (Frenet[2][i - 1][j] + Frenet[2][i][j]) / 2;
			p.coord[j][i] = p.coord[j][i - 1] + Frenet[0][i - 1][j] * ds;
		}
	}
	return p;
}
int main() {
	Para p;
	double s = 50, kappa[p.N] {}, k0 = 0.9, dk0 = 0.9, tau0 = 1, lambda = 0.75, ks = 2;
	double ppp = DER(p, s, kappa, k0, dk0, tau0, ks, lambda).erfen();
	Para q = DER(p, s, kappa, k0, ppp, tau0, ks, lambda).SolveER();
	//	dk0=DER(p,s,kappa,k0,dk0,tau0,ks,lambda).erfen();
	Para ret = ks2coord(q, s);
	cout << ret.kappa[q.N - 1] << endl;
	ofstream out("out.data");
	for (int i = 0; i < p.N; i++) {
		out << ret.coord[0][i] << "   " << ret.coord[1][i] << "   " << ret.coord[2][i]<<"    "<<ret.kappa[i]<< endl;
	}
	out.close();
	return 0;
}

本文采用CC-BY-SA-3.0协议，转载请注明出处
作者: 得意喵~