maoyj 发表于 2008-5-18 10:01:18

[原创]ADINA轴对称模型的反力换算

ADINA的轴对称模型显示反力（接触反力和支座反力）时，由于半径的影响，不同半径的单元面积是不一样的，这就造成了反力的不直观，不直接。如果是压力，特别是均布压力，无法直接从反力结果判断出压力的数值或者趋势。本文讨论了不同节点数单元的反力到压力的换算方法，并对2D、3D模型及不同单元进行了比较。
完整版请见附录PDF文件，无图版在2楼以后
PDF文件1.8M，共4个包。

maoyj 发表于 2008-5-18 10:07:58

ADINA轴对称模型的反力换算


1
前言在有限元模型中,轴对称模型被广泛应用。通用有限元软ADINA的轴对称模型显示反力（接触反力和支座反力）时，由于半径的影响，不同半径的单元面积是不一样的，这就造成了反力的不直观，不直接。如果是压力，特别是均布压力，无法直接从反力结果判断出压力的数值或者趋势。本文讨论了不同节点数单元的反力到压力的换算方法，并对2D、3D模型及不同单元进行了比较。

2
均布荷载下的反力在各节点的分配1图为均布荷载（大小1000）作用于四个完全相同的固定于XY平面的壳单元（Z方向约束），支座反力图。（命令流文件见附录）

左下为4节点面单元，右下为8节点单元，左上为9节点单元，右上为16节点单元。
通过对反力数值提取后的分析，可以得出：4节点单元各节点分配1/4。8节点单元角节点分配-1/12，边节点分配1/3。9节点单元角节点分配1/36，边节点分配1/9，中心节点分配4/9。16节点单元角节点分配1/64，边节点分配3/64，中心节点分配9/64。在边界处互相叠加。（见下图）

2D的平面的线型单元（iso-beam的轴对称壳）：
2节点：（1/2，1/2）。3节点：（1/6，2/3，1/6）。4节点：（1/8，3/8，3/8，1/8）。

3
平面轴对称问题的反力换算
通过对例子（IN文件见附录）中的不同节点的反力比较，可以按以下方法在EXCEL中将反力还原成均布压力：
3.1
二节点单元
单元编号
z向反力
换算后反力
换算公式
计算半径
p=换算反力/计算半径
1
166.667
333.334
=反力*2
0.333333
1000.002
2
1000
1000
=反力
1
1000
3
2000
2000
=反力
2
1000
4
3000
3000
=反力
3
1000
5
4000
4000
=反力
4
1000
6
5000
5000
=反力
5
1000
7
6000
6000
=反力
6
1000
8
7000
7000
=反力
7
1000
9
8000
8000
=反力
8
1000
10
9000
9000
=反力
9
1000
11
4833.33333
9666.6667
=反力*2
9.666667
1000
说明：1、施加的均布荷载大小为1000。
2、边界节点由于分配系数是1/2，所以要*2换算。
3、单元公用点为1/2+1/2=1，不用换算。
4、对于2D问题，计算中已经包含了2pi项，单元面积=单元径向长度*环向长度=2pi*r*1（本算例中单元长度1，10个单元共10）。
5、交界处的面积可以看作左单元的右半部分和右单元的左半部分，计算半径为节点所在位置。但是模型边界节点，轴处计算长度为1/3个单元长，最远边界处为节点位置-1/3个单元长=10-1/3。

3.2
三节点单元节点号
z向反力
换算反力
换算公式
计算半径r
p=换算反力/计算半径
12
0
反力*6
0.0
#DIV/0!
23
3.33E+02
499.9995
反力*3/2
0.5
999.999
13
3.33E+02
999.999
反力*3
1.0
999.999
24
1.00E+03
1500
反力*3/2
1.5
1000.000
14
6.67E+02
2000.001
反力*3
2.0
1000.001
25
1.67E+03
2500.005
反力*3/2
2.5
1000.002
15
1.00E+03
3000
反力*3
3.0
1000.000
26
2.33E+03
3499.995
反力*3/2
3.5
999.999
16
1.33E+03
3999.99
反力*3
4.0
999.998
27
3.00E+03
4500
反力*3/2
4.5
1000.000
17
1.67E+03
5000.01
反力*3
5.0
1000.002
28
3.67E+03
5500.005
反力*3/2
5.5
1000.001
18
2.00E+03
6000
反力*3
6.0
1000.000
29
4.33E+03
6499.995
反力*3/2
6.5
999.999
19
2.33E+03
6999.99
反力*3
7.0
999.999
30
5.00E+03
7500
反力*3/2
7.5
1000.000
20
2.67E+03
8000.01
反力*3
8.0
1000.001
31
5.67E+03
8500.005
反力*3/2
8.5
1000.001
21
3.00E+03
9000
反力*3
9.0
1000.000
32
6.33E+03
9499.995
反力*3/2
9.5
999.999
22
1.67E+03
1.00E+04
反力*6
10.0
1000.002
说明：1、施加均布荷载1000。
2、对于3节点单元，轴节点（12）没有反力数值。
3、由于ADINA软件本身在创建节点时会根据位置依次创建，单元中心的节点编号较大，可在EXCEL中按下公式输入：
= IF(A89>=23,B89*3/2,B89*3)
其中A89为节点11节点号所在单元格，B89位节点11的反力所在单元格
4、上公式在模型边界点不适用。轴心节点和最远端边界节点由于分配系数1/6，应*6。中间的单元交界节点分配系数1/6+1/6=1/3，应*3。单元内节点分配系数2/3，应*1.5
5、三节点单元的计算半径就是节点所在位置。
3.3
四节点单元节点号
z向反力
换算反力
换算公式
计算半径
p=换算反力/计算半径
33
1.67E+01
133.3336
反力*8
0.133
1000.002
44
7.50E+01
200
反力*8/3
0.200
1000
45
3.00E+02
800
反力*8/3
0.800
1000
34
2.50E+02
1000
反力*4
1.000
1000
46
4.50E+02
1200
反力*8/3
1.200
1000
47
6.75E+02
1800
反力*8/3
1.800
1000
35
5.00E+02
2000
反力*4
2.000
1000
48
8.25E+02
2200
反力*8/3
2.200
1000
49
1.05E+03
2800
反力*8/3
2.800
1000
36
7.50E+02
3000
反力*4
3.000
1000
50
1.20E+03
3200
反力*8/3
3.200
1000
51
1.43E+03
3800
反力*8/3
3.800
1000
37
1.00E+03
4000
反力*4
4.000
1000
52
1.58E+03
4200
反力*8/3
4.200
1000
53
1.80E+03
4800
反力*8/3
4.800
1000
38
1.25E+03
5000
反力*4
5.000
1000
54
1.95E+03
5200
反力*8/3
5.200
1000
55
2.18E+03
5800
反力*8/3
5.800
1000
39
1.50E+03
6000
反力*4
6.000
1000
56
2.33E+03
6200
反力*8/3
6.200
1000
57
2.55E+03
6800
反力*8/3
6.800
1000
40
1.75E+03
7000
反力*4
7.000
1000
58
2.70E+03
7200
反力*8/3
7.200
1000
59
2.93E+03
7800
反力*8/3
7.800
1000
41
2.00E+03
8000
反力*4
8.000
1000
60
3.08E+03
8200
反力*8/3
8.200
1000
61
3.30E+03
8800
反力*8/3
8.800
1000
42
2.25E+03
9000
反力*4
9.000
1000
62
3.45E+03
9200
反力*8/3
9.200
1000
63
3.68E+03
9800
反力*8/3
9.800
1000
43
1.23E+03
9866.64
反力*8
9.867
999.9972973
说明：1、施加荷载1000
2、单元中心节点编号较大，换算式可按下公式在EXCEL中输入：
=IF(A124>=44,B124/3*8,B124*4)
中A124为节点33节点号所在单元格，B124为节点33的反力在单元格
3、上公式在模型边界点不适用。轴心节点和最远端边界节点由于分配系数1/8，应*8。中间的单元交界节点分配系数1/8+1/8=2/8，应*4。单元内节点分配系数3/8，应*8/3
4、交界处的面积可以看作左单元的右半部分和右单元的左半部分，计算半径为节点所在位置。但是模型边界节点，轴处计算长度为2/15个单元长，最远边界处为节点位置-2/15个单元长=10-2/15。单元的内节点应按以下方法采用：左侧内节点为节点位置-2/15=*.3333-2/15=*.2；右侧内节点为节点位置+2/15=*.6667+2/15=*.8。其中*为各节点位置的整数部分。




[ 本帖最后由 maoyj 于 2008-5-18 10:09 编辑 ]

maoyj 发表于 2008-5-18 10:13:08

4
空间轴对称问题的反力换算通过对例子（IN文件见附录）中的不同节点的反力比较，可以按以下方法在EXCEL中将反力还原成均布压力：

4.1
四节点单元

节点号z向反力换算反力换算公式计算半径单元面积压力计算结果4065.21E+042.89E+03反力/36*20.33330.05817849746.58314058.68E+038682.41=反力10.17453349746.54494041.74E+0417364.8=反力20.34906649746.48764032.60E+0426047.2=反力30.52359949746.48764023.47E+0434729.6=反力40.69813249746.48764014.34E+0443412=反力50.87266549746.48764005.21E+0452094.5=反力61.04719849746.58313996.08E+0460776.9=反力71.2217349746.569463986.95E+0469459.3=反力81.39626349746.559223977.81E+0478141.7=反力91.57079649746.551273614.20E+0483930反力*29.66671.68715249746.56663
说明：1、施加荷载为50000均布压力
2、环向单元划分36；轴向10个单元，每个单元长度1
3、边界节点由于分配系数1/4+1/5=1/2，应*2，四个单元交界节点分配系数1/4+1/4+1/4+1/4=1不做变动。轴心节点由于是环向36个单元的节点都退化成的点，应/36份。
4、交界处的面积可以看作内单元的外半部分和外侧单元的内半部分，计算半径为节点所在位置。但是模型边界节点，轴处计算长度为1/3个单元长，最远边界处为节点位置-1/3个单元长=10-1/3。
5、单元面积为PI()*计算半径*2/36。3D模型的计算中没有简化2PI，必须包含。
6、最终计算结果和实际值偏差0.005094

4.2
九节点单元

节点号z向反力换算反力换算公式计算半径单元面积压力计算结果1522-3.98E+01-6.63E+00反力/36*180.0 0#DIV/0!18472.91E+0313083.3反力*9/20.5 0.087266149923.574415214.86E+024377.69反力*91.0 0.17453325082.316118462.91E+0313103.19反力*9/21.5 0.26179950050.4990115201.94E+0317470.98反力*92.0 0.34906650050.670918454.85E+0321838.68反力*9/22.5 0.43633250050.5677715192.91E+0326206.38反力*93.0 0.52359950050.4990118446.79E+0330574.17反力*9/23.5 0.61086550050.5972315183.88E+0334941.87反力*94.0 0.69813250050.5419818438.74E+0339309.615反力*9/24.5 0.78539850050.5563115174.85E+0343677.36反力*95.0 0.87266550050.5677718421.07E+0448045.15反力*9/25.5 0.95993150050.6240215165.82E+0352412.85反力*96.0 1.04719850050.5849518411.26E+0456780.55反力*9/26.5 1.13446450050.551915156.79E+0361148.34反力*97.0 1.2217350050.5972318401.46E+0465515.95反力*9/27.5 1.30899750050.4990115147.76E+0369883.83反力*98.0 1.39626350050.6064418391.65E+0474251.35反力*9/28.5 1.4835350050.4585715138.74E+0378619.23反力*99.0 1.57079650050.5563118381.84E+0482987.2反力*9/29.5 1.65806350050.698043614.85E+038.74E+04反力*1810.0 1.74532950050.56777说明：1、施加荷载为50000均布压力
2、环向单元划分36；轴向10个单元，每个单元长度1
3、边界节点由于分配系数1/36+1/36=1/18，应*18，四个单元交界节点分配系数1/36+1/36+1/36+1/36=1/9应*9。两个单元交界的边缘节点分配系数1/9+1/9=2/9，应*9/2。轴心节点由于是环向36个单元的节点都退化成的点，应/36份。单元中心节点编号较大，换算式可按下公式在EXCEL中输入：
=IF(A56>1800,B56*9/2,B56*9)
中A56为节点1522节点号所在单元格，B56为节点1522的反力在单元格
4、交界处的面积可以看作内单元的外半部分和外侧单元的内半部分，计算半径为节点所在位置。
5、单元面积为PI()*计算半径*2/36。3D模型的计算中没有简化2PI，必须包含。
6、最内侧单元反力存在畸变，应抛弃其数据
7、最终计算结果和实际值偏差-0.00101
4.3
十六节点单元

节点号z向反力换算反力换算公式计算半径单元面积压力计算结果33585.24E+031163.5533反力/36*80.1333330.02327150000.0214336936.54E+021745.3307反力*16/60.20.03490750000.0405336922.62E+036981.3067反力*16/60.80.13962649999.925943357-9.27E+028726.6353见说明510.17453349999.9371736919.82E+0210471.989反力*64/61.20.2094450000.0659936901.47E+0315707.947反力*64/61.80.31415949999.9471633561.09E+0317453.28反力*64/420.34906649999.96413……中间略……33504.36E+0369813.28反力*64/481.39626350000.0787236776.71E+0371558.613反力*64/68.21.4311750000.0796636767.20E+0376794.56反力*64/68.81.5358950000.0474733494.91E+0378540反力*64/491.57079650000.1169236757.53E+0380285.227反力*64/69.21.60570350000.0504936748.02E+0385521.28反力*64/69.81.71042350000.085743612.69E+038.61E+04反力*329.8666671.72205850017.543112节点单元

-1/8

1/16

-1/8

-1/4

3/16

3/16

3/8

3/8

3/8

3/8

3/16

3/16

说明：1、此处仅列出内侧两单元节点和最外侧两单元节点
2、施加荷载为50000均布压力
3、环向单元划分36；轴向10个单元，每个单元长度1
4、边界节点由于分配系数1/64+1/64=1/32，应*32，四个单元交界节点分配系数1/64+1/64=4/64应*64/4。两个单元交界的边缘节点分配系数3/64+3/64=6/64，应*64/6。轴心节点由于是环向36个单元的节点都退化成的点，应/36份。单元中心节点编号较大，换算式可按下公式在EXCEL中输入：
=IF(A150>=3600,B150/6*64,B150/4*64)
中A150为节点3358节点号所在单元格，B150为节点3358的反力单元格
5、四个单元交界处的面积可以看作内单元的外半部分和外侧单元的内半部分，计算半径为节点所在位置。两个单元交界处的边缘节点，和2D情况类似，内、外侧的边缘节点分别在节点位置-+2/15得到*.2和*.8。最外层节点计算半径10-2/15，最内层单元为12节点的矩形单元退化而成。其分配系数为



故轴心节点的分配系数为3/16+3/16-1/8-1/8=1/8,应*8
最内层单元的两单元交界面的边缘节点,分配系数3/8,应*8/3
最内层单元和次内层单元的交界节点,为12节点单元和16节点单元交界,其计算系数本来应该为(1-2/15)*(-1/16)+(1+2/15)*(1/32)=-3/32+1/3/6(对于16节点内外单元交界可理解为(1-2/15)*(1/32)+(1+2/15)*(1/32)=1/32+1/32=1/16),但是实际发现使用(1+2/15)*(-1/16)+(1+2/15)*(1/32)=-3/32-1/5/6可以得到正确的结果，原因不明。该节点的反力换算公式为：反力/(-3/32-1/5/16)

6、单元面积为PI()*计算半径*2/36。3D模型的计算中没有简化2PI，必须包含。
8、最终计算结果和实际值偏差1.32E-6


5
计算结果误差比较通过实际的模型，可以看出：
2D轴对称模型（YZ平面）的误差明显比3D模型小，因为3D模型的环向单元划分数少（划分越密结果越好，但不经济）。2D模型的平均误差在2E-6级别。其中轴心节点和最边缘节点的误差较大，其他节点的误差非常小。如果不是特殊要求，适用2D模型可以得到很好的结果
3D轴对称模型的误差明显比2D大。在环向36个单元的情况下，4节点单元误差在0.005级别，9节点单元的误差在0.001级别，16节点单元的误差在1E-6级别。可见，随着节点数的增加，计算精度迅速增加
无中心节点的8节点单元（3D）和12节点单元（3D），由于分配系数存在负值，应尽量不使用。

武汉游侠 发表于 2008-5-19 09:30:11

支持一下！

skyhawkwu 发表于 2008-5-19 14:00:04

回去研究一下，太好了……

minsen83712 发表于 2008-5-20 11:32:21

不好意思，我简单看了一下，没有看明白，好像挺难的。好好研究一下才行啊

肖伟 发表于 2008-5-27 15:26:42

似乎很难哦

jiyafeng 发表于 2009-5-9 21:47:19

头有点晕晕的感觉，还得静下心来仔细推敲
页: [1]
查看完整版本: [原创]ADINA轴对称模型的反力换算