FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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!-------------------------------------------------------------------------------
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
44contains
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
363end 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