邓肯张模型FORTRAN子程序源代码
SUBROUTINE UMA T(STRESS,STA TEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MA TERL,NDI,NSHR,NTENS,
3 NSTA TV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 MA TERL
DIMENSION STRESS(NTENS),STA TEV(NSTA TV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),
4 DFGRD0(3,3),DFGRD1(3,3)
C
DIMENSION PS(3),DSTRESS(NTENS),TDSTRESS(NTENS),TSTRESS(NTENS) PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
K=PROPS(1)
N=PROPS(2)
RF=PROPS(3)
C=PROPS(4)
FAI=PROPS(5)/180.0*3.1415926
G=PROPS(6)
D=PROPS(7)
F=PROPS(8)
KUR=PROPS(9)
PA=PROPS(10)
DFAI=PROPS(11)/180.0*3.1415926
S1S3O=STA TEV(1)
S3O=STA TEV(2)
SSS=STA TEV(3)
CALL GETPS(STRESS,PS,NTENS)
FAI=FAI-DFAI*LOG10(S3O/PA)
CALL GETEMOD(PS,K,N,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O,G,D,F
1 ,SSS,S1S3O)
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
TDSTRESS=0.0
CALL GETSTRESS(DDSDDE,TDSTRESS,DSTRAN,NTENS)
DO 701 I1=1,NTENS
TSTRESS(I1)=STRESS(I1)+TDSTRESS(I1)*0.5
701 CONTINUE
CALL GETPS(TSTRESS,PS,NTENS)
CALL GETEMOD(PS,K,N,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O,G,D,F,
1 SSS,S1S3O)
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
DSTRESS=0.0
CALL GETSTRESS(DDSDDE,DSTRESS,DSTRAN,NTENS)
DO 702 I1=1,NTENS
STRESS(I1)=STRESS(I1)+DSTRESS(I1)
702 CONTINUE
CALL GETPS(STRESS,PS,NTENS)
CALL GETEMOD(PS,K,N,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O,G,D,F,
1 SSS,S1S3O)
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
IF(PS(3).GT.S3O)S3O=PS(3)
IF((PS(1)-PS(3)).GT.S1S3O)S1S3O=PS(1)-PS(3)
IF(S.GT.SSS)SSS=S
STA TEV(1)=S1S3O
STA TEV(2)=S3O
STA TEV(3)=SSS
END
SUBROUTINE GETPS(STRESS,PS,NTENS)
INCLUDE 'ABA_PARAM.INC'
DIMENSION PS(3),STRESS(NTENS)
CALL SPRINC(STRESS,PS,1,3,3)
DO 310 I=1,2
DO 320 J=I+1,3
IF(PS(I).GT.PS(J))THEN
PPS=PS(I)
PS(I)=PS(J)
PS(J)=PPS
END IF
320 CONTINUE
310 CONTINUE
DO 330 K1=1,3
PS(K1)=-PS(K1)
330 CONTINUE
RETURN
END
SUBROUTINE GETEMOD(PS,K,EN,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O
1 ,G,D,F,SSS,S1S3O)
INCLUDE 'ABA_PARAM.INC'
DIMENSION PS(3)
S=(1-SIN(FAI))*(PS(1)-PS(3))
IF(PS(3).L T.0) THEN
PSFEI=0.1
ELSE
PSFEI=PS(3)
END IF
S=S/(2*C*COS(FAI)+2*PSFEI*SIN(FAI))
IF(S.GE.0.99) THEN
S=0.99
END IF
AA=D*(PS(1)-PS(3))
AA=AA/(K*PA*((S3O/PA)**N))
AA=AA/(1-RF*S)
ENU=G-F*LOG10(S3O/PA)
ENU=ENU/(1-AA)/(1-AA)
IF(ENU.GT.0.49)ENU=0.49
EMOD=K*PA*((S3O/PA)**N)*((1-RF*S)**2)
IF(S.LT.SSS.AND.(PS(1)-PS(3)).LT.S1S3O)THEN
EMOD=KUR*PA*((S3O/PA)**N)
END IF
END
SUBROUTINE GETFEI(STRESS,SMISES,NTENS,NDI)
INCLUDE 'ABA_PARAM.INC'
DIMENSION STRESS(NTENS)
SMISES=STRESS(1)**2+STRESS(2)**2+STRESS(3)**2+
1STRESS(4)**2+STRESS(5)**2+STRESS(6)**2
SMISES=SQRT(SMISES)
RETURN
END
SUBROUTINE GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
INCLUDE 'ABA_PARAM.INC'
DIMENSION DDSDDE(NTENS,NTENS)
DO 20 K1=1,NTENS
DO 10 K2=1,NTENS
DDSDDE(K2,K1)=0.0
10 CONTINUE
20 CONTINUE
DO 40 K1=1,NDI
DO 30 K2=1,NDI
DDSDDE(K2,K1)=ELAM
30 CONTINUE
DDSDDE(K1,K1)=EG2+ELAM
40 CONTINUE
DO 50 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EG
50 CONTINUE
RETURN
END
SUBROUTINE GETSTRESS(DDSDDE,STRESS,DSTRAN,NTENS)
INCLUDE 'ABA_PARAM.INC'
DIMENSION DDSDDE(NTENS,NTENS),STRESS(NTENS),DSTRAN(NTENS) DO 70 K1=1,NTENS
DO 60 K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
60 CONTINUE
70 CONTINUE
RETURN
END
第一种方法： / o/ J5 @6 U/ ^- o$ 1. 建立工作目录/ ]" 2. 将Abaqus安装目录\6.4-pr11\site下的aba_param_dp.inc或aba_param_sp.inc拷贝到工作目录，并改名为aba_param.inc； # ~/ |0 I0 E6 {, @4 X3 q: W3. 将编译的fortran程序拷贝到工作目录； 4. 将.obj文件拷贝到工作目录； 5. 建立好输入文件.inp； 6. 运行abaqusjob=inp_name user=fortran name即可。 第二种方法： 在Job模块里，创建工作，在EditJob对话框中选择General选项卡，在Usersubroutine file中点击Select 按钮，从弹出对话框中选择你要调用的子程序文件（后缀为.for或.f）。 , D8 i7 d/r c6 @" | 以下是网上摘录的资料，供参考：. |$ t/ }$W7 Y6 m4 h6 D6 j 用户进行二次开发时，要在命令行窗口执行下面的命令： 4 O. R+ ^,@( ? abaqus job=job_name user=sub_name ABAQUS会把用户的源程序编译成obj文件，然后临时生成一个静态库standardU.lib和动态库standardU.dll，还有其它一些临时文件，而它的主程序（如standard.exe和explicit.exe等）则没有任何改变，由此看来ABAQUS是通过加载上述2个库文件来实现对用户程序的连接，而一旦运行结束则删除所有的临时文件。这种运行机制与ANSYS、LS-DYNA、marc等都不同。 : j6 g' R-o( {0 [* N2 J3 X这些生成的临时文件要到文件夹C:\Documentsand Settings\Administrator\Local Settings\Temp\中才能找到，这也是6楼所说的藏了一些工作吧，大家不妨试一下。 1子程序格式(程序后缀是.f; .f90; .for;.obj??) 答：我试过，.for格是应该是不可以的，至少6.2和6.3版本应该是不行，其他的没用过，没有发言权。在Abaqus中，运行abaqusj=jobname user=username时，默认的用户子程序后缀名是.for（.f,.f90应该都不行的，手册上也有讲过），只有在username.for文件没有找到的情况下，才会去搜索username.obj，如果两者都没有，就会报错误信息。 如果username包括扩展名for或obj，那么就根据各自的扩展名ABAQUS会自动选择进行操作。 2CAE中如何调用？Command下如何调用？ 答：CAE中在creat job的jobmanager中的general中可以指定子程序； Command下用命令:abaqus j=jobnameuser=userfilename （无后缀）； 3若有多个子程序同时存在，如何处理 答：将其写在一个文件中即可，然后用一个总的子程序调用（具体参见手册） 4我对VF不是很熟，是否可以用VC，C＋＋编写子程序？ A: 若要在vf中调试，那么应该根据需要把SITE文件夹中的ABA_PARAM_DP.INC（双精度）或ABA_PARAM_SP.INC（单精度）拷到相应的位置，并改名为ABA_PARAM.INC即可。 据说6.4的将可以，6.3的你可以尝试着将VC，C＋＋程序编译为obj文件，没试过。在你的工作目录下应该已经存在ufield.obj和uvarm.obj这两个文件（这两个文件应该是你分别单独调试ufield.FOR和uvarm.FOR时自动编译生成的，你可以将他们删掉试试看），但是由于你的FOR文件中已经有了UV ARM 和UFIELD这两个subroutine，显然会造成重复定义，请查实。 用户子程序的使用 假设你的输入文件为：a.inp b.for 那么在ABAQUS Command 中的命令应该是这样的： abaqusjob=a user=b
塑性力学读书报告 邓肯-张模型研究认识 学院：建设工程 姓名：王吉亮 学号：2006631011 专业：地质工程
教师：金英玉
邓肯-张模型研究认识 王吉亮（83分） 摘 要:从邓肯-张模型的本源开始，分析研究了邓肯-张模型与E-B 模型的建立过程和模型中参数如何确定的问题，结合对该模型的认识，提出该模型具有的缺点与不足。 关键词:邓肯-张模型；E-B 模型；参数确定 CONGNITION ON THE STUDY OF DUNCAN-CHANG MODEL Wang Jiliang Abstract: rom the parent of Duncan-Chang model, studing the establish procedure of Duncan-Chang model and E-B model, introducing the problem of how to define the indexes in the model. Associate the congnition on this model, present the shortcomings. Keywords: Duncan-Chang model; E-B model; indexes define 1 引言 邓肯－张模型是一个非线性本构模型，既然是一个本构模型，可想而之他反应的是应力与应变之间的关系。说它是非线性的，那么反映应力应变关系的模量就不是一个常数Ｅ那么简单。在介绍该模型之前，先要介绍一个概念，就是反映非线性关系的增量广义胡克定律： 1123()t t t v d d d d E E σεσσ= -+ （1） 1963年，康纳（Kondner ）根据大量土的三 轴试验的应力应变关系曲线，提出可以用双曲线拟合出一般土的三轴试验13()~a σσε-曲线，即： 13a a a b εσσε-= + （2） 其中，a 、b 为试验常数。对于常规三轴压缩试验，1a εε=。邓肯等人根据这一双曲线应力应变关系提出了一种目前被广泛的增量弹性模型， 一般被称为邓肯－张（Duncan-Chang ）模型。 在常规三轴压缩试验中，13a a a b εσσε-=+可以写成： 1113 a b εεσσ=+- （3） 将常规三轴压缩试验的结果按 11 13 ~εεσσ-的关系进行整理，则二者近似成线性关系(见图1)。其中，a 为直线的截距；b 为直线的斜率。 在常规三轴压缩试验中，由于 230d d σσ==，所以切线模量为 ε1/(σ1 -σ3 ) -σ3 )ult 图1 1113 ~εεσσ-线性关系图 132 11()() t d a E d a b σσεε-= =+ （4） 在试验的起始点，10ε=，t i E E =，则： 1 i E a = ，这表明a 表示的是在这个试验中的起始变形模量E i 的倒数。如果1ε→∞，则： 131 ()ult b σσ-= （5）
FORTRAN 90 程序编程规范 Fortran 90 编程规范，使程序代码高度组织化，更加易读、易懂、易于维护，程序更加高效。使编出的程序更易懂、易于维护。 1 语言选择 数值预报创新系统软件开发应避免使用Fortran77 的某些过时特征以Fortran 90不一致的特征。选择Fortran 90 作为开发语言，并采用Fortran 90 的新功能，如动态内存的分配(dynamic memory allocation)、递归（recursion ）, 模块(modules)、POINTER 、长变量名、自由格式等。 Fortran 77其中某些只是一些冗余的功能，这些功能已经过时，另外，还有一些在Fortran90 中被证明是不好的用法，建议不要使用。 2 Fortran 90 的新特性 2.1.1 建议使用的Fortran 90 新特性 建议使用Fortran 90 提供的模块（module ），并用Use ONLY 指定module 中哪些变量或派生类型定义可用于调用程序。 尽量使用数组下标三元组，这样可优化并减少所需的代码行数。为提高可读性，要在括号内表明数组的维数，例如： 1dArrayA(:) = 1dArrayB(:) + 1dArrayC(:) 2dArray(: , :) = scalar * Another2dArray(: , :) 当访问数组的子集时，例如在有限差分等式中，可以通过使用下标三元组实现。例如：2dArray(: , 2:len2) = scalar *（ & Another2dArray(:, 1:len2 -1) & - Another2dArray(:, 2:len2) & ) 对程序单元（program units ）命名，并使用End program ，End subroutine ，End interface ，End module 等结构再次指定“program unit ”的名称。 在逻辑表达式中使用>、 >=、 ==、 <、 <=、 /=，它们分别代 替.gt.、.ge.、.eq.、.lt.、.le.、.ne. 。新的表示方法更接近标准的数学符号 在变量定义中始终使用“：：”；始终用“DIMENSION ”定义数组形状；始终用（len=）的语法格式声明字符变量的长度。
FORTRAN程序设计复习题 一、选择题 B (1)下列各FORTRAN表达式中合法的是 A) S+T*2P >= B) .NOT. (A*B+C) C) A2+B2/(C+D) <= D) (A+B).NOT．A*B．GT．(.NOT.只跟一个表达式) C (2)数学式(3/5)ex+y的FORTRAN表达式是 A) 3*EXP(X+Y)/5 B) 3*E* *(X+Y)/ C) (3/5)*EXP(X+Y)D) EXP(X+Y) D (3)下列FORTRAN77表达式中不合法的是 A) A．GT．B．EQV．C．GT．D B) A．AND．B．AND．C．AND．D C) .NOT．(X．LE．D) A．LT．B．LT．C．LT．D D(4)下列叙述中不正确的是 A) FORTRAN子程序可以单独编译 B) 对一个FORTRAN源程序进行编译和连接无误后可生成可执行文件 C) 即使编译和连接都正确无误，FORTRAN程序运行时仍可能出错 D) FORTRAN连接的主要任务是把函数库中的函数翻译成机器指令(正确描述：主要任务为连接目标文件) B (5)在下列FORTRAN77运算符中，优先级最高的是 A) .AND. B) .NOT. C) .OR. D) .EQ. B (6)FORTRAN表达式"6/5+9/2**3/2"的值为 A) 33 B) 1 C) 5 D) 3 A (7)下列FORTRAN77表达式中，合法的是： A) .AND.. B) 10.0 C) D) 提示：A）相当于 .AND.（.NOT.（）） D (8)关于编译一个FORTRAN源程序文件，下列说法中错误的是 A) 允许编译只有一个主程序而没有子程序的源文件 B) 允许编译有多个子程序的源文件 C) 允许编译只有一个子程序而没有主程序的源文件 D) 允许编译有多个主程序的源文件 C (9)在FORTRAN77源程序中，续行标志符必须放在 A) 第1列 B) 第1-6列C) 第6列D) 第5列 D (10)下列关于"SUBROUTIN E MAP(X，Y)"语句行的叙述中，不正确的是 A) 这是子程序的第一个语句 B) 字符串"MAP"是子程序名 C) 变量X是子程序的形参D) 子程序执行后，MAP将返回整型数据 提示：子程序无返回值，自定义函数才有） A (11)FORTRAN表达式"2/4+"的值是 A) B) 1 C) D) 0 提示：2/4默认等于整型，=》 D (12)FORTRAN表达式"MOD，"的值是 A) B)0.0 C) D) A (13下列FORTRAN运算符中，优先级最低的是 A）逻辑运算符.AND. B）算术运算符*
1、找出100-150之间和400-450之间能被9整除的数 2、找出100-999之间的回文数 3、找出水仙花数
4、输出“*”图形 5、输入n个数，找出大于平均值的数和最小数
6、把一个数值型数组的相同数删除到只剩一个 7、形成一个5*5矩阵，对角线元素为“i”，其余为“j”program juzhen implicit none integer I,J integer,parameter::size=5 integer::a(size,size) forall(I=1:size,J=1:size,I>J) a(I,J)=j forall(I=1:size,J=1:size,I==J) a(I,J)=i forall(I=1:size,J=1:size,I integer,allocatable::c(:) write(*,*) "输入数组A的数据个数" read(*,*) n write(*,*) "输入数组B的数据个数" read(*,*) m l=m+n allocate (a(n)) allocate (b(m)) allocate (c(l)) write(*,*) "从小到大输入A的元素" do i=1,n read(*,*) a(i) end do write(*,*) "从小到大输入B的元素" do i=1,m read(*,*) b(i) end do do i=1,n c(i)=a(i) end do i=1 do while(i 5.370569 应力差(б1-б3/100kPa轴向变形ε1 体应变εv ε3ε1/(б1-б3 0.5080.002250.00074-0.0007554.42913E-051.0020.004490.0013- 0.0015954.48104E-051.4630.006740.00176-0.002494.60697E-051.8490.008980.00223-0.0033754.85668E-052.1490.011230.00223-0.00455.22569E-052.3310.013480.00214-0.005675.78293E-052.4770.015720.00176-0.006986.34639E-052.5870.017970.00158-0.0081956.94627E-052.6650.020210.00111-0.009557.58349E-052.730.022460.00056-0.010958.22711E-052.7770.024710.00009-0.012318.89809E-052.8120.026950.00003-0.013469.58393E-052.8450.03032-0.00012-0.015220.0001065732.8780.03369-0.00195-0.017820.000117062.8890.03706-0.00288-0.019970.000128282.8940.04043-0.00381-0.022120.0001397032.8890.04492-0.00474-0.024830.0001554862.8790.04941-0.00567-0.027540.0001716222.8690.0539-0.0065-0.03020.000187872.8530.05839-0.00752- 0.0329550.0002046622.8360.06289-0.00827-0.035580.0002217562.809 0.06738 -0.00891 -0.038145 0.000239872 应力差(б1-б3/100kPa轴向变形ε1 体应变εv ε3ε1/(б1-б3 0.9090.001250.00075-0.000250.1375137512.2860.00350.00151- 0.0009950.1531058623.5330.005750.00236-0.0016950.1627512034.3410.0080.00302- 1)实例3—求多个半径下的圆周长 ! z3.f90 --Fortran95 ! FUNCTIONS: ! z3 - Entry point of console application. !************************************************************************* ! PROGRAM: z3 ! PURPOSE: Entry point for the console application. !************************************************************************ program z3 ! 求多个半径下的圆周长 ! 主程序 ! PROGRAM Z3 PRINT *, 'R=',1.2,'C=',C(1.2) PRINT *, 'R=',3.4,'C=',C(3.4) PRINT *, 'R=',15.6,'C=',C(15.6) PRINT *, 'R=',567.3,'C=',C(567.3) END program z3 !子程序 FUNCTION C(R) PI=3.1415926 C=2*PI*R RETURN ! Body of z3 end 2)实例4—键盘与显示器输入/输出 a)Fortran 基本操作 b)程序指令 ! ZXZ_I_O.f90 ! FUNCTIONS: ! ZXZ_I_O - Entry point of console application. ! PROGRAM: ZXZ_I_O ! PURPOSE: Entry point for the console application. !***************输入、输出样式种种************************** program ZXZ_I_O implicit none !变量声明的位置 INTEGER(2) i; INTEGER(4) j; INTEGER(4) m; REAL n INTEGER A,B ! Variables PRINT*,'输入整数A'; READ*, A PRINT*,'输入整数B'; READ*, B B=A+B PRINT*,'B=A+B=',B WRITE(*,*) 'A*B=',A*B PRINT* ,'以上为计算机的计算结果，注意B的值' 邓肯－张模型是一个非线性本构模型，既然是一个本构模型，可想而之他反应的是应力与应变之间的关系。说它是非线性的，那么反映应力应变关系的模量就不是一个常数Ｅ那么简单。在介绍该模型之前，先要介绍一个概念，就是反映非线性关系的增量广义胡克定律： 1123()t t t v d d d d E E σεσσ= -+ （1） 1963年，康纳（Kondner ）根据大量土的三 轴试验的应力应变关系曲线，提出可以用双曲线拟合出一般土的三轴试验13()~a σσε-曲线，即： 13a a a b εσσε-= + （2） 其中，a 、b 为试验常数。对于常规三轴压缩试验，1a εε=。邓肯等人根据这一双曲线应力应变关系提出了一种目前被广泛的增量弹性模型，一般被称为邓肯－张（Duncan-Chang ）模型。 在常规三轴压缩试验中，13a a a b εσσε-=+可以写成： 1113 a b εεσσ=+- （3） 将常规三轴压缩试验的结果按 11 13 ~εεσσ-的关系进行整理，则二者近似成线性关系(见图 1)。其中，a 为直线的截距；b 为直线的斜率。 在常规三轴压缩试验中，由于 230d d σσ==，所以切线模量为 ε1 /(σ1 -σ3 ) 1 b=1/(σ1 -σ3 )ult a =1/E i 图1 1113 ~εεσσ-线性关系图 132 11()() t d a E d a b σσεε-= =+ （4） 在试验的起始点，10ε=，t i E E =，则： 1i E a = ，这表明a 表示的是在这个试验中的起始变形模量E i 的倒数。如果1ε→∞，则： 131 ()ult b σσ-= （5） 由此可以看出b 代表的是双曲线的渐近线所对应的极限偏差应力13()ult σσ-的倒数。 在土的试样中，如果应力应变曲线近似于双曲线关系，则往往是根据一定的应变值（如 115%ε=）来确定土的强度13()f σσ-，而不可 能在试验中使1ε无限大，求取13()ult σσ-；对于有峰值点的情况，取1313()()f σσσσ-=-峰， 这样1313()()f σσσσ--ult <。定义破坏比R f 为： 1313()()f f R σσσσ-=-ult （6） 而 13131 ()()f f R b σσσσ== --ult （7） 将上式与1 i E a = 代入 132 11()() t d a E d a b σσεε-==+ （8） 得到： 邓肯张模型FORTRAN子程序源代码 SUBROUTINE UMA T(STRESS,STA TEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MA TERL,NDI,NSHR,NTENS, 3 NSTA TV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 MA TERL DIMENSION STRESS(NTENS),STA TEV(NSTA TV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PROPS(NPROPS),COORDS(3),DROT(3,3), 4 DFGRD0(3,3),DFGRD1(3,3) C DIMENSION PS(3),DSTRESS(NTENS),TDSTRESS(NTENS),TSTRESS(NTENS) PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0) K=PROPS(1) N=PROPS(2) RF=PROPS(3) C=PROPS(4) FAI=PROPS(5)/180.0*3.1415926 G=PROPS(6) D=PROPS(7) F=PROPS(8) KUR=PROPS(9) PA=PROPS(10) DFAI=PROPS(11)/180.0*3.1415926 S1S3O=STA TEV(1) S3O=STA TEV(2) SSS=STA TEV(3) CALL GETPS(STRESS,PS,NTENS) FAI=FAI-DFAI*LOG10(S3O/PA) CALL GETEMOD(PS,K,N,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O,G,D,F 1 ,SSS,S1S3O) EBULK3=EMOD/(ONE-TWO*ENU) EG2=EMOD/(ONE+ENU) EG=EG2/TWO EG3=THREE*EG ELAM=(EBULK3-EG2)/THREE CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG) Aspen plus fortran 子程序的编译、连接与运行 Aspen Plus可以用来模拟反应精馏，但提供的动力学方程形式比较死板，很多类型的动力学方程都不能直接输入，不得不外接一个fortran 子程序。本来用到fortran 的机会就不多，所以会的人也不多。何况编写好fortran 子程序后，还要用fortran 编译器来编译生成obj 文件，而且每次运行的时候都得连接一次，比较麻烦。我的电脑以前是装了编译器这些软件 的，但最近出了问题，只能编译fortran 子程序，运行aspen外接fortran 子程序的时候就会出错。恰好这几天要帮师兄写个fortran 子程序，写完后需要调试，但又运行不了，重装 呢也没把握，想起网上有兄弟说过可以用asplink命令生成一个dll文件，这样就不用运行 aspen的时候每次都连接编译器了。于是找来资料学习学习，居然成功了，欣喜之余，分享给大家。第一步当然是编写动力学子程序了（此处省略若干字） 第二步编译写好的动力学子程序（以网友soaty的动力学子程序为例） 打开Aspen Plus Simulation Engine JS M... 打开后是这样的 接下来改变文件夹路径，转到存放动力学子程序的那个文件夹 比如我的文件路径是C:\ProgramData\Aspe nTech\Aspen Plus \IPP 于是输入cd C:\ProgramData\AspenTech\Aspen Plus \IPP 回车后就转到这个文件夹了，这时候就可以进行编译了：编译命令：aspcomp 输入命令aspcomp ,回车 到文件夹一看你会发现已经生成了obj文件了, 计算圆周率 REAL R,R1,R2,PI ISEED=RTC() N0=0 N=300000 DO I=1,N R1=RAN(ISEED) R2=RAN(ISEED) R=SQRT(R1*R1+R2*R2) IF(R<1.0)N0=N0+1 END DO PI=4.0*N0/N WRITE(*,*)PI END 一）蒙特卡洛计算生日问题 假设有N个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50％的几率人群中至少有2个人生日相同。 INTEGER M(1:10000), NUMBER1(0:364), NUMBER2 REAL X,Y ISEED=RTC() DO J=1, 10000 NUMBER1=0 X=RAN(ISEED) NUMBER1(0)=INT(365*X+1) JJJ=1 DO I=1,364 Y=RAN(ISEED) NUMBER2=INT(365*Y+1) ETR=COUNT(NUMBER1.EQ.NUMBER2) IF (ETR= =1) THEN EXIT ELSE JJJ=JJJ+1 M(J)=JJJ NUMBER1(I)=NUMBER2 END IF END DO END DO DO I=1,10000 IF(M(I).LE.23) SUM=SUM+1 END DO PRINT *,SUM/10000 END 二）MONTE CARLO SIMULATION OF ONE DIMENSIONAL DIFFUSION 蒙特卡罗计算一维扩散问题 INTEGER X,XX(1:1000,1:1000) REAL XXM(1:1000) ! X:INSTANTANEOUS POSITION OF ATOM ! XX(J,I)：X*X ,J:第几天实验,I:第几步跳跃 ! XXM(I): THE MEAN OF XX WRITE(*,*) "实验天数JMAX，实验次数IMAX" READ(*,*) JMAX,IMAX ISEED=RTC() DO J=1,JMAX !第几天实验 X=0 !!! DO I=1,IMAX !第几步跳跃 RN=RAN(ISEED) IF(RN<0.5)THEN X=X+1 ELSE X=X-1 END IF XX(J,I)=X*X END DO END DO OPEN(1,FILE="C:\DIF1.DAT") DO I=1,IMAX XXM=0.0 XXM(I)=1.0*SUM(XX(1:JMAX,I))/JMAX !! WRITE(1,*) I, XXM(I) END DO CLOSE(1) END 三维的！ 三）通过该程序了解FORTRAN语言如何画图（通过像素画图）USE MSFLIB INTEGER XR，YR !在的区域中画一个圆 PARAMETER XR=400,YR=400 INTEGER R, S(1:XR,1:YR) X0=XR/2 ! 圆心位置X0，YO Y0=YR/2 R=MIN(X0-10,Y0-10) !圆半径 S=0 !像素的初始状态（颜色） DO I=1,XR DO J=1,YR IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10 IER=SETCOLOR(S(I,J)) IER=SETPIXEL(I,J) END DO END DO END 程序1 PROGRAM chaper1272 INTEGER :: A(10)=(/5,7,4,8,12,2,10,3,9,11/) INTEGER :: sum=0,ave !打开一数据文件,设置一个有格式直接存取文件,将10个数分2个记录写入文件。 !数据文件生成2个记录,每个记录行长度相同,记录长度为25 OPEN(1,FILE='input21.dat',FORM='FORMATTED',ACCESS='D IRECT',RECL=25) WRITE(1,"(5I5)",REC=1)(A(I)+10,I=1,5) !按格式说明将头5个数写入第1个记录 WRITE(1,"(5I5)",REC=2)(A(I)+10,I=6,10)!按格式说明将后5个数写入第2个记录 READ(1,"(5I5)",REC=2)(A(I),I=6,10) !按格式说明从第2个记录中读取后5个数 READ(1,"(5I5)",REC=1)(A(I),I=1,5) !按格式说明从第1个记录中读取头5个数 DO I=1,10 sum=sum+A(I) ENDDO ave=sum/10 !打开一个最大记录长度为22的有格式顺序存取文件 OPEN(2,FILE='input22.dat',FORM='FORMATTED',ACCESS='D IRECT',RECL=22) WRITE(2,"('10个数之和为：',I5)",REC=1) sum !输出1记录行,记录长度为22 WRITE(2,"('10个数平均值为：',I5)",REC=2)ave !输出1记录行,记录长度为22 WRITE(2,"(A)",REC=3) '程序运行正常结束。' !输出一个记录行,记录长度为22 END 程序2 PROGRAM average REAL sum,ave INTEGER n OPEN(1,file='score.dat') PRINT*,'正在统计平均成绩,请等待。' sum=0.0;n=0 DO READ(1,*,END=100) S sum=sum+s n=n+1 ENDDO 100 ave=sum/n 研究生课程作业邓肯张模型参数计算 学生姓名李俊 学科专业岩土工程 学号201420105614 任课教师周小文教授 作业提交日期2014年12月 1．计算轴向应变 c h h ?∑= 1ε 式中 1ε－轴向应变； h ?∑－固结下沉量，由轴向位移计测得 0h －土样初始高度 c h —按实测固结下沉的试样高度 c h ?—试样固结下沉量 2．计算按实测固结下沉的试样高度，面积： 式中 Ac －按实测固结下沉的试样面积 0V －土样初始体积 3．计算剪切过程中试样的平均面积： 式中 a A －剪切过程中平均断面积 c V －按实测固结下沉的试样的体积 i V ?－排水剪中剪切时的试样体积变化 按体变管或排水管读数求得 1h ?－固结下沉量，由轴向位移计测得 3. 计算主应力差 c i c h V V A ?-= 01 h h V V A c i c a ?-?-= C c c A h V ?= 103 1?=-a A CR σσ 式中 31σσ- － 主应力差 1σ―大主应力 3σ－小主应力 C －测力计率定系数 R －测力计读数 2 数据处理 2.1 3σ=100kPa 数据初步计算 当3σ=100kPa 时，各数据初步计算如表1所示。 围压100kPa 数据初步计算表 表1 2.1.1 由切线模量计算数据 对公式 ) (311 σσε-=a +b 1ε进行直线拟合，如图1所示。 图1 1131 /()~εσσε-拟合曲线 a =0.0002，1 i E a = =5000kPa b ==0.0028，()131 ult b σσ-= =263.16kPa ()13f σσ-=204.26kPa ，()()1313f f ult R σσσσ-= -=0.7762 2.1.2 由泊松比计算数据 对公式()313/f D εεε-=+-进行直线拟合，如图2所示。 第一种方法：在Job模块里，创建工作，在Edit Job对话框中选择General选项卡，在User subroutine file中点击Select按钮，从弹出对话框中选择你要调用的子程序文件（后缀为.for或.f）。 第二种方法： 1. 建立工作目录?崠 2. 将Abaqus安装目录\6.4-pr11\site下的aba_param_dp.inc 或aba_param_sp.inc 拷贝到工作目录，并改名为aba_param.inc； 3. 将编译的fortran程序拷贝到工作目录； 4.将.obj文件拷贝到工作目录； 5. 建立好输入文件.inp； 6. 运行abaqus job=inp_name user=fortran name即可。 以下是网上摘录的资料，供参考： 用户进行二次开发时，要在命令行窗口执行下面的命令： abaqus job=job_name user=sub_name ABAQUS会把用户的源程序编译成obj文件，然后临时生成一个静态库standardU.lib和动态库standardU.dll，还有其它一些临时文件，而它的主程序（如standard.exe和explicit.exe等）则没有任何改变，由此看来ABAQUS是通过加载上述2个库文件来实现对用户程序的连接，而一旦运行结束则删除所有的临时文件。这种运行机制与ANSYS、LS-DYNA、marc等都不同。 这些生成的临时文件要到文件夹C:\Documents and Settings\Administrator\Local Settings\Temp\中才能找到，这也是6楼所说的藏了一些工作吧，大家不妨试一下。1子程序格式(程序后缀是.f; .f90; .for;.obj??) 答：我试过，.for格是应该是不可以的，至少6.2和6.3版本应该是不行，其他的没用过，没有发言权。 在Abaqus中，运行abaqus j=jobname user=username时，默认的用户子程序后缀名是.for（.f,.f90应该都不行的，手册上也有讲过），只有在username.for文件没有找到的情况下，才会去搜索username.obj，如果两者都没有，就会报错误信息。如果username包括扩展名for或obj，那么就根据各自的扩展名ABAQUS会自动选择进行操作。 2CAE中如何调用？Command下如何调用？ 答：CAE中在creat job的job manager中的general中可以指定子程序；Command下用命令:abaqus j=jobname user=userfilename （无后缀）； 3若有多个子程序同时存在，如何处理 答：将其写在一个文件中即可，然后用一个总的子程序调用（具体参见手册）＋＋编写子程序？C，VC不是很熟，是否可以用VF我对4． A: 若要在vf中调试，那么应该根据需要把SITE文件夹中的 ABA_PARAM_DP.INC（双精度）或ABA_PARAM_SP.INC（单精度）拷到相应 第四章 1.program main implicit none write(*,*) "Have a good time." write(*,*) "That's not bad." write(*,*) '"Mary" isn''t my name.' end program 2.program main real, parameter :: PI=3 implicit none.14159 real radius write(*,*) "请输入半径长" read(*,*) radius write(*,"(' 面积='f8. 3)") radius*radius*PI end program 3.program main implicit none real grades write(*,*) "请输入成绩" read(*,*) grades write(*,"(' 调整后成绩为'f8.3)") SQRT(grades)*10.0 end program 4.integer a,b real ra,rb a=2 b=3 ra=2.0 rb=3.0 write(*,*) b/a ! 输出1, 因为使用整数计算, 小数部分会无条件舍去write(*,*) rb/ra ! 输出1.5 5.p rogram main implicit none type distance real meter, inch, cm end type type(distance) :: d write(*,*) "请输入长度:" read(*,*) d%meter d%cm = d%meter*100 d%inch = d%cm/2.54 write(*,"(f8.3'米='f8.3'厘米='f8.3'英寸')") d%meter, d%cm, d%inch end program 第五章 1.program main implicit none integer money real tax write(*,*) "请输入月收入" read(*,*) money if ( money<1000 ) then tax = 0.03 else if ( money<5000) then tax = 0.1 else tax = 0.15 end if write(*,"(' 税金为'I8)") nint(money*tax) end program 2.program main implicit none integer day character(len=20) :: tv write(*,*) "请输入星期几" read(*,*) day select case(day) case(1,4) tv = "新闻" case(2,5) tv = "电视剧" case(3,6) tv = "卡通" case(7) tv = "电影" case default write(*,*) "错误的输入" stop end select write(*,*) tv end program 3.program main implicit none integer age, money real tax write(*,*) "请输入年龄" read(*,*) age write(*,*) "请输入月收入" read(*,*) money if ( age<50 ) then if ( money<1000 ) then tax = 0.03 else if ( money<5000 )then tax = 0.10 else tax = 0.15 end if else if ( money<1000 ) then tax = 0.5 else if ( money<5000 )then tax = 0.7 else tax = 0.10 end if end if write(*,"(' 税金为'I8)") nint(money*tax) end program 4.program main implicit none integer year, days logical mod_4, mod_100, mod_400 write(*,*) "请输入年份" read(*,*) year mod_4 = ( MOD(year,4) == 0 ) mod_100 = ( MOD(year,100) == 0 ) mod_400 = ( MOD(year,400) == 0 ) if ( (mod_4 .NEQV. mod_100) .or. mod_400 ) then days = 366 else days = 365 end if write(*,"('这一年有'I3'天')") days stop end program 第六章 1.program main implicit none integer i do i=1,5 write(*,*) "Fortran" end do stop end program 2.program main implicit none integer i,sum sum = 0 do i=1,99,2 sum = sum+i end do write(*,*) sum stop end program 3.program main implicit none integer, parameter :: answer = 45 integer, parameter :: max = 5 integer weight, i do i=1,max write(*,*) "请输入体重" read(*,*) weight if ( weight==answer ) exit end do if ( i<=max ) then write(*,*) "猜对了" else write(*,*) "猜错了" end if stop end program 4.program main implicit none integer, parameter :: max=10 integer i real item real ans ans = 1.0 item = 1.0 do i=2,max item = item/real(i) ans = ans+item邓肯张模型参数(精)
fortran程序实例
邓肯-张模型公式推导 高土
邓肯张模型FORTRAN子程序源代码
aspenfortran子程序的编译、连接、运行
几个fortran程序
fortran程序30个
邓肯张模型模拟
完整word版ABAQUS中Fortran子程序调用方法—自己总结
Fortran95程序设计课后习题答案(word版方便)
相关主题
文本预览