二维FEM建模及求解（入门级：数理方程到代码编写）

前言：好久没有更新了，最近一直在看关于脑电正问题建模的问题，重点是关于FEM头模型搭建的问题。网上好多资料要不就是纯理论讲解，要不就是直接上打断代码的。所以趁着现在的理解还是在由少到多的阶段，记录一下，面向群体是跟我一样基础比较薄弱，想从理论到代码编写捋一遍FEM方法的同志们。

先把要解决的问题摆出来：

1）相关数理方程化简步骤；

2）离散化，用什么原理？合理性？误差来源；

3）代码实现。

注：1）数理方程就是将物理现象搭建为数学模型，很有趣的一门课，感兴趣的可以了解；2）关于基函数的理解还不够到位，如果有人想了解基函数相关内容请移步。

1. 数学模型搭建和相关处理

数学模型为：

将公式弱化（试函数为v）：

函数弱化思想见相应参考文章。

采用分部积分有：

注：公式弱化还能很好反映原微分公式的前提是积分的面积足够小，可以参考高数积分定理的证明思想，这也是需要进行离散化的原因。

2. 离散化（三角形剖分）

需要求解的u用uh表示，j 表示第j个element，每个element包含3个顶点。 $\varphi_j$ 为插值基函数，为3*1的关于坐标x、y的线性函数矩阵， $\hat{u}_j$ 为相应的权值（常数），为1*3的矩阵。

注：1）线性函数：是因为我们采用的是线性插值；

2）3*1的函数矩阵是因为每个element由3个节点组成。第 j 个element可由如下公式表示：

其中，N的推理和计算见下个文档。

对离散化的公式应用弱形式有：（取 $\varphi_j$ 的转置为试函数）

3. 代码实现

代码实际上就是对 $K_j$ 和 $f_j$ 进行迭代积分。

具体见如下：（已经进行了详细注释）

function u = FEM_solver(nds,els,bcs)nnd = size(nds,1);u = zeros(nnd,1);K =zeros (nnd,nnd);f = zeros(nnd,1);
for j = 1:size(els,1)%sti = stima(nds(els(j,:),:))K(els(j,:),els(j,:)) = K(els(j,:),els(j,:)) + stima(nds(els(j,:),:));    % 仅对角线上有值,stima为第j个Kf(els(j,:)) = f(els(j,:)) + ones(3,1) * det([1,1,1;nds(els(j,:),:)'])/6;
end
freends = setdiff(1:nnd,bcs);     % 去除边界上的值，即边界条件u=0的描述。
u(freends) = K(freends,freends)\f(freends);function sti = stima(vertices)   % 计算第j个Kndim = size(vertices,2);J = [ones(1,ndim+1);vertices'];   % J可逆 ，构建目的：1）用于求偏导；2）计算element的面积B = J\[zeros(1,ndim);eye(ndim)];   % 方程组求解：求试函数的偏导sti = det(J).*B*B'/prod(1:ndim);   % 试函数的转置与试函数的乘积再乘以element的面积就是Kend
end

调用主程序为：

其中网格划分采用的是distmesh工具包进行的。

相应链接为：DistMesh - A Simple Mesh Generator in MATLABhttp://persson.berkeley.edu/distmesh/

% 采用distmesh工具包进行网格生成
clear
figure
fd=@(p) sqrt(sum(p.^2,2))-1;
[p,t]=distmesh2d(fd,@huniform,0.2,[-1,-1;1,1],[]);
e=boundedges(p,t);
% 采用FEM_solver求解器求解
figure
nodes = p;
element = t;
boundaries = e;
x = FEM_solver(nodes,element,boundaries);
pp = {'FaceColor','Interp'};
trisurf(element,nodes(:,1),nodes(:,2),x,pp{:});
axis square;
colorbar

结果为：

1）网格划分结果：

2）FEM结果：

最后，回答一下一开始提出的问题：

1）数学推理主要用到采用积分的形式对微分函数进行弱化，简化计算；然后就是采用分部积分将二阶导数转化为一阶。

2）离散化是采用基函数和相应系数的方式进行的，离散化是为了配合函数弱化进行的。误差来源是采用线性插值使得小范围内的u出现了平面。合理性在于当网格足够小时，误差将越来越小。

最后的最后，我也是刚入门，若是有疑问希望不吝赐教。

参考的文章：（这三篇文章给了我很大的启示，在此谢谢他们的分享）

你的第一个有限元求解器：仅十行 MATLAB - 知乎

什么是弱形式？ - 知乎

有限元分析简介及伽辽金法_feyily的博客-CSDN博客_伽辽金有限元方法

二维FEM建模及求解（入门级：数理方程到代码编写）相关推荐

doa的matlab算法,基于MATLAB的DOA估计算法的二维仿真建模
基本描述: DOA算法使用T形或L形天线接收微波信号,计算到达角,然后对二维坐标建模. 要解决的问题: 实际环境中存在多个相干源(例如多径效应,信号反射等),并给出了最佳DOA估计算法. 功能描述: ...
二维稳态热传导基本方程的有限元方程
二维稳态热传导基本方程的有限元方程 - 知乎二维稳态热传导基本方程的有限元求解(2) - 知乎关于学习拉格朗日矩形单元和serendipity四边形单元形函数的构造方法 - 知乎 CFD理论学习- ...
matlab二维谐振子,基于有限差分法求解的二维谐振子的MATLAB程序如下。哪位大神能帮我做个注明啊，完全看不懂啊，，急...
基于有限差分法求解的二维谐振子的MATLAB程序如下.哪位大神能帮我做个注明啊,完全看不懂啊,,急0 ____丿呆呆丶2017.04.15浏览20次分享举报 tic clc clear L=20; W ...
pythoncad标注教程_CAD 2014二维三维建模渲染标注基础与提升视频教程
CAD 2014基础与提升视频教程课程描述 CAD 2014基础与提升视频教程,偏机械类学习人群零基础,室内都可以:CAD制图者,特别是机械制图:机械.建筑.电子行业从业者主图需求:准备转行CA ...
自动驾驶 4-3 二维动态建模Dynamic Modeling in 2D
sprung mass: [车辆] 簧上质量 spring: 弹簧 damper: [车辆] 减震器 unsprung mass: 非簧载质量 tire: 轮胎 vehicle shock absor ...
求马鞍点java_二维数组马鞍点求解算法
若在矩阵A 中存在一个元素ai,j(0≤i≤n-1,0≤j≤m-1),该元素是第i行元素中最小值且又是第j 列元素中最大值,则称此元素为该矩阵的一个马鞍点.假设以二维数组存储矩阵A,试设计一个求该矩阵 ...
Python对二维矩阵沿主对角线（次对角线）翻转变换代码实现
Python对二维数组(矩阵)沿主对角线(次对角线)翻转变换代码实现目录 Python对二维数组(矩阵)沿主对角线(次对角线)翻转变换代码实现 1. 原始数据以及图示 2. 主对角线翻转及图示 3. ...
ZPL 打印条码、二维码及小票（中文/汉字），生成条码、二维码图片【Asp.Net】-含示例代码
生成条码(图1).二维码(图2)图片及打印出二维码标签(图3)效果图1 图2 ...
二维图像中的Hessian矩阵（及MATLAB代码）
文章目录一.图像中Hessian矩阵的定义及公式推导二.MATLAB代码一.图像中Hessian矩阵的定义及公式推导对于二维图像 f ( x , y ) f(x,y) f(x,y),在点 x ...

二维FEM建模及求解（入门级：数理方程到代码编写）

二维FEM建模及求解（入门级：数理方程到代码编写）相关推荐

最新文章

热门文章