1 !=====================================================================
2 ! DISCRETE LIMIT ANALYSIS PROGRAM !
3 ! USING KAWAI'S MODEL !
4 ! (RBSM : Rigid Bodies-Spring Model) !
5 ! PC-Ver2.00 !
6 ! (2015.7) !
7 ! PROGRAM BY Prof. N.TAKEUCHI (Hosei University) !
8 !=====================================================================
9 PROGRAM RBSM !
10 !=====================================================================
11 IMPLICIT NONE
12 !< Node type >--------------------------
13 TYPE :: typeNode !
14 INTEGER :: no !
15 INTEGER :: dim !
16 REAL(8), POINTER :: coord(:,:) !
17 END TYPE !
18 !---------------------------------------
19 !< Element type >-----------------------
20 TYPE :: typeElement !
21 INTEGER :: no !
22 INTEGER :: dof !
23 INTEGER :: type !
24 INTEGER, POINTER :: node(:,:) !
25 INTEGER, POINTER :: material(:) !
26 REAL(8), POINTER :: area(:) !
27 REAL(8), POINTER :: center(:,:) !
28 REAL(8), POINTER :: load(:) !
29 REAL(8), POINTER :: disp(:) !
30 REAL(8), POINTER :: ddisp(:) !
31 END TYPE !
32 !---------------------------------------
33 !< Material type >-----------------------
34 TYPE :: typeMaterial !
35 INTEGER :: no !
36 REAL(8), POINTER :: young(:) !
37 REAL(8), POINTER :: poisson(:) !
38 REAL(8), POINTER :: weight(:) !
39 REAL(8), POINTER :: c(:) !
40 REAL(8), POINTER :: phi(:) !
41 REAL(8), POINTER :: thick(:) !
42 END TYPE !
43 !---------------------------------------
44 !< Spring type >------------------------
45 TYPE :: typeSpring !
46 INTEGER :: no !
47 INTEGER, POINTER :: element(:,:) !
48 INTEGER, POINTER :: node(:,:) !
49 REAL(8), POINTER :: length(:) !
50 REAL(8), POINTER :: hline(:,:) !
51 INTEGER, POINTER :: flag(:) !
52 REAL(8), POINTER :: strain(:,:) !
53 REAL(8), POINTER :: dstrain(:,:) !
54 REAL(8), POINTER :: stress(:,:) !
55 REAL(8), POINTER :: dstress(:,:) !
56 END TYPE !
57 !---------------------------------------
58 !< Load type >--------------------------
59 TYPE :: typeLoad !
60 INTEGER :: no !
61 INTEGER, POINTER :: element(:) !
62 REAL(8), POINTER :: force(:,:) !
63 END TYPE !
64 !---------------------------------------
65 !< Boundary(fix) type >-----------------
66 TYPE :: typeBoundary !
67 INTEGER :: no !
68 INTEGER, POINTER :: element(:) !
69 INTEGER, POINTER :: direction(:) !
70 END TYPE !
71 !---------------------------------------
72 !< Solver(band) type >------------------
73 TYPE :: typeBand !
74 INTEGER :: no !
75 INTEGER :: width !
76 REAL(8), POINTER :: stiff(:) !
77 REAL(8), POINTER :: load(:) !
78 END TYPE !
79 !---------------------------------------
80 !< File type >--------------------------
81 TYPE :: typeFile !
82 INTEGER :: io !
83 INTEGER :: lp !
84 END TYPE !
85 !---------------------------------------
86 !< Nonlinear type >---------------------
87 TYPE :: typeNonLinear !
88 INTEGER :: iteration !
89 REAL(8) :: maxdisp !
90 REAL(8) :: rmin !
91 REAL(8) :: step !
92 END TYPE !
93 !---------------------------------------
94 !=============================
95 TYPE(typeNode) :: Node !
96 TYPE(typeElement) :: Element !
97 TYPE(typeMaterial) :: Material !
98 TYPE(typeSpring) :: Spring !
99 TYPE(typeLoad) :: Load !
100 TYPE(typeBoundary) :: Boundary !
101 TYPE(typeNonLinear):: NonLinear !
102 TYPE(typeBand) :: Solver !
103 TYPE(typeFile) :: File !
104 !=======================================
105 CHARACTER*80 Title
106 File%io = 5
107 File%lp = 6
108
109 CALL Prepocessing(Title,Node,Element,Material,Spring,Load,Boundary, &
110 Nonlinear,Solver,File)
111 CALL Analysis(Node,Element,Material,Spring,Load,Boundary,Nonlinear, &
112 Solver,File)
113 !-------
114 CONTAINS
115 !-------
116 !===============================================================================
117 SUBROUTINE Prepocessing( Title, Node, Element, Material, Spring, Load, &
118 Boundary, Nonlinear, Solver, File)
119
120 CHARACTER*80, INTENT(OUT) :: Title
121 TYPE(typeNode), INTENT(OUT) :: Node
122 TYPE(typeElement), INTENT(OUT) :: Element
123 TYPE(typeMaterial), INTENT(OUT) :: Material
124 TYPE(typeSpring), INTENT(OUT) :: Spring
125 TYPE(typeLoad), INTENT(OUT) :: Load
126 TYPE(typeBoundary), INTENT(OUT) :: Boundary
127 TYPE(typeNonLinear),INTENT(OUT) :: NonLinear
128 TYPE(typeBand), INTENT(OUT) :: Solver
129 TYPE(typeFile), INTENT(IN) :: File
130
131 CALL InputData(Title,Node,Element,Material,Load,Boundary,Nonlinear,Solver,File)
132 Spring%no = noSpring(Element)
133 CALL SpringData(Element,Spring)
134 CALL Area(Node,Element)
135 CALL gCenter(Node,Element)
136 CALL Perpendicular(Node,Element,Spring)
137 Solver%width = BandWidth(Element,Spring)
138 CALL PrintData(Title,Node,Element,Material,Spring,Load,Boundary,Nonlinear, &
139 Solver,File)
140 END SUBROUTINE Prepocessing
141 !-------------------------------------------------------------------------------
142 SUBROUTINE InputData( Title, Node, Element, Material, Load, Boundary, &
143 Nonlinear, Solver, File)
144
145 CHARACTER*80, INTENT(OUT) :: Title
146 TYPE(typeNode), INTENT(OUT) :: Node
147 TYPE(typeElement), INTENT(OUT) :: Element
148 TYPE(typeMaterial), INTENT(OUT) :: Material
149 TYPE(typeLoad), INTENT(OUT) :: Load
150 TYPE(typeBoundary), INTENT(OUT) :: Boundary
151 TYPE(typeNonLinear),INTENT(OUT) :: NonLinear
152 TYPE(typeBand), INTENT(OUT) :: Solver
153 TYPE(typeFile), INTENT(IN) :: File
154
155 INTEGER :: in,ie,im,il,ib,i,j
156
157 Node%dim = 2
158 Element%dof = 3
159 NonLinear%maxdisp = 1000.D0
160
161 READ(File%io,'(A80)') Title
162 READ(File%io,*) NonLinear%iteration
163 READ(File%io,*) Node%no
164 ALLOCATE( Node%coord(Node%dim,Node%no) )
165 DO i = 1, Node%no
166 READ(File%io,*) in,(Node%coord(j,in),j=1,Node%dim)
167 END DO
168 READ(File%io,*) Element%no,Element%type
169 ALLOCATE( Element%node(3,Element%no) )
170 ALLOCATE( Element%material(Element%no) )
171 DO i = 1, Element%no
172 READ(File%io,*) ie,Element%material(ie),(Element%node(j,ie),j=1,3)
173 END DO
174 READ(File%io,*) Material%no
175 ALLOCATE( Material%young(Material%no) )
176 ALLOCATE( Material%poisson(Material%no) )
177 ALLOCATE( Material%weight(Material%no) )
178 ALLOCATE( Material%c(Material%no) )
179 ALLOCATE( Material%phi(Material%no) )
180 ALLOCATE( Material%thick(Material%no) )
181 DO i = 1, Material%no
182 READ(File%io,*) im,Material%Young(im), Material%poisson(im), &
183 Material%weight(im),Material%c(im), &
184 Material%phi(im), Material%thick(im)
185 IF(Material%thick(im) == 0.D0) THEN
186 Material%thick(im) = 1.D0
187 END IF
188 END DO
189 READ(File%io,*) Load%no
190 IF(Load%no /= 0) THEN
191 ALLOCATE( Load%element(Load%no) )
192 ALLOCATE( Load%force(Element%dof,Load%no) )
193 DO i = 1, Load%no
194 READ(File%io,*) il,Load%element(il), &
195 (Load%force(j,il),j=1,Element%dof)
196 END DO
197 END IF
198 READ(File%io,*) Boundary%no
199 ALLOCATE( Boundary%element(Boundary%no) )
200 ALLOCATE( Boundary%direction(Boundary%no) )
201 DO i = 1, Boundary%no
202 READ(File%io,*) ib,Boundary%element(ib),Boundary%direction(ib)
203 END DO
204 ALLOCATE( Element%area(Element%no) )
205 ALLOCATE( Element%center(2,Element%no) )
206 Solver%no = Element%no*Element%dof
207 ALLOCATE( Element%load(Solver%no) )
208 ALLOCATE( Element%disp(Solver%no) )
209 ALLOCATE( Element%ddisp(Solver%no) )
210 END SUBROUTINE InputData
211 !-------------------------------------------------------------------------------
212 FUNCTION noSpring( Element ) RESULT( no )
213
214 TYPE(typeElement), INTENT(IN) :: Element
215 INTEGER :: no
216
217 INTEGER :: ie,je,same,edge,error,i,j
218
219 no = 0
220 error = 0
221 DO ie = 1, Element%no - 1
222 edge = 0
223 DO je = ie+1, Element%no
224 same = 0
225 DO i = 1, 3
226 IF(Element%node(i,ie) /= 0) THEN
227 DO j = 1, 3
228 IF(Element%node(i,ie) == Element%node(j,je)) THEN
229 same = same + 1
230 END IF
231 END DO
232 END IF
233 END DO
234 IF(same > 2) THEN
235 error = error + 1
236 PRINT 6000, ie,je
237 6000 FORMAT(/,' *** ERROR *** SAME NODE ELEM. NO.(',I3,',',I3,')')
238 ELSE IF(same == 2) THEN
239 no = no + 1
240 edge = edge + 1
241 END IF
242 END DO
243 IF(edge > 3) THEN
244 error = error + 1
245 PRINT 6010, ie,(Element%node(i,ie),i=1,3)
246 6010 FORMAT(/,' *** ERROR *** NODE DATA ELEM. NO.',I3, &
247 ' NODE NO.(',I3,',',I3,',',I3,')')
248 END IF
249 END DO
250 IF(error == 0) RETURN
251 STOP 100
252 END FUNCTION noSpring
253 !-------------------------------------------------------------------------------
254 SUBROUTINE SpringData( Element, Spring )
255
256 TYPE(typeElement), INTENT(IN) :: Element
257 TYPE(typeSpring), INTENT(INOUT) :: Spring
258
259 INTEGER :: no,is,ie,je,in,jn,n1,n2,same,i,j,ibit
260
261 ALLOCATE( Spring%node(2,Spring%no) )
262 ALLOCATE( Spring%element(2,Spring%no) )
263 ALLOCATE( Spring%length(Spring%no) )
264 ALLOCATE( Spring%hline(2,Spring%no) )
265 ALLOCATE( Spring%flag(Spring%no) )
266 ALLOCATE( Spring%strain(2,Spring%no) )
267 ALLOCATE( Spring%dstrain(2,Spring%no) )
268 ALLOCATE( Spring%stress(2,Spring%no) )
269 ALLOCATE( Spring%dstress(2,Spring%no) )
270
271 no = 0
272 DO ie = 1, Element%no - 1
273 DO je = ie+1, Element%no
274 same = 0
275 DO i = 1, 3
276 IF(Element%node(i,ie) /= 0) THEN
277 DO j = 1, 3
278 IF(Element%node(i,ie) == Element%node(j,je)) THEN
279 same = same + 1
280 END IF
281 END DO
282 END IF
283 END DO
284 IF(same == 2) THEN
285 no = no + 1
286 Spring%element(1,no) = ie
287 Spring%element(2,no) = je
288 END IF
289 END DO
290 END DO
291 Loop : DO is = 1, Spring%no
292 ie = Spring%element(1,is)
293 je = Spring%element(2,is)
294 DO i = 1, 3
295 in = Element%node(i,ie)
296 ibit = 0
297 DO j = 1, 3
298 jn = Element%node(j,je)
299 IF(in == jn) THEN
300 ibit = 1
301 END IF
302 END DO
303 IF(ibit == 0) THEN
304 n1 = i + 1
305 n2 = i + 2
306 IF(n1 > 3) n1 = n1 - 3
307 IF(n2 > 3) n2 = n2 - 3
308 Spring%node(1,is) = Element%node(n1,ie)
309 Spring%node(2,is) = Element%node(n2,ie)
310 CYCLE Loop
311 END IF
312 END DO
313 END DO Loop
314 END SUBROUTINE SpringData
315 !-------------------------------------------------------------------------------
316 SUBROUTINE Area( Node, Element )
317
318 TYPE(typeNode), INTENT(IN) :: Node
319 TYPE(typeElement), INTENT(INOUT) :: Element
320
321 INTEGER :: ie,n1,n2,n3,error
322 REAL(8) :: x1,x2,x3,y1,y2,y3
323
324 error = 0
325 DO ie = 1, Element%no
326 n1 = Element%node(1,ie)
327 n2 = Element%node(2,ie)
328 n3 = Element%node(3,ie)
329 IF(n3 /= 0) THEN
330 x1 = Node%coord(1,n1)
331 x2 = Node%coord(1,n2)
332 x3 = Node%coord(1,n3)
333 y1 = Node%coord(2,n1)
334 y2 = Node%coord(2,n2)
335 y3 = Node%coord(2,n3)
336 Element%area(ie) = ((y2-y3)*(x1-x3)-(y3-y1)*(x3-x2))/2.D0
337 IF(Element%area(ie) <= 0.D0) THEN
338 error = error + 1
339 PRINT 6000, ie,Element%area(ie)
340 6000 FORMAT(/,' *** ERROR *** AREA =< 0 ELEM. NO.=',I4,' AREA =',D15.5)
341 END IF
342 ELSE
343 Element%area(ie) = 0.D0
344 END IF
345 END DO
346 IF(error == 0) RETURN
347 STOP 100
348 END SUBROUTINE Area
349 !-------------------------------------------------------------------------------
350 SUBROUTINE gCenter( Node, Element )
351
352 TYPE(typeNode), INTENT(IN) :: Node
353 TYPE(typeElement), INTENT(INOUT) :: Element
354
355 INTEGER :: ie,n,i
356 REAL(8) :: ww,wx,wy
357
358 DO ie = 1, Element%no
359 ww = 0.0
360 wx = 0.0
361 wy = 0.0
362 DO i = 1, 3
363 n = Element%node(i,ie)
364 IF(n /= 0) THEN
365 ww = ww + 1
366 wx = wx + Node%coord(1,n)
367 wy = wy + Node%coord(2,n)
368 END IF
369 END DO
370 Element%center(1,ie) = wx/ww
371 Element%center(2,ie) = wy/ww
372 END DO
373 END SUBROUTINE gCenter
374 !-------------------------------------------------------------------------------
375 SUBROUTINE Perpendicular( Node, Element, Spring )
376
377 TYPE(typeNode), INTENT(IN) :: Node
378 TYPE(typeElement), INTENT(IN) :: Element
379 TYPE(typeSpring), INTENT(INOUT) :: Spring
380
381 INTEGER :: is,n1,n2,e1,e2
382 REAL(8) :: sl,sm,xm,ym,h1,h2,x,y,xl,yl
383
384 DO is = 1, Spring%no
385 n1 = Spring%node(1,is)
386 n2 = Spring%node(2,is)
387 x = Node%coord(1,n2) - Node%coord(1,n1)
388 y = Node%coord(2,n2) - Node%coord(2,n1)
389 Spring%length(is) = DSQRT(x*x+y*y)
390 sl = y/Spring%length(is)
391 sm = x/Spring%length(is)
392 xm = (Node%coord(1,n1) + Node%coord(1,n2)) /2.D0
393 ym = (Node%coord(2,n1) + Node%coord(2,n2)) /2.D0
394 e1 = Spring%element(1,is)
395 e2 = Spring%element(2,is)
396 xl = xm - Element%center(1,e1)
397 yl = ym - Element%center(2,e1)
398 h1 = sl*xl - sm*yl
399 xl = Element%center(1,e2) - xm
400 yl = Element%center(2,e2) - ym
401 h2 = sl*xl - sm*yl
402 Spring%hline(1,is) = ABS(h1)
403 Spring%hline(2,is) = ABS(h2)
404 END DO
405 END SUBROUTINE Perpendicular
406 !-------------------------------------------------------------------------------
407 FUNCTION BandWidth( Element, Spring ) RESULT( width )
408
409 TYPE(typeElement), INTENT(IN) :: Element
410 TYPE(typeSpring), INTENT(IN) :: Spring
411
412 INTEGER :: width,is,wd
413
414 width = 0
415 DO is = 1, Spring%no
416 wd = Spring%element(2,is) - Spring%element(1,is)
417 IF(wd > width) THEN
418 width = wd
419 END IF
420 END DO
421 width = (width+1)*Element%dof
422 ALLOCATE( Solver%stiff(Solver%no*width) )
423 ALLOCATE( Solver%load(Solver%no) )
424 END FUNCTION BandWidth
425 !-------------------------------------------------------------------------------
426 SUBROUTINE PrintData( Title, Node, Element, Material, Spring, Load, &
427 Boundary, Nonlinear, Solver, File)
428
429 CHARACTER*80, INTENT(IN) :: Title
430 TYPE(typeNode), INTENT(IN) :: Node
431 TYPE(typeElement), INTENT(IN) :: Element
432 TYPE(typeMaterial), INTENT(IN) :: Material
433 TYPE(typeSpring), INTENT(IN) :: Spring
434 TYPE(typeLoad), INTENT(IN) :: Load
435 TYPE(typeBoundary), INTENT(IN) :: Boundary
436 TYPE(typeNonLinear),INTENT(IN) :: NonLinear
437 TYPE(typeBand), INTENT(IN) :: Solver
438 TYPE(typeFile), INTENT(IN) :: File
439
440 INTEGER :: in,ie,is,im,il,ib,i
441
442 WRITE(File%lp,6000) Title
443 IF(Element%type == 0) THEN
444 WRITE(File%lp,6010)
445 ELSE
446 WRITE(File%lp,6020)
447 END IF
448 WRITE(File%lp,6030) Node%no,Element%no,Material%no,Spring%no, &
449 Load%no,Boundary%no,Solver%width, &
450 Nonlinear%iteration,Nonlinear%maxdisp
451 WRITE(File%lp,6040)
452 DO in = 1, Node%no
453 WRITE(File%lp,6050) in,Node%coord(1,in),Node%coord(2,in)
454 END DO
455 WRITE(File%lp,6060)
456 DO ie = 1, Element%no
457 WRITE(File%lp,6070) ie,(Element%node(i,ie),i=1,3),Element%material(ie), &
458 Element%area(ie),(Element%center(i,ie),i=1,2)
459 END DO
460 WRITE(File%lp,6080)
461 DO is = 1, Spring%no
462 WRITE(File%lp,6090) is,(Spring%element(i,is),i=1,2), &
463 (Spring%node(i,is),i=1,2), &
464 Spring%length(is),(Spring%hline(i,is),i=1,2)
465 END DO
466 WRITE(File%lp,6100)
467 DO im = 1, Material%no
468 WRITE(File%lp,6110) im,Material%young(im),Material%poisson(im), &
469 Material%weight(im),Material%c(im), &
470 Material%phi(im),Material%thick(im)
471 END DO
472 IF(Load%no /= 0) THEN
473 WRITE(File%lp,6120)
474 DO il = 1, Load%no
475 WRITE(File%lp,6130) il,Load%element(il),(Load%force(i,il),i=1,3)
476 END DO
477 END IF
478 WRITE(File%lp,6140)
479 DO ib = 1, Boundary%no
480 WRITE(File%lp,6150) ib,Boundary%element(ib),Boundary%direction(ib)
481 END DO
482 RETURN
483 6000 FORMAT(' TITLE : ',A80,//)
484 6010 FORMAT(' INFORMATION :',//, &
485 ' ELEMET TYPE ------------------> PLANE STRAIN')
486 6020 FORMAT(' INFORMATION :',//, &
487 ' ELEMET TYPE ------------------> PLANE STRESS')
488 6030 FORMAT(' NUMBER OF NODE --------------->',I5,/, &
489 ' NUMBER OF ELEMENTS ----------->',I5,/, &
490 ' NUMBER OF MATERIALS ---------->',I5,/, &
491 ' NUMBER OF SPRINGS ------------>',I5,/, &
492 ' NUMBER OF LOADS -------------->',I5,/, &
493 ' NUMBER OF SUPPORT D.O.F. ----->',I5,/, &
494 ' MAXIMUM BAND WIDTH ----------->',I5,/, &
495 ' NUMBER OF ITARATIONS --------->',I5,/, &
496 ' MAXIMUM DISPLACEMENT --------->',E10.1,//)
497 6040 FORMAT(' COORDINATE DATA :',//, &
498 ' NODE NO. X-COORDINATE Y-COORDINATE',/)
499 6050 FORMAT(3X,I5,2E15.5)
500 6060 FORMAT(//, &
501 ' ELEMENT DATA :',//, &
502 ' ELEM.NO. (1) (2) (3) MAT.NO. AREA', &
503 ' XG YG',/)
504 6070 FORMAT(3X,4I5,3X,I5,3E15.5)
505 6080 FORMAT(//, &
506 ' SPRING DATA :',//, &
507 ' SPRG.NO. E-1 E-2 N-1 N-2 L', &
508 ' H1 H2',/)
509 6090 FORMAT(3X,5I5,3E15.5)
510 6100 FORMAT(//, &
511 ' MATERIAL DATA :',//, &
512 ' MAT.NO. E P BODY FORCE', &
513 ' C PHI t',/)
514 6110 FORMAT(2X,I5,7E12.5)
515 6120 FORMAT(//, &
516 ' LOAD DATA :',//, &
517 ' NO. ELEM.NO. X Y M',/)
518 6130 FORMAT(I6,5X,I5,3E12.3)
519 6140 FORMAT(//, &
520 ' SUPPORT D.O.F.:',//, &
521 ' NO. ELEM.NO. DIREC.',/)
522 6150 FORMAT(I6,5X,I5,I7)
523 END SUBROUTINE PrintData
524 !===============================================================================
525 SUBROUTINE Analysis( Node, Element, Material, Spring, Load, Boundary, &
526 Nonlinear, Solver, File)
527
528 TYPE(typeNode), INTENT(IN) :: Node
529 TYPE(typeElement), INTENT(INOUT) :: Element
530 TYPE(typeMaterial), INTENT(IN) :: Material
531 TYPE(typeSpring), INTENT(INOUT) :: Spring
532 TYPE(typeLoad), INTENT(IN) :: Load
533 TYPE(typeBoundary), INTENT(IN) :: Boundary
534 TYPE(typeNonLinear),INTENT(INOUT) :: NonLinear
535 TYPE(typeBand), INTENT(INOUT) :: Solver
536 TYPE(typeFile), INTENT(IN) :: File
537
538 INTEGER :: iste,jramc,jramd,error
539
540 Nonlinear%rmin = 0.D0
541 CALL ClearArray(Element,Spring)
542 CALL formLoad(Element,Material,Load)
543 iste = 0
544 Iteration : DO WHILE(.TRUE.)
545 iste = iste + 1
546 jramc = 0
547 jramd = 1
548 ! unloading & loading again iteration ===>
549 Unload : DO WHILE( jramd == 1)
550 Solver%stiff = 0.D0
551 Solver%load = Element%load
552 CALL formStiff(Node,Element,Material,Spring,Solver)
553 CALL setBoundary(Element,Boundary,Solver)
554 CALL Solve(Solver%no,Solver%width,Solver%stiff,Solver%load,error)
555 CALL setDisplacement(Element,Boundary,Solver)
556 IF(.not. MechanismCheck(Element,File,error)) RETURN
557 CALL SpringStress(Node,Element,Material,Spring,jramc,jramd)
558 IF(jramd == 1) CYCLE Unload
559 CALL incRatio(Element,Material,Spring,Nonlinear)
560 CALL totalValue(Element,Spring,Nonlinear)
561 CALL YieldCheck(Element,Material,Spring)
562 CALL OUTPUT(iste,Element,Spring,Nonlinear,File)
563 END DO Unload
564 IF(NonLinear%rmin > 0.998D0) EXIT
565 IF(iste == NonLinear%iteration) EXIT
566 END DO Iteration
567 END SUBROUTINE Analysis
568 !-------------------------------------------------------------------------------
569 SUBROUTINE ClearArray( Element, Spring )
570
571 TYPE(typeElement), INTENT(INOUT) :: Element
572 TYPE(typeSpring), INTENT(INOUT) :: Spring
573
574 Element%disp = 0.D0
575 Spring%flag = 0
576 Spring%strain = 0.D0
577 Spring%stress = 0.D0
578 END SUBROUTINE ClearArray
579 !-------------------------------------------------------------------------------
580 SUBROUTINE formLoad( Element, Material, Load )
581
582 TYPE(typeElement), INTENT(INOUT) :: Element
583 TYPE(typeMaterial), INTENT(IN) :: Material
584 TYPE(typeLoad), INTENT(IN) :: Load
585
586 INTEGER :: ie,im,in,il,i
587 REAL(8) :: weight
588
589 Element%load = 0.D0
590 DO ie =1, Element%no
591 IF(Element%node(3,ie) /= 0) THEN
592 im = Element%material(ie)
593 IF(Material%weight(im) /= 0.D0) THEN
594 in = 2 + (ie-1)*Element%dof
595 weight = Material%weight(im)*Element%area(ie)*Material%thick(im)
596 Element%load(in) = Element%load(in) - weight
597 END IF
598 END IF
599 END DO
600 DO il = 1, Load%no
601 ie = Load%element(il)
602 DO i = 1, 3
603 in = i + (ie-1)*Element%dof
604 Element%load(in) = Element%load(in) + Load%force(i,il)
605 END DO
606 END DO
607 END SUBROUTINE formLoad
608 !-------------------------------------------------------------------------------
609 SUBROUTINE formStiff( Node, Element, Material, Spring, Solver)
610
611 TYPE(typeNode), INTENT(IN) :: Node
612 TYPE(typeElement), INTENT(IN) :: Element
613 TYPE(typeMaterial), INTENT(IN) :: Material
614 TYPE(typeSpring), INTENT(IN) :: Spring
615 TYPE(typeBand), INTENT(INOUT) :: Solver
616
617 REAL(8) :: sfm(6,6),coef(6),D(2,2)
618 INTEGER :: is
619
620 DO is = 1, Spring%no
621 CALL formSpring(is,Element,Material,Spring,D)
622 CALL coefStiff(is,Node,Spring,coef)
623 IF(Spring%flag(is) /= 1) THEN
624 CALL elasticStiff(is,Element,Material,Spring,D,coef,sfm)
625 ELSE
626 CALL plasticStiff(is,Element,Material,Spring,D,coef,sfm)
627 END IF
628 CALL Assemble(is,sfm,Spring,Solver)
629 END DO
630 END SUBROUTINE formStiff
631 !-------------------------------------------------------------------------------
632 SUBROUTINE formSpring( is, Element, Material, Spring, D )
633
634 INTEGER, INTENT(IN) :: is
635 TYPE(typeElement), INTENT(IN) :: Element
636 TYPE(typeMaterial), INTENT(IN) :: Material
637 TYPE(typeSpring), INTENT(IN) :: Spring
638 REAL(8), INTENT(INOUT) :: D(:,:)
639
640 REAL(8) :: Dp(2,2),E,P,phi,h
641 INTEGER :: ie1,ie2,im1,im2
642
643 D = 0.D0
644 h = Spring%hline(1,is) + Spring%hline(2,is)
645 ie1 = Spring%element(1,is)
646 ie2 = Spring%element(2,is)
647 im1 = Element%material(ie1)
648 im2 = Element%material(ie2)
649 E = (Material%young(im1)*Spring%hline(1,is) &
650 + Material%young(im2)*Spring%hline(2,is))/h
651 P = (Material%poisson(im1)*Spring%hline(1,is) &
652 + Material%poisson(im2)*Spring%hline(2,is))/h
653 IF(Element%type == 0) THEN
654 D(1,1) = E*(1.0-P)/(1.0+P)/(1.0-2.0*P)
655 ELSE
656 D(1,1) = E/(1.0-P**2)
657 END IF
658 D(2,2) = E/(1.0+P)
659
660 IF(Spring%flag(is) == 1) THEN
661 IF(Element%type == 0) THEN
662 phi = (Material%phi(im1)*Spring%hline(1,is) &
663 + Material%phi(im2)*Spring%hline(2,is))/h
664 CALL MohrCoulomb(is,phi,Spring,D,Dp)
665 ELSE
666 CALL vonMises(is,Spring,D,Dp)
667 END IF
668 D = D - Dp
669 END IF
670 END SUBROUTINE formSpring
671 !-------------------------------------------------------------------------------
672 SUBROUTINE MohrCoulomb( is, phi, Spring, D, Dp )
673
674 INTEGER, INTENT(IN) :: is
675 REAL(8), INTENT(IN) :: phi
676 TYPE(typeSpring), INTENT(IN) :: Spring
677 REAL(8), INTENT(IN) :: D(:,:)
678 REAL(8), INTENT(OUT) :: Dp(:,:)
679
680 REAL(8) :: pt,f,D1,D2
681 INTEGER :: i,j
682 REAL(8) :: PAI = 3.141592654D0
683
684 pt = DTAN(phi*PAI/180.D0)
685 D1 = D(1,1)
686 D2 = D(2,2)
687 f = 1.D0/(D2+D1*pt**2)
688 Dp(1,1) = (D1*pt)**2*f
689 Dp(2,2) = D2*D2*f
690 Dp(1,2) = D1*D2*pt*f
691 IF(Spring%stress(2,is) < 0.D0) THEN
692 Dp(1,2) = -Dp(1,2)
693 END IF
694 Dp(2,1) = Dp(1,2)
695 END SUBROUTINE MohrCoulomb
696 !-------------------------------------------------------------------------------
697 SUBROUTINE vonMises( is, Spring, D, Dp )
698
699 INTEGER, INTENT(IN) :: is
700 TYPE(typeSpring), INTENT(IN) :: Spring
701 REAL(8), INTENT(IN) :: D(:,:)
702 REAL(8), INTENT(OUT) :: Dp(:,:)
703
704 REAL(8) :: sig,tau,f,D1,D2
705 INTEGER :: i,j
706
707 D1 = D(1,1)
708 D2 = D(2,2)
709 sig = Spring%stress(1,is)
710 tau = Spring%Stress(2,is)
711 f = 1.D0/(D1*sig*sig/4.D0 + 4.D0*D2*tau*tau)
712 Dp(1,1) = D1*D1*sig*sig/4.D0*f
713 Dp(2,2) = 4.D0*D2*D2*tau*tau*f
714 Dp(1,2) = D1*D2*sig*tau*f
715 Dp(2,1) = Dp(1,2)
716 END SUBROUTINE vonMises
717 !-------------------------------------------------------------------------------
718 SUBROUTINE coefStiff( is, Node, Spring, coef )
719
720 INTEGER, INTENT(IN) :: is
721 TYPE(typeNode), INTENT(IN) :: Node
722 TYPE(typeSpring), INTENT(IN) :: Spring
723 REAL(8), INTENT(OUT) :: coef(:)
724
725 INTEGER :: ie1,ie2,in1,in2
726 REAL(8) :: x31,x32,x51,x52,y31,y32,y51,y52
727
728 ie1 = Spring%element(1,is)
729 ie2 = Spring%element(2,is)
730 in1 = Spring%node(1,is)
731 in2 = Spring%node(2,is)
732 x31 = Node%coord(1,in2) - Element%center(1,ie1)
733 x32 = Node%coord(1,in2) - Element%center(1,ie2)
734 x51 = Node%coord(1,in1) - Element%center(1,ie1)
735 x52 = Node%coord(1,in1) - Element%center(1,ie2)
736 y31 = Node%coord(2,in2) - Element%center(2,ie1)
737 y32 = Node%coord(2,in2) - Element%center(2,ie2)
738 y51 = Node%coord(2,in1) - Element%center(2,ie1)
739 y52 = Node%coord(2,in1) - Element%center(2,ie2)
740 coef(1) = Node%coord(1,in2) - Node%coord(1,in1)
741 coef(2) = Node%coord(2,in2) - Node%coord(2,in1)
742 coef(3) = 0.5*( coef(1)*(x31+x51) + coef(2)*(y31+y51))
743 coef(4) = 0.5*( coef(1)*(y32+y52) - coef(2)*(x32+x52))
744 coef(5) = 0.5*(-coef(1)*(y31+y51) + coef(2)*(x31+x51))
745 coef(6) = 0.5*(-coef(1)*(x32+x52) - coef(2)*(y32+y52))
746 END SUBROUTINE coefStiff
747 !-------------------------------------------------------------------------------
748 SUBROUTINE elasticStiff( is, Element, Material, Spring, D, coef, sfm )
749
750 INTEGER, INTENT(IN) :: is
751 TYPE(typeElement), INTENT(IN) :: Element
752 TYPE(typeMaterial), INTENT(IN) :: Material
753 TYPE(typeSpring), INTENT(IN) :: Spring
754 REAL(8), INTENT(IN) :: D(:,:)
755 REAL(8), INTENT(IN) :: coef(:)
756 REAL(8), INTENT(OUT) :: sfm(:,:)
757
758 REAL(8) :: h,thick,Kd,Ks,Kr
759 REAL(8) :: x35,y35,D11,D12,D21,D22
760 INTEGER :: ie1,ie2,im1,im2,i,j
761
762 h = Spring%hline(1,is) + Spring%hline(2,is)
763 ie1 = Spring%element(1,is)
764 ie2 = Spring%element(2,is)
765 im1 = Element%material(ie1)
766 im2 = Element%material(ie2)
767 thick = (Material%thick(im1)*Spring%hline(1,is) &
768 + Material%thick(im2)*Spring%hline(2,is))/H
769 Kd = D(1,1)/Spring%length(is)/h*thick
770 Ks = D(2,2)/Spring%length(is)/h*thick
771 Kr = Kd*Spring%length(is)**4/12.D0
772 x35 = coef(1)
773 y35 = coef(2)
774 D11 = coef(3)
775 D12 = coef(4)
776 D21 = coef(5)
777 D22 = coef(6)
778 sfm(1,2) = -(Kd-Ks)*x35*y35
779 sfm(1,3) = Kd*y35*D11 - Ks*x35*D21
780 sfm(1,4) = -(Kd*y35*y35 + Ks*x35*x35)
781 sfm(1,5) = -sfm(1,2)
782 sfm(1,6) = Kd*y35*D22 - Ks*x35*D12
783 sfm(2,3) = -(Kd*x35*D11 + Ks*y35*D21)
784 sfm(2,4) = -sfm(1,2)
785 sfm(2,5) = -(Kd*x35*x35 + Ks*y35*y35)
786 sfm(2,6) = -(Kd*x35*D22 + Ks*y35*D12)
787 sfm(3,4) = -sfm(1,3)
788 sfm(3,5) = -sfm(2,3)
789 sfm(3,6) = Kd*D11*D22 + Ks*D21*D12 - Kr
790 sfm(4,5) = sfm(1,2)
791 sfm(4,6) = -sfm(1,6)
792 sfm(5,6) = -sfm(2,6)
793 DO i = 1, 6
794 DO j = i, 6
795 sfm(j,i) = sfm(i,j)
796 END DO
797 END DO
798 sfm(1,1) = -sfm(1,4)
799 sfm(2,2) = -sfm(2,5)
800 sfm(3,3) = Kd*D11*D11 + Ks*D21*D21 + Kr
801 sfm(4,4) = sfm(1,1)
802 sfm(5,5) = sfm(2,2)
803 sfm(6,6) = Kd*D22*D22 + Ks*D12*D12 + Kr
804 END SUBROUTINE elasticStiff
805 !-------------------------------------------------------------------------------
806 SUBROUTINE plasticStiff( is, Element, Material, Spring, D, coef, sfm )
807
808 INTEGER, INTENT(IN) :: is ! spring number
809 TYPE(typeElement), INTENT(IN) :: Element ! Element Data
810 TYPE(typeMaterial), INTENT(IN) :: Material ! Material Data
811 TYPE(typeSpring), INTENT(IN) :: Spring ! Spring Data
812 REAL(8), INTENT(IN) :: D(:,:) ! Sprin matrix
813 REAL(8), INTENT(IN) :: coef(:) ! coefficient of spring stiffness
814 REAL(8), INTENT(OUT) :: sfm(:,:) ! spring matrix
815
816 REAL(8) :: h,thick,Kd,Ks,Kds,Kr
817 REAL(8) :: x35,y35,D11,D12,D21,D22
818 INTEGER :: ie1,ie2,im1,im2,i,j
819
820 h = Spring%hline(1,is) + Spring%hline(2,is)
821 ie1 = Spring%element(1,is)
822 ie2 = Spring%element(2,is)
823 im1 = Element%material(ie1)
824 im2 = Element%material(ie2)
825 thick = (Material%thick(im1)*Spring%hline(1,is) &
826 + Material%thick(im2)*Spring%hline(2,is))/H
827 Kd = D(1,1)/Spring%length(is)/h*thick
828 Ks = D(2,2)/Spring%length(is)/h*thick
829 Kds = D(1,2)/Spring%length(is)/h*thick
830 Kr = Kd*Spring%length(is)**4/12.D0
831 x35 = coef(1)
832 y35 = coef(2)
833 D11 = coef(3)
834 D12 = coef(4)
835 D21 = coef(5)
836 D22 = coef(6)
837 sfm(1,2) = -(Kd-Ks)*x35*y35 + (y35*y35-x35*x35)*Kds
838 sfm(1,3) = Kd*y35*D11 - Ks*x35*D21 + (x35*D11-y35*D21)*Kds
839 sfm(1,4) = -(Kd*y35*y35 + Ks*x35*x35) -2.0*x35*y35*Kds
840 sfm(1,5) = -sfm(1,2)
841 sfm(1,6) = Kd*y35*D22 - Ks*x35*D12 + (x35*D22-y35*D12)*Kds
842 sfm(2,3) = -(Kd*x35*D11 + Ks*y35*D21) + (x35*D21+y35*D11)*Kds
843 sfm(2,4) = -sfm(1,2)
844 sfm(2,5) = -(Kd*x35*x35 + Ks*y35*y35) + 2.0*x35*y35*Kds
845 sfm(2,6) = -(Kd*x35*D22 + Ks*Y35*D12) + (x35*D12+y35*D22)*Kds
846 sfm(3,4) = -sfm(1,3)
847 sfm(3,5) = -sfm(2,3)
848 sfm(3,6) = Kd*D11*D22 + Ks*D21*D12 - Kr - (D11*D12 + D21*D22)*Kds
849 sfm(4,5) = sfm(1,2)
850 sfm(4,6) = -sfm(1,6)
851 sfm(5,6) = -sfm(2,6)
852 DO i = 1, 6
853 DO j = i, 6
854 sfm(j,i) = sfm(i,j)
855 END DO
856 END DO
857 sfm(1,1) = -sfm(1,4)
858 sfm(2,2) = -sfm(2,5)
859 sfm(3,3) = Kd*D11*D11 + Ks*D21*D21 + Kr - 2.D0*D11*D21*Kds
860 sfm(4,4) = sfm(1,1)
861 sfm(5,5) = sfm(2,2)
862 sfm(6,6) = Kd*D22*D22 + Ks*D12*D12 + Kr - 2.D0*D22*D12*Kds
863 END SUBROUTINE plasticStiff
864 !-------------------------------------------------------------------------------
865 SUBROUTINE Assemble( is, sfm, Spring, Solver)
866
867 INTEGER, INTENT(IN) :: is
868 REAL(8), INTENT(IN) :: sfm(:,:)
869 TYPE(typeSpring), INTENT(IN) :: Spring
870 TYPE(typeBand), INTENT(INOUT) :: Solver
871
872 INTEGER :: ie1,ie2,i,j,k,l,k1,ka,kb,la,lb
873
874 DO i = 1, 2
875 ie1 = Spring%element(i,is)
876 DO j = i, 2
877 ie2 = Spring%element(j,is)
878 DO k = 1, 3
879 k1 = k
880 IF(i /= j) THEN
881 k1 = 1
882 END IF
883 DO l = k1, 3
884 ka = k + (i-1)*3
885 kb = l + (j-1)*3
886 la = k + (ie1-1)*3
887 lb = l + (ie2-1)*3
888 IF(lb < la) THEN
889 kb = k + (i-1)*3
890 ka = l + (j-1)*3
891 lb = k + (ie1-1)*3
892 la = l + (ie2-1)*3
893 END IF
894 lb = lb + (la-1)*Solver%width - la + 1
895 Solver%stiff(lb) = Solver%stiff(lb) + sfm(ka,kb)
896 END DO
897 END DO
898 END DO
899 END DO
900 END SUBROUTINE Assemble
901 !-------------------------------------------------------------------------------
902 SUBROUTINE setBoundary( Element, Boundary, Solver)
903
904 TYPE(typeElement), INTENT(IN) :: Element
905 TYPE(typeBoundary), INTENT(IN) :: Boundary
906 TYPE(typeBand), INTENT(INOUT) :: Solver
907
908 INTEGER :: ib,ie,id,in,jn,kn,i
909
910 DO ib = 1, Boundary%no
911 ie = Boundary%element(ib)
912 id = Boundary%direction(ib)
913 in = id + (ie-1)*Element%dof
914 DO i = 2, Solver%width
915 jn = i + (in-1)*Solver%width
916 Solver%stiff(jn) = 0.D0
917 kn = in - i + 1
918 IF(kn > 0) THEN
919 jn = i + (kn-1)*Solver%width
920 Solver%stiff(jn) = 0.0
921 END IF
922 END DO
923 jn = 1 + (in-1)*Solver%width
924 Solver%stiff(jn) = 1.D0
925 END DO
926 END SUBROUTINE setBoundary
927 !-------------------------------------------------------------------------------
928 SUBROUTINE setDisplacement( Element, Boundary, Solver )
929
930 TYPE(typeElement), INTENT(INOUT) :: Element
931 TYPE(typeBoundary), INTENT(IN) :: Boundary
932 TYPE(typeBand), INTENT(INOUT) :: Solver
933
934 INTEGER :: ib,ie,id,in
935
936 DO ib = 1, Boundary%no
937 ie = Boundary%element(ib)
938 id = Boundary%direction(ib)
939 in = id + (ie-1)*Element%dof
940 Solver%load(in) = 0.D0
941 END DO
942 Element%ddisp = Solver%load
943 END SUBROUTINE setDisplacement
944 !-------------------------------------------------------------------------------
945 FUNCTION MechanismCheck( Element, File, error ) RESULT( check )
946
947 TYPE(typeElement), INTENT(IN) :: Element
948 TYPE(typeFile), INTENT(IN) :: File
949 INTEGER, INTENT(IN) :: error
950 LOGICAL :: check
951
952 REAL(8) :: dmax
953
954 check = .TRUE.
955
956 IF(error == 1) THEN
957 check = .FALSE.
958 WRITE(File%lp,6000)
959 6000 FORMAT(' ',//////////,' ',12('*'),/, &
960 ' ',' COLLAPS !!',/, &
961 ' ',12('*'))
962 RETURN
963 END IF
964 dmax = MAXVAL(ABS(Element%ddisp))
965 IF(dmax > NonLinear%maxdisp) THEN
966 check = .FALSE.
967 WRITE(File%lp,6010) dmax,NonLinear%maxdisp
968 6010 FORMAT(' ',////////,' ',44('*'),/, &
969 ' COLLAPS !! DMAX =',E15.5,' (',E15.5,' )',/, &
970 ' ',44('*'))
971 RETURN
972 END IF
973 END FUNCTION MechanismCheck
974 !-------------------------------------------------------------------------------
975 SUBROUTINE SpringStress( Node, Element, Material, Spring, jramc, jramd )
976
977 TYPE(typeNode), INTENT(IN) :: Node
978 TYPE(typeElement), INTENT(IN) :: Element
979 TYPE(typeMaterial), INTENT(IN) :: Material
980 TYPE(typeSpring), INTENT(INOUT) :: Spring
981 INTEGER, INTENT(INOUT) :: jramc
982 INTEGER, INTENT(INOUT) :: jramd
983
984 REAL(8) :: D(2,2),ww
985 INTEGER :: is,jbit,i,j
986
987 jramd = 0
988 DO is = 1, Spring%no
989 IF(jramd == 1 .AND. Spring%flag(is) == 0) CYCLE
990 IF(jramd == 1 .AND. Spring%flag(is) == -2) CYCLE
991 CALL Strain(is,Node,Element,Spring)
992 CALL formSpring(is,Element,Material,Spring,D)
993 DO i = 1, 2
994 ww = 0.0
995 DO j = 1, 2
996 ww = ww + D(i,j)*Spring%dstrain(j,is)
997 END DO
998 Spring%dstress(i,is) = ww
999 END DO
1000 IF(jramc < 2) THEN
1001 IF(Spring%flag(is) == 1 .OR. Spring%flag(is) == 2) THEN
1002 jbit = 0
1003 IF(Spring%flag(is) == 2) jbit = 1
1004 CALL Unload(is,Element,Material,Spring,D,jramd,jbit)
1005 END IF
1006 END IF
1007 END DO
1008 jramc = jramc + 1
1009 END SUBROUTINE SpringStress
1010 !-------------------------------------------------------------------------------
1011 SUBROUTINE Strain( is, Node, Element, Spring )
1012
1013 INTEGER, INTENT(IN) :: is
1014 TYPE(typeNode), INTENT(IN) :: Node
1015 TYPE(typeElement), INTENT(IN) :: Element
1016 TYPE(typeSpring), INTENT(INOUT) :: Spring
1017
1018 REAL(8) :: uu(2*Element%dof),BB(2,2*Element%dof),h,ww
1019 INTEGER :: ie,i,j,ii,jj
1020
1021 ii = 0
1022 DO i = 1, 2
1023 ie = Spring%element(i,is)
1024 DO j = 1, 3
1025 ii = ii + 1
1026 jj = j + (ie-1)*Element%dof
1027 uu(ii) = Element%ddisp(jj)
1028 END DO
1029 END DO
1030 CALL Bmatrix(is,Node,Element,Spring,BB)
1031 h = Spring%hline(1,is) + Spring%hline(2,is)
1032 DO i = 1, 2
1033 ww = 0.D0
1034 DO j = 1, 2*Element%dof
1035 ww = ww + BB(i,j)*uu(j)
1036 END DO
1037 Spring%dstrain(i,is) = ww/h
1038 END DO
1039 END SUBROUTINE Strain
1040 !-------------------------------------------------------------------------------
1041 SUBROUTINE Bmatrix( is, Node, Element, Spring, BB )
1042
1043 INTEGER, INTENT(IN) :: is
1044 TYPE(typeNode), INTENT(IN) :: Node
1045 TYPE(typeElement), INTENT(IN) :: Element
1046 TYPE(typeSpring), INTENT(IN) :: Spring
1047 REAL(8), INTENT(OUT) :: BB(:,:)
1048
1049 REAL(8) :: yl1,yl2,ym1,ym2,x,y
1050 INTEGER :: in1,in2,ie1,ie2
1051
1052 in1 = Spring%node(1,is)
1053 in2 = Spring%node(2,is)
1054 ie1 = Spring%element(1,is)
1055 ie2 = Spring%element(2,is)
1056 yl1 = (Node%coord(2,in2)-Node%coord(2,in1))/Spring%length(is)
1057 yl2 = (Node%coord(1,in2)-Node%coord(1,in1))/Spring%length(is)
1058 ym1 =-yl2
1059 ym2 = yl1
1060 x = (Node%coord(1,in1)+Node%coord(1,in2))/2.D0
1061 y = (Node%coord(2,in1)+Node%coord(2,in2))/2.D0
1062 BB(1,1) = -yl1
1063 BB(1,2) = -ym1
1064 BB(1,3) = -yl1*(y-Element%center(2,ie1)) + ym1*(x-Element%center(1,ie1))
1065 BB(1,4) = -BB(1,1)
1066 BB(1,5) = -BB(1,2)
1067 BB(1,6) = yl1*(y-Element%center(2,ie2)) - ym1*(x-Element%center(1,ie2))
1068 BB(2,1) = -yl2
1069 BB(2,2) = -ym2
1070 BB(2,3) = -yl2*(y-Element%center(2,ie1)) + ym2*(x-Element%center(1,ie1))
1071 BB(2,4) = -BB(2,1)
1072 BB(2,5) = -BB(2,2)
1073 BB(2,6) = yl2*(y-Element%center(2,ie2)) - ym2*(x-Element%center(1,ie2))
1074 END SUBROUTINE Bmatrix
1075 !-------------------------------------------------------------------------------
1076 SUBROUTINE Unload( is, Element, Material, Spring, D, jramd, jbit)
1077
1078 INTEGER, INTENT(IN) :: is
1079 TYPE(typeElement), INTENT(IN) :: Element
1080 TYPE(typeMaterial), INTENT(IN) :: Material
1081 TYPE(typeSpring), INTENT(INOUT) :: Spring
1082 REAL(8), INTENT(IN) :: D(:,:)
1083 INTEGER, INTENT(INOUT) :: jramd
1084 INTEGER, INTENT(IN) :: jbit
1085
1086 REAL(8) :: lam ! lambda
1087
1088 lam = Lambda(is,Element,Material,Spring,D)
1089 IF(jbit /= 1) THEN
1090 IF(lam < 0.D0) THEN
1091 Spring%flag(is) = 2
1092 jramd = 1
1093 END IF
1094 ELSE
1095 IF(lam > 0.D0) THEN
1096 Spring%flag(is) = 1
1097 jramd = 1
1098 END IF
1099 END IF
1100 END SUBROUTINE Unload
1101 !-------------------------------------------------------------------------------
1102 FUNCTION Lambda( is, Element, Material, Spring, D ) RESULT( lam )
1103
1104 INTEGER, INTENT(IN) :: is
1105 TYPE(typeElement), INTENT(IN) :: Element
1106 TYPE(typeMaterial), INTENT(IN) :: Material
1107 TYPE(typeSpring), INTENT(IN) :: Spring
1108 REAL(8), INTENT(IN) :: D(:,:)
1109 REAL(8) :: lam
1110
1111 REAL(8) :: cm,phi,pt
1112 REAL(8) :: PAI = 3.141592654D0
1113
1114 lam = 0.0
1115 IF(Element%type == 0) THEN
1116 CALL strength(is,Element,Material,Spring,cm,phi)
1117 pt = DTAN(phi*PAI/180.D0)
1118 lam = D(2,2)*Spring%stress(2,is)*Spring%dstrain(2,is) &
1119 + D(1,1)*(cm-Spring%stress(1,is)*pt)*pt*Spring%dstrain(1,is)
1120 ELSE
1121 lam = 2.0*D(2,2)*Spring%stress(2,is)*Spring%dstrain(2,is) &
1122 + 0.5*D(1,1)*Spring%stress(1,is)*Spring%dstrain(1,is)
1123 END IF
1124 END FUNCTION Lambda
1125 !-------------------------------------------------------------------------------
1126 SUBROUTINE strength( is, Element, Material, Spring, cm, phi )
1127
1128 INTEGER, INTENT(IN) :: is
1129 TYPE(typeElement), INTENT(IN) :: Element
1130 TYPE(typeMaterial), INTENT(IN) :: Material
1131 TYPE(typeSpring), INTENT(IN) :: Spring
1132 REAL(8), INTENT(OUT) :: cm
1133 REAL(8), INTENT(OUT) :: phi
1134
1135 REAL(8) :: h
1136 INTEGER :: ie1,ie2,im1,im2
1137
1138 h = Spring%hline(1,is) + Spring%hline(2,is)
1139 ie1 = Spring%element(1,is)
1140 ie2 = Spring%element(2,is)
1141 im1 = Element%material(ie1)
1142 im2 = Element%material(ie2)
1143 cm = (Material%c(im1)*Spring%hline(1,is) &
1144 + Material%c(im2)*Spring%hline(2,is))/h
1145 phi = (Material%phi(im1)*Spring%hline(1,is) &
1146 + Material%phi(im2)*Spring%hline(2,is))/h
1147 END SUBROUTINE strength
1148 !-------------------------------------------------------------------------------
1149 SUBROUTINE incRatio( Element, Material, Spring, Nonlinear )
1150
1151 TYPE(typeElement), INTENT(IN) :: Element
1152 TYPE(typeMaterial), INTENT(IN) :: Material
1153 TYPE(typeSpring), INTENT(IN) :: Spring
1154 TYPE(typeNonLinear),INTENT(INOUT) :: NonLinear
1155
1156 REAL(8) :: minStep,step,cm,phi,pt,a,b,c
1157 INTEGER :: is
1158 REAL(8) :: PAI = 3.141592654D0
1159
1160 minStep = 9999.D0
1161 DO is = 1, Spring%no
1162 CALL strength(is,Element,Material,Spring,cm,phi)
1163 IF(Element%type == 0) THEN
1164 pt = DTAN(phi*PAI/180.D0)
1165 a = Spring%dstress(2,is)**2 - (Spring%dstress(1,is)*pt)**2
1166 b = 2.D0*(Spring%stress(2,is)*Spring%dstress(2,is) &
1167 + (cm-Spring%stress(1,is)*pt)*Spring%dstress(1,is)*pt)
1168 c = Spring%stress(2,is)**2 - (cm-Spring%stress(1,is)*pt)**2
1169 ELSE
1170 a = (Spring%dstress(1,is)**2)/4.D0 + Spring%dstress(2,is)**2
1171 b = Spring%stress(1,is)*Spring%dstress(1,is)/2.D0 &
1172 + 2.D0*Spring%stress(2,is)*Spring%dstress(2,is)
1173 c = (Spring%stress(1,is)**2)/4.D0 + Spring%stress(2,is)**2 &
1174 - cm**2
1175 END IF
1176 step = ratio(a,b,c)
1177 IF(minStep > step) THEN
1178 minStep = step
1179 END IF
1180 END DO
1181 NonLinear%step = minStep
1182 NonLinear%rmin = NonLinear%rmin + minStep
1183 IF(NonLinear%rmin > 0.9998D0) THEN
1184 NonLinear%step = NonLinear%step - NonLinear%rmin + 1.D0
1185 NonLinear%rmin = 1.D0
1186 END IF
1187 END SUBROUTINE incRatio
1188 !-------------------------------------------------------------------------------
1189 FUNCTION ratio( a, b, c ) RESULT( step )
1190
1191 REAL(8), INTENT(IN) :: a,b,c
1192 REAL(8) :: step
1193
1194 REAL(8) :: r1,r2,b1
1195 REAL(8) :: eps = 1.0D-8
1196
1197 step = 9999.D0
1198 IF(DABS(a) <= eps) THEN
1199 IF(DABS(b)>= eps) THEN
1200 r2 = -c/b
1201 IF(r2 >0.D0) THEN
1202 step = r2
1203 END IF
1204 END IF
1205 ELSE
1206 b1 = b**2 - 4.D0*a*c
1207 IF(b1 >= eps) THEN
1208 r1 = (-b+DSQRT(b1))/(2.D0*a)
1209 r2 = (-b-DSQRT(b1))/(2.D0*a)
1210 IF(r1 < 0.D0) THEN
1211 IF(r2 > 0.D0) THEN
1212 step = r2
1213 END IF
1214 ELSE
1215 IF(r2 < 0.D0) THEN
1216 step = r1
1217 ELSE
1218 IF(r1 < r2) THEN
1219 step = r1
1220 ELSE
1221 step = r2
1222 END IF
1223 END IF
1224 END IF
1225 END IF
1226 END IF
1227 END FUNCTION ratio
1228 !-------------------------------------------------------------------------------
1229 SUBROUTINE totalValue( Element, Spring, Nonlinear )
1230
1231 TYPE(typeElement), INTENT(INOUT) :: Element
1232 TYPE(typeSpring), INTENT(INOUT) :: Spring
1233 TYPE(typeNonLinear),INTENT(IN) :: NonLinear
1234
1235 Element%disp = Element%disp + Element%ddisp*NonLinear%step
1236 Spring%strain = Spring%strain + Spring%dstrain*NonLinear%step
1237 Spring%stress = Spring%stress + Spring%dstress*NonLinear%step
1238 END SUBROUTINE totalValue
1239 !-------------------------------------------------------------------------------
1240 SUBROUTINE YieldCheck( Element, Material, Spring )
1241
1242 TYPE(typeElement), INTENT(IN) :: Element ! Element Data
1243 TYPE(typeMaterial), INTENT(IN) :: Material ! Material Data
1244 TYPE(typeSpring), INTENT(IN) :: Spring ! Spring Data
1245
1246 INTEGER :: is
1247 REAL(8) :: cm,phi,pt,tau
1248 REAL(8) :: PAI = 3.141592654D0
1249
1250 DO is = 1, Spring%no
1251 IF(Spring%flag(is) /= 1) THEN
1252 IF(Spring%flag(is) ==2) THEN
1253 Spring%flag(is) = -2
1254 ELSE
1255 CALL strength(is,Element,Material,Spring,cm,phi)
1256 IF(Element%type == 0) THEN
1257 cm = cm*0.998D0
1258 pt = DTAN(phi*PAI/180.0)
1259 tau = ABS(Spring%stress(2,is))+Spring%stress(1,is)*pt
1260 IF(tau > cm) THEN
1261 Spring%flag(is) = 1
1262 END IF
1263 ELSE
1264 cm = cm**2*0.998D0
1265 tau = Spring%stress(1,is)**2/4.D0 + Spring%stress(2,is)**2
1266 IF(tau > cm) THEN
1267 Spring%flag(is) = 1
1268 END IF
1269 END IF
1270 END IF
1271 END IF
1272 END DO
1273 END SUBROUTINE YieldCheck
1274 !-------------------------------------------------------------------------------
1275 SUBROUTINE Output( iste, Element, Spring, Nonlinear, File)
1276
1277 INTEGER, INTENT(IN) :: iste
1278 TYPE(typeElement), INTENT(IN) :: Element
1279 TYPE(typeSpring), INTENT(IN) :: Spring
1280 TYPE(typeNonLinear),INTENT(IN) :: NonLinear
1281 TYPE(typeFile), INTENT(IN) :: File
1282
1283 INTEGER :: ie,is,i,i1,i2,i3
1284
1285 WRITE(File%lp,6000) iste,NonLinear%step,NonLinear%rmin
1286 WRITE(File%lp,6010)
1287 DO ie = 1, Element%no
1288 i1 = 1 + (ie-1)*Element%dof
1289 i2 = i1 + 1
1290 i3 = i2 + 1
1291 WRITE(File%lp,6020) ie,Element%disp(i1), &
1292 Element%disp(i2), &
1293 Element%disp(i3)
1294 END DO
1295 WRITE(File%lp,6030)
1296 DO is = 1, Spring%no
1297 WRITE(File%lp,6040) is,Spring%flag(is), &
1298 (Spring%stress(i,is),i=1,2), &
1299 (Spring%strain(i,is),i=1,2)
1300 END DO
1301 6000 FORMAT(//, &
1302 ' OUTPUT DATA :',//, &
1303 ' STEP NO.=',I3,' STEP =',E15.5,' RMIN =',E15.5)
1304 6010 FORMAT(//, &
1305 ' < DISPLACEMENT >',//, &
1306 ' ELEM.NO. U V', &
1307 ' T',/)
1308 6020 FORMAT(3X,I5,3E15.5)
1309 6030 FORMAT(//, &
1310 ' < STRESS & STRAIN >',//, &
1311 ' SPRING.NO. YIELD SIG-N TAU-S', &
1312 ' EPS-N GAM-S',/)
1313 6040 FORMAT(5X,I5,2X,I5,4E15.5)
1314 END SUBROUTINE OUTPUT
1315 !===============================================================================
1316 SUBROUTINE Solve( nfree, band, Tsf, Uv, error )
1317
1318 INTEGER, INTENT(IN) :: nfree
1319 INTEGER, INTENT(IN) :: band
1320 REAL(8), INTENT(INOUT) :: Tsf(:)
1321 REAL(8), INTENT(INOUT) :: Uv(:)
1322 INTEGER, INTENT(OUT) :: error
1323
1324 INTEGER :: ik,jk,lk,mk,in,jn,kn,nb1,nfre1,ibet,i1,ii,ij,ialp
1325 REAL(8) :: eps = 1.0D-8
1326
1327 error = 0
1328 nb1 = band - 1
1329 DO in = 2, nfree
1330 i1 = in - 1
1331 DO jn = 1, nb1
1332 ij = jn + (in-1)*band
1333 ialp = in + jn - band
1334 IF(ialp < 1) THEN
1335 ialp = 1
1336 END IF
1337 DO ik = ialp, i1
1338 jk = in - ik + 1 + (ik-1)*band
1339 lk = in + jn - ik + (ik-1)*band
1340 mk = 1 + (ik-1)*band
1341 IF(ABS(Tsf(mk)) <= eps) THEN
1342 error = 1
1343 RETURN
1344 END IF
1345 Tsf(ij) = Tsf(ij) - Tsf(jk)*Tsf(lk)/Tsf(mk)
1346 END DO
1347 END DO
1348 END DO
1349 DO in = 2, nfree
1350 i1 = in - 1
1351 ialp = in - band + 1
1352 IF(ialp < 1) THEN
1353 ialp = 1
1354 END IF
1355 DO ik = ialp, i1
1356 jk = in - ik + 1 + (ik-1)*band
1357 mk = 1 + (ik-1)*band
1358 Uv(in) = Uv(in) - Tsf(jk)/Tsf(mk)*Uv(ik)
1359 END DO
1360 END DO
1361 ii = 1 + (nfree-1)*band
1362 Uv(nfree) = Uv(nfree)/Tsf(ii)
1363 nfre1 = nfree - 1
1364 DO in = 1, nfre1
1365 jn = nfree - in
1366 kn = 1 + (jn-1)*band
1367 ibet = nfree - jn + 1
1368 IF(ibet > band) THEN
1369 ibet = band
1370 END IF
1371 DO ik = 2, ibet
1372 jk = ik + (jn-1)*band
1373 lk = ik + jn - 1
1374 Uv(jn) = Uv(jn) - Tsf(jk)*Uv(lk)
1375 END DO
1376 Uv(jn) = Uv(jn)/Tsf(kn)
1377 END DO
1378 END SUBROUTINE Solve
1379 END PROGRAM RBSM
図8.1 プログラムリスト
表8.1 副プログラム一覧表
| 内部副プログラム |
説 明 |
|---|
| Analysis | 変位・表面力計算などの解析コントロール |
| Area | 要素面積の計算．ただし，境界用要素の面積は0 |
| Assemble | ばね剛性行列の全体剛性行列への組み込み（重ね合わせ） |
| BandWidth | バンド幅の計算 |
| Bmatrix | 行列（変位－相対変位関係）の作成 |
| ClearArray | 全変位, 全荷重, 全ひずみ, 全表面力, 降伏判定フラッグのクリア |
| coefStiff | ばね剛性行列を作成するための付録1に示された係数の計算 |
| elasticStiff | 弾性時のばね剛性行列の作成 |
| formLoad | 荷重項の作成 |
| formSpring | 弾性ばね行列と非線形ばね行列作成のコントロール |
| formStiff | ばね剛性行列を作成するためのコントロール |
| gCenter | 要素図心座標の計算 |
| incRatio | 荷重増分率の計算 |
| InputData | 入力データの読み込み |
| Lambda | 除荷・負荷判定のためのの計算 |
| MechanismCheck | 最大変位を基にしたメカニズムのチェック |
| MohrCoulomb | モール・クーロン条件による塑性化後のばね行列の作成 |
| noSpring | ばね数（境界辺の数）のカウント |
| Output | 解析結果の出力 |
| Perpendicular | 垂線の計算 |
| plasticStiff | 塑性化後のばね剛性行列の作成 |
| Prepocessing | データの入力，ばね関連データや幾何データの作成，出力 |
| PrintData | 入力情報の出力 |
| ratio | 2次方程式を解いて荷重増分率を計算 |
| setBoundary | 拘束条件の処理 |
| setDisplacement | 解の配列に対する強制変位（拘束条件）のセット |
| Solve | 変形コレスキー法により連立方程式を解き，変位を求める |
| SpringData | ばね構成要素番号と節点番号の作成 |
| SpringStress | 表面力の計算 |
| Strain | ひずみ（相対変位）の計算 |
| strength | 境界辺上の平均的な強度定数の計算 |
| totalValue | 全変位，全ひずみ，全応力（表面力）の計算 |
| Unload | 除荷・負荷判定 |
| vonMises | ミーゼス条件による塑性化後のばね行列の作成 |
| YieldCheck | 降伏判定 |