吾生有涯 学海无涯
析模有界 知识无界

利用OpenFOAM求解常微分方程

内容纲要

源地址:http://hassankassem.me/posts/ode/。作者:Hassan Kassem


本文描述如何利用OpenFOAM求解常微分方程。

本文并非要对每个常微分方程(ODE)求解器进行全面的回顾或解释,而是一步一步地展示如何定义一个ODE,以及如何编写一个小程序来求解此ODE。

1 数学背景

为了利用数值方法求解高阶常微分方程,可以其简化为一阶常微分方程组进行求解。降阶产生的常微分方程中方程的数量等于原始常微分方程的阶数。

一般的常微分方程可以表示为:

式中,为导数阶数。

可以将式(1)转化为:

最后,常微分方程可以表示为以下的形式:

转换过程可概括为两个步骤;变量转换及微分。

在OpenFOAM中求解ODE,刚性求解器额外需要雅可比矩阵和时间偏导数。

雅克比矩阵如下:

其中是简化系统中的方程式。及其关于自变量的偏导数:

2 例子:运动方程

这里以单自由度系统(简单质量-弹簧系统)的运动方程为例,其为一个二阶常微分方程:

式中,为质量,为刚度,为位移。

该方程有一个精确解,表示系统的谐波振荡为(假设初始速度为零);

式中,为初始位移,为系统的固有频率,

求解此方程的分为三步:

  1. 转化为
  1. 微分(先替换)

则有:

  1. Jacobian

雅克比矩阵仅用于刚性求解器。对于这里的微分方程,其雅克比矩阵为矩阵:

时间导数为:

3 程序设计

在OpenFOAM中,典型的做法是搜索类似的实现作为起点。幸运的是,OpenFOAM源文件中包含了一个小测试程序。所以将其复制到运行目录,然后编译并运行它。(RKCK45是OpenFOAM ODE解算器之一)。

3.1 准备文件

运行一下命令:

run
mkdir myODE
cd myODE
cp -r $FOAM_APP/test/ODE .
cd ODE

ODE文件夹中包含一个Make文件夹以及一个名为Test-ODE.C的源文件。

编译源文件并测试运行:

wmake
Test-ODE RKCK45

如下图所示表示编译且运行成功。

注:这里的RKCK45为四阶-五阶Cash-Karp Runge-Kutta法。除该算法外,OpenFOAM还提供了  SIBS, Rosenbrock34, RKF45, RKDP45, rodas34, rodas23, Rosenbrock23, Rosenbrock12, EulerSI, Euler, 及 Trapezoid等算法,有兴趣的话可以自行尝试。

3.2 修改源文件

修改Test-ODE.C文件,如下所示。

#include "argList.H"
#include "IOmanip.H"
#include "ODESystem.H"
#include "ODESolver.H"
using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * //
class testODE
:
public ODESystem
{

public:
// 质量
const scalar m_;
// 刚度系数
const scalar k_;
testODE(const scalar &mass, const scalar &stiff) : ODESystem(),
m_(mass),
k_(stiff)
{
}

label nEqns() const
{
return 2;
}

// 必须实现的虚函数,计算导数,与原文有差异
void derivatives(
const scalar x,
const scalarField &y,
const label li,
scalarField &dydx)
const
{
dydx[0] = y[1];
dydx[1] = -(k_ / m_ + VSMALL) * y[0];
}

// 必须实现的虚函数,计算雅克比矩阵,与原文有差异
void jacobian(
const scalar x,
const scalarField &y,
const label li,
scalarField &dfdx,
scalarSquareMatrix &dfdy)
const
{
dfdx[0] = 0.0;
dfdx[1] = 0;

dfdy(0, 0) = 0.0;
dfdy(0, 1) = 1.0;

dfdy(1, 0) = -(k_ / m_ + VSMALL);
dfdy(1, 1) = 0;
}
};

// * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
argList::validArgs.append("ODESolver");
argList args(argc, argv);

// 指定质量与刚度系数
const scalar mass = 2.0;
const scalar stiff = 2048.0;

// 指定初始条件
const scalar y0 = 0.2;
const label n = 200;
const scalar endTime = 1.0;

// 创建ODE对象
testODE ode(mass, stiff);

dictionary dict;
dict.add("solver", args[1]);

autoPtr odeSolver = ODESolver::New(ode, dict);

scalar xStart = 0.0;
const scalar dx = endTime / n;
scalarField yStart(ode.nEqns());
yStart[0] = y0;
yStart[1] = 0.0;
scalar dxEst = 0.1;
scalar xEnd = 0.0;
label li =1;
scalarField dyStart(ode.nEqns());

for (label i = 0; i < n; i++)
{
xEnd = xStart + dx;
ode.derivatives(xStart, yStart, li, dyStart);
odeSolver->solve(xStart, xEnd, yStart, li, dxEst);
xStart = xEnd;
Info << xStart << " " << yStart[0] << endl;
}
Info << "n Endn"
<< endl;
return 0;
}
  • 运行命令编译源代码:
wmake

如下图所示表示编译成功。

  • 运行以下命令执行程序
Test-ODE RKCK45 > log
  • 启动gnuplot执行以下命令
mass = 2.0
stiff = 2048.0
x0 = 0.2

set grid
set samples 200
set xlabel 'Time (s)'
set ylabel 'Displacement (m)'
set yrange [-0.3:0.3]

f(x) = x0*cos(sqrt(stiff/mass)*x)
plot 'log' u 1:2 w p lw 2 t "OpenFOAM", f(x) w l lw 2 t "Exact"

结果如下图所示。数值解与解析解吻合良好。


利用OpenFOAM中自带的ODE求解器还是蛮麻烦的,像上面简单的二阶常微分方程,利用matlab只需要几行程序就能搞定。


(本文结束)

本篇文章来源于微信公众号: CFD之道

赞(0) 打赏
版权声明:未经允许,请勿随意用于商业用途。
文章名称:《利用OpenFOAM求解常微分方程》
文章链接:https://www.topcfd.cn/18788/
本站资源仅供个人学习交流,请于下载后24小时内删除,不允许用于商业用途,否则法律问题自行承担。
分享到

说两句 抢沙发

评论前必须登录!

 

觉得文章有用就打赏一下文章作者吧

非常感谢你的打赏,我们将继续给力更多优质内容,让我们一起创建更加美好的网络世界!

支付宝扫一扫

微信扫一扫

登录

找回密码

注册