1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine springdamp_n2f(xl,elas,voldl,s,imat,elcon, 20 & ncmat_,ntmat_,nope,iperturb,springarea,nmethod,mi, 21 & reltime,nasym) 22! 23! calculates the contact damping matrix (node-to-face penalty) 24! (User's manual -> theory -> boundary conditions -> 25! node-to-face penalty contact) 26! 27 implicit none 28! 29 integer i,j,imat,ncmat_,ntmat_,k,l,nope,nterms,iflag,i1, 30 & iperturb(*),nmethod,mi(*),nasym 31! 32 real*8 xl(3,10),elas(21),ratio(9),pproj(3),val,shp2(7,9), 33 & al(3),s(60,60),voldl(0:mi(2),10),pl(3,10),xn(3),dm, 34 & c1,c2,c3,c4,elcon(0:ncmat_,ntmat_,*),xm(3),xmu(3,3,10), 35 & dxmu(3,10),dval(3,10),fpu(3,3,10),xi,et,xs2(3,7),xk, 36 & a11,a12,a22,b1(3,10),b2(3,10),dal(3,3,10),qxxy(3),fnl(3), 37 & qxyy(3),dxi(3,10),det(3,10),determinant,c11,c12,c22, 38 & qxyx(3),qyxy(3),springarea(2),dist,tu(3,3,10), 39 & clear,reltime,alnew(3) 40! 41 data iflag /4/ 42! 43! actual positions of the nodes belonging to the contact spring 44! (otherwise no contact force) 45! 46 do i=1,nope 47 do j=1,3 48 pl(j,i)=xl(j,i)+voldl(j,i) 49 enddo 50 enddo 51! 52! contact springs 53! 54 nterms=nope-1 55! 56! vector al connects the actual position of the slave node 57! with its projection on the master face = vec_r (User's 58! manual -> theory -> boundary conditions -> node-to-face 59! penalty contact) 60! 61 do i=1,3 62 pproj(i)=pl(i,nope) 63 enddo 64 call attach_2d(pl,pproj,nterms,ratio,dist,xi,et) 65 do i=1,3 66 al(i)=pl(i,nope)-pproj(i) 67 enddo 68! 69! determining the jacobian vector on the surface 70! 71 if(nterms.eq.8) then 72 call shape8q(xi,et,pl,xm,xs2,shp2,iflag) 73 elseif(nterms.eq.4) then 74 call shape4q(xi,et,pl,xm,xs2,shp2,iflag) 75 elseif(nterms.eq.6) then 76 call shape6tri(xi,et,pl,xm,xs2,shp2,iflag) 77 else 78 call shape3tri(xi,et,pl,xm,xs2,shp2,iflag) 79 endif 80! 81! d xi / d vec_u_j 82! d eta / d vec_u_j 83! 84! dxi and det are determined from the orthogonality 85! condition 86! 87 a11=-(xs2(1,1)*xs2(1,1)+xs2(2,1)*xs2(2,1)+xs2(3,1)*xs2(3,1)) 88 & +al(1)*xs2(1,5)+al(2)*xs2(2,5)+al(3)*xs2(3,5) 89 a12=-(xs2(1,1)*xs2(1,2)+xs2(2,1)*xs2(2,2)+xs2(3,1)*xs2(3,2)) 90 & +al(1)*xs2(1,6)+al(2)*xs2(2,6)+al(3)*xs2(3,6) 91 a22=-(xs2(1,2)*xs2(1,2)+xs2(2,2)*xs2(2,2)+xs2(3,2)*xs2(3,2)) 92 & +al(1)*xs2(1,7)+al(2)*xs2(2,7)+al(3)*xs2(3,7) 93! 94 do i=1,3 95 do j=1,nterms 96 b1(i,j)=shp2(4,j)*xs2(i,1)-shp2(1,j)*al(i) 97 b2(i,j)=shp2(4,j)*xs2(i,2)-shp2(2,j)*al(i) 98 enddo 99 b1(i,nope)=-xs2(i,1) 100 b2(i,nope)=-xs2(i,2) 101 enddo 102! 103 determinant=a11*a22-a12*a12 104 c11=a22/determinant 105 c12=-a12/determinant 106 c22=a11/determinant 107! 108 do i=1,3 109 do j=1,nope 110 dxi(i,j)=c11*b1(i,j)+c12*b2(i,j) 111 det(i,j)=c12*b1(i,j)+c22*b2(i,j) 112 enddo 113 enddo 114! 115! d vec_r / d vec_u_k 116! 117 do i=1,nope 118 do j=1,3 119 do k=1,3 120 dal(j,k,i)=-xs2(j,1)*dxi(k,i)-xs2(j,2)*det(k,i) 121 enddo 122 enddo 123 enddo 124 do i=1,nterms 125 do j=1,3 126 dal(j,j,i)=dal(j,j,i)-shp2(4,i) 127 enddo 128 enddo 129 do j=1,3 130 dal(j,j,nope)=dal(j,j,nope)+1.d0 131 enddo 132! 133! d2q/dxx x dq/dy 134! 135 qxxy(1)=xs2(2,5)*xs2(3,2)-xs2(3,5)*xs2(2,2) 136 qxxy(2)=xs2(3,5)*xs2(1,2)-xs2(1,5)*xs2(3,2) 137 qxxy(3)=xs2(1,5)*xs2(2,2)-xs2(2,5)*xs2(1,2) 138! 139! dq/dx x d2q/dyy 140! 141 qxyy(1)=xs2(2,1)*xs2(3,7)-xs2(3,1)*xs2(2,7) 142 qxyy(2)=xs2(3,1)*xs2(1,7)-xs2(1,1)*xs2(3,7) 143 qxyy(3)=xs2(1,1)*xs2(2,7)-xs2(2,1)*xs2(1,7) 144! 145! Modified by Stefan Sicklinger 146! 147! dq/dx x d2q/dxy 148! 149 qxyx(1)=xs2(2,1)*xs2(3,6)-xs2(3,1)*xs2(2,6) 150 qxyx(2)=xs2(3,1)*xs2(1,6)-xs2(1,1)*xs2(3,6) 151 qxyx(3)=xs2(1,1)*xs2(2,6)-xs2(2,1)*xs2(1,6) 152! 153! d2q/dxy x dq/dy 154! 155 qyxy(1)=xs2(2,6)*xs2(3,2)-xs2(3,6)*xs2(2,2) 156 qyxy(2)=xs2(3,6)*xs2(1,2)-xs2(1,6)*xs2(3,2) 157 qyxy(3)=xs2(1,6)*xs2(2,2)-xs2(2,6)*xs2(1,2) 158! 159! 160! End modifications 161! 162! normal on the surface 163! 164 dm=dsqrt(xm(1)*xm(1)+xm(2)*xm(2)+xm(3)*xm(3)) 165 do i=1,3 166 xn(i)=xm(i)/dm 167 enddo 168! 169! distance from surface along normal (= clearance) 170! 171 clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3) 172 if(nmethod.eq.1) then 173 clear=clear-springarea(2)*(1.d0-reltime) 174 endif 175! 176! representative area: usually the slave surface stored in 177! springarea; however, if no area was assigned because the 178! node does not belong to any element, the master surface 179! is used 180! 181 if(springarea(1).le.0.d0) then 182 if(nterms.eq.3) then 183 springarea(1)=dm/2.d0 184 else 185 springarea(1)=dm*4.d0 186 endif 187 endif 188! 189 elas(2)=-springarea(1)*elcon(5,1,imat) 190 elas(1)=elas(2)*clear 191! 192! contact force 193! 194 do i=1,3 195 fnl(i)=-elas(1)*xn(i) 196 enddo 197! 198! derivatives of the jacobian vector w.r.t. the displacement 199! vectors (d m / d u_k) 200! 201 do k=1,nterms 202 xmu(1,1,k)=0.d0 203 xmu(2,2,k)=0.d0 204 xmu(3,3,k)=0.d0 205 xmu(1,2,k)=shp2(1,k)*xs2(3,2)-shp2(2,k)*xs2(3,1) 206 xmu(2,3,k)=shp2(1,k)*xs2(1,2)-shp2(2,k)*xs2(1,1) 207 xmu(3,1,k)=shp2(1,k)*xs2(2,2)-shp2(2,k)*xs2(2,1) 208 xmu(1,3,k)=-xmu(3,1,k) 209 xmu(2,1,k)=-xmu(1,2,k) 210 xmu(3,2,k)=-xmu(2,3,k) 211 enddo 212 do i=1,3 213 do j=1,3 214 xmu(i,j,nope)=0.d0 215 enddo 216 enddo 217! 218! correction due to change of xi and eta 219! 220 do k=1,nope 221 do j=1,3 222 do i=1,3 223! 224! modified by Stefan Sicklinger 225! 226 xmu(i,j,k)=xmu(i,j,k)+(qxxy(i)+qxyx(i))*dxi(j,k) 227 & +(qxyy(i)+qyxy(i))*det(j,k) 228 enddo 229 enddo 230 enddo 231! 232! derivatives of the size of the jacobian vector w.r.t. the 233! displacement vectors (d ||m||/d u_k) 234! 235 do k=1,nope 236 do i=1,3 237 dxmu(i,k)=xn(1)*xmu(1,i,k)+xn(2)*xmu(2,i,k)+ 238 & xn(3)*xmu(3,i,k) 239 enddo 240! 241! auxiliary variable: (d clear d u_k)*||m|| 242! 243 do i=1,3 244 dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k)+ 245 & al(3)*xmu(3,i,k)-clear*dxmu(i,k)+ 246 & xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k)+xm(3)*dal(3,i,k) 247 enddo 248! 249 enddo 250! 251 c1=1.d0/dm 252 c2=c1*c1 253 c3=elas(2)*c2 254 c4=elas(1)*c1 255! 256! derivatives of the forces w.r.t. the displacement vectors 257! 258 do k=1,nope 259 do j=1,3 260 do i=1,3 261 fpu(i,j,k)=-c3*xm(i)*dval(j,k) 262 & +c4*(xn(i)*dxmu(j,k)-xmu(i,j,k)) 263 enddo 264 enddo 265 enddo 266! 267! tangential damping 268! 269 if(elcon(8,1,imat).gt.0.d0) then 270! 271! stiffness of shear stress versus slip curve 272! 273 xk=elcon(8,1,imat)*elcon(5,1,imat)*springarea(1) 274! 275! calculating the relative displacement between the slave node 276! and its projection on the master surface 277! 278 do i=1,3 279 alnew(i)=voldl(i,nope) 280 do j=1,nterms 281 alnew(i)=alnew(i)-ratio(j)*voldl(i,j) 282 enddo 283 enddo 284! 285! calculating the difference in relative displacement since 286! the start of the increment = vec_s 287! 288 do i=1,3 289 al(i)=alnew(i) 290 enddo 291! 292! s=||vec_s|| 293! 294 val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3) 295! 296! d vec_s/ d vec_u_k 297! notice: xi & et are const. 298! 299 do k=1,nope 300 do i=1,3 301 do j=1,3 302 dal(i,j,k)=0.d0 303 enddo 304 enddo 305 enddo 306! 307 do i=1,nterms 308 do j=1,3 309 dal(j,j,i)=-shp2(4,i) 310 enddo 311 enddo 312! 313 do j=1,3 314 dal(j,j,nope)=1.d0 315 enddo 316! 317! d s/ dvec_u_k.||m|| 318! 319 do k=1,nope 320 do i=1,3 321 dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k) 322 & +al(3)*xmu(3,i,k)-val*dxmu(i,k) 323 & +xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k) 324 & +xm(3)*dal(3,i,k) 325 enddo 326 enddo 327! 328! d vec_t/d vec_u_k 329! 330 do k=1,nope 331 do j=1,3 332 do i=1,3 333 tu(i,j,k)=dal(i,j,k) 334 & -c1*(xn(i)*(dval(j,k)-val*dxmu(j,k)) 335 & +val*xmu(i,j,k)) 336 enddo 337 enddo 338 enddo 339! 340! damping matrix 341! 342 do k=1,nope 343 do j=1,3 344 do i=1,3 345 fpu(i,j,k)=fpu(i,j,k)+xk*tu(i,j,k) 346 enddo 347 enddo 348 enddo 349 endif 350! 351! determining the damping matrix contributions 352! 353! complete field shp2 354! 355 shp2(1,nope)=0.d0 356 shp2(2,nope)=0.d0 357 shp2(4,nope)=-1.d0 358! 359 do k=1,nope 360 do l=1,nope 361 do i=1,3 362 i1=i+(k-1)*3 363 do j=1,3 364 s(i1,j+(l-1)*3)=-shp2(4,k)*fpu(i,j,l) 365 & -shp2(1,k)*fnl(i)*dxi(j,l) 366 & -shp2(2,k)*fnl(i)*det(j,l) 367 enddo 368 enddo 369 enddo 370 enddo 371! 372! symmetrizing the matrix 373! this is done in the absence of friction or for modal dynamic 374! calculations 375! 376 if((nasym.eq.0).or.((nmethod.eq.4).and.(iperturb(1).le.1))) then 377 do j=1,3*nope 378 do i=1,j-1 379 s(i,j)=(s(i,j)+s(j,i))/2.d0 380 enddo 381 enddo 382 endif 383! 384 return 385 end 386 387