分享 fortran调用ansys做桁架结构优化例子

sphinix

ANSYS命令流部分：
/clear
B=1000                   !基本尺寸
!进入PREP7并建模
/PREP7

*dim,A,,3,1  ！定义一个3*1的数组
*VREAD,A,E:\PROJECTS\TRUSS\trussdata_1,txt,,IJK,3,1   !IJK表示读入的顺序
(F7.2)
A1=A(1,1)
A2=A(2,1)
A3=A(3,1)
*cfopen,E:\PROJECTS\TRUSS\trussdata_3,txt
*vwrite,A1,A2,A3
(F7.2,/F7.2,/F7.2)
*cfclos

ET,1,LINK1               !二维杆单元
R,1,A1                   !以参数形式的实参
R,2,A2
R,3,A3
MP,EX,1,2.1E6            !杨氏模量
N,1,-B,0,0                !定义结点
N,2,0,0,0
N,3,B,0,0
N,4,0,-B,0
E,1,4                    !定义单元
REAL,2
E,2,4
REAL,3
E,3,4
FINISH
!
!进入求解器，定义载荷和求解
/SOLU
D,1,ALL,0,,3                 !结点UX=UY=0
F,4,FX,200000                !结点4上的X方向载荷分量(工况1)

F,4,FY,-200000               !结点4上的Y方向载荷分量(工况2)

SOLVE
FINISH

!进入POST1并读出状态变量数值
/POST1
SET,LAST
ETABLE,EVOL,VOLU         !将每个单元的体积放入ETABLE
SSUM                       !将单元表格内数据求和
*GET,VTOT,SSUM,,ITEM,EVOL  !VTOT=总体积
RHO=2.85E-4
WT=RHO*VTOT             !计算总体积

ETABLE,SIG,LS,1            !将轴向应力放入ETABLE
*GET,SIG1,ELEM,1,ETAB,SIG  !SIG1=第一个单元的轴向应力
*GET,SIG2,ELEM,2,ETAB,SIG  !SIG2=第二个单元的轴向应力


SIG1=ABS(SIG1)            !计算轴向应力的绝对值
SIG2=ABS(SIG2)

!/ESHAPE,2                !以实体单元模式显示壳单元
!/VIEW,1,1,1,1              !轴测视图
!EPLOT       　          !画单元
FINISH

*cfopen,E:\PROJECTS\TRUSS\trussdata_2,txt
*vwrite,SIG1,SIG2,WT
(F7.2,/F7.2,/F7.2)
*cfclos



FORTRAN程序部分：
PROGRAM TRUSS
USE DFLIB
IMPLICIT NONE

character(len=20) :: filename1="trussdata_2.txt"
character(len=20) :: filename2="trussdata_1.txt"
REAL,parameter :: e=0.01  !e为计算精度
REAL,parameter :: xu=400   !xu为许用应力
INTEGER :: fileid,k

INTEGER i
LOGICAL(4) result
REAL u1(50),u2(50),ww(50),umax    !u1,u2为应力比,umax为最大应力
REAL sig1,sig2,W,ap  !sig1为杆1的最大应力,sig2为杆2的最大应力

REAL x1(50),x2(50),xx1,xx2  !x为杆截面面积

i=1
k=1
x1(1)=707  !定义杆的初始值
x2(1)=707

ap=10


DO        WHILE( ap>=e )

!调用ANSYS做结构分析
result=SYSTEMQQ('C:\Ansys81\v81\ANSYS\bin\intel\ANSYS81 -b -p &
& ane3fl -i E:\ANSYSOBJECT\truss.txt -o E:\ANSYSOBJECT\trussanswer.txt')

fileid=10
open(fileid,file=filename1)  !从ansys写出的文件中读入数据(应力和求得的重量)
read(fileid,*) sig1,sig2,W
close(fileid)

ww(k)=w

u1(i)=sig1/xu
u2(i)=sig2/xu

if ( u1(i)>=u2(i) ) then  !判断最大应力
        umax=u1(i)
else
        umax=u2(i)
end if

!射线步
x1(i+1)=umax*x1(i)
x2(i+1)=umax*x2(i)

u1(i+1)=u1(i)/umax  !求出新的应力比
u2(i+1)=u2(i)/umax

!调整步
x1(i+2)=u1(i+1)*x1(i+1)
x2(i+2)=u2(i+1)*x2(i+1)

xx1=x1(i+2)/(500.0*1.414)
xx2=x2(i+2)/(500.0*1.414)

fileid=20
open(fileid,file=filename2)  !
write(fileid,"(A3,F7.2,/A3,F7.2,/A3,F7.2)") "A1=",x1(i+2),"A2=",x2(i+2),"A3=",x1(i+2)
close(fileid)

if ( i>=3) then
        if ( ww(k)-ww(k-1)>0 ) then
                write(*,"(A3,F7.2,/A3,F7.2,/A3,F7.2)") "A1=",x1(i),"A2=",x2(i),"A3=",x1(i)
                write(*,*) ww(k),ww(k-1)
                pause
        end if
end if
               
i=i+2

if ( i>=5 ) then        !ap为精度控制
        ap=sqrt((x1(i+1)-x1(i-1))**2+(x2(i+1)-x2(i-1))**2)
end if

k=k+1
END DO

STOP

END

icepoint

这种二次开发功能偶也想做，只是很无奈

zzzxxxx

不错，只是用满应力法是不是太落伍了啊？有没有求敏度的啊？

AaronSpark

原帖由 zzzxxxx 于 2006-8-20 16:33 发表
不错，只是用满应力法是不是太落伍了啊？有没有求敏度的啊？

效果好就行，不管用什么方法，否则算法什么的再漂亮那也是虚的

eongoo

fortran不会啊
有限元编程好像大多用fortran、c、c++
有没有用java的啊

xinyuxf

java进行矩阵运算太麻烦了吧，得自己编出来啊

mjhzhjg

楼上的说得很有道理，java面向工程要差一些，fortran就是针对工程运算的，就目前来说，matlab 应该是最佳选择，他有许多内置程序，直接调用，矩阵运算更为简单了。当然，每个软件都有优缺点，还是需要自己多学习，掌握更多的工具，才能跟上时代

AaronSpark

科学计算主要用的还是fortran，不过现在用c的也不少，主要还是从计算效率的角度考虑的比较多

maltab处理一些小规模计算问题还行，大型计算和fortran差很远，上次看到有人说几乎差了3倍左右

xingyuan085

我对ansys不熟，看其它软件可以不，做前处理的。

schultz

感谢楼主的例子，让我大开眼界，打消了一些疑虑。
BOW

jli

思路不错，是一个集成优化的例子。鼓励原创。

fuchengwang

但是FORTRAN编程没有MATLAB要好一些。

mjhzhjg

你的话怎么看不明白呢？
fortran对于大型计算确实要好一些，因为它从来就是针对工程开发的，matlab几乎涉及了各个方面，所谓少者专，多者泛，杂或许是它慢的原因之一