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

maoyj

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

maoyj

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

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

4.1
四节点单元

节点号	z向反力	换算反力	换算公式	计算半径	单元面积	压力计算结果
406	5.21E+04	2.89E+03	反力/36*2	0.3333	0.058178	49746.5831
405	8.68E+03	8682.41	=反力	1	0.174533	49746.5449
404	1.74E+04	17364.8	=反力	2	0.349066	49746.4876
403	2.60E+04	26047.2	=反力	3	0.523599	49746.4876
402	3.47E+04	34729.6	=反力	4	0.698132	49746.4876
401	4.34E+04	43412	=反力	5	0.872665	49746.4876
400	5.21E+04	52094.5	=反力	6	1.047198	49746.5831
399	6.08E+04	60776.9	=反力	7	1.22173	49746.56946
398	6.95E+04	69459.3	=反力	8	1.396263	49746.55922
397	7.81E+04	78141.7	=反力	9	1.570796	49746.55127
361	4.20E+04	83930	反力*2	9.6667	1.687152	49746.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*18	0.0	0	#DIV/0!
1847	2.91E+03	13083.3	反力*9/2	0.5	0.087266	149923.5744
1521	4.86E+02	4377.69	反力*9	1.0	0.174533	25082.3161
1846	2.91E+03	13103.19	反力*9/2	1.5	0.261799	50050.49901
1520	1.94E+03	17470.98	反力*9	2.0	0.349066	50050.6709
1845	4.85E+03	21838.68	反力*9/2	2.5	0.436332	50050.56777
1519	2.91E+03	26206.38	反力*9	3.0	0.523599	50050.49901
1844	6.79E+03	30574.17	反力*9/2	3.5	0.610865	50050.59723
1518	3.88E+03	34941.87	反力*9	4.0	0.698132	50050.54198
1843	8.74E+03	39309.615	反力*9/2	4.5	0.785398	50050.55631
1517	4.85E+03	43677.36	反力*9	5.0	0.872665	50050.56777
1842	1.07E+04	48045.15	反力*9/2	5.5	0.959931	50050.62402
1516	5.82E+03	52412.85	反力*9	6.0	1.047198	50050.58495
1841	1.26E+04	56780.55	反力*9/2	6.5	1.134464	50050.5519
1515	6.79E+03	61148.34	反力*9	7.0	1.22173	50050.59723
1840	1.46E+04	65515.95	反力*9/2	7.5	1.308997	50050.49901
1514	7.76E+03	69883.83	反力*9	8.0	1.396263	50050.60644
1839	1.65E+04	74251.35	反力*9/2	8.5	1.48353	50050.45857
1513	8.74E+03	78619.23	反力*9	9.0	1.570796	50050.55631
1838	1.84E+04	82987.2	反力*9/2	9.5	1.658063	50050.69804
361	4.85E+03	8.74E+04	反力*18	10.0	1.745329	50050.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向反力	换算反力	换算公式	计算半径	单元面积	压力计算结果
3358	5.24E+03	1163.5533	反力/36*8	0.133333	0.023271	50000.02143
3693	6.54E+02	1745.3307	反力*16/6	0.2	0.034907	50000.04053
3692	2.62E+03	6981.3067	反力*16/6	0.8	0.139626	49999.92594
3357	-9.27E+02	8726.6353	见说明5	1	0.174533	49999.93717
3691	9.82E+02	10471.989	反力*64/6	1.2	0.20944	50000.06599
3690	1.47E+03	15707.947	反力*64/6	1.8	0.314159	49999.94716
3356	1.09E+03	17453.28	反力*64/4	2	0.349066	49999.96413
……中间略……
3350	4.36E+03	69813.28	反力*64/4	8	1.396263	50000.07872
3677	6.71E+03	71558.613	反力*64/6	8.2	1.43117	50000.07966
3676	7.20E+03	76794.56	反力*64/6	8.8	1.53589	50000.04747
3349	4.91E+03	78540	反力*64/4	9	1.570796	50000.11692
3675	7.53E+03	80285.227	反力*64/6	9.2	1.605703	50000.05049
3674	8.02E+03	85521.28	反力*64/6	9.8	1.710423	50000.08574
361	2.69E+03	8.61E+04	反力*32	9.866667	1.722058	50017.5431

12节点单元

-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），由于分配系数存在负值，应尽量不使用。

武汉游侠

支持一下！

skyhawkwu

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

minsen83712

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

肖伟

似乎很难哦

jiyafeng

头有点晕晕的感觉，还得静下心来仔细推敲