计算实习报告.docx

上传人:b****8 文档编号:9358703 上传时间:2023-02-04 格式:DOCX 页数:11 大小:604.94KB
下载 相关 举报
计算实习报告.docx_第1页
第1页 / 共11页
计算实习报告.docx_第2页
第2页 / 共11页
计算实习报告.docx_第3页
第3页 / 共11页
计算实习报告.docx_第4页
第4页 / 共11页
计算实习报告.docx_第5页
第5页 / 共11页
点击查看更多>>
下载资源
资源描述

计算实习报告.docx

《计算实习报告.docx》由会员分享,可在线阅读,更多相关《计算实习报告.docx(11页珍藏版)》请在冰豆网上搜索。

计算实习报告.docx

计算实习报告

计算实习报告

部门:

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、本次计算深化了对有限元方法原理的理解,且所得结果比较理想。

取得了预期的效果。

申明:

所有资料为本人收集整理,仅限个人学习使用,勿做商业用途。

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 工程科技 > 机械仪表

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1