粘弹性人工边界在ANSYS中的实现
粘弹性人工边界的有限元分析
粘弹性人工边界的有限元分析徐浩【摘要】采用数值方法集中比较了工程中广泛应用的几种人工边界在SSI体系分析中的动力反应,研究结果表明:粘性边界能较好的模拟土的边界但计算的位移会发生整体的飘逸,精度也不如粘弹性边界;粘弹性边界能很好的模拟土一结构相互作用体系中土的边界问题,但实现过程比较麻烦.【期刊名称】《山西建筑》【年(卷),期】2011(037)009【总页数】2页(P72-73)【关键词】人工边界;SSI体系;粘弹性人工边界;有限元法【作者】徐浩【作者单位】同济大学土木工程学院建筑工程系,上海,200092【正文语种】中文【中图分类】TU411土木工程中有许多涉及弹性半无限地基的振动及波动问题,诸如土—结构动力相互作用(SSI)问题、地震波的传播问题、动力机器基础的振动问题、打桩及道路交通引起的振动问题等。
对这类问题用有限元法求解与用其他计算方法相比具有可灵活地适用于地基物性的不均匀性并可考虑地基的非线性特性等优点。
但由于有限元法必须对求解对象的全领域进行离散,故在分析弹性半无限地基问题时必须截取一定范围的计算模型,这就要求在切取的边界上建立人工边界,用于模拟切除的无限域影响。
目前人工边界主要分为两类：一类是全局人工边界条件,如边界元法等;另一类是局部人工边界条件,如旁轴近似人工边界、透射人工边界、粘弹性人工边界等。
局部人工边界具有时空解耦的特点,得到了广泛的应用。
其中粘弹性人工边界是通过沿人工边界设置一系列由弹簧和阻尼器组成的简单物理元件来吸收射向人工边界的波动能量和反射波的散射,其模型简单,物理意义清晰,便于在大型通用有限元计算软件中得到实现。
ANSYS是一功能强大的有限元计算软件,其中的 ANSYS中的Combin14单元和 LSDYNA中的Combin165单元,是弹簧与阻尼器的元件,易于实现粘弹性人工边界。
本文在 ANSYS中实现了粘弹性边界并进行了验证和对比。
粘弹性人工边界作为一种应力边界条件,该应力是边界结点位移和速度的函数,一般形式写为：由式(1)可以看出,粘弹性边界相当于在边界结点每个方向施加一个一段固定的单向弹簧—阻尼元件,而且该弹簧—阻尼元件的刚度和阻尼系数仅与该边界结点在该方向该时刻的反应相关,从而通过粘性阻尼的吸能作用和弹簧的刚性恢复作用模拟无限域对广义结构的影响。
CAE-ANSYS Mechanical介绍(超弹性、粘弹性)
应用多种超弹性本构的案例：膝盖韧带分析
Cartilage(软骨) Isotropic Mooney-Rivlin
Ligament(韧带) Transversely isotropic Neo-Hookean
Meniscus(半月板)
Fung
W
2
I1
3
F2 (I4 )
Orthotropic Hyperelastic:
i 1
2i
3 i
3
第三类分子统计学本构模型，Neo-Hookean、Gent。
W
2
I1
3
5
超弹性模型 介绍
6
超弹性模型 选择
• 2阶项的 Mooney-Rivlin模型（ANSYS提供2,3,5,9参数的 模型）
W C1I1 3 C2 I2 3
W C1(I1 3) C2 (I2 3) C3 (I1 3)(I2 3)
Cauchy应力 剪切松弛模量 Green应变
t
0
2G t
de
d
d
t
I 0
K
t
d
d
d
式中σ——Cauchy应力(真实应力)
G(t)——为剪切松弛核函数
K(t)——为体积松弛核函数(由于体积几乎不变，可不考虑)
e——为应变偏量部分(剪切变形)
Δ——为应变体积部分(体积变形)
t——当前时间
τ——过去时间
粘弹性模型 静态 时温等效
• 时温等效原理 升高温度与延长时间对分子运动或高聚物的粘弹行为都
是等效的，这个等效性可以借助转换因子aT，将在某一温度 下测定的力学数据转换成另一温度下的数据.
静态下，升高温度与延长时间等效 降低温度与缩短时间等效
粘弹性人工边界在ANSYS中实现
从半空间无限域取一4X2的矩形平面结构，顶部中间一定范围内受随时间变化的均布荷载，荷载如下p(t)=t 当0< DIV>p(t)=2-t 当1<=t<=2时p(t)=0 当t>2时材料弹性模量E=2.5，泊松比0.25，密度1网格尺寸0.1X0.1，在网格边界上所有结点加法向和切向combin14号单元用以模拟粘弹性人工边界（有关理论可参考刘晶波老师的相关文章）。
combine14单元的两个结点，其中一个与实体单元相连，另一个结点固定。
网格图如图1所示时程分析的时间步长为0.02秒，共计算16秒。
计算得到四个控制点位移时程图如图2所示，控制点坐标A(0,2)、B(0,1)、C(0,0)、D(2,2).计算所用命令流如下：/PREP7L=4 !水平长度H=2 !竖起深度E=2.5 !弹性模量density=1 !密度nu=0.25 !泊松比dxyz=0.1 !网格尺寸G = E/(2.*(1.+nu)) !剪切模量alfa = E*(1-nu)/((1.+nu)*(1.-2.*nu)) !若计算平面应力，此式需要修改Cp=sqrt(alfa/density) !压缩波速Cs=sqrt(g/density) !剪切波速R=sqrt(L*L/4.+H*H/4.) !波源到边界点等效长度KbT=0.5*G/R*dxyzKbN=1.0*G/R*dxyzCbT=density*Cs*dxyzCbN=density*Cp*dxyzET, 1, plane42,,,2 !按平面应变计算et, 2, combin14, ,, 2 !切向et, 3, combin14, ,, 2 !法向r, 2, KbT, CbTr, 3, KbN, CbNMP, EX, 1, EMP, PRXY, 1, nuMP, DENS, 1, densityrectng,-L/2.,L/2,0.,Hasel, allaesize, all, dxyzmshape,0,2Dmshkey,1amesh, all!以下建立底边界法向和切向弹簧阻尼单元nsel,s,loc,y,0.*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,1,npnpnum=node((ip-1)*dxyz-L/2.,0.,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x.,y-dxyz/2,z !定义底边界法向结点以便与边界点形成法向单元type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x-dxyz/2.,y,z !定义底边界切向结点以便与边界点形成切向单元type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddo!以下建立左边界法向和切向弹簧阻尼单元nsel,s,loc,x,-L/2*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,2,np !侧边界最下面一个点按底边界上处理npnum=node(-L/2,(ip-1)*dxyz,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x-dxyz/2.,y,z !定义左边界法向结点以便与边界点形成法向单元type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x,y-dxyz/2.,z !定义左边界切向结点以便与边界点形成切向单元type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddo!以下建立右边界法向和切向弹簧阻尼单元nsel,s,loc,x,L/2*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,2,np !侧边界最下面一个点按底边界上处理npnum=node(L/2,(ip-1)*dxyz,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x+dxyz/2.,y,z !定义右边界法向结点以便与边界点形成法向单元type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x,y-dxyz/2.,z !定义右边界切向结点以便与边界点形成切向单元type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddoallsel,all/pnum,type,1/number,1eplotfinish/soluANTYPE,trans!*TRNOPT,FULLLUMPM,0btime=0.02etime=16.00dtime=0.02*DO,itime,btime,etime,dtimeTIME,itimensel,s,loc,y,H !选中需要加荷载的点nsel,r,loc,x,-L/4,L/4*if,itime,lt,1.,thenf,all,fy,1*itime*elseif,itime,ge,1.0,and,itime,le,2.0f,all,fy,1*(2-itime)*elsef,all,fy,0.0*endifallsel,allSOLVE*ENDDO另外，还用自己编写的有限元程序计算了一下这个例子，并与ANSYS得到的结果进行了比较，结果非常吻合，这里给出A点的比较结果。
ansys高级非线性分析五粘弹性
粘塑性
A. 粘塑性背景
• 与前述的率无关塑性相反, 率相关塑性与应变速率或时间有关。从 材料的观点看, 粘塑性 和蠕变是相同的。
– 通常, 对工程用途, 蠕变用于描述在恒载荷下应变的变化。通常当温度 达到材料熔点的30-60%时, 蠕变变形就很重要，而且, 时间更长，蠕 变和塑性应变解耦。
– ANSYS 中, Perzyna 和 Peirce 模型(TB,RATE) 意味着高应变率(即冲 击)载荷状态，非弹性应变是不解耦的。 – Anand 模型与 ANSYS 蠕变法则类似, 主要区别是使用了一个内部变 量，变形抗力, 来表示对非弹性材料流动的各向同性抗力。
1 o g in
n
Peirce 模型 n 1 o in g
n
g 1 in o in g
Strain Rates Time scale Temperature Effects
Supported Element Types 3
BISO, MISO, NLISO, BKIN, HILL BISO, MISO, NLISO, HILL None Suitable for small strain rates Suitable for large strain rates Suitable for small strain rates Long periods, creep and plasticity Short periods, usually for impactShort/medium periods have different time scales type problems Temperature effects included as Can input temperature-dependent Anand's equation considers part of equation (or material material constants, but equations do temperature effects directly. No constants can be temperaturenot consider temp effects directly need to input temperature2 dependent material constants dependent) Implicit - core and 18x Core and 18x VISCO106-108 Explicit - core and misc
ansys高级非线性分析五粘弹性
& σ = SR (εin )σo
等效有效应力 应变率强化 [应力比率] (TB,RATE) 静态屈服应力 (TB,BISO/MISO/NLISO)
3
应力比率
ε &in m & SR (εin ) = 1+ γ & &in ) = 1+ εin SR (ε γ
m
PERZYNA
Perzyna 模型
2
1
PEIRCE
Peirce 模型
0 0 1 2 3 4 5
正则化应变率
September 30, 2001 Inventory #001491 5-8
粘塑性
... RATE 粘塑性选项
• 这意味着应力比率作为应变率的函数而改变。 这意味着应力比率作为应变率的函数而改变。
用 Perzyna 模型的例子
September 30, 2001 Inventory #001491 5-9
粘塑性
... RATE 粘塑性选项
• 关于 关于RATE模型的一些假设 模型的一些假设: 模型的一些假设
Training Manual
Advanced Structural Nonlinearities 6.0
September 30, 2001 Inventory #001491 5-3
粘塑性
... 粘塑性理论的背景
• ANSYS中率相关塑性模型的总结 中率相关塑性模型的总结: 中率相关塑性模型的总结
CREEP Behavior Isotropic or anisotropic creep (see HILL below) No explicitly defined yield surface RATE 1 Strain rate- or time-dependent Isotropic or anisotropic viscoplasticity (see HILL below) Includes yield surface ANAND Isotropic
粘弹性人工边界及地震动输入在通用有限元软件中的实现
防灾减灾工程学报第２７卷簧一阻尼器系统，图１为三维粘弹性人工边界示意图，人工边界上法向与切向的弹簧刚度和阻尼系数按照公式（１）和（２）取值‘８’９］。
，结构图１人工边界等效弹簧一阻尼器系统Ｋ删＝口Ⅳ号，ＣＢ．Ⅳ＝ｐｃｐ（１）广Ｋ盯＝奸百ｋＪ＂，Ｃ盯＝ｐｃ，（２）式中ＫｓＮ、Ｋ胛分别为法向与切向弹簧刚度；ＧⅣ、Ｃ胛分别为法向与切向阻尼器的阻尼系数；Ｒ为波源至人工边界点的距离；ｏ和以分别为介质的Ｐ波和Ｓ波波速；Ｇ为介质剪切模量；Ｐ为介质质量密度；ａＮ与即分别为法向与切向粘弹性人工边界修正系数。
大量数值计算表明，粘弹性人工边界具有良好的鲁棒性，人工边界参数口Ⅳ与奸在一定范围内取值均可以给出良好的计算结果，经过大量算例分析，推荐使用表１中的数据。
表ｌ粘弹性人工边界参数蛳与昕的取值范围粘弹性人工边界所模拟的是人工边界上的应力条件，因此这是１种连续分布的人工边界条件。
当采用有限元法或其它离散方法将人工边界所包围的计算区离散化时，人工边界面也将随之离散化，此时可以采用有限元法的形函数将连续分布的人工边界物理元件化为耦联的人工边界，可称为一致粘弹性人工边界；也可以简单地采用集中处理方法，形成解耦的人工边界，称为集中粘弹性人工边界。
三维粘弹性人工边界有限元离散系统如图２（ａ）所示，由三维单元边界面上的节点１、２、３、４所确定的平面上均匀分布着法向和切向粘弹性人工边界的阻尼器和弹簧。
（ａ）一致粘弹性人工边界（ｂ）一致粘弹性人工边界单元图２三维一致粘弹性人Ｔ边界及粘弹性人工边界单元示意１．２一致粘弹性人工边界单元为了更为简便地应用粘弹性人工边界，在有限元中使用等效实体单元来替换空间分布的弹簧一阻尼器元件，即在已有有限元模型的边界上沿边界面法向延伸一层一定厚度的实体单元。
并将外层边界固定，如图２（ｂ）所示‘１“。
如果用与计算区域相同的单元来模拟一致粘弹性边界，单元的等效剪切模量０和弹性模量雷为ｆｏ—ｈＫｖｒ＝的＾簧｛卺＝ｇ二≮｝铲＾Ｋ肼＝ｃ３）【蛳ＪＩｌ页Ｇ·型导％型式中Ｋ删、Ｋ朋．为人工边界弹簧法向与切向刚度；ｈ为边界单元的厚度；孝为等效泊松比，按下式取值：移一』者各口≥２移一＿２（口一１）一７。
ansys高级非线性分析五粘弹性
粘塑性
A. 粘塑性背景
Training Manual
• 与前述的率无关塑性相反, 率相关塑性与应变速率或时间有关。从
材料的观点看, 粘塑性 和蠕变是相同的。
– 通常, 对工程用途, 蠕变用于描述在恒载荷下应变的变化。通常当温度 达到材料熔点的30-60%时, 蠕变变形就很重要，而且, 时间更长，蠕 变和塑性应变解耦。
Training Manual
记住首先定义线弹性材 料(EX 和 PRXY)。
否则, 当定义粘塑性时, 材料GUI将提示该信息。
第112页/共34页
Advanced Structural Nonlinearities 6.0
粘塑性
... 定义 RATE 材料特性
• 本例中选择了双线性等向强化的各向异性粘塑性。
• 当输入RATE常数后, 会出现率无关塑性参数对话框。
– 该例中, 输入双线性等向强化模型( BISO)的值, 即屈服应力和剪切模 量。
第143页/共34页
粘塑性
... 定义 RATE 材料特性
• 最后, 需要输入HILL参数来指定各向异性。
– 下面对话框中需要输入6个应力比率。
Training Manual
have different time scales
type problems
Temperature effects included as Can input temperature-dependent Anand's equation considers
part of equation (or material material constants, but equations do temperature effects directly. No
二维黏弹性人工边界单元及地震波输入在ANSYS中的实现
水利与建筑工程学报
ＪｏｕｒｎａｌｏｆＷａｔｅｒＲｅｓｏｕｒｃｅｓａｎｄＡｒｃｈｉｔｅｃｔｕｒａｌＥｎｇｉｎｅｅｒｉｎｇ
ＤＯＩ：１０．３９６９／ｊ．ｉｓｓｎ．１６７２－１１４４．２０１９．０５．０１０
Ｖｏｌ．１７Ｎｏ．５ Ｏｃｔ．，２０１９
二维黏弹性人工边界单元及地震波输入 在 ＡＮＳＹＳ中的实现
对无限域或半无限域进行地震响应分析时，往 往采用有限元法。该方法通常采取的处理措施是在 无限域或半无限域空间中人为地截取一定范围的有 限域，进而可以将无限区域转化为有限区域，然后在 有限域的基础上，对其进行离散法处理，进而无限域
或半无限域最终被转化为有限个单元。同时为了达 到近似处理的效果，通常在截取的有限域的边界上 人为的施加约束，但是这一处理的缺点是在处理过 程中，不仅没有考虑远场无限地基的辐射阻尼效应， 同时半空间无限地基弹卷
截取的有限域的范围存在不确定性，随着有限域截 取范围的增加，自由度也随之增加。基于此离散范 围也更加宽广，进而网格划分的单元数及结点数成 倍增加，因此，在地震动的动力计算中工作量也进一 步被增加，对计算机的要求也更高，进而，在动力学 中非线性问题的求解难度进一步增大。目前，针对 该关键性问题处理的诸多措施中，较好的处理方法 是在截取的有限域边界处设置局部人工边界。目前 在工程中较常用的局部边界主要有黏性边界［１］、透 射边界［２］、黏弹性边界［３］等。黏性边界的优点是物 理概念清晰，在实际动力计算程序中实现起来相对 比较容易，对于 １阶动力响应精确度较高。基于此， 黏性边界单元被嵌入到了很多大型商用软件中以便 于工程的实际应用，如 ＬＳＤＹＮＡ、ＦＩＡＣ、ＡＢＡＱＵＳ［４］ 等。它的缺点是在施加人工边界过程中只着重考虑 了人工边界辐射阻尼效应，忽略了人工边界处半无 限空间介质的弹性恢复性能。因此在实际工程应用 中容易发生整体漂移的现象，进而使计算结果精度 不高。透射边界的优点是对于 ２阶动力响应计算结 果存在较高的精度，其缺点是动力计算分析程序的 代码编写较复杂、在实际工程应用中可能存在高频 失稳的现象。黏弹性边界不仅能有效的模拟半无限 空间介质的辐射阻尼效应，同时还能较好的模拟半 无限地基的弹性恢复性能。此外还具有良好的频率 稳定性和较高的精度。因此在实际工程中得到了广 泛的应用［５－１０］。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
粘弹性人工边界在ANSYS中的实现
(2007-11-07 00:25:58)
标签：
分类：FEM软件
知识/探索
ansys
粘弹性人工边界
动力边界条件
粘弹性人工边界在ANSYS中的实现
从半空间无限域取一4X2的矩形平面结构，顶部中间一定范围内受随时间变化的均布荷载，荷载如下
p(t)=t 当0< DIV>
p(t)=2-t 当1<=t<=2时
p(t)=0 当t>2时
材料弹性模量E=2.5，泊松比0.25，密度1
网格尺寸0.1X0.1，在网格边界上所有结点加法向和切向combin14号单元用以模拟粘弹性人工边界（有关理论可参考刘晶波老师的相关文章）。
combine14单元的两个结点，其中一个与实体单元相连，另一个结点固定。
网格图如图1所示
时程分析的时间步长为0.02秒，共计算16秒。
计算得到四个控制点位移时程图如图2所示，控制点坐标A(0,2)、B(0,1)、C(0,0)、D(2,2).
计算所用命令流如下：
/PREP7
L=4 !水平长度
H=2 !竖起深度
E=2.5 !弹性模量
density=1 !密度
nu=0.25 !泊松比
dxyz=0.1 !网格尺寸
G = E/(2.*(1.+nu)) !剪切模量
alfa = E*(1-nu)/((1.+nu)*(1.-2.*nu)) !若计算平面应力，此式需要修改Cp=sqrt(alfa/density) !压缩波速
Cs=sqrt(g/density) !剪切波速
R=sqrt(L*L/4.+H*H/4.) !波源到边界点等效长度
KbT=0.5*G/R*dxyz
KbN=1.0*G/R*dxyz
CbT=density*Cs*dxyz
CbN=density*Cp*dxyz
ET, 1, plane42,,,2 !按平面应变计算
et, 2, combin14, ,, 2 !切向
et, 3, combin14, ,, 2 !法向
r, 2, KbT, CbT
r, 3, KbN, CbN
MP, EX, 1, E
MP, PRXY, 1, nu
MP, DENS, 1, density
rectng,-L/2.,L/2,0.,H
asel, all
aesize, all, dxyz
mshape,0,2D
mshkey,1
amesh, all
!以下建立底边界法向和切向弹簧阻尼单元
nsel,s,loc,y,0.
*get,np,node,,count !得到选中的结点数，存入np
*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax *do,ip,1,np
npnum=node((ip-1)*dxyz-L/2.,0.,0.)
x=nx(npnum)
y=ny(npnum)
z=nz(npnum)
npmax=npmax+1
n,npmax,x.,y-dxyz/2,z !定义底边界法向结点以便与边界点形成法向单元
type,3
real,3
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
npmax=npmax+1
n,npmax,x-dxyz/2.,y,z !定义底边界切向结点以便与边界点形成切向单元
type,2
real,2
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
*enddo
!以下建立左边界法向和切向弹簧阻尼单元
nsel,s,loc,x,-L/2
*get,np,node,,count !得到选中的结点数，存入np
*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax *do,ip,2,np !侧边界最下面一个点按底边界上处理
npnum=node(-L/2,(ip-1)*dxyz,0.)
x=nx(npnum)
y=ny(npnum)
z=nz(npnum)
npmax=npmax+1
n,npmax,x-dxyz/2.,y,z !定义左边界法向结点以便与边界点形成法向单元
type,3
real,3
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
npmax=npmax+1
n,npmax,x,y-dxyz/2.,z !定义左边界切向结点以便与边界点形成切向单元
type,2
real,2
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
*enddo
!以下建立右边界法向和切向弹簧阻尼单元
nsel,s,loc,x,L/2
*get,np,node,,count !得到选中的结点数，存入np
*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax *do,ip,2,np !侧边界最下面一个点按底边界上处理
npnum=node(L/2,(ip-1)*dxyz,0.)
x=nx(npnum)
y=ny(npnum)
z=nz(npnum)
npmax=npmax+1
n,npmax,x+dxyz/2.,y,z !定义右边界法向结点以便与边界点形成法向单元
type,3
real,3
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
npmax=npmax+1
n,npmax,x,y-dxyz/2.,z !定义右边界切向结点以便与边界点形成切向单元
type,2
real,2
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
*enddo
allsel,all
/pnum,type,1
/number,1
eplot
finish
/solu
ANTYPE,trans
!*
TRNOPT,FULL
LUMPM,0
btime=0.02
etime=16.00
dtime=0.02
*DO,itime,btime,etime,dtime
TIME,itime
nsel,s,loc,y,H !选中需要加荷载的点
nsel,r,loc,x,-L/4,L/4
*if,itime,lt,1.,then
f,all,fy,1*itime
*elseif,itime,ge,1.0,and,itime,le,2.0
f,all,fy,1*(2-itime)
*else
f,all,fy,0.0
*endif
allsel,all
SOLVE
*ENDDO
另外，还用自己编写的有限元程序计算了一下这个例子，并与ANSYS得到的结果进行了比较，结果非常吻合，这里给出A点的比较结果。