我想对岩石力学中的应变和应力变化做一次迭代,但我坚持错误:
实::应变(1:N),西格玛(1:N),DSIGMA (1:N),应变(1:n)=0.0
和
实数,维数(6)::CEL(1:n,1:n)!刚度矩阵
!program elastic_plastic
implicit none
!define all parameter
integer :: i = 1.0,j,K,M,N,inc
real::STRAIN (1:N), SIGMA (1:N),DSIGMA (1:N),DSTRAIN (1:N)=0.0
real,Dimension(6)::CEL(1:N,1:N)!stiffness matrix
real:: YOUNG, NU, COHESION !rock properties
real::ALPHA, KAPPA! cohesion and frictional angle
real::F !function
real::FRICTION_DEG, FRICTION_RAD !friction angle
real::VARJ2 ,VARI1 !stress invariants (I1 and J2)(MPa)
real:: LAMBDA,GMODU !lames constant and shear modulus
real::SIGMA_1,SIGMA_2,SIGMA_3 !principle stresses(MPa)
real::SHEAR_4,SHEAR_5,SHEAR_6 !shear stresses
real,parameter::DEG_2_RAD = 0.01745329
!INPUT
NU = 0.25
COHESION = 15 ! in MPa
YOUNG = 20 ! in GPa
FRICTION_DEG = 30.0d0
FRICTION_RAD = FRICTION_DEG *(DEG_2_RAD)
!perform calculations
KAPPA=6.*COHESION *cos(FRICTION_DEG*DEG_2_RAD)/(sqrt(3.)*(3.-sin(FRICTION_DEG*DEG_2_RAD)))
ALPHA=2.*sin(FRICTION_DEG*DEG_2_RAD)/(sqrt(3.)*(3.-sin(FRICTION_DEG*DEG_2_RAD)))
GMODU=YOUNG/2.*(1.+NU)
LAMBDA=NU*YOUNG/((1.+NU)*(1.-(2.*NU)))
!Set up elastic stiffness matrix (CEL)
CEL(1:N,1:N)=0.0
CEL (1,1)= LAMBDA-(2.*GMODU)
CEL (2,2)= LAMBDA-(2.*GMODU)
CEL (3,3)= LAMBDA-(2.*GMODU)
CEL (4,4)= 2.*GMODU
CEL (5,5)= 2.*GMODU
CEL (6,6)= 2.*GMODU
DO
inc = inc + 1
DSTRAIN(1)=0.00002
DSIGMA = matmul (CEL(1:N,1:N), DSTRAIN)
SIGMA =SIGMA +DSIGMA
STRAIN=STRAIN+DSTRAIN
!calculate I1 AND J2
VARI1=SIGMA_1+SIGMA_2+SIGMA_3
VARJ2=1./6.*((SIGMA_1-SIGMA_2)**2+(SIGMA_2-SIGMA_3)**2+(SIGMA_3- SIGMA_1)**2+SHEAR_4**2+SHEAR_5**2+SHEAR_6**2)
!Yield function (Drucker-prager)
F= ALPHA*VARI1+(sqrt(VARJ2)-KAPPA)
IF (F.LE.0.0d0)then !Elastic step (exit)
SIGMA =SIGMA
STRAIN=STRAIN
exit
endif
if (F.GT.0.0d0)then !Plastic step (continue)
goto 20
end if
20 continue
write(11,*)STRAIN,SIGMA,inc
END DO
end发布于 2017-12-21 01:02:58
不能静态地定义带有变量的数组。你必须使用常量。例如,以下内容将起作用:
real::STRAIN (1:5), SIGMA (1:5),DSIGMA (1:5),DSTRAIN (1:5)=0.0
real,Dimension(6)::CEL(1:5,1:5)!stiffness matrix如果在代码时不知道数组的大小,可以使用“分配”语句。这就是所谓的“动态存储分配”。来自Fortran 90/95中的'Arrays和并行编程
声明可分配数组的方法如下:
integer Nparticles ! number of particles
integer, parameter :: dim=3 ! dimensionality of space
...
real, allocatable :: charge(:) ! defines an array containing the charge of
! each particle
integer, allocatable :: xyz(:,:) ! coordinates of each particle一旦读取了模拟中粒子的实际数量,我们就可以分配以下数组:
read(*,*) Nparticles
allocate (charge(Nparticles),xyz(dim,Nparticles))https://stackoverflow.com/questions/47916291
复制相似问题