JiaHao Li

JiaHao Li

Ph.D. student, Department of Modern Mechanics, University of Science and Techonology of China.

弹性线问题与数值求解

5 分钟阅读时长

发布时间：

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

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

预备工作

分部积分法的妙用

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

\[\int gf^{\alpha}dx+(-1)^{\alpha+1}\int fg^{\alpha}dx=\sum_{i=0}^{\alpha-1}(-1)^{i}f^{\alpha-1-i}g^{i}\]

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

含有高阶导数的Lagrange方程为：

\[\frac{d^2}{dx^2}\left( \frac{\partial L}{\partial y_{xx}} \right) -\frac{d}{dx}\left( \frac{\partial L}{\partial y_x} \right) +\frac{\partial L}{\partial y}=0\]

对应守恒量：

\[\epsilon \left( \frac{\partial L}{\partial y_x}-\frac{d}{dx}\frac{\partial L}{\partial y_{xx}} \right) +\dot{\epsilon}\frac{\partial L}{\partial y_{xx}}=const\]

几个椭圆函数

一些符号变换而已

建立Euler弹性线方程

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

\[E=\frac{1}{2}\int{\kappa ^2\left( s \right) +\tau ^2\left( s \right) ds}\] \[||r'(s)||=1\]

这个不考虑伸长有一定的物理意义，并不意味着弹性线的弹性系数很大。比如说一根松弛状态弯曲的橡皮筋，虽然橡皮筋很容易伸缩变形，但是由于绳子并没有绷紧，因此也可以看成不考虑弯曲能的情况。

在不考虑转动时，能量仅仅由曲率进行贡献。

曲率与矢径的关系为：

\[\tau=||r''(s)||\]

利用Lagrange乘数法，得到能量泛函：

\[E=\int{||r''\left( s \right) ||^2+\Lambda \left( ||r'\left( s \right) ||^2-1 \right)}ds\]

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

\[r''''\left( s \right) -(\Lambda r'\left( s \right))' =0\]

守恒量为：

\[\epsilon \left( s \right) \cdot \left( \Lambda -r'''\left( s \right) \right) +\epsilon'\left( s \right) \cdot r''\left( s \right) =const=J\]

上面得到的只是矢径需要满足的关系，要想把矢径和曲线的内禀量（什么是内禀量？内禀量是系统本身的性质，不随坐标系选取而改变，这里指的是曲率和挠率。）联系起来，需要把矢径投影到Frenet标架（不一定是Frenet标架，只要是曲线的随体标架都可以）中，并利用Frenet关系进行化简。

（Add：如何理解内禀坐标？ 首先，对于一个空间曲线来说，选取不同位置为基点进行描述，矢径参数方程肯定是不同的，但是$r’(s)$这个量与选取基点的方式无关。这也说明了为什么描述曲线自身性质的Frenet标架中不会含有$r$，而含有$r$的各阶导数。）

满足关系式：

\[\begin{cases} r'=T \\ T'=\kappa N \\ N'=-\kappa T+\tau B \\ B'=-\tau N \end{cases}\]

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

\[\begin{aligned} &\left( \Lambda'\left( s \right) +3\kappa \left( s \right) \kappa'\left( s \right) \right) T \\ &\quad -\left( -\Lambda \left( s \right) \kappa \left( s \right) +\kappa'' \left( s \right) -\kappa \left( s \right) \tau ^2\left( s \right) -\kappa ^3\left( s \right) \right) N \\ &\quad -\left( 2\tau \left( s \right) \kappa '\left( s \right) +\kappa \left( s \right) \tau'\left( s \right) \right) B=0 \end{aligned}\]

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

\[\begin{cases} \Lambda' \left( s \right) +3\kappa \left( s \right) \kappa' \left( s \right) =0 \\ \kappa'' \left( s \right) -\Lambda \left( s \right) \kappa \left( s \right) -\kappa \left( s \right) \tau ^2\left( s \right) -\kappa ^3\left( s \right) =0 \\ 2\tau \left( s \right) \kappa'\left( s \right) +\kappa \left( s \right) \tau'\left( s \right) =0 \end{cases}\]

上式积分得到：

\[\begin{cases} \Lambda \left( s \right) =-\frac{3}{2}\kappa ^2\left( s \right) +\frac{\lambda}{2} \\ \kappa''\left( s \right) +\frac{1}{2}\kappa ^3\left( s \right) -\frac{\lambda}{2}\kappa \left( s \right) -\frac{c^2}{k^3\left( s \right)}=0 \\ \kappa ^2\left( s \right) \tau \left( s \right) =const=c \end{cases}\]

后续只需要通过边界条件求解方程，然后把参数空间$(\kappa,\tau,s)$转换到欧氏空间$(x,y,z)$中。

分析求解

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

数值求解

弹性线问题中，为了方便起见，总是以内禀参数进行求解的。但是要想进行可视化，必须进行坐标的变换，简而言之就是从内禀参数空间变换到欧式空间。这个变换过程就是求解Frenet方程。

数值求解的关键： 求解Frenet方程，把参数空间转化为欧氏空间。

首先会想到二维情况如下图：

1

猛一看这不是很简单嘛。啪！很快啊！于是通过几何关系很快就写出了参数空间与欧氏空间的变换关系。

\[\begin{cases} dx=\cos \theta ds \\ dy=\sin \theta ds \\ \kappa =\theta _s \end{cases}\]

离散化：

\[\begin{cases} \theta _k=\theta _{k-1}+\kappa \left( \left( k-1 \right) \Delta s \right) \Delta s \\ x_k=x_{k-1}+\cos \theta _{k-1}\Delta s \\ y_k=y_{k-1}+\sin \theta _{k-1}\Delta s \end{cases}\]

但是，如果考虑三维空间的情况这个几何关系是什么呢？

也就是说，当从曲线上的一点$(x_0,y_0,z_0)$向前前进了$ds$的距离，这个前进的过程中既有弯曲又有扭转，到达的位置怎么表达呢？

要想和二维情况一样通过几何直觉来得到这个关系就很难了。

于是，可以发现，上面的通过所谓的“几何直觉”得到转换关系的过程只不过是参数化的过程，借助中间参数$\theta$进行转换。

\[ds=\sqrt{dx^2+dy^2} \\ \begin{cases} dx=\cos \theta ds \\ dy=\sin \theta ds \end{cases}\]

三维可不可以照猫画虎来一遍？

于是想到，同样把弧长进行参数化如下：

\[ds=\sqrt{dx^2+dy^2+dz^2} \\ \begin{cases} dx=\cos \theta \sin \varphi ds \\ dy=\sin \theta \sin \varphi ds \\ dz=\cos \varphi ds \end{cases}\]

$(\theta,\tau)$与内禀参数的关系怎么写？

通过Frenet方程求解得到：

\[\kappa ^2=\left\| \frac{d\left( \frac{dx}{ds},\frac{dy}{ds},\frac{dz}{ds} \right)}{ds} \right\| ^2 ={\theta _s}^2\sin ^2\varphi +{\varphi _s}^2\] \[\begin{aligned} \tau &=||\dot{B}||=||\dot{T}\times N+T\times \dot{N}|| \\ &=\left\| \frac{d\left( \frac{T \times \dot{T}}{\kappa} \right)}{ds} \right\| \\ &=\theta_s\cos\varphi+\frac{\theta_{ss}\varphi_s\sin\varphi+\theta_s(\varphi_s^2\cos\varphi-\varphi_{ss}\sin\varphi)}{\theta_s^2+\sin^2\varphi+\varphi_s^2} \end{aligned}\]

到这一步，只要求解出$(\kappa,\tau)$和$(\theta,\phi)$的关系就能把参数空间转化为欧氏空间。但是可以发现上面是一个复杂的非线性方程，离散起来需要花费很大功夫，因此这个方法行不通。

绕了一个大弯，发现绕进了胡同里。

既然参数空间的转化与Frenet方程有关，为何不直接求解Frenet方程？

那么如何把Frenet方程转化到欧氏空间呢？

通过$T(s)=r’(s)$这个关系式。

至此已经有了求解思路，感觉这个思考过程很有意思。

也许前面的弯路并没有白走，二维情况下的几何直觉引导我们想到三维情况。而三维情况下几何直觉不管用了，又得回去归纳前面通过“几何直觉”求解二维问题的流程。好不容易归纳好了，却发现同样的方法三维问题根本没法解。而这个时候发现了三维求解过程中与二维情况的不同：利用了Frenet方程。然后想到可能与Frenet方程有关，最终去求解Frenet方程，得到了参数空间向欧氏空间转换的一般方法。

离散Frenet方程以及$T=r’(s)$得到：

\[\begin{cases} T_k=T_{k-1}+\Delta s \kappa ( (k-1)\Delta s ) N_{k-1} \\ N_k=N_{k-1}-\Delta s \kappa ( (k-1)\Delta s ) T_{k-1}+\Delta s\tau ( (k-1)\Delta s ) B_{k-1} \\ B_k=B_{k-1}-\Delta s \tau ( (k-1)\Delta s ) N_{k-1} \\ r_k=r_{k-1}+\Delta s*T_{k-1} \end{cases}\]

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

  1. 求解弹性线方程 \(\kappa'' \left( s \right) +\frac{1}{2}\kappa ^3\left( s \right) -\frac{\lambda}{2}\kappa \left( s \right) -\frac{c^2}{k^3\left( s \right)}=0\)
  2. 得到内禀参数$(\kappa(s),\tau(s))$的演化规律
  3. 通过Frenet方程求解$(T,N,B)$标架的演化规律
  4. 通过$r’(s)=T$把参数空间映射到欧氏空间

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

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

$\kappa(0)$$\tau(0)$$\kappa(s)$$\lambda$s
212150
1.510.70.850
0.7120.7550

2 3 4

黎曼空间中的弹性线

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

Cpp代码

主函数(main.cpp)：

#############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;
}

你可能感兴趣的

如何理解板、壳、杆理论？

少于 1 分钟阅读时长

发布时间：

板、壳、杆理论无论对于工程实践还是理解复杂物理现象都具有重要应用，但是如何很好的理解这套基础理论并不容易。

Non-Euclidean Euler Buckling

6 分钟阅读时长

发布时间：

This notes presents a basic theory for curved surface constrained kirchhoff rods theory.