- Email: [email protected]
An adaptive procedure for fracture simulation in extensive lattice networks
~
Engineering Fracture Mechanics Vol. 54, No. 3, pp. 325-334, 1996
Pergamon
0013-7944(95)00200-6
AN ADAPTIVE SIMULATION
PROCEDURE
IN EXTENSIVE
Copyright !~ 1996 Elsevier Science Ltd. Printed in Great Britain. All rights reserved 0013-7944/96 $15.00+0.00
FOR FRACTURE LATTICE
NETWORKS
J. E. BOLANDER, JR Department of Civil and Environmental Engineering, University of California, Davis, California 95616, U.S.A. T. SHIRA1SHI Kyushu Electric Power Corporation, Fukuoka 810, Japan and Y. ISOGAWA River Bureau, Japanese Ministry of Construction, Tokyo 100, Japan Abstract--The potential for simulating fracture processes using lattice networks remains largely unrealized due to the computational expense associated with such approaches. The range of specimen sizes and mesh densities which can be practically explored is quite limited. This paper shows that, for localized fracture processes, both computing time and storage requirements can be greatly reduced by modeling only the fracture process zone and its immediate vicinity with a lattice network. The material surrounding the lattice network is represented using boundary elements and is assumed to be linear elastic. The model layout is updated as the fracture process evolves; updates are guided by a conventional fuzzy control scheme. Effectiveness of this adaptive procedure is shown through examples involving concrete fracture. Copyright © 1996 Elsevier Science Ltd.
1. INTRODUCTION LATTICEMODELSare useful for studying the basic features of fracture in a variety of materials [1]. When applied to modeling cement-based composites, lattice models have: (a) illustrated the major role aggregate-matrix interfacial strength has in deciding overall fracture properties [2, 3], (b) clarified the process of crack-face bridging in the wake of a macrocrack tip [2, 3] and (c) simulated local variances in fracture energy due to both material heterogeneity and the effects of strain gradient [4]. These studies indicate a great potential for simulating concrete fracture and point to the use of lattice networks in engineering advanced cement-based composites. Several concerns must be addressed on route to achieving such goals. For example, two-dimensional lattices may misrepresent important aspects of three-dimensional fracture processes [4]. There are also questions regarding the effects of mesh orientation and density. Without some means for reducing the computational expense, however, both resolving these issues and the practical application of lattice networks are beyond the reach of most computing facilities. Typically, only small changes are made to a lattice model per solution cycle. Iterative solvers, notably conjugate gradient techniques, are therefore attractive since a converged solution is generally a good initial estimate of the next solution. In this paper, however, a direct solution technique (skyline solver) is employed and we work towards reducing the number of computational degrees of freedom. We show that both computing time and storage requirements can be drastically reduced by modeling only the fracture process zone (FPZ) and its immediate vicinity with a lattice network; a less computationally demanding technique (e.g. boundary elements) is used to represent the surrounding material. The model design is periodically updated so as to adapt to the FPZ as it proceeds along its course. Such a procedure may be applied provided the FPZ remains rather localized during the loading history. Fortunately, this is the case in many practical instances, including fracture driven by flexural deformations. The effectiveness of this procedure is shown through examples involving concrete fracture. Program results also provide insight into the abilities of lattice models to simulate fracture. For instance, there is a notable absence of chaotic behavior in the numerical fracture process, at least for the parameter values adopted here. 325
326
J. E. BOLANDER, JR et al.
Test Specimen
Analysis Model
,.j I E
E~ E~
thickness b = 120 mm 90
mm]---
---I
360 mm--
Fig. 1. Concrete fracture specimen and analysis model.
2. FRACTURE SIMULATION 2.1. Analysis model Concrete fracture specimens tested by Wittmann et al. [5] are modeled as shown in Fig. 1. Loading is applied via displacement control at the points indicated in the figure. A two-dimensional lattice network of beam elements composes the central region where fracture is likely to occur. Here the concrete mesoscale structure is represented by three components: aggregate (a), matrix (m) and aggregate-matrix interface (i). The distribution of aggregate particles within each specimen is generated randomly, the overall statistical distribution agreeing with the actual test material. Beam properties are assigned according to their location relative to the three components [2], as indicated in Table 1, where E and ar represent Young's modulus and fracture strength, respectively. Super elements constructed from plane stress boundary element equations [6] are used to model the regions surrounding the lattice network. As shown in Fig. 2, compatibility is maintained along the lattice region borders by constraining the displacements of the perimeter lattice nodes to follow: ui= ~b~u~ i = 1,2,
(1)
where u~ are the displacement components associated with boundary element node l. ~b~are the shape functions associated with element node l. ~b,= - ½ ~ ( l - ~ ) ~62= (1 - 42)
(2)
q~3= ½~(1 -4- 4). 2.2. Lattice elements dependent on variable damage measures
In local coordinates, beam (lattice) element equilibrium equations are (1 - fl) E I 4~
21 41 !J[00~ ] = [~F,:1 , 0
(3)
where L is beam length, A is beam cross-section area and I is beam moment of inertia; L, A and I are the same for all lattice elements. 0~, 02 and A are the end rotations and relative axial Table 1. Assumed component properties Component Matrix (m) Interface (i) Aggregate (a)
E (GPa)
tr,. (MPa)
25 25 70
5 1.25 10
An adaptive procedure for fracture simulation in extensive lattice networks 1=2
327 l=3
Boundary element in parent domain
Fig. 2. Displacement compatibility along lattice region border.
deformation, respectively; M,, M2 and F are the corresponding end moments and axial force. is a damage parameter which is initially set to zero for all elements. During computations the beam stiffness matrix i,.~expanded and transformed to global coordinates in the usual manner. At each load stage, the effective stress, tr, acting in each element is computed according to:
~r = (1 - f~)
+ ~
'
,
(4)
where W is the beam section modulus, ~ is a parameter to control the influence flexure has on fracture and fl is a parameter for scaling beam effective stress to global stress levels. The significance of parameters ct and fl is discussed in refs [2--4]. The lattice element with the highest ratio R = tr/trr >_ 1 experiences a "fracture event". Element internal forces associated with this state change are released, and both the element and system stiffness matrices are reformed. If R < 1 for all elements, loading is incremented so that the next most critical element just meets the condition R = 1. The calculations proceed in this manner modifying the elastic properties of one lattice element at a time. Conventional algorithms reduce element stiffness to zero upon violation of the fracture criterion. In order to better model the three-dimensioned nature of the material and fracture process, the stiffness of fracturing elements is gradually reduced via damage parameter fZ, which is dependent on component type [4]. 2.3. Typical fracture response The relation between opening displacement (measured at the mouth of the prenotch) and load-point reactive force is shown in Fig. 3. This response is typical of other simulations and will be compared to experimental results later. Fracture, or partial fracture, of a lattice element causes a release of energy which can be computed from the changes in load-point reactions under the 10
4 Z
5
8. 6. 4.
O) re
2. O" 0
"L _---~ •
'
, t
. . . .
0.1
i
. . . .
0.2
i
. . . .
0.3
J
. . . .
0.4
i
.
0.5
.
.
.
.
.
.
.
0.6
Opening displacement (ram) Fig. 3. Load-openingresponse.
0.7
328
J.E. BOLANDER, JR et al.
imposed load-point displacements. By storing and later processing such information, the extent of the active FPZ between any two load-point displacements can be visualized. Figure 4 shows the distributions of fracture energy consumed over the intervals 1-5 indicated in Fig. 3. The distribution of energy consumed over the entire loading history is also shown. These contours are plotted using a log scale (i.e. consecutive contour energies differ by a factor of ten) with darker levels indicating higher energies. Calculations indicate that 97-99% of the fracture energy is being consumed in the formation of the discrete crack; outlying events consume relatively little energy. 3. ADAPTIVE REMESHING SCHEME 3.1. Overview Provided the fracture process is rather localized, as it is for this problem, only the active fracture region and its immediate vicinity need be modeled using the lattice. The boundary element constructs are modified as the fracture process travels along the ligament length (Fig. 5). While the entire lattice network is initialized at the start of computations, lattice elements are simply activated/deactivated as they enter/leave the "window" created by this surrounding elastic region. Again, constraint equations are used to insure compatibility along the window borders. As the fracture process advances, additional nodes are added in the wake of the fracture zone. Only the active mesh region and the nodes along the yet undamaged ligament are modified as the solution proceeds; the positions of peripheral nodes, and of those along the formed macrocrack, do not vary. This remeshing activity is shown in more detail in Fig. 6, where displacements have been scaled to zero. Beam elements which have undergone partial damage are plotted using thinner lines; those which have completely fractured, i.e. f~ = 1, are not plotted. It is not necessary to retain all the boundary elements in the wake of the fracture zone, yet this shows the crack trajectory quite clearly. 3.2. Remeshing control One challenge lies in determining when and to what extent the active lattice region should be updated. Here, each side of the rectangular lattice region is allowed to move independently, guided by a fuzzy control scheme. The fuzzy controller mimics the actions of a user who adjusts the lattice region (i.e. the position of each side) based on the number of and rate at which fracture events occur near each corresponding side. Events are presented to the controller sequentially and in equal size samples over the loading history. One such sample is shown in Fig. 7(a), where each small dot represents a fracture event. For example, we focus on controlling the coordinate of the leading side of the lattice region, yz. To express whether an individual event is near this side, any number of appropriate membership functions could be employed. Referring to Fig. 7(b), we have used: =
1,
0 ~ Oi ~< 6c
6o - 6;
m(3i)-- (~o~ ~ ' =
3c <_ 3i<_ Oo 6; >_ 60,
O,
(5)
where 31 is the distance of event i from the side in question; m(6i) indicates the degree to which event i is near the side. For a sample of size n, the number of events near the side is: n
N = Y'm(6,).
(6)
i=1
The degree to which N belongs to fuzzy subsets representing linguistic values (such as none, few and many), together with a fuzzy measure of the rate at which N is changing, provides the fuzzy input set U to a conventional controller similar to that outlined by Mamdani [7] and recently clarified by Kosko [8]. Correlation-minimum inference is used to define the mapping from fuzzy sets U to fuzzy output sets V, i.e. S ( U ) = V. The fuzzy associated memory (FAM) rules, or "expert opinions", necessary for control are formed during training using event data provided by the complete mesh model (i.e. from simulations using the entire lattice network). The centroidal
A n a d a p t i v e p r o c e d u r e f o r f r a c t u r e s i m u l a t i o n in extensive lattice n e t w o r k s
o ~A o "7. '-o
¢ t~ o~
O0
¢q O0
0
v--4
329
J. E. BOLANDER,JR et
330
al.
I
G I!
©
©
Fig. 5. Adaptive remeshing.
method is used to de-fuzzify output set V and produce a real value for movement, Aye. Coordinate values for all sides of the lattice region are updated in this manner after each sampling of events. When remeshing is performed, however, boundaries are adjusted to meet the lattice nodes, whose coordinates remain fixed over the course of the analysis. Because crack-face bridges may be present well behind the fracture front, special consideration is necessary to ensure the trailing side of the lattice region is not advanced beyond the macro-crack tip. Figure 8 shows the performance of the fuzzy controller during a typical training run; only intermediate updates are shown in the figure. Some outlying events are missed during the remeshing process. FAM rules could be decided to capture all the events. However, making the lattice region larger would increase computational expense without any significant gain in performance. As will be seen in the next section, missing some of these low energy events has virtually no effect on the important aspects of the fracture response.
I
i
Fig. 6. Remeshingto contain fracture process.
An adaptiveprocedurefor fracturesimulationin extensivelatticenetworks
331
(a)
event i ~
i ............... (b) •
"
5i
m(5) 1
Yzl .................. ::t
Yt_.
Yl
±
5
5c 50 27 Fig. 7. Controllinglattice region boundaries. (a) Fracture events withinlattice region. (b) Membership function representing"near".
Note that this remeshing process could be performed using a binary logic control scheme. However, the fuzzy approach leads to simple, compact coding and would be easier to extend to more complex situations, such as when adapting to three-dimensional fracture processes, or when curved fracture paths or crack branching might occur. The abilities of fuzzy systems to generalize when given incomplete or partly erroneous data make them better suited to handling more complex situations. 3.3 Computational! effectiveness and accuracy Results from the complete lattice model serve as a benchmark for testing the accuracy of the adaptive procedure. Figure 9 shows the damage patterns obtained from the complete mesh and adaptive mesh models, as well as the load-opening response according to each model. Both the damage patterns and the global response appear identical (there are minor differences in the damage patterns, as will be shown later). The experimental scatter from five specimens is also shown for comparison. Rather than reproducing the experimental response, our goal here has been to study the performance of the adaptive procedure. The boundary element regions are therefore given an effective shear modulus and effective Poisson ratio obtained from numerical simulations of uniform stretching of undamaged lattice material. This accounts for the difference in initi~,l slopes of the numerical and experimental load-opening response. Inability of the numerical model to reproduce certain other aspects of the experimental response is discussed in ref. [4]. Table 2 highlights the main advantages of the adaptive remeshing procedure. By reducing the number of computational degrees of freedom, there are other associated benefits such as reduced program output and post-processing requirements. The relative degree of computer time and
-LL .;..: ".:
Fig. 8. Performanceof fuzzycontroller.
332
J.E. BOLANDER, JR et al.
!
(a)
(b)
9
_
8-
i¸ >~
o
i
~
Experimental scatter
.>_ n"
> complete mesh
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
adaptive remeshing
Opening displacement(mm)
Fig. 9. Verification of adaptive procedure. (a) Damage patterns. (b) Load-opening response.
storage savings is highly problem dependent and would, in most cases, be even greater When simulating fracture in larger specimens. It is instructive to have a quantitative measure of how well the adaptive procedure reproduces the complete mesh results. Let the sequence of fracture events define a time scale (i.e. fracture events 1, 2, 3... correspond to times t = 1, 2, 3...)• The effective distance between two damage configurations at time t can then be represented using: l
D(0
=
K
I
~ ZIn,
- f~l,
(7)
where K is the total number of lattice elements; f~) and f~ are damage measures of element i in configurations 1 and 2, respectively• D(t) ranges from 0 to 1, the lower bound value indicating both configurations are identical at time t. The upper bound value will not be approached when comparing ordinary fracture processes. Figure 10(a) shows the variance of D during crack propagation in fracture specimens A, B and C, which differ only in particle layout. Results presented in Fig. 9 have been obtained using specimen A. All comparisons show an initial time period where configuration 1 (complete mesh) and configuration 2 (adaptive mesh) evolve identically• At some point, however, D becomes nonzero and then tends to increase over time. Rather than increasing monotonically, there are numerous brief intervals over which the two configurations diverge and then reunite. These small peaks arise whenever the local sequence of events is different for the two configurations. Slight differences in the effective stress field do occur between the complete mesh and the adaptive procedure. Since, in general, many elements are very close to fracturing at any given instant, it is not surprising that the local sequence of events may differ, even though the corresponding local damages eventually become identical. A more interesting, meaningful feature of these diagrams is the presence of relatively long time intervals over which D remains constant, except for the small variations just described. During these intervals, some of which span more than 100 fracture events, the two configurations are evolving almost identically. The magnitude of D is actually not a good performance indicator, since missing low energy peripheral events will surely increase D, though this seems to have little impact on the Table 2. Computational effectiveness of adaptive remeshing procedure Method Complete mesh Adaptive remeshing
Computing timet (hr)
Degrees of freedom
Mean half bandwidth
Required memory++ (MB)
31 1.9
5727 2201§
406 233§
18.6 4.1§
tComputations performed on Macintosh PowerPC Series computers. ,+Memory for storing global stiffness matrix (in double precision). §Indicates maximum value during solution process.
An adaptive procedure for fracture simulation in extensive lattice networks (a)
0.008 o
•
~
~- -~--
P°~';o~ ~,
o
~o
333
00°6: o oo,
o.oo3-.
A
B
O
,
,
f
~
c
(
o.oc,2t
'
0.001 t, I
1'
'
'
0
(b)
'
I
.
.
.
100
.
i
'
'
'
200
'
I
.
.
.
300
.
i
.
.
.
400
.
I
.
.
.
500
.
1
.
600
.
.
.
I
.
700
.
.
.
800
t
0.12 0.1. 0.08.
f " 0.06. 0.04 • 0.02-
A,B,C .
0
.
.
.
i
.
1O0
.
.
.
,
.
200
.
.
.
I
.
300
.
.
.
i
.
.
400
.
.
i
500
.
.
.
.
i
i
600
700
800
t Fig. 10. Effective distance between damage configurations, (a) for same particle layouts, (b) for differing particle layout.
future course of ~he fracture process. Rather, it is the presence of these plateaus which indicates the stability of the fracture process and viability of the adaptive procedure. Amid the disorder seen in the random material structure, there is an order which keeps damage configurations from diverging despite differences in the effective stress fields and event sequences. To put these trends into better perspective, this same information is plotted in Fig. 10(b), along with a relation representing the distance between damage configurations evolving in different particle layouts. This new relation is labeled AB, since the complete mesh results from specimens A and B are being compared. From the onset of fracture, these two configurations diverge at essentially a con,;tant rate. The effective distances arising from trends A, B and C are virtually nothing by comparison. At least for the parameter values used here, the numerical fracture process is not chaotic. However, weak interfaciai strength (relative to matrix strength) establishes preferred pathways for crack propagation. Additional computations are being performed to determine whether the adaptive procedure remains effective for interfacial strengths approaching the matrix strength. Furthermore, bizs due to mesh density and orientation may be suppressing chaotic behavior. It is therefore unclear whether these observations extend to actual fracture processes, or are merely inherent to this type of lattice modeling. 4. SUMMARY Lattice models are useful for investigating the fracture behavior of many materials, including cement-based composites. Because lattice simulations are computationally intensive, however, the range of specimen sizes and mesh densities which can be practically solved is quite limited. For
334
J.E. BOLANDER, JR et al.
this main reason, a procedure for reducing computing time and storage requirements has been developed in this study. Effectiveness of this procedure is shown through examples involving concrete fracture. Program results also reflect on the ability of lattice models to simulate certain aspects of fracture. Conclusions are as follows: 1. Provided the fracture process is localized, only the fracture process zone and its immediate vicinity need be modeled with a lattice network; the surrounding material may be assumed linear elastic and represented by some appropriate technique such as boundary element methods. 2. The active lattice region and surrounding constructs must adapt to the evolving fracture process. This updating is accomplished using a conventional fuzzy control scheme. The fuzzy approach is effective, leads to compact coding and would be easier to extend to more complex situations, such as when following curved crack trajectories or when adapting to three-dimensional fracture processes. 3. Savings in computing time and storage requirements are highly problem dependent. Even for the rather small scale simulations presented here, however, more than an order of magnitude savings in computer time is realized. Greater savings are anticipated when addressing larger scale problems, three-dimensional simulations in particular. 4. Virtually identical solutions were provided by the adaptive meshing procedure and a model which utilizes the lattice mesh over the entire fracture domain. This suggests that the numerical fracture process is not chaotic. However, the remarkable stability of the fracture process may have been promoted by the material parameter values adopted here, particularly the relatively weak interfacial strengths. Furthermore, peculiarities of the lattice model, such as mesh bias on fracture direction, may be suppressing chaotic behavior. It is, therefore, not clear whether the observations concerning fracture stability apply to actual fracture processes as well.
REFERENCES [1] Statistical Models for the Fracture of Disordered Media (Edited by H. J. Herrmann and S. Roux). Elsevier Science (1990). [2] E. Schlangen and J. G. M. van Mier, Experimental and numerical analysis of micromechanisms of fracture of cement-based composites. Cem. Conc. Composites 14, 105-118 (1992). [3] E. Schlangen and J. G. M. van Mier, Simplelattice model for numerical simulation of concrete materialsand structures. Mater. Structures 25, 534-542 (1992). [4] J. E. Bolander, Jr, Mesoscopicanalyses of fracture in cement-basedcomposites. Univ. California, Davis, Dept. of Civil and Envir. Engng Report (1996). [5] F. H. Wittmann, H. Mihashi and N. Nomura, Size effect on fracture energy of concrete. Engng Fracture Mech. 35, 107-115 (1990). [6] C. A. Brebbia and J. Dominguez, Boundary Elements--An Introductory Course. Computational Mechanics Publications, Southampton, U.K. (1989). [7] E. H. Mamdani, Application of fuzzy algorithms for control of simple dynamical plant. Proc. lEE 121, 1585-1588 (1974). [8] B. Kosko, Neural Networks and Fuzzy Systems. Prentice-Hall, Englewood Cliffs, NJ, U.S.A. (1992). (Received 30 December 1994)
Engineering Fracture Mechanics Vol. 54, No. 3, pp. 325-334, 1996
Pergamon
0013-7944(95)00200-6
AN ADAPTIVE SIMULATION
PROCEDURE
IN EXTENSIVE
Copyright !~ 1996 Elsevier Science Ltd. Printed in Great Britain. All rights reserved 0013-7944/96 $15.00+0.00
FOR FRACTURE LATTICE
NETWORKS
J. E. BOLANDER, JR Department of Civil and Environmental Engineering, University of California, Davis, California 95616, U.S.A. T. SHIRA1SHI Kyushu Electric Power Corporation, Fukuoka 810, Japan and Y. ISOGAWA River Bureau, Japanese Ministry of Construction, Tokyo 100, Japan Abstract--The potential for simulating fracture processes using lattice networks remains largely unrealized due to the computational expense associated with such approaches. The range of specimen sizes and mesh densities which can be practically explored is quite limited. This paper shows that, for localized fracture processes, both computing time and storage requirements can be greatly reduced by modeling only the fracture process zone and its immediate vicinity with a lattice network. The material surrounding the lattice network is represented using boundary elements and is assumed to be linear elastic. The model layout is updated as the fracture process evolves; updates are guided by a conventional fuzzy control scheme. Effectiveness of this adaptive procedure is shown through examples involving concrete fracture. Copyright © 1996 Elsevier Science Ltd.
1. INTRODUCTION LATTICEMODELSare useful for studying the basic features of fracture in a variety of materials [1]. When applied to modeling cement-based composites, lattice models have: (a) illustrated the major role aggregate-matrix interfacial strength has in deciding overall fracture properties [2, 3], (b) clarified the process of crack-face bridging in the wake of a macrocrack tip [2, 3] and (c) simulated local variances in fracture energy due to both material heterogeneity and the effects of strain gradient [4]. These studies indicate a great potential for simulating concrete fracture and point to the use of lattice networks in engineering advanced cement-based composites. Several concerns must be addressed on route to achieving such goals. For example, two-dimensional lattices may misrepresent important aspects of three-dimensional fracture processes [4]. There are also questions regarding the effects of mesh orientation and density. Without some means for reducing the computational expense, however, both resolving these issues and the practical application of lattice networks are beyond the reach of most computing facilities. Typically, only small changes are made to a lattice model per solution cycle. Iterative solvers, notably conjugate gradient techniques, are therefore attractive since a converged solution is generally a good initial estimate of the next solution. In this paper, however, a direct solution technique (skyline solver) is employed and we work towards reducing the number of computational degrees of freedom. We show that both computing time and storage requirements can be drastically reduced by modeling only the fracture process zone (FPZ) and its immediate vicinity with a lattice network; a less computationally demanding technique (e.g. boundary elements) is used to represent the surrounding material. The model design is periodically updated so as to adapt to the FPZ as it proceeds along its course. Such a procedure may be applied provided the FPZ remains rather localized during the loading history. Fortunately, this is the case in many practical instances, including fracture driven by flexural deformations. The effectiveness of this procedure is shown through examples involving concrete fracture. Program results also provide insight into the abilities of lattice models to simulate fracture. For instance, there is a notable absence of chaotic behavior in the numerical fracture process, at least for the parameter values adopted here. 325
326
J. E. BOLANDER, JR et al.
Test Specimen
Analysis Model
,.j I E
E~ E~
thickness b = 120 mm 90
mm]---
---I
360 mm--
Fig. 1. Concrete fracture specimen and analysis model.
2. FRACTURE SIMULATION 2.1. Analysis model Concrete fracture specimens tested by Wittmann et al. [5] are modeled as shown in Fig. 1. Loading is applied via displacement control at the points indicated in the figure. A two-dimensional lattice network of beam elements composes the central region where fracture is likely to occur. Here the concrete mesoscale structure is represented by three components: aggregate (a), matrix (m) and aggregate-matrix interface (i). The distribution of aggregate particles within each specimen is generated randomly, the overall statistical distribution agreeing with the actual test material. Beam properties are assigned according to their location relative to the three components [2], as indicated in Table 1, where E and ar represent Young's modulus and fracture strength, respectively. Super elements constructed from plane stress boundary element equations [6] are used to model the regions surrounding the lattice network. As shown in Fig. 2, compatibility is maintained along the lattice region borders by constraining the displacements of the perimeter lattice nodes to follow: ui= ~b~u~ i = 1,2,
(1)
where u~ are the displacement components associated with boundary element node l. ~b~are the shape functions associated with element node l. ~b,= - ½ ~ ( l - ~ ) ~62= (1 - 42)
(2)
q~3= ½~(1 -4- 4). 2.2. Lattice elements dependent on variable damage measures
In local coordinates, beam (lattice) element equilibrium equations are (1 - fl) E I 4~
21 41 !J[00~ ] = [~F,:1 , 0
(3)
where L is beam length, A is beam cross-section area and I is beam moment of inertia; L, A and I are the same for all lattice elements. 0~, 02 and A are the end rotations and relative axial Table 1. Assumed component properties Component Matrix (m) Interface (i) Aggregate (a)
E (GPa)
tr,. (MPa)
25 25 70
5 1.25 10
An adaptive procedure for fracture simulation in extensive lattice networks 1=2
327 l=3
Boundary element in parent domain
Fig. 2. Displacement compatibility along lattice region border.
deformation, respectively; M,, M2 and F are the corresponding end moments and axial force. is a damage parameter which is initially set to zero for all elements. During computations the beam stiffness matrix i,.~expanded and transformed to global coordinates in the usual manner. At each load stage, the effective stress, tr, acting in each element is computed according to:
~r = (1 - f~)
+ ~
'
,
(4)
where W is the beam section modulus, ~ is a parameter to control the influence flexure has on fracture and fl is a parameter for scaling beam effective stress to global stress levels. The significance of parameters ct and fl is discussed in refs [2--4]. The lattice element with the highest ratio R = tr/trr >_ 1 experiences a "fracture event". Element internal forces associated with this state change are released, and both the element and system stiffness matrices are reformed. If R < 1 for all elements, loading is incremented so that the next most critical element just meets the condition R = 1. The calculations proceed in this manner modifying the elastic properties of one lattice element at a time. Conventional algorithms reduce element stiffness to zero upon violation of the fracture criterion. In order to better model the three-dimensioned nature of the material and fracture process, the stiffness of fracturing elements is gradually reduced via damage parameter fZ, which is dependent on component type [4]. 2.3. Typical fracture response The relation between opening displacement (measured at the mouth of the prenotch) and load-point reactive force is shown in Fig. 3. This response is typical of other simulations and will be compared to experimental results later. Fracture, or partial fracture, of a lattice element causes a release of energy which can be computed from the changes in load-point reactions under the 10
4 Z
5
8. 6. 4.
O) re
2. O" 0
"L _---~ •
'
, t
. . . .
0.1
i
. . . .
0.2
i
. . . .
0.3
J
. . . .
0.4
i
.
0.5
.
.
.
.
.
.
.
0.6
Opening displacement (ram) Fig. 3. Load-openingresponse.
0.7
328
J.E. BOLANDER, JR et al.
imposed load-point displacements. By storing and later processing such information, the extent of the active FPZ between any two load-point displacements can be visualized. Figure 4 shows the distributions of fracture energy consumed over the intervals 1-5 indicated in Fig. 3. The distribution of energy consumed over the entire loading history is also shown. These contours are plotted using a log scale (i.e. consecutive contour energies differ by a factor of ten) with darker levels indicating higher energies. Calculations indicate that 97-99% of the fracture energy is being consumed in the formation of the discrete crack; outlying events consume relatively little energy. 3. ADAPTIVE REMESHING SCHEME 3.1. Overview Provided the fracture process is rather localized, as it is for this problem, only the active fracture region and its immediate vicinity need be modeled using the lattice. The boundary element constructs are modified as the fracture process travels along the ligament length (Fig. 5). While the entire lattice network is initialized at the start of computations, lattice elements are simply activated/deactivated as they enter/leave the "window" created by this surrounding elastic region. Again, constraint equations are used to insure compatibility along the window borders. As the fracture process advances, additional nodes are added in the wake of the fracture zone. Only the active mesh region and the nodes along the yet undamaged ligament are modified as the solution proceeds; the positions of peripheral nodes, and of those along the formed macrocrack, do not vary. This remeshing activity is shown in more detail in Fig. 6, where displacements have been scaled to zero. Beam elements which have undergone partial damage are plotted using thinner lines; those which have completely fractured, i.e. f~ = 1, are not plotted. It is not necessary to retain all the boundary elements in the wake of the fracture zone, yet this shows the crack trajectory quite clearly. 3.2. Remeshing control One challenge lies in determining when and to what extent the active lattice region should be updated. Here, each side of the rectangular lattice region is allowed to move independently, guided by a fuzzy control scheme. The fuzzy controller mimics the actions of a user who adjusts the lattice region (i.e. the position of each side) based on the number of and rate at which fracture events occur near each corresponding side. Events are presented to the controller sequentially and in equal size samples over the loading history. One such sample is shown in Fig. 7(a), where each small dot represents a fracture event. For example, we focus on controlling the coordinate of the leading side of the lattice region, yz. To express whether an individual event is near this side, any number of appropriate membership functions could be employed. Referring to Fig. 7(b), we have used: =
1,
0 ~ Oi ~< 6c
6o - 6;
m(3i)-- (~o~ ~ ' =
3c <_ 3i<_ Oo 6; >_ 60,
O,
(5)
where 31 is the distance of event i from the side in question; m(6i) indicates the degree to which event i is near the side. For a sample of size n, the number of events near the side is: n
N = Y'm(6,).
(6)
i=1
The degree to which N belongs to fuzzy subsets representing linguistic values (such as none, few and many), together with a fuzzy measure of the rate at which N is changing, provides the fuzzy input set U to a conventional controller similar to that outlined by Mamdani [7] and recently clarified by Kosko [8]. Correlation-minimum inference is used to define the mapping from fuzzy sets U to fuzzy output sets V, i.e. S ( U ) = V. The fuzzy associated memory (FAM) rules, or "expert opinions", necessary for control are formed during training using event data provided by the complete mesh model (i.e. from simulations using the entire lattice network). The centroidal
A n a d a p t i v e p r o c e d u r e f o r f r a c t u r e s i m u l a t i o n in extensive lattice n e t w o r k s
o ~A o "7. '-o
¢ t~ o~
O0
¢q O0
0
v--4
329
J. E. BOLANDER,JR et
330
al.
I
G I!
©
©
Fig. 5. Adaptive remeshing.
method is used to de-fuzzify output set V and produce a real value for movement, Aye. Coordinate values for all sides of the lattice region are updated in this manner after each sampling of events. When remeshing is performed, however, boundaries are adjusted to meet the lattice nodes, whose coordinates remain fixed over the course of the analysis. Because crack-face bridges may be present well behind the fracture front, special consideration is necessary to ensure the trailing side of the lattice region is not advanced beyond the macro-crack tip. Figure 8 shows the performance of the fuzzy controller during a typical training run; only intermediate updates are shown in the figure. Some outlying events are missed during the remeshing process. FAM rules could be decided to capture all the events. However, making the lattice region larger would increase computational expense without any significant gain in performance. As will be seen in the next section, missing some of these low energy events has virtually no effect on the important aspects of the fracture response.
I
i
Fig. 6. Remeshingto contain fracture process.
An adaptiveprocedurefor fracturesimulationin extensivelatticenetworks
331
(a)
event i ~
i ............... (b) •
"
5i
m(5) 1
Yzl .................. ::t
Yt_.
Yl
±
5
5c 50 27 Fig. 7. Controllinglattice region boundaries. (a) Fracture events withinlattice region. (b) Membership function representing"near".
Note that this remeshing process could be performed using a binary logic control scheme. However, the fuzzy approach leads to simple, compact coding and would be easier to extend to more complex situations, such as when adapting to three-dimensional fracture processes, or when curved fracture paths or crack branching might occur. The abilities of fuzzy systems to generalize when given incomplete or partly erroneous data make them better suited to handling more complex situations. 3.3 Computational! effectiveness and accuracy Results from the complete lattice model serve as a benchmark for testing the accuracy of the adaptive procedure. Figure 9 shows the damage patterns obtained from the complete mesh and adaptive mesh models, as well as the load-opening response according to each model. Both the damage patterns and the global response appear identical (there are minor differences in the damage patterns, as will be shown later). The experimental scatter from five specimens is also shown for comparison. Rather than reproducing the experimental response, our goal here has been to study the performance of the adaptive procedure. The boundary element regions are therefore given an effective shear modulus and effective Poisson ratio obtained from numerical simulations of uniform stretching of undamaged lattice material. This accounts for the difference in initi~,l slopes of the numerical and experimental load-opening response. Inability of the numerical model to reproduce certain other aspects of the experimental response is discussed in ref. [4]. Table 2 highlights the main advantages of the adaptive remeshing procedure. By reducing the number of computational degrees of freedom, there are other associated benefits such as reduced program output and post-processing requirements. The relative degree of computer time and
-LL .;..: ".:
Fig. 8. Performanceof fuzzycontroller.
332
J.E. BOLANDER, JR et al.
!
(a)
(b)
9
_
8-
i¸ >~
o
i
~
Experimental scatter
.>_ n"
> complete mesh
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
adaptive remeshing
Opening displacement(mm)
Fig. 9. Verification of adaptive procedure. (a) Damage patterns. (b) Load-opening response.
storage savings is highly problem dependent and would, in most cases, be even greater When simulating fracture in larger specimens. It is instructive to have a quantitative measure of how well the adaptive procedure reproduces the complete mesh results. Let the sequence of fracture events define a time scale (i.e. fracture events 1, 2, 3... correspond to times t = 1, 2, 3...)• The effective distance between two damage configurations at time t can then be represented using: l
D(0
=
K
I
~ ZIn,
- f~l,
(7)
where K is the total number of lattice elements; f~) and f~ are damage measures of element i in configurations 1 and 2, respectively• D(t) ranges from 0 to 1, the lower bound value indicating both configurations are identical at time t. The upper bound value will not be approached when comparing ordinary fracture processes. Figure 10(a) shows the variance of D during crack propagation in fracture specimens A, B and C, which differ only in particle layout. Results presented in Fig. 9 have been obtained using specimen A. All comparisons show an initial time period where configuration 1 (complete mesh) and configuration 2 (adaptive mesh) evolve identically• At some point, however, D becomes nonzero and then tends to increase over time. Rather than increasing monotonically, there are numerous brief intervals over which the two configurations diverge and then reunite. These small peaks arise whenever the local sequence of events is different for the two configurations. Slight differences in the effective stress field do occur between the complete mesh and the adaptive procedure. Since, in general, many elements are very close to fracturing at any given instant, it is not surprising that the local sequence of events may differ, even though the corresponding local damages eventually become identical. A more interesting, meaningful feature of these diagrams is the presence of relatively long time intervals over which D remains constant, except for the small variations just described. During these intervals, some of which span more than 100 fracture events, the two configurations are evolving almost identically. The magnitude of D is actually not a good performance indicator, since missing low energy peripheral events will surely increase D, though this seems to have little impact on the Table 2. Computational effectiveness of adaptive remeshing procedure Method Complete mesh Adaptive remeshing
Computing timet (hr)
Degrees of freedom
Mean half bandwidth
Required memory++ (MB)
31 1.9
5727 2201§
406 233§
18.6 4.1§
tComputations performed on Macintosh PowerPC Series computers. ,+Memory for storing global stiffness matrix (in double precision). §Indicates maximum value during solution process.
An adaptive procedure for fracture simulation in extensive lattice networks (a)
0.008 o
•
~
~- -~--
P°~';o~ ~,
o
~o
333
00°6: o oo,
o.oo3-.
A
B
O
,
,
f
~
c
(
o.oc,2t
'
0.001 t, I
1'
'
'
0
(b)
'
I
.
.
.
100
.
i
'
'
'
200
'
I
.
.
.
300
.
i
.
.
.
400
.
I
.
.
.
500
.
1
.
600
.
.
.
I
.
700
.
.
.
800
t
0.12 0.1. 0.08.
f " 0.06. 0.04 • 0.02-
A,B,C .
0
.
.
.
i
.
1O0
.
.
.
,
.
200
.
.
.
I
.
300
.
.
.
i
.
.
400
.
.
i
500
.
.
.
.
i
i
600
700
800
t Fig. 10. Effective distance between damage configurations, (a) for same particle layouts, (b) for differing particle layout.
future course of ~he fracture process. Rather, it is the presence of these plateaus which indicates the stability of the fracture process and viability of the adaptive procedure. Amid the disorder seen in the random material structure, there is an order which keeps damage configurations from diverging despite differences in the effective stress fields and event sequences. To put these trends into better perspective, this same information is plotted in Fig. 10(b), along with a relation representing the distance between damage configurations evolving in different particle layouts. This new relation is labeled AB, since the complete mesh results from specimens A and B are being compared. From the onset of fracture, these two configurations diverge at essentially a con,;tant rate. The effective distances arising from trends A, B and C are virtually nothing by comparison. At least for the parameter values used here, the numerical fracture process is not chaotic. However, weak interfaciai strength (relative to matrix strength) establishes preferred pathways for crack propagation. Additional computations are being performed to determine whether the adaptive procedure remains effective for interfacial strengths approaching the matrix strength. Furthermore, bizs due to mesh density and orientation may be suppressing chaotic behavior. It is therefore unclear whether these observations extend to actual fracture processes, or are merely inherent to this type of lattice modeling. 4. SUMMARY Lattice models are useful for investigating the fracture behavior of many materials, including cement-based composites. Because lattice simulations are computationally intensive, however, the range of specimen sizes and mesh densities which can be practically solved is quite limited. For
334
J.E. BOLANDER, JR et al.
this main reason, a procedure for reducing computing time and storage requirements has been developed in this study. Effectiveness of this procedure is shown through examples involving concrete fracture. Program results also reflect on the ability of lattice models to simulate certain aspects of fracture. Conclusions are as follows: 1. Provided the fracture process is localized, only the fracture process zone and its immediate vicinity need be modeled with a lattice network; the surrounding material may be assumed linear elastic and represented by some appropriate technique such as boundary element methods. 2. The active lattice region and surrounding constructs must adapt to the evolving fracture process. This updating is accomplished using a conventional fuzzy control scheme. The fuzzy approach is effective, leads to compact coding and would be easier to extend to more complex situations, such as when following curved crack trajectories or when adapting to three-dimensional fracture processes. 3. Savings in computing time and storage requirements are highly problem dependent. Even for the rather small scale simulations presented here, however, more than an order of magnitude savings in computer time is realized. Greater savings are anticipated when addressing larger scale problems, three-dimensional simulations in particular. 4. Virtually identical solutions were provided by the adaptive meshing procedure and a model which utilizes the lattice mesh over the entire fracture domain. This suggests that the numerical fracture process is not chaotic. However, the remarkable stability of the fracture process may have been promoted by the material parameter values adopted here, particularly the relatively weak interfacial strengths. Furthermore, peculiarities of the lattice model, such as mesh bias on fracture direction, may be suppressing chaotic behavior. It is, therefore, not clear whether the observations concerning fracture stability apply to actual fracture processes as well.
REFERENCES [1] Statistical Models for the Fracture of Disordered Media (Edited by H. J. Herrmann and S. Roux). Elsevier Science (1990). [2] E. Schlangen and J. G. M. van Mier, Experimental and numerical analysis of micromechanisms of fracture of cement-based composites. Cem. Conc. Composites 14, 105-118 (1992). [3] E. Schlangen and J. G. M. van Mier, Simplelattice model for numerical simulation of concrete materialsand structures. Mater. Structures 25, 534-542 (1992). [4] J. E. Bolander, Jr, Mesoscopicanalyses of fracture in cement-basedcomposites. Univ. California, Davis, Dept. of Civil and Envir. Engng Report (1996). [5] F. H. Wittmann, H. Mihashi and N. Nomura, Size effect on fracture energy of concrete. Engng Fracture Mech. 35, 107-115 (1990). [6] C. A. Brebbia and J. Dominguez, Boundary Elements--An Introductory Course. Computational Mechanics Publications, Southampton, U.K. (1989). [7] E. H. Mamdani, Application of fuzzy algorithms for control of simple dynamical plant. Proc. lEE 121, 1585-1588 (1974). [8] B. Kosko, Neural Networks and Fuzzy Systems. Prentice-Hall, Englewood Cliffs, NJ, U.S.A. (1992). (Received 30 December 1994)