ABAQUS-二次开发资料-UMAT - 图文

各个楼层及内容索引

2-------------------------------------什么是UMAT 3-------------------------------------UMAT功能简介

4-------------------------------------UMAT开始的变量声明

5-------------------------------------UMAT中各个变量的详细解释 6-------------------------------------关于沙漏和横向剪切刚度

7-------------------------------------UMAT流程和参数表格实例展示

8-------------------------------------FORTRAN语言中的接口程序Interface

9-------------------------------------关于UMAT是否可以用Fortran90编写的问题 10-17--------------------------------Fortran77的一些有用的知识简介 20-25\\30-32-----------------------弹塑性力学相关知识简介

34-37--------------------------------用户材料子程序实例JOhn-cook模型压缩包下载 38-------------------------------------JOhn-cook模型本构简介图

40-------------------------------------用户材料子程序实例JOhn-cook模型完整程序+david详细注

解[欢迎大家来看看,并提供意见,完全是自己的diy的,不保证完全正确,希望共同探讨,以便更正,带\部分,还望各位大师\\同仁指教]

1什么是UMAT???

1.1 UMAT功能简介!!![-摘自庄茁老师的书

UMAT子程序具有强大的功能，使用UMAT子程序：

(1)可以定义材料的本构关系，使用ABAQUS材料库中没有包含的材料进行计算，扩充程序 功能。ABAQUS软件2003年度用户年会论文集

(2)几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予ABAQUS中 的任何单元；

(3)必须在UMAT中提供材料本构模型的雅可比（Jacobian）矩阵，即应力增量对应变增量 的变化率。

(4)可以和用户子程序“USDFLD”联合使用，通过“USDFLD”重新定义单元每一物质点上传 递到UMAT中场变量的数值。

1.2 UMAT开始的变量声明

由于主程序与UMAT之间存在数据传递，甚至共用一些变量，因此必须遵守有关书写格式，UMAT中常用的变量在文件开头予以定义，通常格式为：

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,

2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) C

INCLUDE'ABA_PARAM.INC'-----此处是将ABAQUS本身自带的参量精度定义的文件包含进来[后面详说]

C

CHARACTER*80 CMNAME

DIMENSION STRESS(NTENS),STATEV(NSTATV),

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),DFGRD0(3,3),DFGRD1(3,3) -------------------------------以上是变量声明吧!

user coding to define DDSDDE,STRESS,STATEV,SSE,SPD,SCD and,if necessary,RPL,DDSDDT,DRPLDE,DRPLDT,PNEWDT

-------------------------------此处,看来是将用户定义材料属性的fortran程序编入 RETURN------------------这是返回值 END------------------------结束

UMAT中各个变量的详细解释[凡是-注明david的,都是我自己猜的,仅供参考] DDSDDE (NTENS ,NTENS)

是一个NTENS[Number of the Tensions----david]维的方阵，称作雅可比矩阵，应力增量/应变增量的偏导数，DDSDDE (I ,J)表示增量步结束时第J个应变分量的改变引起的第I个应力增量的变化!雅可比是一个对称矩阵，除非在“*USER MATERIAL”语句中加\参数 STRESS (NTENS)

应力张量矩阵，对应NDI[Number of the Direct Components--david]个直接分量和NSHR[Number of the shear Components-david]个剪切分量.在增量步的开始,应力张量矩阵中的数值通过UMAT和主程序之间的接口传递到UMAT中,在增量步的结束,UMAT将对应力张量矩阵更新,即[return].对于包含刚体转动的有限应变问题，一个增量步条用UMAT之前就已经对应力张量进行了刚体转动,因此在UMAT中只需处理应力张量的共旋部分-------这部分我没看明白,敬请高手指点.UMAT中应力张量的度量为柯西（真实）应力

STATEV (NSTATEV)[STATE VARIABLES (Number of the State Variables)]

用于存储状态变量的矩阵，在增量步开始时将数值传递到UMAT中。也可在子程序USDFLD或UEXPAN中先更新数据，然后增量步开始时将更新后的数据传递到UMAT中。在增量步的结束必须更新状态变量矩阵中的数据。

和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的数据更新以外，状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。-----不懂,请高手指教

状态变量矩阵的维数NATATEV，等于关键字“*DEPVAR”定义的数值。状态变量矩阵的维数通过ABAQUS输入文件中的关键字“*DEPVAR”定义，关键字下面数据行的数值即为状态变量矩阵的维数。 材料常数的个数，等于关键字“*USER MATERIAL”中“CONSTANTS”常数设定的值。 PROPS (NPROPS)

材料常数矩阵，矩阵中元素的数值对应于关键字“*USER MATERIAL”下面的数据行。 SSE，SPD，SCD

分别定义每一增量步的弹性应变能[Elastic Strain Energy]，塑性耗散[Plastic Dissipation]和蠕变耗散[Creep Dissipation]。它们对计算结果没有影响，仅仅作为能量输出。 STRAN (NTENS)：应变矩阵；

DSTRAN (NTENS)：[D--大抵代表Deta,增量的意思-david]应变增量矩阵； DTIME：增量步的时间增量； NDI：直接应力分量的个数； NSHR：剪切应力分量的个数；

NTENS：总应力分量的个数，NTENS =NDI +NSHR。

1.3关于沙漏刚度控制和横向剪切刚度

使用UMAT时需要注意单元的沙漏控制刚度和横向剪切刚度。通常减缩积分单元的沙漏控制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定义的。这些刚度基于材料初始剪切模量的值，通常在材料定义中通过“*ELASTIC”选项定义。但是使用UMAT的时候，ABAQUS对程序输入文件进行预处理的时候得不到剪切模量的数值。所以这时候用户必须使用“*HOURGLASS STIFFNESS”选项来定义具有沙漏模式的单元的沙漏控制刚度，使用“*TRANSVERSE SHEAR STIFFNESS”选项来定义板、壳、梁单元的横向剪切刚度。 ------------

一个问题:得不到剪切模量的数值 和

解决方案:用户必须使用“*HOURGLASS STIFFNESS”选项来定义具有沙漏模式的单元的沙漏控制刚度，使用“*TRANSVERSE SHEAR STIFFNESS”选项来定义板、壳、梁单元的横向剪切刚度

1.4关于UMAT的流程图和参数表格实例

跟大家说说所谓的接口程序Interface--FORTRAN的知识

在Fortan语言中,主调程序和被调程序是分别编译的.由于Fortran90对过程的许多功能做了扩充,有些功能单靠简单的调用语句已经无法反应,因而系统也就无法进行正确的编译,这时需要在主调程序中加入interface接口块,通过它为主调程序和被调程序指明一个显示的接口.如果被调用中哑元含有假定形

状[assumed-shape]数组,或可选变元,或含键盘输入的参数,就需要interface接口块说明.一般来讲,在Fortran90程序之间需要提供interface块有三种方法:

1.将interface接口块直接写入调用程序,并复制被调用程序的参数列表这种方法简单易用,但也相应增加了维护代码的工作量,因为只要被调用程序的参数列表发生变化,就必须相应改变interface接口块和调用[call]语句.

2.可以将一个软件包中所有程序的interface块写入一个模块中,该模块被软件包中的所有程序使用.这样做的优点是只需一个模块来检查interface定义,缺点是仍需对此模块和调用语句进行维护.

3.Fortran90编译器可在contains语句后自动提供过程之间的interface块,这种interface块可用于使用模块的任何程序.

建议在同一个软件包中使用2\\3的形式,在调用软件包的入口程序时采用1\\2的形式!

[是不是在UMAT中,我们所编译的带接口的Fortran程序为调用程序,原ABAQUS主程序为被调用程序,调用程序中的第一部分我们先复制被调入程序的参数列表????????似乎和ABAQUS主程序调用UMAT有些相反了???????不过个人认为interface作为一个接口块,在Fortran语法中应该放在主调程序中,且复制被调程序的参数列表.而UMAT的参数变量的声明,只不过是为了和ABAQUS使用一致的变量格式,当Fortran程序处理完材料的本构定义之后再以这一致的形式将变量值返回到ABAQUS主程序中吧]

2 关于UMAT用Fortran90可不可以的问题

在论坛上搜索了关于这个问题,某位高手说是可以的,只要你自己装的Fortran编译器能成功编译你的Fortran90程序即可,个人认为也是如此,不过我还没有证明过!

2.1关于Fortran77的一些有用的简介[1]:

Fortran77的基本结构:

1. 一个Fortran源程序由一个或者多个程序单位组成,每个独立的程序单位以'end'语句结束

2. 每个程序单位包括若干行[不能一行写多条语句,但是可以一个语句写成行];分为语句行[执行语

句行和非执行语句]和非语句行[注释行]

源程序的书写格式:

1. 不区分大小写:每行只能80列以内,并把80列分为四个区 2. 1-5列:标号区[1-5位整数;第一列为'*'或者'C'时,为注释行] 3. 第6列,续行标志区[非空格或非0;最多19个续行] 4. 7-72列:语句区[书写语句:一行只能写一个语句] 5. 73-80列:语句注释区[一般做编号注释] 程序的编辑与运行:

1. 创建源程序文件并编写源程序 2. 编译并连接源文件

3. 运行程序编译生成的可执行文件 常量: 六种

1. 整型常量[Integer]4位:正\\负整数和0 2. 实型常量[Real]4位:小数和指数形式 3. 双精度常量[Double precision]8位 4. 复型常量[Complex]8位 5. 逻辑性常量[Logical]4位 6. 字符型常量[Character]1位 变量

? 变量名

第一个字符是字母第一个字符后可跟1-5个字母或者数字 不区分大小写 空格不起作用

允许变量名与语言中具有特定含义的字同名,但尽量不适用 尽量\见名知义\? 变量类型

不同的变量类型用来存放不同类型的常量数据.变量相应的也有六种;变量在使用前必须首先对其进行类型说明,三种说明方法:

按优先级别排列

1. 类型说明语句:类型变量名列表,多个变量名之间用逗号隔开,如 REAL A,B,C

DOUBLE PRECISION X,Y,Z[或者REAL*8 X,Y,Z] CHARACTER*5 [缺省字符长度5 ] STR1,STR2*8,STR3*19

[此处,STR1未指明长度,则默认使用缺省长度5;STR2的长度为 8;STR3的长度为19] 2. IMPLCIT语句:IMPLICIT 类型[字母表],类型[字母表],....

如: IMPLICIT REAL(A-D),INTERGER(I-M),DOUBLE PRECISION (X,Z) 3. I-N规则:Fortran规定,不加说明的情况下,I-N为整型,其他都为实型 几点说明

类型语句说明位于所有可执行语句的前面;IMPLICIT语句位于最前面;

IMPLICIT NONE取消IMPLICIT说明和I-N规则,所有的变量必须显式说明;只在本程序单位有效.

1. 2. 3. 4. 5.

2.2关于Fortran77的一些有用的简介[2]:

算术表达式:

1. 运算符: +, -, *, /, ** 2. 优先级: ( ), ** , *\\ / , +\\- 3. 书写问题

1. * 不能省略

2. 括号不分大小写,成对出现

3. 多次乘方,按'先右边后左边'处理

4. 运算符不能连续出现,要用小括号隔开

5. 运算顺序

)---->函数---->**----->*,/----->+,-

6. 运算中的类型问题:运算类型相同:结果仍为原类型;不同,则自动转换成同一类型 7. 误差问题:

1. 溢出:超出有效数字范围-------->解决:很大或者很小的数用实型的指数表示

2. 误差:由于有效数字的位数限制,实型数运算存在误差-------->解决:转换成双精度

型或者避免因为书写不当造成有效数字的丢失

简单输出\\输入语句:

输入\\输出三要素:对象[哪些数据];格式;设备. 输出语句

输出语句的分类:格式输出'表控格式输出[系统标准格式];无格式输出[二进制]

1. 表控输出语句:按计算机系统规定的格式输出:系统自动为每种类型的数据规定了列数

1. 整数的表控输出格式[与具体的计算机系统有关]:规定每个整数占13列,打印在右端,左

补空格;

2. 实数的表控输出格式:规定每个实数占17列,打印在右端,左补空格,小数部分占7列;[当实数的绝对值>=10**7或<1时,按标准的指数格式输出,共占15列,指数4列,小数6列

2. 表控格式输出语句:

1. print*,<输出表列>: print*,56.8,125 2. write(*,*)<输出表列>

输入语句

输入\\输出三要素:对象;格式;设备. 分类:同上

1. 表控输入语句

1. 自由格式输入-->语句:read*,<输入表列>;read(*,*)<输入表列> 2. 输入数据以逗号或者空格作为间隔 3. 变量名称为输入表

4. 输入的数据应和输入表的变量个数\\类型\\次序严格地一一对应;少了,程序停止,等待继续输入;多了,程序继续进行,多余的不起作用;较多的数据可以几个一组,回车,再输入几个一组,回车...

5. 重复数据,可以7*3---->7,7,7

6. 每一个read(*,*)和write(*,*)语句从一个新的记录[以回车结束的一批输入\\输出数据]开始读数\\输出

1. 例如:read(*,*) A,B,C 2. read(*,*) D,I,J

3. 输入: 2.3,-63.5[回车] 4. 6.4,91.0[回车] 5. 5,8[回车]

6. 结果: A=2.3,B=-63.5,C=6.4, 7. [从新记录开始读数] 8. D=5.0,I=8,J未被赋值 PARAMETER语句

作用:将程序中经常用到的参数或字符串定义成一个符号常量,其值不可改变. 语句:parameter(p1=c1,p2=c2,...,pn=cn) 注意:

1. 符号常量的命名规则与变量名相同,但在程序中其值不可改变,也不能赋值; 2. 符号变量也有类型,可用前面的三种类型说明方法说明类型;

3. 参数语句是非执行语句,位于所有可执行语句的前面,单位与类型说明语句的后面; 4. 一条语句可以定义多个符号常量; 5. 优点:方便修改程序 END,STOP,PAUSE语句

END语句:结束标志,有且仅有一条

PAUSE[n]语句:暂定执行;用于调试程序,n可以是一个字符串或不超过5位的数 STOP[n]语句:停止运行语句;用于调试程序,n可以是一个字符串或不超过5位的数

2.3关于Fortran77的一些有用的简介[3]:

逻辑运算和选择结构 ? 关系表达式

1. 构成选择判断的基本式子

2. 关系运算符:

1. .GT.[greater than] > 2. .GE.[greater than or equal to] >= 3. .LT.[limiter than] < 4. .LE.[limiter than or equal to] <= 5. .EQ.[equal to] = 6. .NE.[not equal to] ≠

3. 一般形式:<算术量或者算术表达式><关系运算符><算术量或者算术表达式> 4. 运算结果:逻辑值:真[.TRUE.]\\假[.FALSE.] 5. 运算顺序:算术运算>关系运算

? 逻辑表达式

1. 运算符:

1. .and. 2. .or. 3. .not.

4. .eqv.逻辑等 5. .neqv.逻辑不等

2. 一般形式:<逻辑变量\\逻辑常量\\关系表达式><逻辑运算符><逻辑变量\\逻辑常量\\关系表达式>

3. 结果:逻辑值:真[.TRUE.]\\假[.FALSE.]

4. 运算顺序:算术运算--->关系运算--->逻辑运算

5. 逻辑运算优先级:.not.--->.and.--->.or.--->.eqv.--->.neqv.

关于Fortran77的一些有用的简介[4]: IF类选择结构

? 用块IF实现选择结构:三种典型形式

1. 基本形式

1. IF(条件) THEN (块IF语句) 2. 块1 (THEN块) 3. ELSE (ELSE语句) 4. 块2 (ELSE块) 5. ENDIF (ENDIF语句)

6. 说明:IF...THEN语句为块IF结构的入口语句;ENDIF语句为出口语句,必须一一对应,

配对使用

2. 简单结构

1. IF(条件) THEN 2. 块 3. ENDIF

4. 说明:没有else块 3. 嵌套结构

1. IF( ) THEN 2. 块1

3. ELSE IF( ) THEN 4. 块2 5. ...

6. ELSE IF( ) THEN 7. 块n

8. [ELSE 块n+1]

9. ENDIF ? 逻辑IF语句 只用一行表示一种选择结构,当且仅当条件成立时执行,并且只执行一条语句; IF(条件) 语句 ? 算术IF语句 IF<算术表达式> N1,N2,N3 当算术表达式的值 <0执行标号为N1的语句; =0执行标号为N2的语句; >0执行标号为N1的语句; 关于Fortran77的一些有用的简介[5]: 循环结构 ? 结构形式:循环体[由一些可执行的语句组成]+循环控制语句[控制循环的开始和结束] ? 分类:条件型循环和计数型循环[DO循环] GOTO语句实现循环 ? 一般形式:GOTO 其中:S1为语句标号 ? 功能:程序执行到此语句时,无条件的转向标号为S1的语句 DO语句实现循环 ? 当循环的初值\\终值\\循环次数都已知时,可用; ? 组成:一个DO语句和循环体组成 ? 一般形式: DO S1 I=E1,E2 [,E3] ...... S1 <终端语句> 例如 DO 10 I=1,19,2 SUM=SUM+1 S1 CONTINUE DO I=1,19,2 SUM=SUM+1 ENDDO 说明

1. I为循环变量,S1为语句标号,是本程序单位中另一可执行语句的标号; 2. 步长可以省略,缺省值=1;

3. 循环初值[E1],终值[E2]和步长[E3]都可以是常量\\变量\\表达式; 4. 由于实数在内存中存储的误差,I,E1,E2,E3尽可能用振型量 5. E1,E2,E3都可正可负,E1,E2,可为0,但是E3不能为0.

? 具体执行过程

1. 执行DO语句,首先计算表达式E1,E2,E3的值,若他们的类型与循环变量I不一致,则自动转换成循环变量的类型

2. 将E1的值赋予循环变量I,及执行赋值语句:I=E1;

?

DO I=E1,E2 [,E3] ......ENDDO DO 10 I=E1,19 ,2 10 SUM=SUM+1 计算循环次数:R=MAX0(E2-E1+E3)/E3,MAX0表示从多个整型变量中取最大的一个;

检查循环次数:若R=0则不执行循环体内的语句,跳出循环;R≠0 则执行循环体内的语句 执行循环终端语句:I=I+E3,即是循环变量获得一个新值,而循环次数R自动减1; 返回步骤4,继续执行,直到R=0. ? CONTINUE语句

循环终端语句必须是可执行语句;那么,这种作为循环终端的语句具有双重作用:一是作为循环终端的标志;而是要完成自身的功能.因此影响了程序的可读性.FORTRAN用一个专门的语句作为DO循环的终端语句,即CONTINUE语句.它自身没有任何功能.

? 一些规定

1. 循环变量在循环体内只能被引用,不能被赋值;

2. 在执行DO循环体期间,E1,E2,E3的值不能被改变,因为他们决定了循环的次数 3. 离开DO循环后,循环变量可以在循环体外被引用,它的值为脱离循环时最后一次

被赋的值;

4. 程序中用到转移语句,规定:只允许从循环体内----->体外;反之不行;

5. 循环终端语句必须是除GOTO,块IF,ENDIF,END和STOP语句外的任何可执行语句 ? DO循环的嵌套

在一个DO循环中还可以包含一个或者多个完整的DO循环,这就是DO循环的嵌套. 一般形式:

DO 10 I=1,10 . . .

DO 20 J=1,10 . . . 20 CONTINUE . . . 10 CONTINUE 说明:

? 嵌套要完整,不能交叉

? 循环变量的名字,规定:并列的循环:循环变量的名字可以相同;嵌套的循环:循环变量的名字不可

以相同

? 若多层循环的结束语句在同一个地方,可以共用一条CONTINUE语句 ? 控制转向语句的使用[体内----->体外] 当型循环的实现

在无法确定循环次数的情况下可以使用当型循环.当型循环是指执行循环体要依据实现给定的条件:当条件成立时执行循环,否则不执行. ? 用DO-WHILE语句实现当型循环 一般形式:

DO S1 [,] WHILE(条件) ... S1 <终端语句>

? 用块IF和GOTO语句实现循环 一般形式:

3. 4. 5. 6.

S1 IF(条件) THEN 块

GOTO S1 ENDIF 直到型循环的实现

所谓直到型循环,是指先执行循环体,再判断条件.如果条件为'假',继续执行循环体,直到条件为'真'时终止循环.

? 用逻辑IF语句实现: S1 循环体

IF(条件) GOTO S1 几种循环形式的关系和比较

? DO循环适用于已知循环次数的情况 ? 几种循环可以互换

DO循环:条件型循环[可用次数作为条件] 当型循环:直到型循环

当型:块IF语句(单边)+GOTO语句(先判断后执行) 直到型:逻辑IF语句+GOTO语句(先执行后判断) ? 各种循环可以相互嵌套

2.4关于Fortran77的一些有用的简介[6]:

数据的输入和输出

数据输入\\输出需要确定的三个基本要素: 输入\\输出的设备 输入\\输出的格式 输入\\输出的数据

系统中隐含的输入\\输出的设备为:键盘\\显示器和打印机

[说明:####(设备,格式)数据列表,当设备显示为*,为默认设备输出,好像是显示器\\或默认设备输入,键盘吧;格式为*,默认格式输出\\输入-david]

有格式的输出

输出语句的一般形式:

WRITE (*,S1) <输出列表> S1 FORMAT(格式说明) 或者

PRINT S1,<输出列表> S1 FORMAT(格式说明)

格式说明符: 主要介绍:I,F,E,D,G,L,A,'(撇号),H,X,r(重复系数),/(斜杠)

I 编辑符(Integer)

作用:用于整型数据的输出.一般形式:Iw或Iw.m

其中:I表示整型输出,w为字段宽度,m表示输出数据的最少数字位数

注意:数据输出时,在指定的区域内向右靠齐;如果数据的实际位数大于指定的字段宽度w,则不输出数据,而在该区域内充满\号;当m大于数据的实际位数时,前面添0,若小于数据实际位数,则不起作用

F 编辑符(Fixed point number)

作用:用于实数的小数形式输出,一般形式:Fw.d

其中:F表示实数的小数形式输出;w为字段宽度;d为输出数据的小数位数

E 编辑符(IExponent)

作用:用于实数的指数形式输出,一般形式:Ew.d

其中:E表示实数的指数形式输出;w为字段宽度;d为数字部分的小数位数 注意:指数部分占4列,负号占1列,小数点前为0.如123.45--->0.12345E+03

D 编辑符(Double precision)

作用:用于双精度的指数形式输出,用法和E 编辑符相仿.一般形式w.d

G 编辑符

作用:由系统根据实际数据的大小来决定使用F编辑符还是E编辑符.一般形式:Gw.d

L 编辑符

作用:用于逻辑型数据的输出,一般形式w 其中表示整型输出,w为字段宽度

A 编辑符

作用:用于字符型数据的输出,一般形式:Aw或A

其中:A表示整型输出,w为字段宽度;若不指定,则表示按实际长度输出

' (撇号) 编辑符

作用:用于输出字符常量,即把撇号内的字符串原样输出.

注意:如果输出的字符中包含撇号,则用两个连续的撇号代表一个要输出的撇号

H 编辑符

作用:用于输出字符常量.一般形式:nH

其中:H表示输出字符常量;n为输出字符个数;str为输出的字符串 (较少使用)

X 编辑符

作用:用于输出空格.一般形式:nX

其中:X表示输出空格;n表示输出的空格数

重复系数r

在format语句中,如果出现几个(或者几组)相同的格式编辑符,则可以利用重复系数而只写一个(或者一组)编辑符.

如FORMAT('A=',/,4(5(1X,F4,0),/))

反斜杠/编辑符

作用:结束本记录的输出,开始下一个记录的输出,通常指换行.

WRITE语句和FORMAT语句的相互作用

WRITE语句的输出变量个数与FORMAT语句的编辑符(不含撇号,H和X)个数可以相等,也可以不等;如果编辑符个数多,则剩余的编辑符不起作用;如果变量的个数多,则当编辑符用完后,重新使用该格式说明,当如果格式说明含带重复系数的编辑符组,则格式说明用完后,只有最右面一个带重复系数的编辑符组(包含重复系数)及其右面的编辑符被重复使用. 可以有空格式说明,如FORMAT(),用于输出一个空行.

有格式的输入

有格式的输入语句

一般形式:

READ(*,S1) <输入列表>

S1 FORMAT(格式说明[由各种格式编辑符构成]) 例如:

READ(*,100) A,B,C 100 FORMAT(F5.1,E12.2,F7.2) END

键盘输入:_15.7_2345.67E+04_705.83enter

在PRINT\\WRITE\\READ语句中包含格式说明

例如: PRINT 100,K,Y 100 FORMAT(18,F7.2) 也可以写成: PRINT'(18,F7.2)',K,Y 注意写法: '(...)' 关于Fortran77的一些有用的简介[7]: 数组

使用原则:\先声明,后使用\ 说明方法:

? 用类型说明语句(显式说明) ? 用DIMENSION语句(隐式说明)

数组的说明和数组元素的引用 ? 用类型语句说明数组

1. 一般形式: 类型说明 数组说明符

2. 其中:数组说明符的一般形式为:数组名(维数说明符,...)[维数说明符,由\下标下界:下标上界\组成

3. 例如:REAL X(1:10),W(1:2,1:3),K(10:20)或者INTEGER B(1:100),PY(0:2,0:3,0:5)

? 用DIMENSION语句说明数组

1. 一般形式:DIMENSION 数组说明符,...

? 说明:

1. 在数组说明符中,维数说明符(下标)的个数成为数组的维数

2. 维数说明符只能使用整型常量或者整型符号常量表达式:如 PARANETER(I=1,J=10) REAL KX(I:J+5)

3. 维数说明符的下标下界为1时,可以省略.如:REAL X(1:10)---->REAL X(10) 4. 数组说明语句必须写在所有可执行的语句之前.[属于非执行语句]

? 数组元素的引用

1. 一般形式:数组名(下标,...)

2. 即要有确定的数组名和下标值,如XN(5),W(1,3),KW(1,2,3)

3. 引用数组元素时,下标可用算术表达式,如果算术表达式的值为实行,则自动取整.

数组的逻辑结构和存储结构

逻辑结构:数组所表示的实际数据结构

存储结构:数组在机器内存储时的排列结构 ? 一维数组

逻辑结构:依次排列的一串数据

存储结构:一组连续存放的一列数据块 ? 二维数组

逻辑结构:一张二维数据表

存储结构:一组按列连续存放的数据块

? 三维数组

逻辑结构:若干张二维数据表

存储结构:一组按页连续存放的数据块 数组的输入和输出

三种方式:用DO循环\\用隐含DO循环\\用数组名

? 用DO循环实现数组的输入输出

1. 优点:数组元素的输入输出次序可由用户控制 2. 缺点:做一次循环就换行输入或输出 ? 用隐含DO循环实现数组的输入输出

1. 优点:既能控制数组元素的输入输出顺序,又能控制一行内输入输出数据的个数 2. 例如:READ(*,*) ((G(I,J),J=1,3),I=1,2),由于是一个READ语句,所以既可以一行输入,

也可以多行输入,键盘输入如下:86,75,72[enter]87,70,83[enter]

3. 注意:一个READ语句可以多行输入;但是多个READ语句时,每一个READ语句必须

从心的一行读数.

用数据名进行数组的输入输出

使用时,其顺序要与数组元素在机器内的存储顺序一致 例如:

DIMENSION K(5)

READ *,K----------对数组进行整体操纵 等价于:READ*,K(1),K(2),K(3),K(4),K(5) 也等价于:READ*,(K(I),I=1,5) 使用DATA语句给数组赋初值

? 一般形式: DATA 变量列表\\初值表\\,变量列表\\初值表......

? 功能:在程序编译期间给变量或者数组赋初值.其中,变量列表可以是变量名\\数组名\\数组元素\\

隐DO循法;初值表只能是常量,不允许出现任何形式的表达式

? 例如: DATA A, B/7.85,9.1/[代表赋初值A=7.85,B=9.1--david], I,J, K /5,10,15/[代表赋初值

I=5,J=10,K=15--david], ? 例如: DIMENSION K(2,3)

DATA((K(I,J),J=1,3),I=1,2) /90,23,20,42,14,32/--------初值列表[2维3列, I=1: 90 23 20

I=2: 42 14 32 或

DATA K/90,42,23,14,20,32/----排列为按列排,排满一列之后,再排下一列; 90 23 20 --------david

42 14 32 ? 例如:

DIMENSION A(10)

DATA A/10*1.0/(表示'10个1.0')

注意:DATA语句属于说明语句,但是它可以放在END语句之前的任意行;当程序中有多个DATA语句给同一个变量赋初值时,以最后一条为准;程序在编译期间给变量赋予初值,在程序执行期间,DATA语句不起任何作用!

2.5关于Fortran77的一些有用的简介[8]

子程序 FORTRAN子程序:包括函数子程序,子例行程序,数据块子程序 执行:从主程序开始执行,遇到调用语句再执行相应的子程序. 不同类型的子程序,关键字不同,调用方法也不同 ? 函数子程序:一种可以作为函数来调用的子程序(\外部函数\1. 定义:一般形式,由FUNCTION语句和子程序体组成 类型说明 FUNCTION 函数名([虚参表]) .........(子程序体) END 2. 1. FUNCTION语句:是函数子程序的第一条语句,标志着该函数子程序的开始

1. 类型说明 FUNCTION 函数名([虚参表]) 2. 注意: 1. 虚元也有类型,需在子程序体中说明 1. 例如:REAL FUNCTION INTEP(X1,X2,X3) 2. INTEGER X1,X2,X3 2. 函数名的命名规则和类型都和变量相同 3. 虚参可以是变量名\\数组名\\子程序名,但不允许用常量和数组元素,它表示了函数自变量的个数\\顺序和类型. 2. 子程序体:完成一个具体任务的程序段 3. 注意: 1. 若无虚参时,括号不能省 2. 函数子程序中所有变量和标号(除函数名和虚参外),与其他程序单位无任何关系 FUNCTION 函数名([虚参表]) 类型说明 函数名 ......(子程序体) END 3. 函数体的说明部分包括对虚参和本函数体内所用变量和数组的说明

4. 函数体中可设置一条或者多条RETURN语句,表示执行到此语句时返回调用程序.

1. 当RETURN语句和END语句紧挨着的时候,可省略RETURN语句 2. 也可以不设RETURN语句,但需从中间返回时,必须设置RETURN语句

5. 函数名的作用:函数名在函数体中一定要被赋值,因为函数名把函数值带回调用程序.

? 函数子程序的调用

1. 一般形式:调用方式和内部函数相似: 函数名(实参数) 或函数名( ) 2. 说明:

1. 调用程序中函数名必须与函数子程序中定义的函数名相同

2. 实参与虚参在个数\\类型\\位置上必须一一对应,但名字可以不同

3. 当虚参是变量名的时候,实参可以是常量\\变量\\数组元素或者表达式;但是当虚参

要在函数体中被赋予初值的时候,则实参不可以是常量或者表达式[因为两者共用一个存储单元]

4. 函数子程序是一个独立的程序单位,可以作为一个单独的源程序进行存贮和编译,

并与调用程序连编后才能起作用.

4 材料本构的相关力学知识

虽然目前要做的工作就是用UMAT把一个新的材料的本构加进ABAQUS中调用,但是我并不想仅仅局限于此.在这个过程中,我会把自己遇到的所有相关问题学懂弄透再去编这个程序.所以估计我这个帖子会很长滴!我力学知识不是很好,所以这两天要好好学习一下力学知识,主要是弹塑性力学吧!

说说弹塑性力学----1

弹性力学\\塑性力学\\弹塑性力学

弹性力学和塑性力学时固体力学的两个重要分支.

1. 固体力学:研究固体材料及其构成的物体结构在外部干扰(载荷\\温度\\变化等)下的力学响应的

科学.按不同的研究对象区分为不同的学科分支.

2. 弹性力学:研究固体材料及由其构成的物体结构在弹性变形阶段的力学行为,包括外部干扰下弹

性物体的内力[应力\\,变形[应变]和位移的很不,以及与之相关的原理\\理论和方法. 3. 塑性力学:则研究他们在塑性变形阶段的力学响应.

4. 弹性和塑性的区别与联系:大多数材料都同时具有弹性和塑性性质,当外载较小时,材料呈现为

弹性的或者基本弹性的;当荷载渐渐增加时,材料将进入塑性变形阶段,即材料的行为呈现塑性的.所谓弹性和塑性,只是材料力学性质的流变学分类法中两个典型性质或理想模型;同意种材料在不同条件下可以主要表现为弹性的或塑性的.因此,所谓弹性材料或弹性物体是指在一定条件下主要呈现弹性性质的材料或物体.塑性材料或者塑性无私的含义与此相类.

5. 弹塑性材料:大多数材料往往都同时具有弹性和塑性性质,特别是在塑性变形阶段,变形中既有

可恢复的弹性变形,又有不可恢复的塑性变形;因此有时又称弹塑性材料

6. 弹性设计方法:是以弹性分析为基础的结构设计,假定材料为理想弹性地,相应地这种设计观点

便以分析结果的实际使用范围作为设计的失效准则,即认为应力[严格地说是应力的某一函数值]达到一定限值[弹性界限],将进入塑性变形阶段时,材料将破坏.

7. 塑性设计方法:结构中如果有一处或一部分材料\破坏\则认为结构失效(丧失所规定的效用).由

于一般的结构都处于非均匀受力状态。当高应力点或高应力区的材料到达弹性界限时、结构的大部分材料仍处于弹性界限之内;而实际材料在应力超过弹性界限以后并不实际发生破坏，仍具有一定的继续承受应力(载荷)的能力，只不过刚度相对地降低。因此弹性设计方法不能充分发挥材料的潜力，导致材料的某种浪费。实际上，当结构内的局部材料进入塑性变形阶段，在

继续增加外载时，结构的达力(应力)分布规律与弹性阶段不同，即所谓内力(应力)重分布;这种重分布总的是使内力(应力)分布更趋均匀，使原来处于低应力区的材料承受更大的应力，从而更好地发挥材料的潜力，提高结构的承载能力。显然，以塑性分析为基础的设计比弹性设计更为优越。但是，塑性设计允许结构有更大的变形，以及完全卸载后结构将存在残余变形。因此，对于刚度要求较高及不允许出现残余变形的场合、这种设计方法不适用。

8. 弹塑性力学的研究对象和方法:是研究结构的强度、刚度和稳定性问题(有时统称为强度问题),

以及结构的“破坏”准则或失效准则.在方法上是在一定的边界条件(或再加上初始条件)下求解三类基本方程:平衡(运动)方程、几何方程和本构〔物理)方程。以实验结果为依据，所得结果由实验来检验.

4.1说说弹塑性力学----2

力学模型的相关知识

'模型'是'原型'的近似描述或表示。建立模型的原则，一是科学性----能尽可能地近似表示原型;二是实用性----能方便地应用。显然，一种科学(力学)模型的建立，要受到科学技术水平的制约。总的来说，力学模型大致有三个层次:材料构造模塑，材料力学性质模型，以及结构计算模型。第一类模型属于基本的，它们属于科学假设范畴。因此，往往以“假设”的形式出现。'模型'有时还与一种理论相对应;因而在有些情况下，'模型'、'假设'和'理论'可以是等义的。 1. 材料构造模型:

1. 连续性假设

2. 假定固体材料是连续介质，即组成物体的质点之间不存在

3. 任何空隙，连续紧密地分布于物体所占的整个空间。由此，我们可以认为，一些物理量如应力，应变和位移等可以表示为坐标的连续函数，从而在作数学推导时可方便地运用连续和极限的概念，事实上，一切物体都是由微粒组成的，都不可能符合这个假设。但可以想象，当微粒尺寸及各微粒之间的距离远比物体的几何尺寸小时。运用这个假设不会引起显著的误差. 4. 均匀及各向同性假设

5. 假设物体由同一类型的均匀材料组成，即物体内各点与各 6. 方向上的物理性质相同(各向同性);物体各部分具有相同的物 7. 理性质.不会随坐标的改变而变化(均匀性)

2. 材料力学性质模型

1. 均弹性材料

2. 弹性材料是对实际固体材料的一种抽象.它构成一个近似于真实材料的理想模型。弹性材料的特征是:物体在变形过程中，对应于一定的温度，应力与应变之间呈一一对应的关系，它和载荷的持续时间及变形历史无关;卸载后，其变形可以完全恢复。在变形过程中，应力与应变之间呈线性规律，即服从胡克(Hooke R)规律的弹性材料，称为线性弹性材抖;而某些金属和塑料等，其应力与应变之间呈非线性性质，称为非线性弹性材料。材料弹性规律的应用，就成为弹性力学区别于其它固体力学分支学科的本质特征。 3. 塑性材转

4. 塑性材料也是固体材料的一种理想模型。塑性材料的特征

5. 是:在变形过程中，应力和应变不再具有一一对应的关系，应变的大小与加载的历史有关但与时间无关;卸载过程中，应力与应变之间按材料固有的弹性规律变化，完全卸载后。物体保持一个永久变形，或称残余变形。变形的不可恢复性是塑性材料的基本特征。 6. 粘性材料

7. 当材料的力学性质具有时间效应，即材料的力学性质与载 8. 荷的待续时间和加载速率相关时，称为粘性材料。实际材料都具有不同程度的枯性性质，

只不过有时可以略去不计。

3. 结构计算模型

1. 小变形假设

2. 假定物体在外部因素作用下所产生的位移远小于物体原来

3. 的尺寸。应用这条假设，可使计算模型大为简化。例如，在研究物体的平衡时，可不考虑由于变形所引起的物体尺寸位置的变化;在建立几何方程和物理方程时，可以略去其中的二次及更高次项，使得到的基本方程是线性偏微分方程组。与之相对立的是大变形情况，这时必须考虑几何关系中的二阶或高阶非线性项，导致变形与载荷之间为非线性关系.得到的基本方程是更难求解的非线性偏微分方程组。 4. 无初应力假设

5. 假定物体原来是处于一种无应力的自然状态。即在外力作

6. 用以前，物体内各点应力均为零。我们的分析计算是从这种状态出发的。 7. (3)荷载分类

8. 作用于物体的外力可以分为体积力和表面力，两者分别简 9. 称为体力和面力。

10. 所谓体力,是分布在物体体积内的力，例如重力和惯性力二物体内各点所受的体力一般是不同的.所谓面力，是分布在物体表面上的力，如风力、流体压力、固体间的接触力等二物体上各点所受的面力一般也是不同的。作用在物体表面

11. 上的力都占有一定的面积;当作用面很小或呈狭长形时.可分别理想化为集中力或线集中力。

4.2说说弹塑性力学----3

1.弹塑性材料

固体材料在受力后产生变形，从变形开始到破坏一般要经历弹性变形和塑性变形这两个阶段。根据材料力学性质的不同，有的弹性阶段较明显，而塑性阶段很不明显，象铸铁等脆性材料，往往经历弹性阶段后就破坏。有的则弹性阶段很不明显，从开始变形就伴随着塑性变形，弹塑性变形总是耦联产生，象混凝土材料就是这洋。而大部分固体材料都呈现出明显的弹性变形阶段和塑性变形阶段。今后我们主要是讨论这种有弹性与塑性变形阶段的固体材料，并统称为弹塑性材料。 2.鲍辛格效应

由于预加塑性拉伸荷载而使压缩屈服应力降低的现象称为Bauschinger效应.正是由于这种效应,塑性变形时一种各向异性的过程,Bauschinger效应是一种由塑性应变引起的特殊的方向各向异性的形式,因为在后继逆向荷载作用下,一个方向的初始塑性变形会减小其反方向的屈服一个应力.在多轴应力情况下,与这种现象对应的是具有不同方向屈服应力之间的相互影响和横向效应,某一方向的预加应变达到塑性范围将会改变其所有方向的屈服应力值.因此Bauschinger效应对于多维问题更重要,包括荷载方向有明显改变的复杂应力历史,比如应力改变符号和循环荷载的情况. 3.弹性变形与塑性变形的区别:

1. 卸除载荷后。变形可以完全恢复，是弹性变形的基本特征，而变形的不可恢复性是塑性变形的

基本特征。弹性与塑性的基本区别不在于它们的应力一应变关系是否线性，例如，在比例极限与弹性极限之间的AB曲线段,应力与应变不再成比例，进入了非线性阶段，但在B点以前卸除载荷，变形仍将完成恢复，属于弹性变形阶段。因此，弹性和塑性的基本区别在于卸载后，是否保留一个永久变形(塑性应变〕.

2. 在弹性变形阶段，应力与应变之间呈一一对应的关系。而在塑性变形阶段，应力与应变之间不

再是单值关系，对应于同一个应力状态，如果加载的历史不同，所又寸应的应变就不同.这并不是说塑性应力和应变状态就不能唯一地确定。为了描述材料在塑性变形阶段的应力一应变关系，我们需要知道:

1. 材料的屈服应力或加载应力。它是用来区别材料是处于弹性阶段还是已进入塑性阶段的特证值。在屈服应力之前，应力一应变服从胡克定律; 2. 加载准则。在材料进入塑性变形阶段后，应力和应变在加载和卸载的情况下服从两个不同的规律，需要有一个判别材料是加载还是卸载的准则，称为加载或卸载准则。在应力等于屈服应力或加载应力时，应力的变化有两种可能、它可写成加载和卸载两种不同的公式形式.

3. 从某个初始状态到现时的全部变形(或加载)历史。对某 4. 现时来说，我们知道了应力增量和应变增量之间的关系如(1. 6)式所示。明确了变形或加载的历史，就可以对增量积分，求得应力全量与应变全量的关系，从而确定该现时材料中的应力和应变。

3. 弹性变形是可逆的，而塑性变形是不可逆的，由于卸载后永久变形的存在，导致在塑性变形中

所做的塑性功也是不可逆的.塑性功恒大于零，是耗散功。所以说弹性变形储存能量[变形恢复时释放能量,不耗散能量],而塑性变形耗散能量,耗散大小为滞回环的面积.

4.3说说弹塑性力学----4

1.简单的拉伸试验

请参照某本主流教材即可

2.残余应力

所谓残余应力，就是对一个处干自然状态的结构施加载荷，又完全卸去载荷后，在结构内存在的、自我平衡的应力(没有外载时满足平衡条件的应力)。而残余应变则是载荷完全卸去后，结构仍保留的变形.前面己指出，弹性变形是可逆过程，当加上载荷又卸去之后，结构将回到初始状态、不会存在残余应力和残余变形。由此可见，只有当结构内发生塑性变形(即使是结构的一部分)之后，才可能出现残余应力和残余变形。结构内存在残余应力的必要条件是结构已发生塑性变形，并且已发生的塑性变形不能满足几何连续条件。

3 子例行程序[ -------UMAT的类型 ]

1. 与函数子程序的区别

1. 名字的作用不同:函数子程序除了供调用外,还代表函数值;子例行程序只能被调

用;

2. 要求返回值的方式不同:函数子程序是通过函数名将函数值带回调用程序;子例行

程序是通过'虚实结合'将其新值转送回调用程序

3. 子例行程序可以带回一批值(数组),或者完成特定的操作,如交换\\排序等;

4. 调用方式不同:函数子程序的调用出现在表达式中,而子例行程序必须用一条独立

的语句调用.

2. 定义:

1. 一般形式: SUBROUTINE 子例行程序名(虚参表)

1. ...... 2. END

2. 说明:

1. 命名规则同变量名[但不用标示类型] 2. 如果没有虚参,括号可以省略

3. 虚参可以是变量名\\数组名\\子程序名,但不允许是常量和数组元素 4. 它的说明部分包括对虚参和本子例行程序体所有变量和数组的说明

5. 子例行程序的名字只起标识作用,没有值的概念,仅仅为了调用

3. 调用:

1. 一般形式:CALL 子例行程序名 (实参表) 或 CALL 子例行程序名 4. 实参和虚参之间的数据传递

虚参可以是变量名\\数组名\\子程序名\\函数名和星号*

1. 变量名作为虚参时的'虚实结合'

1. 对应的实参:同类型的常量\\变量\\数组元素 2. 结合方式:按地址结合,两者共用一个存储单元 3. 注意:

1. 当实参是数组元素时,虚实结合的过程和变量名相同

2. 如果虚参是字符变量,则其定义长度必须<=对应实参的长度,或用'*',表示长度不定;当调用子程序时,具有不定长度的虚参将自动定义为与对应实参相同的长度

2. 数组名作为虚参时的'虚实结合'

1. 对应的实参:同类型的数组名或者数组元素 2. 结合方式:

1. 实参为数组名时:按地址结合,即实参数组的第一个元素与对应虚参数组的第一个元素结合;实参数组的第二个元素与对应虚参数组的第二个元素结合,...,以此类推;

2. 实参为数组元素时:仍按地址结合,当该数组元素的第一个元素与对应虚参数组的第一个元素结合;该数组元素的下一个元素与对应虚参数组的第二个元素结合,...,以此类推;

3. 注意:虚参数组的最后一个元素必须落在实参数组的范围内,否则会出现'超界错误'

3. 虚参是可调数组

1. 可调数组:只能在子程序中使用,其上\\下界可以是整型虚参数变量,其值通过对应的实参传递. 2. 注意:

1. 可调数组只能出现在子程序中,不能出现在主程序中 2. 可调数组的数组名和界都必须作为虚参出现在虚参表中 3. 虚参数组的最后一维的上届可以用\表示

4. 如果在实参表中出现内部函数名时,必须在条用程序的说明部分用INTRISIC语句说明

5. 在调用程序单位中,如果实参中出现了函数子程序或者子例行程序时,必须在调用程序单位的说明部分用EXTERNAL语句说明这些名字.

4.4说说弹塑性力学----5

应力状态理论

一个在外界因素作用下的物体将产生内力和变形。用以描述物体中任何部位的内力和变形特征的力学量是应力和应变。 1.应力概念

凡提到应力，需指明它是对物体内哪一点并过该点的哪一个微分 面。因为通过物体内同一点可以作无数个方位不同的微分面。显

然，各微分面上的应力一般是不相同的。 2.张量概念

由三个正应力,六个剪应力组成的九个应力分量定义了一个新的量。它描述了一点处的应力状态。数学上，在坐标变换时，服从一定坐标变换式的九个数所定义的量叫二阶张量,应力为二阶张量，它称为柯西(Cauchy A L)应力张量，简称为应力张量。张量中的每一个分量为应力张量在某基矢量的坐标系中的分量，简称为应力分量。应力张量常用矩阵形式表示. 应当指出，物体内各点的应力状态一般是不相同的。应为坐标x的函数，所以，应力张量与给定点的空问位置有关，应力张量总是针对物体中的某一确定点而言的.只要知道了一点的九个应力分量，就可求出通过该点的各个微分面上的应力.应力张量完全确定了一点处的应力状态. 3.转轴时应力分量的变换

坐标系作平移变换时,同一点的各应力分量不会改变;显然，转轴后各应力分量都改变了.但这九个量作为一个“整体”所描述的一点的应力状态是不会改变的，因而又一次证明了应力是二阶张量，在坐标转换时具有不变性。在不计体力偶时应力张量具对称性,为对称张量，其独立的应力分量只有六个。

4.主应力和应力不变量

当坐标系转动时，受力物体内任一确定点的九个应力分量 将随着改变。在坐标系不断转动的过程中，必然能找到一个坐标 系，使得该点在该坐标系中只有正应力分量，而剪应力分量为 零。也就是说，对于任一确定的点，总能找到三个互相垂直的微 分面，其上只有正应力而无剪应力。我们把这祥的微分面称为主 微分平面，简称主平面，其法线方向称为应力主方向，而其上的

应力称为主应力。在应力状态的特征方程中,它的三个根即为主应力,按代数值大小一次成为第一主应力,第二主应力和第三主应力.他们是三个不同截面上的应力矢量的模,而不是某个应力矢量的三个分量.状态特征方程的三个系数分别称为应力张量的第一\\第二和第三不变量。其不变的含义是:当坐标系转动时，虽然每个应力分量都随之改变，但这三个量是不变的。更直观地说，因为方程的根代表的是一点的三个主应力，它们的大小与方向在物体的形状及引起内力的因素确定后是完全确定的，即它们是不会随坐标系的改变而改变的。由于应力状态特征方程的根不变，故方程中的系数一定为不变量.以三个主应力为坐标曲线的坐标系称为主坐标系，也称为主向空间一般地说，主坐标系是正交曲线坐标系。

主应力的几个重要性质: 1. 不变性

2. 由于特征方程的系数是不变量，所以作为特征根的主应力及相应的主方向，都是不变量。

从物理意义可知，它们

3. 都是物体内部受外部确定因素作用时客观存在的量，与人为选择参考坐标无关。 4. 实数性

5. 由于应力张量为对称张量，其各元素均为实数，故必有实特征根，即三个主应力都是实数。

这意味着任何应力状态都存在三个主应力. 6. 正交性

7. 当特征方程无重根时，三个主应力必两两正交;当特征方程有一对重根时,如第一和第二主应力

相等,与第三主应力不等,则与第三主应力垂直的平面内任意两个相互垂直的方向均可作为主方向(如双向等拉或等压应力状态);当特征方程出现三重根，任意三个相互正交的方向都可作为主方向。 8. 极值性

9. 在通过同一点的所有微分面的正应力中。最大和最小的正应力是主应力。 5.最大剪应力和八面体应力

弹性理论的适用范围是由材料的屈服条件来确定的。大量实验证明，剪应力对材料进入塑性屈服阶段起决定性作用，例如第三强度理论，又称特雷斯加(Tresca H )屈服条件，是以最大剪应力为材

1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTEN)--[雅可比矩阵[行数,列数],参照上面]

2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1)--[数组维数说明符的下标下届为1时,可省略]

3 PROPS(NPROPS),COORDS(3),DROT(3,3), 4 DFGRD0(3,3),DFGRD1(3,3)

--------------------------------------------以上是变量类型定义

C

DIMENSION EELAS(6),EPLAS(6),FLOW(6)--[弹性应变,塑性应变,流动?]

PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0,HALF=0.5d0)--[PARAMETER语句,指定某些不变量]

DATA NEWTON,TOLER/40,1.D-6/--[DATA语句,指定程序中的某些变量或数组的初值] C

-------------------------------------------用户自定义变量参数

C----------------------------------------------------------- C UMAT FOR JOHNSON-COOK MODEL

C-----------------------------------------------------------

C PROPS(1)-YANG'S MODULUS--[杨氏模量] C PROPS(2)-POISSON RATIO--[泊松比]

C PROPS(3)-INELASTIC HEAT FRACTION--[塑性耗散比,表示塑性功转化成为热量的比例] C PARAMETERS OF JOHNSON-COOK MODEL:--[JOHNSON-COOK模型中的常量声明] C PROPS(4)-A C PROPS(5)-B C PROPS(6)-n C PROPS(7)-C C PROPS(8)-m

C----------------------------------------------------------- C

IF(NDI.NE.3)THEN

WRITE(6,1)--[输出设备6,输出格式1]

1 FORMAT(//,30X,'***ERROR-THIS UMAT MAY ONLY BE USED FOT'--[输出格式\\内容]

1 'ELEMENTS WITH THREE DIRECT STRESS ENDIF-----------------------[终止条件定义]

-----------------------------用户子程序头文件,说明程序一些相关信息

C

C ELASTIC PROPERTIES--[弹性性质] C

EMOD=PROPS(1)--[将弹性模量赋予EMOD]

ENU=PROPS(2)---[将泊松比赋予ENU]

IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499)----------------[定义泊松比的取值] EBULK3=EMOD/(ONE-TWO*ENU))------------------[定义EBULK3=E/(1-2v),为弹性体积膨胀系数K*3,对应弹性力学公式,[弹性体积膨胀系数]体积模量K=1/3*E/(1-2v)] EG2=EMOD/(ONE+ENU)-----------------[定义EG2=E/(1+v)]

EG=EG2/TWO-------------------[剪切弹性模量定义,对应于弹性力学中的剪切弹性模量[拉梅弹性常数之一μ=]G=E/2(1+v)]

EG3=THREE*EG-------------------[--定义EG3=3*EG=3*E/2(1+v)=3G,不知何意??] ELAM=(EBULK3-EG2)/THREE--------------------[

定

义

ELAM=1/3*{E/(1-2v)-E/(1+v)}=K-2G/3,对应于弹性力学中的拉梅弹性常数入]

------------------------------以上是根据已经给出的弹性模量和泊松比求解拉梅弹性常数μ和入

C

C ELASTIC STIFFNESS-------------------[弹性刚度矩阵] C

DO 20 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;20代表此行循环代号]

DO 10 K2=1,NTENS-------------------[DO循环,初值1,终值为传递过来的NTENS,步长为1,省略;10代表此行循环代号]

DDSDDE(K2,K1)=0.0------------------[循环体,将DDSDDE中的全部值赋值为0.0,K2代表

行,K1代表列,按列赋值]

10 CONTINUE-------------------[代号为10的循环执行结束--continue为终端语句,作为循环终端的标志]

20 CONTINUE-------------------[代号为20的循环执行结束--continue为终端语句,作为循环终端的标志] C

DO 40 K1=1,NDI-------------------[DO循环,初值1,终值为NDI[直接应力分量的个数],步长为默认1,省略;40代表此行循环代号]

DO 30 K2=1,NDI-------------------[DO循环,初值1,终值为NDI[直接应力分量的个数],步长为默认1,省略;30代表此行循环代号]

DDSDDE(K2,K1)=ELAM------------------[循环体,将DDSDDE中NDI行NDI列的直接应力分量[对应于正应力分量]全部值赋值为拉梅弹性常数入,K2代表行,K1代表列,按列赋值] 30 CONTINUE-------------------[代号为30的循环执行结束--continue为终端语句,作为循环终端的标志]

DDSDDE(K1,K1)=EG2+ELAM------------------[循环体,将DDSDDE中对角线上的NDI个直接应力分量[对应于正应力分量]全部值赋值为拉梅弹性常数EG2+入=E/(1+v),K2代表行,K1代表列,按列赋值]

40 CONTINUE-------------------[代号为40的循环执行结束--continue为终端语句,作为循环终端的标志]

DO 50 K1=NDI+1,NTENS-------------------[DO循环,初值NDI+1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;50代表此行循环代号]

DDSDDE(K1,K1)=EG------------------[循环体,将DDSDDE中NDI---->NTENS的对角线元素的值赋为EG=G,,按列赋值]

50 CONTINUE-------------------[代号为50的循环执行结束--continue为终端语句,作为循环

终端的标志]

-------------------以上是为了获得弹性本构方程中的弹性本构\\弹性模量矩阵,参见弹性力学相关知识[陈惠发弹塑性力学书P105]

C

C CALCULATE STRESS FROM ELASTIC STRAINS-------------------[以弹性应变计算应力] C

DO 70 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;70代表此行循环代号]

DO 60 K2=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;60代表此行循环代号]

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)-------------------[循环体,三维矩阵形式的本构关系:{应力}=[弹性模量矩阵]{应变};此处的循环很好的实现了,应力列变量的每一个元素,是弹性模量矩阵中对应的一行元素*应变列阵中的一列元素,不错!]

60 CONTINUE-------------------[代号为60的循环执行结束--continue为终端语句,作为循环终端的标志]

70 CONTINUE-------------------[代号为70的循环执行结束--continue为终端语句,作为循环终端的标志]

-------------------以上是通过弹性阶段的本构方程,由所给的应变计算应力,参见弹性力学相关知识[陈惠发弹塑性力学书P104]

C

C RECOVER ELASTIC AND PLASTIC STRAINS-------------------[弹性和塑性应变覆盖/更新--david自己猜的] C

DO 80 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;80代表此行循环代号]

EELAS(K1)=STATEV(K1)+DSTRAN(K1)-------------------[弹性应变=状态变量+应变增量,状态变量矩阵中1-6存储弹性应变]

EPLAS(K1)=STATEV(K1+NTENS)-------------------[状态变量矩阵中的7-12存储塑性应变] 80 CONTINUE

EQPLAS=STATEV(1+2*NTENS)-------------------[状态变量矩阵中的13存储等效塑性应变] C

-------------------状态变量在增量步开始时将数值传递到UMAT中.在增量步结束时必须更新状态变量矩阵中的数据

C CALCULATE MISES STRESS-------------------[计算MISE屈服准则] C

IF(NPROPS.GT.5.AND.PROPS(4).GT.0.0)THEN-------------------[如果(材料常数的个数>5并且材料常数(4)>0.0),那么]

SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2))+ 1(STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3))+ 1(STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1)) DO 90 K1=NDI+1,NTENS

SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) 90 CONTINUE

SMISES=SQRT(SMISES/TWO)-------------------[MISE屈服应力计算] C

-------------------MISE屈服应力计算

C C C

CALL USERHARD SUBROUTINE,GET HARDENING RATE AND YIELD

STRESS-------------------[调用用户硬化子程序,获取硬化率和屈服应力]

CALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4))-------------------[调用用户硬化子程序(虚参列表,SYIEL0??)]

C DETERMINE IF ACTIVELY YIELDING-------------------[确定是否屈服] C

IF(SMISES.GT.(1.0+TOLER)*SYIEL0)THEN C

C MATERIAL RESPONSE IS PLASTIC,DETERMINE FLOW DIRECTION-------------------[条件成立,则材料进入塑性,定义流动法则]

C

SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE ONESY=ONE/SMISES DO 110 K1=1,NDI

FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO) 110 CONTINUE

DO 120 K1=NDI+1,NTENS

FLOW(K1)=STRESS(K1)*ONESY 120 CONTINUE

-------------------流动法则定义,参见弹塑性力学有关书籍

C C

READ

PARAMETERS

OF

JOHNSON-COOK

MODEL-------------------[读取

JOHNSON-COOK模型参数] C

A=PROPS(4) B=PROPS(5) EN=PROPS(6) C=PROPS(7) EM=PROPS(8) C

-------------------读取模型参数,需要用户输入

C NEWTON ITERATION-------------------[NEWTON迭代法--数值迭代] C

SYIELD=SYIEL0

DEQPL=(SMISES-SYIELD)/EG3 DSTRES=TOLER*SYIEL0/EG3

DEQMIN=HALF*DTIME*EXP(1.0D-4/C) DO 130 KEWTON=1,NEWTON DEQPL=MAX(DEQPL,DEQMIN)

CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4) TVP=C*LOG(DEQPL/DTIME) TVP1=TVP+ONE

HARD1=HARD*TVP1+SYIELD*C/DEQPL SYIELD=SYIELD*TVP1

RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD1) IF(ABS(RHS/EG3).LE.DSTRES)GOTO 140 130 CONTINUE

WRITE(6,2)NEWTON

2 FORMAT(//,30X,'***WARNING-PLASTICITY ALGORITHM DID N 1'CONVERGE AFTER',I3,'ITERATIONS') 140 CONTINUE

EFFHRD=EG3*HARD1/(EG3+HARD1) C

-------------------牛顿迭代法定义,参见弹塑性力学有关书籍

C CALCULATE STRESS AND UPDATE STRAINS-------------------[重新计算应力,并更新应变] C

DO 150 K1=1,NDI

STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO

EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO 150 CONTINUE

DO 160 K1=NDI+1,NTENS

STRESS(K1)=FLOW(K1)*SYIELD

EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL 160 CONTINUE

EQPLAS=EQPLAS+DEQPL

SPD=DEQPL*(SYIEL0+SYIELD)/TWO RPL=PROPS(3)*SPD/DTIME C

-------------------以上皆是塑性力学的知识,参见弹塑性力学有关书籍

C JACOBIAN-------------------[计算雅可比矩阵] C

EFFG=EG*SYIELD/SMISES EFFG2=TWO*EFFG

EFFG3=THREE*EFFG2/TWO

EFFLAM=(EBULK3-EFFG2)/THREE DO 220 K1=1,NDI DO 210 K2=1,NDI

DDSDDE(K2,K1)=EFFLAM 210 CONTINUE

DDSDDE(K1,K1)=EFFG2+EFFLAM

220 CONTINUE

DO 230 K1=NDI+1,NTENS DDSDDE(K1,K1)=EFFG 230 CONTINUE

DO 250 K1=1,NTENS DO 240 K2=1,NTENS

DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1) 1*(EFFHRD-EFFG3) 240 CONTINUE 250 CONTINUE ENDIF ENDIF C

C STORE STRAINS IN STATE VARIABLE ARRAY-------------------[存储应变到状态变量矩阵] C

DO 310 K1=1,NTENS STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=EPLAS(K1) 310 CONTINUE

STATEV(1+2*NTENS)=EQPLAS C

RETURN

END-------------------[UMAT退出,增量步结束]

-------------------结束

C C C

SUBROUTINE USERHARD(SYIELD,HARD,EQPLAS,TABLE)-------------------[用户自定义硬化子程序--与UMAT类似,把相关硬化定义编进去即可,跟所要定义的本构有关] C

INCLUDE'ABA_PARAM.INC' C

DIMENSION TABLE(3) C

C GET PARAMETERS,SET HARDENING TO ZERO-------------------[获取常数,设初始硬化值=0] C

A=TABLE(1) B=TABLE(2) EN=TABLE(3) HARD=0.0 C

C CALSULATE CURRENT YIELD STRESS AND HARDENING RATE-------------------[计算当前屈服应力和硬化率]

C

IF(EQPLAS.EQ.0.0)THEN SYIELD=A ELSE

HARD=EN*B*EQPLAS**(EN-1) SYIELD=A+B*EQPLAS**EN END IF RETURN END

新发现的另外一个版本的ＵＭＡＴ－－－源程序+注解！［比我那个好多了 ］

******************************主程序****************************** SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,

2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) C

INCLUDE 'ABA_PARAM.INC' C

CHARACTER*80 CMNAME

DIMENSION STRESS(NTENS),STATEV(NSTATV),

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),DFGRD0(3,3),DFGRD1(3,3) C

DIMENSION EELAS(6),EPLAS(6),FLOW(6),DVECT(4)

REAL DAMAGE,DAMAGC,DAMAG0,BITA,EPDC,EPD0,EQPLAF,DDAMAG,FFDET PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0) DATA NEWTON,TOLER/10,1.D-6/

C

IF(NDI.NE.3) THEN WRITE(6,1)

1 FORMAT(//,30X,'***ERROR-THIS UMAT MAY ONLY BE USED FOR','ELEMENTS & WITH THREE DIRECT STRESS COMPONENTS') ENDIF

WRITE(6,*) NTENS C C

C ELASTIC PROPERTIES

C 根据弹性模量，泊松比计算体积模量，剪切模量，拉梅常数 EMOD=PROPS(1) ENU=PROPS(2) c//避免体积模量无穷大

IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499 c//3*体积模量

EBULK3=EMOD/(ONE-TWO*ENU) EG2=EMOD/(ONE+ENU)

EG=EG2/TWO EG3=THREE*EG

ELAM=(EBULK3-EG2)/THREE C

C ELASTIC STIFFNESS C 弹性刚度，先置零，后赋值

DO 20 K1=1,NTENS DO 10 K2=1,NTENS DDSDDE(K2,K1)=0.0 10 CONTINUE 20 CONTINUE C

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 C

C CALCULATE STRESS FROM ELASTIC STRAINS

C 弹性应变引起的应力

DO 70 K1=1,NTENS DO 60 K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) 60 CONTINUE 70 CONTINUE C

C RECVER ELASTIC AND PLASTIC STRAINS

C 更新弹性和塑性应变，说明STATEV(1~6)是6个方向弹性应变，7~12是塑性应变 C 等效应变EQPLAS是STATEV(13)

DO 80 K1=1,NTENS

EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTENS) 80 CONTINUE

EQPLAS=STATEV(1+2*NTENS)

WRITE(6,1000) KINC, NOEL, NPT, EQPLAS

1000 FORMAT(//,30X,\

& 2X,I3,2X,'NPT',2X,I1,2X,'CURRENT EQ PLASTIC STRAIN',2X,F12.6) C

IF(NPROPS.GT.2.AND.PROPS(3).GT.0.0) THEN

SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2))+ &(STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3))+ &(STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1)) DO 90 K1=NDI+1,NTENS

SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) 90 CONTINUE

SMISES=SQRT(SMISES/TWO) C 调用子程序得到屈服应力

NVALUE=NPROPS/2-1

CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE) C 判断是否屈服

IF(SMISES.GT.(1.0+TOLER)*SYIEL0)THEN C 屈服后的流动方向

SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE ONESY=ONE/SMISES DO 110 K1=1,NDI

FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO) 110 CONTINUE

DO 120 K1=NDI+1,NTENS FLOW(K1)=STRESS(K1)*ONESY

120 CONTINUE

C 牛顿迭代求等效塑性应变

C

SYIELD=SYIEL0 DEQPL=0.0

DO 130 KEWTON=1,NEWTON RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD)

CALL AHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE) IF (ABS(RHS).LT.TOLER*SYIEL0) THEN

C Output incremental equivalent plastic strain WRITE(6,2000) KINC,NOEL,NPT,DEQPL

2000 FORMAT(//,30X,\ & 'NPT',2X,I1,2X,\ GOTO 140 ENDIF

130 CONTINUE

WRITE(6,2) NEWTON

2 FORMAT(//,30X,'***WARNING-PLASTICITY ALGORITHM DID NOT', &'CONVERGE AFTER',I3,'ITERATIONS') 140 CONTINUE

EFFHRD=EG3*HARD/(EG3+HARD) C

C计算应力更新应变

DO 150 K1=1,NDI

STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO

EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO 150 CONTINUE

EQPLAS=EQPLAS+DEQPL

SPD=DEQPL*(SYIEL0+SYIELD)/TWO C更新J矩阵

EFFG=EG*SYIELD/SMISES EFFG2=TWO*EFFG

EFFG3=THREE*EFFG2/TWO

EFFLAM=(EBULK3-EFFG2)/THREE DO 220 K1=1,NDI DO 210 K2=1,NDI

DDSDDE(K2,K1)=EFFLAM 210 CONTINUE

DDSDDE(K1,K1)=EFFG2+EFFLAM 220 CONTINUE

DO 230 K1=NDI+1,NTENS DDSDDE(K1,K1)=EFFG

230 CONTINUE

DO 250 K1=1,NTENS DO 240 K2=1,NTENS

DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1) & *(EFFHRD-EFFG3) 240 CONTINUE 250 CONTINUE ENDIF ENDIF C储存状态变量

DO 310 K1=1,NTENS STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=EPLAS(K1) 310 CONTINUE

STATEV(1+2*NTENS)=EQPLAS

C print equivalent plastic strain to .dat file

WRITE(6,3000) KINC, NOEL,NPT, STATEV(1+2*NTENS)

3000 FORMAT(//,30X,\ & 2X,I3,2X,'NPT',2X,I1,2X,'TOTAL EQ PLASTIC STRAIN',2X,F12.6) C

RETURN END C

C***子程序****ahard***************

SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE) INCLUDE 'ABA_PARAM.INC' DIMENSION TABLE(2,NVALUE)

C 求EQPLAS在哪段斜率内，然后线性叠加求应力，返回SYIELD和斜率HARD SYIELD=TABLE(1,NVALUE) HARD=0.0 C

C IF MORE THAN ONE ENTRY,SEARCH TABLE C

IF(NVALUE.GT.1) THEN DO 10 K1=1,NVALUE-1 EQPL1=TABLE(2,K1+1) IF(EQPLAS.LT.EQPL1) THEN EQPL0=TABLE(2,K1)

IF(EQPL1.LE.EQPL0) THEN WRITE(6,1)

1 FORMAT(//,30X,'***ERROR-PLASTIC STRAIN MUST BE','ENTERED IN ASCEND &ING ORDER') CALL XIT ENDIF

C

C CURRENT YIELD STRESS AND HARDENING C

DEQPL=EQPL1-EQPL0 SYIEL0=TABLE(1,K1) SYIEL1=TABLE(1,K1+1) DSYIEL=SYIEL1-SYIEL0 HARD=DSYIEL/DEQPL

SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD GOTO 20 ENDIF

10 CONTINUE 20 CONTINUE ENDIF RETURN END

雅可比是一个对称矩阵，除非在“*USER MATERIAL”语句中加\参数

楼主好认真呀，给我们提供这么好的帖子，我想请问下楼主，如果雅克比矩阵是非对称矩阵，但是没有设置\参数，这样算还对吗，对计算影响 ...

fengxing 发表于 2010-1-27 16:15

ddsdde是自己定义的一个矩阵，如果不确定是对称的还是非对称的，程序应该是会出现不可预料的结果

本文来源：https://www.bwwdw.com/article/ked7.html