Hello,
I would like to implement a UEL in order to have a 2-node Timoshenko beam in 3 dimensions with 3 Gauss points for integration.
As a first step, I'm implementing elastic behavior.
So I've taken the example “Example 1: Planar Beam Element with Nonlinear Section Behavior”.
I'm not sure about what I did.
I computed the local e_x, e_y, e_z (with e_x // beam axis)
I have the local B matrix between strain ε(x) = (ε_x(x) = ∂ux/∂x (x), χ_y(x) = ∂θy/∂x (x), χ_z(x) = ∂θz/∂x (x), χ_x(x) = ∂θx/∂x (x), ε_xy(x) = ∂uy/∂x (x) − θ_z(x), ε_xz(x) = ∂uz/∂x (x) + θ_y(x) ) and local nodal displacements u_nod = (u_x1, u_y1, u_z1, θ_x1, θ_y1, θ_z1, u_x1, u_y2, u_z2, θ_x2, θ_y2, θ_z2)
So ε(x) = B(x) u_nod and since u_nod = Lambda * U_nod (global coordinates)
ε(x) = B(x)*Lambda * U_nod
I have then incremental strains ε(x) = B(x)*Lambda * dU_nod
I can call then a constitutive behavior law, in Timoshenko beam we have
Hence the constitutive law matrix
I can then obtain the internal forces and moments at all Gauss points, and then I can assemble RHS and Amatrx.
What do you think about it?
Thanks a lot
c
c simple 3-d linear beam element with generalized section properties
c
SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
1PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,
2KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,
3NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,
4PERIOD)
C
INCLUDE 'ABA_PARAM.INC'
c
DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),
1 SVARS(NSVARS),ENERGY(8),PROPS(*),COORDS(MCRD,NNODE),
2 U(NDOFEL),DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),
3 PARAMS(3),JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),
4 DDLMAG(MDLOAD,*),PREDEF(2,NPREDF,NNODE),LFLAGS(*),
5 JPROPS(*)
C
c
DIMENSION b(6, 12), gauss(3), weight(3)
DIMENSION blocal(6,12)
DIMENSION e_x(3), e_y(3), e_z(3)
DIMENSION R(3,3), Lambda(12,12)
DIMENSION Hooke(6)
DIMENSION Hmatrix(6,6)
DIMENSION bt(12,6)
DIMENSION btHmatrix(12,6)
DIMENSION btHmatrixb(12,12)
DIMENSION strain(6)
DIMENSION dstrain(6)
DIMENSION tcohesion(6)
DIMENSION rhselementaire(12)
double precision :: h, w, S, hw3, Iy, Iz, ky, kz, E, nu
c
PARAMETER (ZERO=0.d0, one=1.d0, two=2.d0, three=3.d0, four=4.d0,
1 six=6.d0, eight=8.d0, twelve=12.d0)
c
c
c
c Gauss points for [-1,1] : xi -sqrt(3/5),0,sqrt(3/5), respective weights 5/9, 8/9 et 5/9
c Sor for [0,L] : (1+ξi)/2 with ξi = xi/L
data gauss/.112701665d0, 0.5d0, .8872983346d0/
data weight/.555555556d0, .888888889d0, .555555556d0/
c
c
c props(1) - section height
c props(2) - section width
c props(3) - Young’s modulus
c props(4) - Poisson's ratio
c props(5) - Shear coefficient ky
c props(6) - Shear coefficient kz
h=props(1) ! h = 1
w=props(2) ! w = 1
E=props(3) ! E = 200 000
nu=props(4) ! nu = 0.3
ky=props(5) ! ky = 5/6 = 0.833333
kz=props(6) ! kz = 5/6 = 0.833333
c Rectangular section
S=h*w
print *, "Section", S
hw3 = h*w**3
print *, "hw3", hw3
c Second moment of area, rectangular section
Iy = hw3 / 12.d0 ! Iy = 1/12 = 0.0833
Iz = (w*h**3)/12.d0 ! Iy = 1/12 = 0.0833
print '(A, F20.10)', "Section S: ", S
print '(A, F20.10)', "hw3: ", hw3
print '(A, F20.10)', "Iy, expected 1/12 = 0,0833 : ", Iy
print '(A, F20.10)', "Iz, expected 1/12 = 0,0833 : ", Iz
c Shear modulus
G = E/(2.d0*(1.d0+nu)) ! G = 200 000 / (2*(1+0.3)) = 100 000 / (1.3) = 76923.07692
print '(A, F20.10)', "G: ", G
c calculate length and direction cosines
c
dX=coords(1, 2)-coords(1, 1) ! dX = 100
dY=coords(2, 2)-coords(2, 1) ! dY = 0
dZ=coords(3, 2)-coords(3, 1) ! dZ = 0
dl2=dX**2+dY**2+dZ**2
dl=sqrt(dl2)
c Phi_y, Phi_z
Phiy = 12.d0*E*Iy/(kz*S*G*dl2) ! Phiy = 3.12\\times10^{-4}
Phiz = 12.d0*E*Iz/(ky*S*G*dl2) ! Phiz = 3.12\\times10^{-4}
print *, "Phiy", Phiy
print *, "Phiz", Phiz
e_x = (/ dX/dl, dY/dl, dZ/dl /)
! Pick e_z, usually (0,0,1) if not aligned with e_x
e_z = (/ 0.0d0, 0.0d0, 1.0d0 /)
if (abs(dot_product(e_x, e_z)) > 0.999999d0) then
e_z = (/ 0.0d0, 1.0d0, 0.0d0 /) ! Change in case of alignment
endif
! Compute e_y with cross product
call cross_product(e_z, e_x,e_y)
! Recomputation of e_z to ensure orthogonality
call cross_product(e_x,e_y,e_z)
! Build rotation matrix R
do i = 1, 3
R(i, 1) = e_x(i)
R(i, 2) = e_y(i)
R(i, 3) = e_z(i)
end do
c Build transform matrix Lambda = [R 0 0 0; 0 R 0 0; 0 0 R 0; 0 0 0 R]
! Call subroutine to build Lambda
call construct_lambda(Lambda, R, 6)
c
c initialize rhs and lhs
c
do k1=1, 12
rhs(k1, 1)= zero
do k2=1, 12
amatrx(k1, k2)= zero
end do
end do
c
nsvint=nsvars/3
c
c loop over integration points
c
do kintk=1, 3
g=gauss(kintk)
print *, "Gauss point", g
c
c make b-matrix
c
c Shape functions
c Shape functions
c xi = x / L
c N1 = 1.0d0 - xi dN1_dx = -1.0d0 / L
c N2 = xi dN2_dx = 1.0d0 / L
c
c N3 = (1.0d0 / (1.0d0 + phi_y)) * (2.0d0*xi**3 - 3.0d0*xi**2 - phi_y*xi +(1+phi_y))
c dN3_dx = (1.0d0 / (1.0d0 + phi_y)) * (6.0d0 * xi**2 - 6.0d0 * xi - phi_y) / L
c
c N4 = (L / (1.0d0 + phi_y)) * (-xi**3 + ((4.do+phi_y)/2.0d)*xi**2 - ((2.0d0 + phi_y)/2.0d0)*xi)
c dN4_dx = (L / (1.0d0 + phi_y)) * (-3.0d0*xi**2 + (4.do+phi_y)*xi - ((2.0d0 + phi_y)/2.0d0)) / L
c
c N5 = (1.0d0 / (1.0d0 + phi_y)) * (-2.0d0*xi**3 + 3.0d0*xi**2 + phi_y*xi)
c dN5_dx = (1.0d0 / (1.0d0 + phi_y)) * (-6.0d0 * xi**2 + 6.0d0 * xi + phi_y) / L
c
c N6 = (L / (1.0d0 + phi_y)) * (-xi**3 + ((2.do-phi_y)/2.0d)*xi**2 + ((phi_y)/2.0d0)*xi)
c dN6_dx = (L / (1.0d0 + phi_y)) * (-3.0d0*xi**2 + (2.do-phi_y)*xi + ((phi_y)/2.0d0)) / L
c
c N7 = 6.0d0 / (L * (1.0d0 + phi_y)) * xi * (1.0d0 - xi)
c dN7_dx = 6.0d0 / (L * (1.0d0 + phi_y)) * (1.0d0 - 2.0d0 * xi)/L
c
c N8 = 1.0d0 / (1.0d0 + phi_y) * (3.0d0 * xi**2 - (4.0d0 + phi_y) * xi + (1.0d0 + phi_y))
c dN8_dx = 1.0d0 / (1.0d0 + phi_y) * (6.0d0 * xi / L - (4.0d0 + phi_y) / L)
c
c N9 = -6.0d0 / (L * (1.0d0 + phi_y)) * xi * (1.0d0 - xi)
c dN9_dx = -6.0d0 / (L * (1.0d0 + phi_y)) * (1.0d0 - 2.0d0 * xi)/L
c
c N10 = 1.0d0 / (1.0d0 + phi_y) * (3.0d0 * xi**2 + (-2.0d0 + phi_y) * xi)
c dN10_dx = 1.0d0 / (1.0d0 + phi_y) * (6.0d0 * xi / L + (-2.0d0 + phi_y) / L)
c
c
c Traction compression
c u_x = N1 u_x1 + N2 u_x2
c theta_x = N1 theta_x1 + N2 theta_x2
c
c Flexion in (xOz) plane
c u_z = N3 u_z1 + N4 theta_y1 + N5 u_z2 + N6 theta_y2
c theta_y = N7 u_z1 + N8 theta_y1 + N9 u_z2 + N10 theta_y2
c
c Flexion in (xOz) plane , replace Phi_y with Phi_z Phi_z = 12.do*E*Iz/(ky*S*G*dl2)
c u_y = N3 u_y1 - N4 theta_z1 + N5 u_y2 - N6 theta_z2
c theta_z = -N7 u_y1 + N8 theta_z1 - N9 u_y2 + N10 theta_z2
c
c
c First row epsilon_x = du_x/dx = d/dx(N1 u_x1 + N2 u_x2)
blocal(1,1) = -1.0d0/dl ! N1,x
blocal(1,2) = 0.d0
blocal(1,3) = 0.d0
blocal(1,4) = 0.d0
blocal(1,5) = 0.d0
blocal(1,6) = 0.d0
blocal(1,7) = 1.0d0/dl ! N2,x
blocal(1,8) = 0.d0
blocal(1,9) = 0.d0
blocal(1,10) = 0.d0
blocal(1,11) = 0.d0
blocal(1,12) = 0.d0
c Second row chi_y = dtheta_y/dx = d/dx(N7 u_z1 + N8 theta_y1 + N9 u_z2 + N10 theta_y2)
blocal(2,1) = 0.d0
blocal(2,2) = 0.d0
blocal(2,3) = 6.0d0 / (dl * (1.0d0 + Phiy)) * (1.0d0 - 2.0d0 * g)/dl ! N7,x
blocal(2,4) = 0.d0
blocal(2,5) = 1.0d0 / (1.0d0 + Phiy) * (6.0d0 * g / dl - (4.0d0 + Phiy) / dl) ! N8,x
blocal(2,6) = 0.d0
blocal(2,7) = 0.d0
blocal(2,8) = 0.d0
blocal(2,9) = -6.0d0 / (dl * (1.0d0 + Phiy)) * (1.0d0 - 2.0d0 * g)/dl ! N9,x
blocal(2,10) = 0.d0
blocal(2,11) = 1.0d0 / (1.0d0 + Phiy) * (6.0d0 * g / dl + (-2.0d0 + Phiy) / dl) ! N10,x
blocal(2,12) = 0.d0
c Third row chi_z = dtheta_z/dx = d/dx(-N7 u_y1 + N8 theta_z1 - N9 u_y2 + N10 theta_z2)
blocal(3,1) = 0.d0
blocal(3,2) = -1.d0 * (6.0d0 / (dl * (1.0d0 + Phiz)) * (1.0d0 - 2.0d0 * g)/dl) ! -N7,x
blocal(3,3) = 0.d0
blocal(3,4) = 0.d0
blocal(3,5) = 0.d0
blocal(3,6) = 1.0d0 / (1.0d0 + Phiz) * (6.0d0 * g / dl - (4.0d0 + Phiz) / dl) ! N8,x
blocal(3,7) = 0.d0
blocal(3,8) = -1.d0 * (-6.0d0 / (dl * (1.0d0 + Phiz)) * (1.0d0 - 2.0d0 * g)/dl) ! -N9,x
blocal(3,9) = 0.d0
blocal(3,10) = 0.d0
blocal(3,11) = 0.d0
blocal(3,12) = 1.0d0 / (1.0d0 + Phiz) * (6.0d0 * g / dl + (-2.0d0 + Phiz) / dl) ! N10,x
c Fourth row chi_x = dtheta_x/dx = d/dx(N1 theta_x1 + N2 theta_x2)
blocal(4,1) = 0.d0
blocal(4,2) = 0.d0
blocal(4,3) = 0.d0
blocal(4,4) = -1.0d0/dl ! N1,x
blocal(4,5) = 0.d0
blocal(4,6) = 0.d0
blocal(4,7) = 0.d0
blocal(4,8) = 0.d0
blocal(4,9) = 0.d0
blocal(4,10) = 1.0d0/dl ! N2,x
blocal(4,11) = 0.d0
blocal(4,12) = 0.d0
c Fifth row epsilon_xy = du_y/dx - theta_z = d/dx(N3 u_y1 - N4 theta_z1 + N5 u_y2 - N6 theta_z2) - (-N7 u_y1 + N8 theta_z1 - N9 u_y2 + N10 theta_z2)
blocal(5,1) = 0.d0
blocal(5,2) = (1.0d0 / (1.0d0 + Phiz)) * (6.0d0 * g**2 - 6.0d0 * g -Phiz) / dl + 6.0d0 / (dl * (1.0d0 + Phiz)) * g * (1.0d0 - g) ! N3,x + N7
blocal(5,3) = 0.d0
blocal(5,4) = 0.d0
blocal(5,5) = 0.d0
blocal(5,6) = -1.d0 * ( (dl / (1.0d0 + Phiz)) * (-3.0d0*g**2 + (4.d0+Phiz)*g - ((2.0d0 + Phiz)/2.0d0)) / dl)
blocal(5,6) = blocal(5,6) + -1.d0 *( 1.0d0 / (1.0d0 + Phiz) * (3.0d0 * g**2 - (4.0d0 + Phiz) * g + (1.0d0 + Phiz)) ) ! -(N4,x + N8)
blocal(5,7) = 0.d0
blocal(5,8) = (1.0d0 / (1.0d0 + Phiz)) * (-6.0d0 * g**2 + 6.0d0 * g + Phiz) / dl
blocal(5,8) = blocal(5,8) + -6.0d0 / (dl * (1.0d0 + Phiz)) * g * (1.0d0 - g) ! N5,x + N9
blocal(5,9) = 0.d0
blocal(5,10) = 0.d0
blocal(5,11) = 0.d0
blocal(5,12) = -1.d0 * ( (dl / (1.0d0 + Phiz)) * (-3.0d0*g**2 + (2.d0-Phiz)*g + ((Phiz)/2.0d0)) / dl )
blocal(5,12) = blocal(5,12) + -1.d0 * ( 1.0d0 / (1.0d0 + Phiz) * (3.0d0 * g**2 + (-2.0d0 + Phiz) * g) )! -(N6,x + N10)
c Sixth row epsilon_xz = du_z/dx + theta_y = d/dx(N3 u_z1 + N4 theta_y1 + N5 u_z2 + N6 theta_y2) + N7 u_z1 + N8 theta_y1 + N9 u_z2 + N10 theta_y2
blocal(6,1) = 0.d0
blocal(6,2) = 0.d0
blocal(6,3) = (1.0d0 / (1.0d0 + Phiy)) * (6.0d0 * g**2 - 6.0d0 * g - Phiy) / dl + 6.0d0 / (dl * (1.0d0 + Phiy)) * g * (1.0d0 - g) ! N3,x + N7
blocal(6,4) = 0.d0
blocal(6,5) = (dl / (1.0d0 + Phiy)) * (-3.0d0*g**2 + (4.d0+Phiy)*g - ((2.0d0 + Phiy)/2.0d0)) / dl
blocal(6,5) = blocal(6,5) + 1.0d0 / (1.0d0 + Phiy) * (3.0d0 * g**2 - (4.0d0 + Phiy) * g + (1.0d0 + Phiy)) ! N4,x + N8
blocal(6,6) = 0.d0
blocal(6,7) = 0.d0
blocal(6,8) = 0.d0
blocal(6,9) = (1.0d0 / (1.0d0 + Phiy)) * (-6.0d0 * g**2 + 6.0d0 * g + Phiy) / dl + -6.0d0 / (dl * (1.0d0 + Phiy)) * g * (1.0d0 - g) ! N5,x + N9
blocal(6,10) = 0.d0
blocal(6,11) = (dl / (1.0d0 + Phiy)) * (-3.0d0*g**2 + (2.d0-Phiy)*g + ((Phiy)/2.0d0)) / dl
blocal(6,11) = blocal(6,11) + 1.0d0 / (1.0d0 + Phiy) * (3.0d0 * g**2 + (-2.0d0 + Phiy) * g) ! N6,x + N10
blocal(6,12) = 0.d0
c
c matmul(A, B, C, m, n, p)
c b = blocal*Lambda
print *, "blocal au point de Gauss", g, "est", blocal
call matmul(blocal,Lambda,b,6,12,12)
c print *, "b au point de Gauss", g, "est", b
c We have blocal between epsilon and u_locaux epsilon = blocal*u_locaux
c We have u_locaux = Lambda*U_globa
c Hence B = blocal*Lambda so epsilon = B*U_globaux
c Epsilon = (εx, χy, χz, χx, εxy, εxz)
c calculate (incremental) strains
c
c epsilon = (epsx, chiy, chiz, chix, epsxy, epsxz)T
c depsilon = (depsx, dchiy, dchiz, dchix, depsxy, depsxz)T
do i=1,6
strain(i) = zero
dstrain(i) = zero
end do
do k=1, 12
strain(1)=strain(1)+b(1, k)*u(k) !epsilon_x = B(1,:)*u
dstrain(1)=dstrain(1)+b(1, k)*du(k, 1)
strain(2)=strain(2)+b(2, k)*u(k) !chi_y = B(2,:)*u
dstrain(2)=dstrain(2)+b(2, k)*du(k, 1)
strain(3)=strain(3)+b(3, k)*u(k) !chi_z = B(3,:)*u
dstrain(3)=dstrain(3)+b(3, k)*du(k, 1)
strain(4)=strain(4)+b(4, k)*u(k) !chi_x = B(4,:)*u
dstrain(4)=dstrain(4)+b(4, k)*du(k, 1)
strain(5)=strain(5)+b(5, k)*u(k) !epsilon_xy = B(5,:)*u
dstrain(5)=dstrain(5)+b(5, k)*du(k, 1)
strain(6)=strain(6)+b(6, k)*u(k) !epsilon_xz = B(6,:)*u
dstrain(6)=dstrain(6)+b(6, k)*du(k, 1)
end do
c print *, "dstrain at Gauss point", g, "est", dstrain
c
c call constitutive routine ugenb
c
isvint=1+(kintk-1)*nsvint
do k1=1, 6
tcohesion(k1) = zero !Fn=zero !Mfy=zero !Mfz=zero !Mtx=zero !FVy=zero !FVz=zero
Hooke(k1) = zero
end do
C Constitutive behavior, returns efforts Fn, Mfy, Mfz, Mtx, FVy, FVz
C Also return elastic operator
c First elastic Timoshenko beam
call ugenb(tcohesion,Hooke,Hmatrix,strain,dstrain,svars(isvint), nsvint, props, nprops)
print *, "tcohesion au point de Gauss", g, "est", tcohesion
c assemble rhs and lhs
c
c RHS = L/2 Σ ω_i BT(L(1+xi)/2). H_i.(ε_i-ε_ip) avec H.(ε_i-ε_ip) = tcohesion_i
c RHS ∈ Matrice 12x1
c B ∈ Matrice 6x12, BT ∈ Matrice 12x6
c tcohesion_i ∈ Matrice 6x1
c H ∈ Matrice 6x6
c K = L/2 Σ ω_i BT(L(1+xi)/2). H_i. B(L(1+xi)/2)
call transpose_matrix(b, bt, 6, 12) ! Compute transpose of B, BT
call matrix_vector_product(bt, tcohesion, rhselementaire, 12,6) ! RHS_i = BT_i x tcohesion_i
call matmul(bt,Hmatrix,btHmatrix,12,6,6) ! btHmatrix_i = BT_i x H_i
call matmul(btHmatrix,b,btHmatrixb,12,6,12) ! K_i = BT_i x H_i x B
do k1=1,12
rhs(k1,1) = rhs(k1,1) - hdl*weight(kintk)*rhselementaire(k1) ! RHS_i(j,1) = BT_i*tcohesion_i (j,1) = Σ BT(j,i)*tcohesion(:,1) = Σ B(i,j)*tcohesion(i,1) i ∈ [[1,6]]
do k2=1,12
amatrx(k1,k2) = amatrx(k1,k2)+hdl*weight(kintk)*btHmatrixb(k1,k2) ! K_i(j,k) = BT_i*H_i*B_i (j,k), Or K_T = L/2 Σ ω_i BT(L(1+xi)/2). H_i. B(L(1+xi)/2)
end do
end do
end do
c
return
end
SUBROUTINE ugenb(tcohesion,Hooke,Hmatrix,strain,dstrain,svint,nsvint,props,nprops)
c
INCLUDE 'ABA_PARAM.INC'
c
parameter(zero=0.d0,twelve=12.d0)
c
dimension svint(*),props(*)
dimension tcohesion(6)
dimension Hooke(6)
dimension Hmatrix(6,6)
dimension strain(6)
dimension dstrain(6)
double precision :: h, w, S, hw3, Iy, Iz, ky, kz, E, nu
c
c variables to be defined by the user
c
c bn - axial force
c bm - bending moment
c daxial - current tangent axial stiffness
c dbend - current tangent bending stiffness
c dcoupl - tangent coupling term
c
c variables that may be updated
c
c svint - state variables for this integration point
c
c variables passed in for information
c
c eps - axial strain
c deps - incremental axial strain
c cap - curvature change
c dcap - incremental curvature change
c props - element properties
c nprops - # element properties
c nsvint - # state variables
c
c current assumption
c
c props(1) - section height
c props(2) - section width
c props(3) - Young’s modulus
c
h=props(1)
w=props(2)
E=props(3)
nu=props(4)
ky=props(5)
kz=props(6)
c Rectangular section
S=h*w
hw3 = h*w**3
c Second moment of area, rectangular section
Iy = hw3 / 12.d0 ! Iy = 1/12 = 0.0833
Iz = (w*h**3)/12.d0 ! Iy = 1/12 = 0.0833
c Shear modulus
G = E/(2.d0*(1.d0+nu)) ! G = 200 000 / (2*(1+0.3)) = 100 000 / (1.3) = 76923.07692
c
c formulate linear stiffness
c
c N = ES * epsilon_x
Hooke(1)=E*S
c Mfy = EIy * chi_y
Hooke(2)=E*Iy
c Mfz = EIz * chi_z
Hooke(3)=E*Iz
c Mtx = GC * chi_x, C = I_y + I_z
Hooke(4)=G*(Iy+Iz)
c Vy = kySG * epsilon_y
Hooke(5)=ky*S*G
c Vz = kzSG * epsilon_z
Hooke(6)=kz*S*G
c Here Hmatrix is diagonal
do i = 1, 6
do j = 1,6
Hmatrix(i,j) = 0.d0
end do
Hmatrix(i,i) = Hooke(i)
end do
c
c calculate axial force and moment
c
tcohesion(1)=svint(1)+Hooke(1)*dstrain(1) ! Fn_n+1 = Fn_n + ES*depsilon_x
tcohesion(2)=svint(2)+Hooke(2)*dstrain(2) ! Mfy_n+1 = Mfy_n + EIy * dchi_y
tcohesion(3)=svint(3)+Hooke(3)*dstrain(3) ! Mfz_n+1 = Mfz_n + EIz * dchi_z
tcohesion(4)=svint(4)+Hooke(4)*dstrain(4) ! Mtx_n+1 = Mtx_n + GC * dchi_x
tcohesion(5)=svint(5)+Hooke(5)*dstrain(5) ! FVy_n+1 = FVy_n + kySG * depsilon_y
tcohesion(6)=svint(6)+Hooke(6)*dstrain(6) ! FVy_n+1 = FVy_n + kzSG * depsilon_z
c
c store internal variables
c
svint(1)=tcohesion(1) ! Fn
svint(2)=tcohesion(2) ! Mfy
svint(3)=tcohesion(3) ! Mfz
svint(4)=tcohesion(4) ! Mtx
svint(5)=tcohesion(5) ! FVy
svint(6)=tcohesion(6) ! FVz
svint(7)=strain(1) ! epsilon_x
svint(8)=strain(2) ! chi_y
svint(9)=strain(3) ! chi_z
svint(10)=strain(4) ! chi_x
svint(11)=strain(5) ! epsilon_y
svint(12)=strain(6) ! epsilon_z
c
return
end
subroutine matrix_vector_product(A, x, y, m, n)
implicit none
integer :: m, n
integer :: i, j
double precision, dimension(m, n) :: A
double precision, dimension(n) :: x
double precision, dimension(m) :: y
! Initialisation du vecteur résultat y
y = 0.0d0
do i = 1, m
do j = 1, n
y(i) = y(i) + A(i, j) * x(j)
end do
end do
end subroutine matrix_vector_product
subroutine cross_product(v1, v2, v_cross)
implicit none
double precision, dimension(3) :: v1, v2, v_cross
! Calcul du produit vectoriel
v_cross(1) = v1(2) * v2(3) - v1(3) * v2(2)
v_cross(2) = v1(3) * v2(1) - v1(1) * v2(3)
v_cross(3) = v1(1) * v2(2) - v1(2) * v2(1)
RETURN ! Explicitement indiquer la fin de la subroutine
end subroutine cross_product
subroutine construct_lambda(Lambda, R, ndof)
implicit none
double precision, dimension(3,3) :: R
double precision, dimension(12,12) :: Lambda
integer :: i, j, k, l, ndof
! Initialisation de Lambda à zéro
Lambda = 0.0d0
! Remplissage de Lambda avec des blocs de R
do k = 1, ndof, 3
do i = 1, 3
do j = 1, 3
Lambda(k+i-1, k+j-1) = R(i, j)
end do
end do
end do
end subroutine construct_lambda
subroutine matmul(A, B, C, m, n, p)
implicit none
integer :: m, n, p
integer :: i, j, k
double precision, dimension(m, n) :: A
double precision, dimension(n, p) :: B
double precision, dimension(m, p) :: C
! Initialisation de C à zéro
C = 0.0d0
! Calcul du produit matriciel
do i = 1, m
do j = 1, p
do k = 1, n
C(i, j) = C(i, j) + A(i, k) * B(k, j)
end do
end do
end do
end subroutine matmul
subroutine transpose_matrix(A, AT, m, n)
implicit none
integer :: m, n
integer :: i, j
double precision, dimension(m, n) :: A
double precision, dimension(n, m) :: AT
do i = 1, m
do j = 1, n
AT(j, i) = A(i, j)
end do
end do
end subroutine transpose_matrix
