EIT正问题求解--利用有限元求解节点电势

为了理解EIT正问题利用有限元求解的方法，自己通过在网上查找程序，结合EIDORS软件包所建立的EIT模型，编写了求解节点电势的matlab程序，并且进行了验证。

1.利用EIDORE建立模型，建立模型后，能够得到模型划分的单元、节点坐标等信息。程序如下：

imdl= mk_common_model('b2c',16);%建立二维模型，该函数生成的是inv_model结构。
fmdl= imdl.fwd_model;%赋值正问题模型
fmdl.stimulation= mk_stim_patterns(16,1,[0 1],[0 1], {},1);%设置激励模式
img.fwd_model= fmdl;
img.elem_data=ones(size(fmdl.elems,1),1);
img.elem_data([22:26,36:40])=3; %对介质的电导率信息进行赋值
node_v= calc_all_node_voltages(img);%求解的节点电位值，方便验证
figure;

show_fem(fmdl,[1 1 2]);%显示模型

2.将第一步中所得到的的相关信息进行提取，写入txt文件，为有限元计算提供数据

程序如下：

fid=fopen('D:\matlabanzhuang\bin\FEM\data.txt','wt');%写入文件路径
%%%%存取单元的个数、节点的个数%%%%%
[m,n]=size(fmdl.elems); %m为单元数，q为节点个数
q=size(fmdl.nodes,1);
fprintf(fid,'%d %d\n',q,m);

%%%%%存取每个单元的介质信息--电导率
for i=1:1:m
fprintf(fid,'%d\n',img.elem_data(i));
end

%%%%%存储单元的节点（此处建立的模型为二维，三角形剖分）
for i=1:1:m
for j=1:1:n
if j==n %如果一行的个数达到n个则换行，否则空格
fprintf(fid,'%d\n',fmdl.elems(i,j));
else
fprintf(fid,'%d ',fmdl.elems(i,j));
end
end
end

%%%%存储节点的坐标
[m,n]=size(fmdl.nodes);
for i=1:1:m
for j=1:1:n
if j==n %如果一行的个数达到n个则换行，否则空格
fprintf(fid,'%d\n',fmdl.nodes(i,j));
else
fprintf(fid,'%d ',fmdl.nodes(i,j));
end
end
end

fclose(fid);

3.数据得到后，即可进行有限元求解

程序如下：

%有限元方法求解各节点电势（正问题求解）
%通过相关的软件，得到模型划分的相应单元数，单元的节点数，节点的坐标以及每个单元所对应的节点
%打开数据文件 FP1数据文件指针
%读入初始数据
clear
FP1=fopen('data.txt','rt');

NPOIN=fscanf(FP1,'%d',1);%总节点数

NELEM=fscanf(FP1,'%d',1);%划分的单元数
NJIEZHI=fscanf(FP1,'%f',[1,NELEM]);%单元的电导率信息
LNODS=fscanf(FP1,'%f',[3,NELEM])';%单元定义数组，代表每个单元所包含的节点
COORD=fscanf(FP1,'%f',[2,NPOIN])';%结点坐标数组每个节点的坐标
ASTIF=ASSEMBLE(NPOIN,NELEM,NJIEZHI,COORD,LNODS);%生成总刚度矩阵
%%%%添加边界条件
m=fix(NPOIN/4);%m输入电流节点，自己设定的
n=fix(NPOIN/3);%n输出电流节点，自己设定的，如果极板（输入电流区域）不为点电极，即相应的极板对应的节点均为I/A
%%因为是进行电势的求解，需要设定一个电势参考点，即设定某一点的电势为0.
c=1;%c表示电势为0的节点，此处设定节点一为0电势点
ASTIF(c,1:c-1)=0;
ASTIF(c,c+1:NPOIN)=0;
ASTIF(c,c)=1;%c参考电势为0
b=zeros(NPOIN,1);%b为激励模式，即边界条件
b(114,1)=-1;%设置激励模式，输入节点为1，输出节点电流激励为-1
b(116,1)=1;
%b(m,1)=-1;
%b(n,1)=1;
[Q1,Q2,Q3,Q4]=gaus(ASTIF,b);%Q4即为所求电势解。
scatter(COORD(:,1),COORD(:,2),5,Q4);%散点图
[X,Y,Z]=griddata(COORD(:,1),COORD(:,2),Q4,linspace(-1,1,100)',linspace(-1,1,100),'v4');%将数据进行拟合并践行插值，x、y坐标，电势大小作为z值
figure(1)
contourf(X,Y,Z,20); %等值线图，等势线
figure(2)
surf(X,Y,Z);%三维曲面
figure(3)
scatter(COORD(:,1),COORD(:,2),5,Q4);%散点图
fclose(FP1);%关闭文件

求得各节点的电势，并画出相应的图。

其图像为：

得到的节点电势的值与用EIDORS软件得到的电势值误差不大，说明计算有效。此处激励采用的是1.2电极激励时，各节点的电势大小。

第三步程序中相关函数，可进行下载。

代码

EIT正问题求解--利用有限元求解节点电势相关推荐

利用python求解节点介数和边介数
利用python求解节点介数和边介数利用networkx里面的函数betweenness_centrality(G)来求解节点介数和函数edge_betweenness_centrality(G)来 ...
python有限元传热求解_二维稳态热传导基本方程的有限元求解(2)
四节点矩形单元在二维稳态热传导基本方程的有限元求解(1)这篇文章中,我们仅仅给出了有限元单元方程的一种比较标准的推导步骤,并未涉及某种具体的单元.且在式(20)中,单元上温度的近似函数表示成节点 ...
利用python求解规划问题
规划问题分为两个大类:线性规划和非线性规划以及下面分支的小类,我们观看这个树状图来粗略的了解一下. 首先我们讲解最简单的线性规划模型,通常线性规划均属于凸规划,通常都是用python中的cvxpy进行 ...
利用树求解算术表达式的值
利用树求解算术表达式的值 #include <stdio.h> #include <string.h> #include <malloc.h> //#include ...
matlab 读取图片后分区域编号_你的第一个有限元求解器——仅十行MATLAB代码
有限元分析话题中有不少讨论有限元求解器的问题,但大都停留在概念层面,未见实际代码.望本文能略起抛砖引玉之作用. 以下代码是基于MATLAB编写. 问题描述考虑一平面有界区域 ,设其边界为 .我们求解 ...
四阶龙格库塔法的基本思想_利用龙格库塔法求解郎之万方程.doc
利用龙格库塔法求解郎之万方程.doc 利用龙格-库塔法求解朗之万方程1. 待解问题布朗颗粒是非常微小的宏观颗粒,其直径的典型大小为10-710-6m.颗粒不断受到液体介质分子的碰撞,在任一瞬间,一个颗 ...
【翻译】利用加速度求解位置的算法——三轴传感器
cposture 一个小白的技术成长之路 [翻译]利用加速度求解位置的算法--三轴传感器 http://www.cnblogs.com/cposture/p/4378922.html 摘要 ...
利用有限元数值模拟技术辅助静电场学习
随着科学研究和工程技术中对于静电场的应用以及人们对于静电应用和静电现象研究的日益深入,近年来对于各类静电场电位和电场强度的分布的了解也是越来越迫切[1]:数值传热学但静电场作为电磁学的重要组成部分,与 ...
利用加速度求解位置的算法——三轴传感器
转载的一篇文章,跟自己做过的一个车载项目类似,也算是标记一下吧. ---------------------------------------分割线------------------------- ...

EIT正问题求解--利用有限元求解节点电势

EIT正问题求解--利用有限元求解节点电势相关推荐

最新文章

热门文章