计算实习报告.docx
- 文档编号:9358703
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:11
- 大小:604.94KB
计算实习报告.docx
《计算实习报告.docx》由会员分享,可在线阅读,更多相关《计算实习报告.docx(11页珍藏版)》请在冰豆网上搜索。
计算实习报告
计算实习报告
部门:
xxx
时间:
xxx
整理范文,仅供参考,可下载自行编辑
一、算例
编写程序来实现本算例,采用有限差分法或有限元法,对圆柱绕流问题进行求解。
图一:
计算实例
二、解决方案
基于Galerkin加权余量法的有限元方法求解,同时借助AutoCAD编制网格和Matlab编程计算。
三、有限元分析
可知此算例所示流动满足由流函数
表示的Laplace方程:
式中:
边界上为强制<本质)边界条件,函数值为给定值
;
边界上为自然边界条件,变量的法向导数值由
确定。
近似解:
,将近似解代入上式,得到余量
由迦辽金加权余量法:
取权函数为插值函数,则有:
其弱形式:
式中:
为边界,
为在边界
上的边界条件。
故可写成:
式中:
式中:
为
;k为单元特征矩阵或刚度矩阵,f为单元荷载列阵。
三角形单元的局部插值函数
为:
式中:
,|D|等于三角形单元的面积A的两倍。
所以可得
四、边界条件处理
由于本算例为对称结构,故可取1/4结构分析。
其边界条件表示如下:
图二:
算例边界条件
<1)第一类边界条件<强制边界条件)
由图可知,用第一类边界条件表示的边界有bc边:
。
ae、ed边:
。
对于第一类边界条件,通过改写总体系数矩阵和总体荷载矩阵得到满足,稍后会在程序中说明。
<2)第二类边界条件<自然边界条件)
由图可知,ab、cd边用第二类边界条件表示:
。
对于边界上,
是给定的,取沿边界的坐标系s-n,它与x-y坐标系的关系如下,
其中l为边界2-3的长度,将上式代入得到:
在边界上的法向坐标为n=0,使用Lagrange插值多项式条件,得到边界上的插值函数为:
,
,
最终的到结果为:
代入此算例的第二类边界条件,可知对应的所需改写的总体荷载矩阵的部分为0。
第二类边界不需要改写总体系数矩阵,所以此算例中的第二类边界条件已经自动满足。
b5E2RGbCAP
五、网格划分
此算例采用AutoCAD手动划分网格,在突变处加密。
共298个节点,516个单元。
图三:
算例网格图
六、程序编写
在编写主程序之前,先将AutoCAD中各节点坐标按编号顺序输出,导入Matlab中的M文件,为
p1EanqFDPw
编写各单元的单元定位向量,为
编写含强制边界的节点信息,为
<由于以上三个文件内容较多,故不在此处列出,请详见相应文件)
以下为求解主程序
clearall
K=zeros(298,298>。
%总体系数矩阵
F=zeros(298,1>。
%总体荷载矩阵
ZB=dlmread('CoordinateMatrix.m'>-ones(298,2>。
%读入节点坐标。
由于编制网格时左小角坐标为<1,1),而不是<0,0)。
为了流线图显示程序,因此将每个坐标减1.DXDiTa9E3d
DW=dlmread('VectorMatrix.m'>。
%读入单元定位向量
QZ=dlmread('ForcedBoundaryConditions.m'>。
%读入强制边界条件RTCrpUDGiT
fori=1:
516
S(i>=abs(1/2*(ZB(DW(i,1>,1>*ZB(DW(i,2>,2>+ZB(DW(i,2>,1>*ZB(DW(i,3>,2>+ZB(DW(i,3>,1>*ZB(DW(i,1>,2>-ZB(DW(i,3>,1>*ZB(DW(i,2>,2>-ZB(DW(i,2>,1>*ZB(DW(i,1>,2>-ZB(DW(i,1>,1>*ZB(DW(i,3>,2>>>。
5PCzVD7HxA
%S(i>为第i个单元的面积
b=zeros(1,3>。
c=zeros(1,3>。
b(1>=(ZB(DW(i,2>,2>-ZB(DW(i,3>,2>>/(2*S(i>>。
c(1>=(ZB(DW(i,3>,1>-ZB(DW(i,2>,1>>/(2*S(i>>。
b(2>=(ZB(DW(i,3>,2>-ZB(DW(i,1>,2>>/(2*S(i>>。
c(2>=(ZB(DW(i,1>,1>-ZB(DW(i,3>,1>>/(2*S(i>>。
b(3>=(ZB(DW(i,1>,2>-ZB(DW(i,2>,2>>/(2*S(i>>。
c(3>=(ZB(DW(i,2>,1>-ZB(DW(i,1>,1>>/(2*S(i>>。
%为每个单元所对应的系数
forj=1:
3
fork=1:
3
K(DW(i,j>,DW(i,k>>=K(DW(i,j>,DW(i,k>>+(b(j>*b(k>+c(j>*c(k>>*S(i>。
jLBHrnAILg
end
end
end
%以下为代入强制边界条件,并改写总体系数矩阵和总体荷载矩阵。
自然边界条件已自行满足
fori=1:
64%64为需要代入强制边界条件的节点个数,即length(QZ>
forj=1:
298
ifQZ(i,1>==j
K(j,j>=1。
F(j>=QZ(i,2>。
continue。
end
K(QZ(i,1>,j>=0。
end
end
Psi=K\F。
%Psi为各结点的流函数值
x=ZB(:
1>。
y=ZB(:
2>。
xx=linspace(min(x>,max(x>,100*(max(x>-min(x>>>。
yy=linspace(min(y>,max(y>,100*(max(y>-min(y>>>。
[X,Y]=meshgrid(xx,yy>。
%作图范围数组
Z=griddata(x,y,Psi,X,Y,'v4'>。
%'v4'为自动处理边界的命令
v=[00.20.40.60.811.21.41.61.82.0]'。
contour(X,Y,Z,v,'k'>%自动插值,画等值线,即流线
axisequaltight%使横纵坐标比例尺一致
运行程序得流线图,如下:
图四:
1/4流函数图
利用所得节点流函数值,做最窄断面流速分布图:
图五:
最窄断面流速分布图
七、总结
1、本次算例采用了有限元法,而没有采用有限差分法。
有限元法在编程上虽然较有限差分法较复杂,但是在边界问题的处理上,有限元法明显优于有限差分法。
而且有限元法精度较高。
故本算例用有限元法计算。
xHAQX74J0X
2、编写程序输出流线图时,采用了Matlab自带的contour命令,自动在求得的节点流函数值中插值得流线图。
同时,使用了命令‘V4’,此命令可自动处理不连续边界,使流线连续。
但是加入此命令后,程序运行时间明显增长。
LDAYtRyKfE
3、给定边界条件时,ab边也可给出第一类边界条件,但是综合各方面考虑后,最后以第二类边界条件给出。
4、本次计算深化了对有限元方法原理的理解,且所得结果比较理想。
取得了预期的效果。
申明:
所有资料为本人收集整理,仅限个人学习使用,勿做商业用途。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算 实习 报告