1!------------------------------------------------------------------------------- 2! Copyright (c) 2019 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5!> \brief This module provide functions of contact stiffness calculation 6module m_contact_lib 7 use elementInfo 8 implicit none 9 10 integer, parameter, private :: kreal = kind(0.0d0) 11 integer, parameter, private :: l_max_surface_node = 20 12 integer, parameter, private :: l_max_elem_node = 100 13 14 integer, parameter :: CONTACTUNKNOWN = -1 15 !> contact state definition 16 integer, parameter :: CONTACTFREE = -1 17 integer, parameter :: CONTACTSTICK = 1 18 integer, parameter :: CONTACTSLIP = 2 19 20 !> contact type or algorithm definition 21 integer, parameter :: CONTACTTIED = 1 22 integer, parameter :: CONTACTGLUED = 2 23 integer, parameter :: CONTACTSSLID = 3 24 integer, parameter :: CONTACTFSLID = 4 25 26 !> This structure records contact status 27 type tContactState 28 integer :: state !< -1:free, 1:in contact, or other needed 29 integer :: surface !< contacting surface number 30 real(kind=kreal) :: distance !< penetration value 31 real(kind=kreal) :: wkdist !< copy of penetration value 32 real(kind=kreal) :: lpos(2) !< contact position(local coordinate) 33 real(kind=kreal) :: gpos(3) !< contact position(global coordinate) 34 real(kind=kreal) :: direction(3) !< contact direction 35 real(kind=kreal) :: multiplier(3) !< Lagrangian multiplier or contact force 36 !< 1: normal 2:tangent component 37 real(kind=kreal) :: tangentForce(3) !< friction force 38 real(kind=kreal) :: tangentForce1(3) !< friction force rotated by element(for trial friction force) 39 real(kind=kreal) :: tangentForce_trial(3) !< trial friction force 40 real(kind=kreal) :: tangentForce_final(3) !< final friction force 41 real(kind=kreal) :: reldisp(3) 42 end type 43 44contains 45 46 !> Initializer 47 subroutine contact_state_init(cstate) 48 type(tContactState), intent(inout) :: cstate !< contact state 49 cstate%state = -1 50 cstate%surface = -1 51 end subroutine 52 53 !> Copy 54 subroutine contact_state_copy(cstate1, cstate2) 55 type(tContactState), intent(in) :: cstate1 !< contact state 56 type(tContactState), intent(inout) :: cstate2 !< contact state 57 cstate2 = cstate1 58 end subroutine 59 60 !> Print out contact state 61 subroutine print_contact_state(fnum, cstate) 62 integer, intent(in) :: fnum !< file number 63 type(tContactState), intent(in) :: cstate !< contact state 64 write(fnum, *) "--Contact state=",cstate%state 65 write(fnum, *) cstate%surface, cstate%distance 66 write(fnum, *) cstate%lpos 67 write(fnum, *) cstate%direction 68 write(fnum, *) cstate%multiplier 69 end subroutine 70 71 !> Transfer contact condition int mpc bundary conditions 72 subroutine contact2mpcval( cstate, etype, nnode, mpcval ) 73 type(tContactState), intent(in) :: cstate !< contact state 74 integer, intent(in) :: etype !< type of contacting surface 75 integer, intent(in) :: nnode !< number of elemental nodes 76 real(kind=kreal), intent(out) :: mpcval(nnode*3 + 4) !< MPC constraint 77 78 integer :: i,j 79 real(kind=kreal) :: shapefunc(nnode) 80 81 call getShapeFunc( etype, cstate%lpos(:), shapefunc ) 82 mpcval(1:3) = cstate%direction(1:3) 83 do i=1,nnode 84 do j=1,3 85 mpcval( i*3+j ) = -cstate%direction(j)*shapefunc(i) 86 enddo 87 enddo 88 mpcval( 3*nnode+4 )=cstate%distance 89 end subroutine 90 91 !> This subroutine calculate contact stiff matrix and contact force 92 subroutine contact2stiff( flag, cstate, etype, nnode, ele, mu, mut, & 93 fcoeff, symm, stiff, force ) 94 integer, intent(in) :: flag !< small slid or finite slide 95 type(tContactState), intent(in) :: cstate !< contact state 96 integer, intent(in) :: etype !< type of contacting surface 97 integer, intent(in) :: nnode !< number of elemental nodes 98 real(kind=kreal), intent(in) :: ele(3,nnode) !< coord of surface element 99 real(kind=kreal), intent(in) :: mu !< penalty 100 real(kind=kreal), intent(in) :: mut !< penalty along tangent 101 real(kind=kreal), intent(in) :: fcoeff !< friction coefficient 102 logical, intent(in) :: symm !< symmtricalize 103 real(kind=kreal), intent(out) :: stiff(:,:) !< contact stiffness 104 real(kind=kreal), intent(out) :: force(:) !< contact force 105 106 integer :: i,j 107 real(kind=kreal) :: shapefunc(nnode) 108 real(kind=kreal) :: N(nnode*3+3), dispmat(2,nnode*3+3) 109 real(kind=kreal) :: metric(2,2) 110 real(kind=kreal) :: det, inverse(2,2), ff(2), cff(2) 111 real(kind=kreal) :: dum11(nnode*3+3,nnode*3+3), dum12(nnode*3+3,nnode*3+3) 112 real(kind=kreal) :: dum21(nnode*3+3,nnode*3+3), dum22(nnode*3+3,nnode*3+3) 113 real(kind=kreal) :: tangent(3,2) 114 115 call getShapeFunc( etype, cstate%lpos(:), shapefunc ) 116 N(1:3) = cstate%direction(1:3) 117 do i=1,nnode 118 N( i*3+1:i*3+3 ) = -shapefunc(i)*cstate%direction(1:3) 119 enddo 120 forall( i=1:nnode*3+3, j=1:nnode*3+3 ) 121 stiff(i,j) = mu* N(i)*N(j) 122 end forall 123 force(1:nnode*3+3) = N(:) 124 125 if( fcoeff/=0.d0 .or. flag==CONTACTFSLID ) & 126 call DispIncreMatrix( cstate%lpos, etype, nnode, ele, tangent, metric, dispmat ) 127 128 ! frictional component 129 if( fcoeff/=0.d0 ) then 130 forall(i=1:nnode*3+3, j=1:nnode*3+3) 131 dum11(i,j) = mut*dispmat(1,i)*dispmat(1,j) 132 dum12(i,j) = mut*dispmat(1,i)*dispmat(2,j) 133 dum21(i,j) = mut*dispmat(2,i)*dispmat(1,j) 134 dum22(i,j) = mut*dispmat(2,i)*dispmat(2,j) 135 end forall 136 stiff(1:nnode*3+3,1:nnode*3+3) & 137 = stiff(1:nnode*3+3,1:nnode*3+3) & 138 + metric(1,1)*dum11 + metric(1,2)*dum12 & 139 + metric(2,1)*dum21 + metric(2,2)*dum22 140 141 if( cstate%state == CONTACTSLIP ) then 142 det = metric(1,1)*metric(2,2)-metric(1,2)*metric(2,1) 143 if( det==0.d0 ) stop "Math error in contact stiff calculation" 144 inverse(1,1) = metric(2,2)/det 145 inverse(2,1) = -metric(2,1)/det 146 inverse(1,2) = -metric(1,2)/det 147 inverse(2,2) = metric(1,1)/det 148 ff(:) = cstate%multiplier(2:3) 149 cff(:) = matmul( inverse, ff ) 150 ff(:) = ff(:)/dsqrt( ff(1)*ff(1)+ff(2)*ff(2) ) 151 cff(:) = cff(:)/dsqrt( cff(1)*cff(1)+cff(2)*cff(2) ) 152 stiff(1:nnode*3+3,1:nnode*3+3) = stiff(1:nnode*3+3,1:nnode*3+3) - & 153 ( cff(1)*ff(1)*metric(1,1)+ cff(2)*ff(1)*metric(1,2) )*dum11 - & 154 ( cff(2)*ff(2)*metric(1,2)+ cff(1)*ff(2)*metric(1,1) )*dum21 - & 155 ( cff(1)*ff(1)*metric(1,2)+ cff(2)*ff(1)*metric(2,2) )*dum12 - & 156 ( cff(2)*ff(2)*metric(2,2)+ cff(1)*ff(2)*metric(1,2) )*dum22 157 endif 158 endif 159 160 end subroutine 161 162 !> This subroutine calculate the metric tensor of a elemental surface 163 subroutine getMetricTensor( pos, etype, ele, tensor ) 164 real(kind=kreal), intent(in) :: pos(2) !< current position(local coordinate) 165 integer, intent(in) :: etype !< surface element type 166 real(kind=kreal), intent(in) :: ele(:,:) !< elemental coordinates 167 real(kind=kreal), intent(out) :: tensor(2,2) !< metric tensor 168 169 integer :: nn 170 real(kind=kreal) :: tangent(3,2) 171 nn= getNumberOfNodes(etype) 172 call TangentBase( etype, nn, pos, ele, tangent ) 173 tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) ) 174 tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) ) 175 tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) ) 176 tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) ) 177 end subroutine 178 179 !> This subroutine calculate the relation between global disp and displacement 180 !> along natural coordinate of master surface supposing penetration is small 181 subroutine DispIncreMatrix( pos, etype, nnode, ele, tangent, tensor, matrix ) 182 real(kind=kreal), intent(in) :: pos(2) !< current position(local coordinate) 183 integer, intent(in) :: etype !< surface element type 184 integer, intent(in) :: nnode !< number of nodes in surface 185 real(kind=kreal), intent(in) :: ele(3,nnode) !< elemental coordinates 186 real(kind=kreal), intent(out) :: tangent(3,2) !< tangent basis 187 real(kind=kreal), intent(out) :: tensor(2,2) !< metric tensor 188 real(kind=kreal), intent(out) :: matrix(2,nnode*3+3) !< relation between local and global disp increment 189 190 integer :: i,j 191 real(kind=kreal) :: det 192 real(kind=kreal) :: shapefunc(nnode), t1(nnode*3+3), t2(nnode*3+3) 193 call TangentBase( etype, nnode, pos, ele, tangent ) 194 tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) ) 195 tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) ) 196 tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) ) 197 tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) ) 198 det = tensor(1,1)*tensor(2,2)-tensor(1,2)*tensor(2,1) 199 if( det==0.d0 ) stop "Error in calculate DispIncreMatrix" 200 ! inverse(1,1) = tensor(2,2)/det 201 ! inverse(1,2) = -tensor(1,2)/det 202 ! inverse(2,1) = -tensor(2,1)/det 203 ! inverse(2,2) = tensor(1,1)/det 204 205 call getShapeFunc( etype, pos(:), shapefunc ) 206 forall( j=1:3 ) 207 t1( j ) = tangent(j,1) 208 t2( j ) = tangent(j,2) 209 end forall 210 forall( i=1:nnode, j=1:3 ) 211 t1( i*3+j ) = -tangent(j,1)*shapefunc(i) 212 t2( i*3+j ) = -tangent(j,2)*shapefunc(i) 213 end forall 214 !matrix( 1:2,: ) = matmul( inverse(:,:), matrix ) 215 matrix(1,:) = (tensor(2,2)*t1(:)-tensor(1,2)*t2(:))/det 216 matrix(2,:) = (tensor(1,1)*t2(:)-tensor(2,1)*t1(:))/det 217 tangent(:,1) = tangent(:,1)/dsqrt(dot_product(tangent(:,1),tangent(:,1))) 218 tangent(:,2) = tangent(:,2)/dsqrt(dot_product(tangent(:,2),tangent(:,2))) 219 end subroutine 220 221 !> This subroutine find the projection of a slave point onto master surface 222 subroutine project_Point2Element(xyz,etype,nn,elemt,reflen,cstate,isin,distclr,ctpos,localclr) 223 real(kind=kreal),intent(in) :: xyz(3) !< Coordinates of a spacial point, whose projecting point is to be computed 224 integer, intent(in) :: etype !< surface element type 225 integer, intent(in) :: nn !< number of elemental nodes 226 real(kind=kreal),intent(in) :: elemt(3,nn) !< nodes coordinates of surface element 227 real(kind=kreal),intent(in) :: reflen !< reference length of surface element 228 type(tContactState),intent(inout) :: cstate !< Recorde of contact information 229 logical, intent(out) :: isin !< in contact or not 230 real(kind=kreal), intent(in) :: distclr !< clearance of contact distance 231 real(kind=kreal), optional :: ctpos(2) !< curr contact position( natural coord ) 232 real(kind=kreal), optional :: localclr !< clearance of contact local coord 233 234 integer :: count,order, initstate 235 real(kind=kreal) :: determ, inverse(2,2) 236 real(kind=kreal) :: sfunc(nn), curv(3,2,2) 237 real(kind=kreal) :: r(2), dr(2), r_tmp(2) ! natural coordinate 238 real(kind=kreal) :: xyz_out(3) ! curr. projection position 239 real(kind=kreal) :: dist_last,dist_now, dxyz(3) ! dist between the point and its projection 240 real(kind=kreal) :: tangent(3,2) ! base vectors in tangent space 241 real(kind=kreal) :: dF(2),d2F(2,2),normal(3) 242 real(kind=kreal),parameter :: eps = 1.0D-8 243 real(kind=kreal) :: clr, tol, factor 244 245 initstate = cstate%state 246 clr = 1.d-4 247 if( present( localclr ) ) clr=localclr 248 if( present( ctpos ) ) then 249 r(:)= ctpos 250 else 251 call getElementCenter( etype, r(:) ) 252 endif 253 254 tol = 1.0D0 255 do count=1,100 256 call getShapeFunc( etype, r, sfunc ) 257 xyz_out = matmul( elemt(:,:), sfunc ) 258 dxyz(1:3) = xyz_out(1:3) - xyz(1:3) 259 dist_last = dot_product( dxyz, dxyz(:) ) 260 261 call TangentBase( etype, nn, r, elemt, tangent ) 262 call Curvature( etype, nn, r, elemt, curv ) 263 264 ! dF(1:2) 265 dF(1:2) = -matmul( dxyz(:), tangent(:,:) ) 266 ! d2F(1:2,1:2) 267 d2F(1,1)= dot_product( tangent(:,1), tangent(:,1) ) - dot_product( dxyz, curv(:,1,1) ) 268 d2F(1,2)= dot_product( tangent(:,1), tangent(:,2) ) - dot_product( dxyz, curv(:,1,2) ) 269 d2F(2,1)= dot_product( tangent(:,2), tangent(:,1) ) - dot_product( dxyz, curv(:,2,1) ) 270 d2F(2,2)= dot_product( tangent(:,2), tangent(:,2) ) - dot_product( dxyz, curv(:,2,2) ) 271 272 ! inverse of d2F 273 determ = d2F(1,1)*d2F(2,2) - d2F(1,2)*d2f(2,1) 274 if( determ==0.d0 ) stop "Math error in contact searching" 275 inverse(1,1) = d2F(2,2) / determ 276 inverse(2,2) = d2F(1,1) / determ 277 inverse(1,2) = -d2F(1,2) / determ 278 inverse(2,1) = -d2F(2,1) / determ 279 dr=matmul(inverse,dF) 280 281 tol=dot_product(dr,dr) 282 if( dsqrt(tol)> 3.d0 ) then ! too far away 283 r= -100.d0; exit 284 endif 285 286 factor = 1.d0 287 do order=1,10 288 r_tmp(1:2) = r(1:2) + factor*dr(1:2) 289 call getShapeFunc( etype, r_tmp, sfunc ) 290 xyz_out(1:3) = matmul( elemt(:,:), sfunc(:) ) 291 dxyz(1:3) = xyz(1:3)-xyz_out(:) 292 dist_now = dot_product( dxyz, dxyz ) 293 if(dist_now <= dist_last) exit 294 factor = factor*0.7D0 295 enddo 296 r(1:2) = r_tmp(1:2) 297 298 if( tol<eps ) exit 299 enddo 300 301 isin = .false. 302 cstate%state = CONTACTFREE 303 if( isInsideElement( etype, r, clr )>=0 ) then 304 dxyz(:)=xyz_out(:)-xyz(:) 305 normal(:) = SurfaceNormal( etype, nn, r, elemt ) 306 normal(:) = normal(:)/dsqrt( dot_product(normal, normal) ) 307 do count = 1,3 308 if( dabs(normal(count))<1.D-10 ) normal(count) =0.d0 309 if( dabs(1.d0-dabs(normal(count)))<1.D-10 ) normal(count) =sign(1.d0, normal(count)) 310 enddo 311 cstate%distance = dot_product( dxyz, normal ) 312 313 if( cstate%distance < distclr*reflen .and. cstate%distance > -5.0d-01*reflen ) isin = .true. 314 315 if( isin ) then 316 if( initstate== CONTACTFREE ) then 317 cstate%state = CONTACTSTICK 318 else 319 cstate%state = initstate 320 endif 321 cstate%gpos(:)=xyz_out(:) 322 cstate%lpos(:)=r(:) 323 cstate%direction(:) = normal(:) 324 cstate%wkdist = cstate%distance 325 endif 326 endif 327 end subroutine project_Point2Element 328 329 !> This subroutine find the projection of a slave point onto master surface 330 subroutine update_TangentForce(etype,nn,elemt0,elemt,cstate) 331 integer, intent(in) :: etype !< surface element type 332 integer, intent(in) :: nn !< number of elemental nodes 333 real(kind=kreal),intent(in) :: elemt0(3,nn) !< nodes coordinates of surface element at t 334 real(kind=kreal),intent(in) :: elemt(3,nn) !< nodes coordinates of surface element at t+dt 335 type(tContactState), intent(inout) :: cstate !< Recorde of contact information 336 337 integer :: i 338 real(kind=kreal) :: tangent0(3,2), tangent(3,2) ! base vectors in tangent space 339 real(kind=kreal) :: coeff(2), norm, norm_tmp 340 real(kind=kreal) :: tangentForce_tmp(3) 341 342 call TangentBase( etype, nn, cstate%lpos, elemt0, tangent0 ) 343 call TangentBase( etype, nn, cstate%lpos, elemt, tangent ) 344 345 !project tangentforce to base vector tangent0 346 do i=1,2 347 coeff(i) = dot_product(cstate%tangentForce(1:3),tangent0(1:3,i)) 348 coeff(i) = coeff(i)/dot_product(tangent0(1:3,i),tangent0(1:3,i)) 349 enddo 350 tangentForce_tmp(1:3) = coeff(1)*tangent0(1:3,1) + coeff(2)*tangent0(1:3,2) 351 norm_tmp = dsqrt(dot_product(tangentForce_tmp,tangentForce_tmp)) 352 !adjust tangent force of slave point which moved over element boundary 353 if( norm_tmp > 1.d-6 ) then 354 norm = dsqrt(dot_product(cstate%tangentForce,cstate%tangentForce)) 355 coeff(1:2) = (norm/norm_tmp)*coeff(1:2) 356 end if 357 358 !set rotated tangentforce to tangentforce1 359 cstate%tangentForce1(1:3) = coeff(1)*tangent(1:3,1) + coeff(2)*tangent(1:3,2) 360 361 end subroutine update_TangentForce 362 363end module m_contact_lib 364 365 366