|
阅读：2220回复：1
三维有限元网格加密方法及其工程应用
<P align=center><FONT face=仿宋_GB2312>强天驰 寇晓东 周维垣</FONT></P>
<P align=center><EM><FONT face=宋体>清华大学水利水电工程系</FONT></EM></P> <P><FONT face=幼圆><SMALL><STRONG>摘 要：</STRONG>本文利用最小势能原理将三维有限元的疏密网格通过界面单元耦合起来，解决了网格突变情况下的位移协调问题，简化了有限元网格的局部加密问题，从而协调了计算精度与效率之间的矛盾.通过算例把本文推荐的方法与精确解和普通有限元解作了对比，验证了方法的合理性.</SMALL></FONT></P> <P><FONT face=幼圆><SMALL><STRONG>关键词：</STRONG>三维有限元，位移协调，网格加密，界面单元，最小势能原理.</SMALL></FONT></P> <P><EM><FONT face=宋体><SMALL>基金项目：国家自然科学资金资助项目(59939190)</SMALL></FONT></EM></P> <P><EM><SMALL><FONT face=宋体>作者简介：强天驰(1972-)，陕西西安人，清华大学水电系博士生，主要从事岩石力学与工 程方面的研究.</FONT></SMALL></EM></P> <P ><FONT face=宋体> 在有限元计算中，人们常常对某一局部区域的位移和应力更感兴趣，或是希望在此范围内能获得更高的精度.同时在对拱坝等结构物进行开裂分析时，人们也希望在裂尖周围能有一套比较致密的网格，这样起裂后就不需要重新剖分网格.但相邻节点密度相差太大而不能保证在其交界面上的位移和应变协调时，则会引起较大的计算误差.目前的解决方法大致有3种：(1)用高次单元与线性单元结合实现疏密网格的过渡，使各个节点满足协调关系；(2)人为的调整有限元基本方程的右端项，通过迭代计算来消除这一影响；(3)通过调整刚度矩阵，在刚度里体现位移不协调而产生的附加能量，然后进行整体求解.在二维情况下，第1种方法是一种十分可行的，但三维情况则会变的非常困难，且由于二次单元与线性单元的形函数分别为二次和一次，使得交界面上仅能保证节点位移协调而非交界面上各点处处位移协调，因此交界面上产生的附加势能导致计算误差在所难免(见算例一).第2种方法计算步骤复杂，而且界面的虚拟刚度如何确定也是一个非常敏感的问题，计算的收敛性不易解决.针对第3种方法，作者提出用界面过渡单元来耦合疏密网格，使其在界面上达到位移协调，以减少计算误差并提高有限元数据前处理的工作效率.</FONT></P> <P ><STRONG><FONT face=宋体>1 界面过渡单元方法基本原理</FONT></STRONG></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="37%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0701.gif"><FONT face=宋体><BR><SMALL>图1 有限元与界面单元的耦合</SMALL></FONT></P></TD> <TD width="63%"> <P ><FONT face=宋体> 设两部分计算域单元密度不同，将两部分边界用界面过渡单元来模拟，两部分计算域内部各自采用疏密不同的有限元网格，计算域通过界面过渡单元耦合起来，用界面过渡单元把两部分计算域的边界刚度耦合到总刚中，那么过渡区就是一个界面过渡单元，两边的单元密度可以随意调整，见图1.</FONT></P> <P ><FONT face=宋体> 各计算域对整体平衡方程的贡献与有限元的讨论相同，不再赘述，这里着重讨论通过界面的耦合原理.</FONT></P></TD></TR></TABLE> <P ><FONT face=宋体> 如图1所示，计算域1，2通过界面元耦合起来，界面元上任一点s对应于1，2块体的两个点，因此s点的变形为1，2块体在相应点变形的差值，以Δu(s)，Δv(s)，Δw(s)表示界面元上s点沿x，y，z 3个方向的变形，上标(1)，(2)表示相应变量在计算域1，2内的取值，则</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="87%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0701.jpg"></P></TD> <TD width="13%"> <P align=right>(1)</P></TD></TR></TABLE> <P ><FONT face=宋体>由有限元的基本理论</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="86%"> <P align=center>d<SUP>(i)</SUP>(s)=N<SUP>(i)</SUP>(s)u,i=1,2</P></TD> <TD width="14%"> <P align=right>(2)</P></TD></TR></TABLE> <P ><FONT face=宋体>其中</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="86%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0703.jpg"></P></TD> <TD width="14%"> <P align=right>(3)</P></TD></TR></TABLE> <P ><FONT face=宋体> u<SUB>i</SUB>=(u<SUB>i</SUB> v<SUB>i</SUB> w<SUB>i</SUB>)T,i=1,2,…n,为s点所在有限单元的第i个节点的位移，d<SUP>(i)</SUP>(s)为第i个块体上s点的位移，N<SUP>(i)</SUP>(s)为有限单元形函数的矩阵，即</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="86%"> <P align=center>N<SUP><SUP>(i)</SUP></SUP>(s),=[N<FONT size=2><SUB>1</SUB><SUP><SUP>(i)</SUP></SUP></FONT>(s),N<FONT size=2><SUB>2</SUB><SUP><SUP>(i)</SUP></SUP></FONT>(s),…,N<FONT size=2><SUB>n</SUB><SUP><SUP>(i)</SUP></SUP></FONT>(s)]</P></TD> <TD width="14%"> <P align=right>(4)</P></TD></TR></TABLE> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="85%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0705.jpg"></P></TD> <TD width="15%"> <P align=right>(5)</P></TD></TR></TABLE> <P ><FONT face=宋体>n表示该单元为n节点单元.</FONT></P> <P ><FONT face=宋体> 该界面过渡单元的作用相当于把两个块体在界面上任意对应的两个点用界面虚拟弹簧固定在一起，如果两个点之间有相对变形则用虚拟弹簧保证协调，当界面虚拟弹簧刚度足够大时这两个点相互之间变形足以忽略以保证位移协调.这相当于s点的X，Y，Z方向有相对变形Δu(s)，Δv(s)，Δw(s)时，用系数分别为：k<SUB>x</SUB>,k<SUB>y</SUB>,k<SUB>z</SUB>的虚拟弹簧将X,Y,Z三个方向分别约束，其系数矩阵为</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="87%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0706.jpg"></P></TD> <TD width="13%"> <P align=right>(6)</P></TD></TR></TABLE> <P ><FONT face=宋体>当取三个方向的虚拟刚度系数为相同数值时，即可将弹簧系数矩阵简化为球形张量，</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="86%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0707.jpg"></P></TD> <TD width="14%"> <P align=right>(7)</P></TD></TR></TABLE> <P ><FONT face=宋体>因此界面产生的新的位能为</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="85%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0708.jpg"></P></TD> <TD width="15%"> <P align=right>(8)</P></TD></TR></TABLE> <P ><FONT face=宋体>利用式(1)，(2)得</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="84%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0709.jpg"></P></TD> <TD width="16%"> <P align=right>(9)</P></TD></TR></TABLE> <P ><FONT face=宋体>其中P见式(7).</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="83%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0710.jpg"></P></TD> <TD width="17%"> <P align=right>(10)</P></TD></TR></TABLE> <P ><FONT face=宋体>故由式(10)可得界面过渡单元对有限元基本方程刚度矩阵的贡献为</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="84%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0711a.jpg"></P></TD> <TD width="16%"> <P align=right>(11a)</P></TD></TR></TABLE> <P ><FONT face=宋体>界面过渡单元对有限元基本方程右端项力矩阵的贡献为</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="84%"> <P align=center>f=(0,0,…,0)<SUP>T</SUP></P></TD> <TD width="16%"> <P align=right>(11b)</P></TD></TR></TABLE> <P ><FONT face=宋体>分解(11a)式有</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="83%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0712a.jpg"></P></TD> <TD width="17%"> <P align=right>(12a)</P></TD></TR></TABLE> <P ><FONT face=宋体>n<SUB>1</SUB>为大单元的界面节点数</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="85%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0712b.jpg"></P></TD> <TD width="15%"> <P align=right>(12b)</P></TD></TR></TABLE> <P ><FONT face=宋体>n为小单元的界面节点数，N<SUB>j</SUB><SUP><SUP>(2)</SUP></SUP>(s)为小单元上的形函数.</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="86%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0712c.jpg"></P></TD> <TD width="14%"> <P align=right>(12c)</P></TD></TR></TABLE> <P ><FONT face=宋体>n<SUB>1</SUB>为大单元的界面节点数，n为小单元的界面节点数.其中N<SUB>i</SUB><SUP><SUP>(1)</SUP></SUP>(s)为s点在大单元上的形函数，N<SUB>j</SUB><SUP><SUP>(2)</SUP></SUP>(s)为s点所在单元上的形函数.</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="72%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0712d.jpg"></P></TD> <TD width="15%"> <P align=right>(12d)</P></TD></TR></TABLE> <P ><FONT face=宋体>该式与上式在总刚矩阵中为对称，n<SUB>1</SUB>,n的意义同上式</FONT></P> <P ><FONT face=宋体> 式(12)即为界面元对整体平衡方程的贡献，利用二维二阶高斯数值积分方法将式(12)集入整体平衡方程，从而完成了两计算域间的耦合.式(12a)与式(12b)在同一套局部坐标下进行二维二阶高斯积分；式(12c)与式(12d)中的各个形函数不在同一套局部坐标下的，须进行坐标转换，然后在同一套坐标下利用二维二阶高斯积分求出其对整体刚度的贡献.</FONT></P> <P ><FONT face=宋体> 以上4个公式可统一写成以下形式，</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="85%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0712e.jpg"></P></TD> <TD width="15%"> <P align=right>(12e)</P></TD></TR></TABLE> <P ><FONT face=宋体>由界面和有限单元各自的刚度集成的整体刚度矩阵K具有对称正定性，可以用通常的求解系数矩阵为对称正定的方程组的方法求解，而且在存储方式上可以利用其对称性来节省存储空间.</FONT></P> <P ><STRONG><FONT face=宋体>2 界面虚拟弹簧刚度系数</FONT></STRONG></P> <P ><FONT face=宋体> 虚拟弹簧刚度系数k的选择非常重要，k值需足够大，以保证界面的变形足够小.但如果k取值太大，则会造成整体刚度矩阵的病态，使解的误差超出允许范围.因此，k必须取合适的值.</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="68%"> <P ><FONT face=宋体> 为了估计k值，需求出弹簧刚度系数k与材料弹性模量E，泊松比μ之间的关系，考虑如图2的一个尺寸为a×b的矩形块体，其上边界受均布面力p，下边界由均匀弹簧垫支撑，弹簧垫刚度为k，为保守起见，块体两侧受x向零位移约束，并考虑平面应变情况.</FONT></P> <P ><FONT face=宋体><SMALL></SMALL> 由平面应变本构关系及已知条件ε<SUB>x</SUB>=0,得</FONT></P></TD> <TD width="32%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0702.gif"><FONT face=宋体><BR><SMALL>图2 虚拟弹簧刚度的估算</SMALL></FONT></P></TD></TR></TABLE> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="100%"> <P align=center>ε<SUB>x</SUB>=1/E[（1-μ<SUP><SUP><FONT size=2>2</FONT></SUP></SUP>）σ<SUB>x</SUB>-μ（1+μ）σ<SUB>y</SUB>]=0</P> <P align=center>ε<SUB>y</SUB>=1/E[-μ（1+μ）σ<SUB>x</SUB>+（1-μ<SUP><SUP><FONT size=2>2</FONT></SUP></SUP>）σ<SUB>y</SUB>]</P></TD></TR></TABLE> <P ><FONT face=宋体>由σ<SUB>y</SUB>=p，利用上面两式可解得</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="100%"> <P align=center>ε<SUB>y</SUB>=[(1+μ)(1-2μ)]/[E(1-μ)]p</P></TD></TR></TABLE> <P ><FONT face=宋体>该d<SUB>b</SUB>为矩形块体的y向变形量，d<SUB>s</SUB>为弹簧的y向变形量.则</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="100%"> <P align=center>d<SUB>b</SUB>=b·ε<SUB>y</SUB>=[b(1+μ)(1-2μ)/[E(1-μ)]p</P> <P align=center>d<SUB>s</SUB>=p/k</P></TD></TR></TABLE> <P ><FONT face=宋体>要想正确模拟约束位移，就要求：d<SUB>b</SUB>>>d<SUB>s</SUB>，</FONT></P> <P ><FONT face=宋体>设：d<SUB>b</SUB>=cd<SUB>s</SUB>，则得</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="86%"> <P align=center>k=[c(1-μ)]/[b(1+μ)(1-2μ)]E</P></TD> <TD width="14%"> <P align=right>(13)</P></TD></TR></TABLE> <P ><FONT face=宋体>其中c为一常数，b为一特征长度，可取为</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="85%"> <P align=center>c=10～100,b=min(d<SUB>n</SUB>)</P></TD> <TD width="15%"> <P align=right>(14)</P></TD></TR></TABLE> <P ><FONT face=宋体>d<SUB>n</SUB>为节点间距.</FONT></P> <P ><STRONG><FONT face=宋体>3 虚拟弹簧刚度系数大小对计算结果的影响</FONT></STRONG></P> <P ><FONT face=宋体> 由于虚拟弹簧的刚度系数为式(13)</FONT>k=[c(1-μ)]/[b(1+μ)(1-2μ)]E<FONT face=宋体>为同一个数量级，如果设b=1m，则弹簧刚度将与cE同一个数量级. 因此，值越大，则界面过渡单元的刚度值越大，这样刚度矩阵就会出现奇异.如果不能精确保证其正定性，则计算结果就可能出现较大的误差.当取c=100时，在总刚中与界面有关的自由度相应的刚度系数上比其他自由度相应的刚度系数要大2个数量级.由于计算界面刚度时(12c)式和(12d)式两套局部坐标必须进行坐标转换，因此会有一定的数值计算误差，当这一计算误差为1%时，相当于界面刚度系数计算误差为1%，但相对于总刚其他自由度相应的刚度系数，则误差就要增大2个数量级，达到100%，因此必须保证坐标转换计算误差足够小.若降低c值，如取c=10，则界面刚度计算误差也可放宽一些.三维情况下疏密单元之间的坐标转换是一个比较复杂的问题，本文作者采用数值方法计算局部坐标转换，通过较少的计算时间，达到较高的计算精度，因此，最终结果的计算误差可以忽略不计.限于篇幅，详细的数值方法本文这里不作介绍.</FONT></P> <P ><STRONG><FONT face=宋体>4 算例及结论</FONT></STRONG></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="100%"> <P align=center><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0703.gif"><FONT face=宋体><BR><SMALL>图3 计算试件</SMALL></FONT></P></TD></TR></TABLE> <P ><FONT face=宋体> 作者根据以上理论编制了三维有限元程序，对如图3所示的试件进行计算.图3(a)为采用二次单元与线性单元过渡的网格其中2号单元为二次单元；图3(b)中所有单元都是线性单元，在疏密网格交界处节点不协调；图3(c)为采用界面过渡单元的网格，所有单元都是线性单元，但疏密网格间由界面过渡单元耦合.该试件上部受Z向均布荷载1MPa，取c=100.试件的材料参数为弹性模量E=2.0×10<SUP>4</SUP>MPa，泊松比μ=0.167.</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD align=middle width="100%" colSpan=5> <P align=center><STRONG>表1 各节点Z向（垂直）位移结果比较（mm）</STRONG></P></TD></TR> <TR> <TD align=middle width="100%" colSpan=5> <HR color=#000000> </TD></TR> <TR> <TD align=middle width="20%">Z向变形</TD> <TD align=middle width="20%">3(a)协调</TD> <TD align=middle width="20%">3(b)不协调</TD> <TD align=middle width="20%">3(c)界面耦合</TD> <TD align=middle width="20%">理论解</TD></TR> <TR> <TD align=middle width="100%" colSpan=5> <HR color=#000000 SIZE=1> </TD></TR> <TR> <TD align=middle width="20%">1</TD> <TD align=middle width="20%">3.1290</TD> <TD align=middle width="20%">3.9719</TD> <TD align=middle width="20%">2.9942</TD> <TD align=middle width="20%">3.0</TD></TR> <TR> <TD align=middle width="20%">2</TD> <TD align=middle width="20%">1.9716</TD> <TD align=middle width="20%">2.9725</TD> <TD align=middle width="20%">1.9942</TD> <TD align=middle width="20%">2.0</TD></TR> <TR> <TD align=middle width="20%">5</TD> <TD align=middle width="20%">1.3059</TD> <TD align=middle width="20%">2.4209</TD> <TD align=middle width="20%">1.4939</TD> <TD align=middle width="20%">1.5</TD></TR> <TR> <TD align=middle width="20%">11</TD> <TD align=middle width="20%"></TD> <TD align=middle width="20%"></TD> <TD align=middle width="20%">1.4873</TD> <TD align=middle width="20%">1.4995</TD></TR> <TR> <TD align=middle width="20%">7</TD> <TD align=middle width="20%">1.3973</TD> <TD align=middle width="20%">0.4958</TD> <TD align=middle width="20%">1.2484</TD> <TD align=middle width="20%">1.25</TD></TR> <TR> <TD align=middle width="20%">6</TD> <TD align=middle width="20%">0.9076</TD> <TD align=middle width="20%">1.5665</TD> <TD align=middle width="20%">1.0094</TD> <TD align=middle width="20%">1.0</TD></TR> <TR> <TD align=middle width="100%" colSpan=5> <HR color=#000000> </TD></TR></TABLE> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="100%"></TD></TR></TABLE> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="100%" colSpan=6> <P align=center><STRONG>表2 应力σ<SUB>z</SUB>计算结果比较(MPa)</STRONG></P></TD></TR> <TR> <TD width="100%" colSpan=6> <HR color=#000000> </TD></TR> <TR> <TD align=middle width="16%">单元号</TD> <TD align=middle width="16%">节点号</TD> <TD align=middle width="20%">3(a)协调</TD> <TD align=middle width="20%">3(b)不协调</TD> <TD align=middle width="20%">3(c)界面耦合</TD> <TD align=middle width="20%">理论解</TD></TR> <TR> <TD width="100%" colSpan=6> <HR color=#000000 SIZE=1> </TD></TR> <TR> <TD align=middle width="16%">单元1</TD> <TD align=middle width="16%">1</TD> <TD align=middle width="17%">-0.9963</TD> <TD align=middle width="17%">-0.9994</TD> <TD align=middle width="17%">-1.0000</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">2</TD> <TD align=middle width="17%">-1.0038</TD> <TD align=middle width="17%">-1.0005</TD> <TD align=middle width="17%">-1.0000</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">3</TD> <TD align=middle width="17%">-1.0034</TD> <TD align=middle width="17%">-1.0009</TD> <TD align=middle width="17%">-1.0000</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">4</TD> <TD align=middle width="17%">-0.9965</TD> <TD align=middle width="17%">-0.9993</TD> <TD align=middle width="17%">-1.0000</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%">单元2</TD> <TD align=middle width="16%">2</TD> <TD align=middle width="17%">-0.9802</TD> <TD align=middle width="17%">-1.0471</TD> <TD align=middle width="17%">-1.0003</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">5</TD> <TD align=middle width="17%">-0.9726</TD> <TD align=middle width="17%">-1.0127</TD> <TD align=middle width="17%">-1.0002</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">6(12)</TD> <TD align=middle width="17%">-1.0209</TD> <TD align=middle width="17%">-0.9496</TD> <TD align=middle width="17%">0.9996</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">3</TD> <TD align=middle width="17%">-1.0290</TD> <TD align=middle width="17%">-0.9712</TD> <TD align=middle width="17%">-0.9998</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%">单元3</TD> <TD align=middle width="16%">5(11)</TD> <TD align=middle width="17%">-0.9106</TD> <TD align=middle width="17%">-1.3970</TD> <TD align=middle width="17%">-0.9926</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">8</TD> <TD align=middle width="17%">-0.9033</TD> <TD align=middle width="17%">-1.4554</TD> <TD align=middle width="17%">-0.9925</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">9</TD> <TD align=middle width="17%">-1.0548</TD> <TD align=middle width="17%">-0.7091</TD> <TD align=middle width="17%">-0.9968</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">7</TD> <TD align=middle width="17%">-1.0621</TD> <TD align=middle width="17%">-0.6469</TD> <TD align=middle width="17%">-0.9969</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%">单元4</TD> <TD align=middle width="16%">7</TD> <TD align=middle width="17%">-1.0823</TD> <TD align=middle width="17%">-0.5605</TD> <TD align=middle width="17%">-1.0006</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">9</TD> <TD align=middle width="17%">-1.0745</TD> <TD align=middle width="17%">-0.6386</TD> <TD align=middle width="17%">-1.0004</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">10</TD> <TD align=middle width="17%">-0.9465</TD> <TD align=middle width="17%">-1.3475</TD> <TD align=middle width="17%">-1.0068</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD align=middle width="16%"></TD> <TD align=middle width="16%">6</TD> <TD align=middle width="17%">-0.9559</TD> <TD align=middle width="17%">-1.2547</TD> <TD align=middle width="17%">-1.0071</TD> <TD align=middle width="17%">-1.0</TD></TR> <TR> <TD width="100%" colSpan=6> <HR color=#000000> </TD></TR></TABLE> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD width="100%"> <P ><FONT face=宋体>注：对于3(c)界面耦合情况单元3的节点5换为节点11，单元2的节点6换为节点12.</FONT></P></TD></TR></TABLE> <P ><FONT face=宋体> 由计算结果可以看出，当节点不协调时，计算误差无法接受，用二次单元过渡时，计算误差较小，但最大计算误差约为13%，这是由于在疏密单元的界面上产生了多余的附加势能，当采用界面过渡单元时，计算误差在1%以内，而且还可以视问题的要求改变参数以达到更高的计算精度.可见界面过渡单元较好的消除了由于疏密网格界面位移不协调所产生的误差.</FONT></P> <P ><FONT face=宋体> 作者用该三维有限元程序对某井壁开裂问题进行了开裂计算，计算区域如图4所示，用界面 协调方法对井壁周围区域进行加密，利用界面过渡单元耦合疏密两部分单元，如图5所示，对图5所示单元进行开裂追踪，裂纹扩展过程见图6(a)、(b)所示.</FONT></P> <P ><FONT face=宋体> 由此可见，利用界面过渡单元可以很好的解决疏密网格的衔接问题，使得有限元网格的局部加密变得易于实现，使三维有限元裂纹扩展的计算效率大大提高.</FONT></P> <TABLE cellSpacing=0 cellPadding=0 width="100%" border=0> <TR> <TD align=middle width="50%"><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0704.gif"><BR><SMALL><FONT face=宋体>图4 计算区域</FONT></SMALL></TD> <TD align=middle width="50%"><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0705.gif"><BR><SMALL><FONT face=宋体>图5 局部加密网格</FONT></SMALL></TD></TR> <TR> <TD align=middle width="100%" colSpan=2><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/0706.gif"><BR><SMALL><FONT face=宋体>图6 裂纹扩展过程</FONT></SMALL></TD></TR></TABLE> <P align=center><STRONG><SMALL><FONT face=宋体>参 考 文 献</FONT></SMALL></STRONG></P> <P><SMALL><FONT face=宋体>[1] 周维垣，杨若琼.混凝土拱坝整体三维有限元非线性分析(TFINE有限元程序系列)[R].北京：清华大学水电系，1988.</FONT></SMALL></P> <P><SMALL><FONT face=宋体>[2] 寇晓东.无单元法追踪结构开裂及拱坝稳定研究[D].北京：清华大学，1998</FONT></SMALL></P> <P><SMALL><FONT face=宋体>[3] 关治，陈景良.数值计算方法[M].北京：清华大学出版社，1990.</FONT></SMALL></P> <P><SMALL><FONT face=宋体>[4] 王勖成，邵敏.有限单元法基本原理与数值方法.[M].北京：清华大学出版社，1988.</FONT></SMALL></P> <P><SMALL><FONT face=宋体>[5] H J Melosh, A Raefsky. A simple and efficient method for introducing faul ts into finite element computations [Ｊ]. Bulletin of the SSA,1981,71:1391～1399.</FONT></SMALL></P> |
|
|
1楼#
发布于：2005-07-31 10:54
会不会太专业了？<img src="images/post/smile/dvbbs/em01.gif" /><img src="images/post/smile/dvbbs/em02.gif" /><img src="images/post/smile/dvbbs/em08.gif" />
|
|
