FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
contact_lib.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
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
16  integer, parameter :: contactfree = -1
17  integer, parameter :: contactstick = 1
18  integer, parameter :: contactslip = 2
19 
21  integer, parameter :: contacttied = 1
22  integer, parameter :: contactglued = 2
23  integer, parameter :: contactsslid = 3
24  integer, parameter :: contactfslid = 4
25 
28  integer :: state
29  integer :: surface
30  real(kind=kreal) :: distance
31  real(kind=kreal) :: wkdist
32  real(kind=kreal) :: lpos(2)
33  real(kind=kreal) :: gpos(3)
34  real(kind=kreal) :: direction(3)
35  real(kind=kreal) :: multiplier(3)
37  real(kind=kreal) :: tangentforce(3)
38  real(kind=kreal) :: tangentforce1(3)
39  real(kind=kreal) :: tangentforce_trial(3)
40  real(kind=kreal) :: tangentforce_final(3)
41  real(kind=kreal) :: reldisp(3)
42  end type
43 
44 contains
45 
47  subroutine contact_state_init(cstate)
48  type(tcontactstate), intent(inout) :: cstate
49  cstate%state = -1
50  cstate%surface = -1
51  end subroutine
52 
54  subroutine contact_state_copy(cstate1, cstate2)
55  type(tcontactstate), intent(in) :: cstate1
56  type(tcontactstate), intent(inout) :: cstate2
57  cstate2 = cstate1
58  end subroutine
59 
61  subroutine print_contact_state(fnum, cstate)
62  integer, intent(in) :: fnum
63  type(tcontactstate), intent(in) :: cstate
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 
72  subroutine contact2mpcval( cstate, etype, nnode, mpcval )
73  type(tcontactstate), intent(in) :: cstate
74  integer, intent(in) :: etype
75  integer, intent(in) :: nnode
76  real(kind=kreal), intent(out) :: mpcval(nnode*3 + 4)
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 
92  subroutine contact2stiff( flag, cstate, etype, nnode, ele, mu, mut, &
93  fcoeff, symm, stiff, force )
94  integer, intent(in) :: flag
95  type(tcontactstate), intent(in) :: cstate
96  integer, intent(in) :: etype
97  integer, intent(in) :: nnode
98  real(kind=kreal), intent(in) :: ele(3,nnode)
99  real(kind=kreal), intent(in) :: mu
100  real(kind=kreal), intent(in) :: mut
101  real(kind=kreal), intent(in) :: fcoeff
102  logical, intent(in) :: symm
103  real(kind=kreal), intent(out) :: stiff(:,:)
104  real(kind=kreal), intent(out) :: 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 
163  subroutine getmetrictensor( pos, etype, ele, tensor )
164  real(kind=kreal), intent(in) :: pos(2)
165  integer, intent(in) :: etype
166  real(kind=kreal), intent(in) :: ele(:,:)
167  real(kind=kreal), intent(out) :: tensor(2,2)
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 
181  subroutine dispincrematrix( pos, etype, nnode, ele, tangent, tensor, matrix )
182  real(kind=kreal), intent(in) :: pos(2)
183  integer, intent(in) :: etype
184  integer, intent(in) :: nnode
185  real(kind=kreal), intent(in) :: ele(3,nnode)
186  real(kind=kreal), intent(out) :: tangent(3,2)
187  real(kind=kreal), intent(out) :: tensor(2,2)
188  real(kind=kreal), intent(out) :: matrix(2,nnode*3+3)
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 
222  subroutine project_point2element(xyz,etype,nn,elemt,reflen,cstate,isin,distclr,ctpos,localclr)
223  real(kind=kreal),intent(in) :: xyz(3)
224  integer, intent(in) :: etype
225  integer, intent(in) :: nn
226  real(kind=kreal),intent(in) :: elemt(3,nn)
227  real(kind=kreal),intent(in) :: reflen
228  type(tcontactstate),intent(inout) :: cstate
229  logical, intent(out) :: isin
230  real(kind=kreal), intent(in) :: distclr
231  real(kind=kreal), optional :: ctpos(2)
232  real(kind=kreal), optional :: localclr
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 
330  subroutine update_tangentforce(etype,nn,elemt0,elemt,cstate)
331  integer, intent(in) :: etype
332  integer, intent(in) :: nn
333  real(kind=kreal),intent(in) :: elemt0(3,nn)
334  real(kind=kreal),intent(in) :: elemt(3,nn)
335  type(tcontactstate), intent(inout) :: cstate
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 
363 end module m_contact_lib
364 
365 
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:126
subroutine getelementcenter(fetype, localcoord)
Return natural coordinate of the center of surface element.
Definition: element.f90:1005
subroutine tangentbase(fetype, nn, localcoord, elecoord, tangent)
Calculate base vector of tangent space of 3d surface.
Definition: element.f90:926
integer function isinsideelement(fetype, localcoord, clearance)
if a point is inside a surface element -1: No; 0: Yes; >0: Node's (vertex) number
Definition: element.f90:1022
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:863
subroutine curvature(fetype, nn, localcoord, elecoord, l2ndderiv, normal, curv)
Calculate curvature tensor at a point along 3d surface.
Definition: element.f90:960
This module provide functions of contact stiffness calculation.
Definition: contact_lib.f90:6
integer, parameter contacttied
contact type or algorithm definition
Definition: contact_lib.f90:21
integer, parameter contactfslid
Definition: contact_lib.f90:24
subroutine print_contact_state(fnum, cstate)
Print out contact state.
Definition: contact_lib.f90:62
subroutine getmetrictensor(pos, etype, ele, tensor)
This subroutine calculate the metric tensor of a elemental surface.
integer, parameter contactglued
Definition: contact_lib.f90:22
integer, parameter contactunknown
Definition: contact_lib.f90:14
subroutine contact2mpcval(cstate, etype, nnode, mpcval)
Transfer contact condition int mpc bundary conditions.
Definition: contact_lib.f90:73
subroutine project_point2element(xyz, etype, nn, elemt, reflen, cstate, isin, distclr, ctpos, localclr)
This subroutine find the projection of a slave point onto master surface.
subroutine update_tangentforce(etype, nn, elemt0, elemt, cstate)
This subroutine find the projection of a slave point onto master surface.
subroutine contact_state_init(cstate)
Initializer.
Definition: contact_lib.f90:48
subroutine contact2stiff(flag, cstate, etype, nnode, ele, mu, mut, fcoeff, symm, stiff, force)
This subroutine calculate contact stiff matrix and contact force.
Definition: contact_lib.f90:94
integer, parameter contactslip
Definition: contact_lib.f90:18
integer, parameter contactsslid
Definition: contact_lib.f90:23
integer, parameter contactfree
contact state definition
Definition: contact_lib.f90:16
subroutine dispincrematrix(pos, etype, nnode, ele, tangent, tensor, matrix)
This subroutine calculate the relation between global disp and displacement along natural coordinate ...
subroutine contact_state_copy(cstate1, cstate2)
Copy.
Definition: contact_lib.f90:55
integer, parameter contactstick
Definition: contact_lib.f90:17
This structure records contact status.
Definition: contact_lib.f90:27