SIMULAÇÃO NUMERICA DO DESMONTE DE ROCHAS POR EXPLOSAO

Prévia do material em texto

Mecânica das Rochas para Recursos Naturais e Infraestrutura 
SBMR 2014 – Conferência Especializada ISRM 09-13 Setembro 2014 
© CBMR/ABMS e ISRM, 2014 
 
SBMR 2014 
Simulação Numérica do Desmonte de Rochas por Explosão 
 
Marko A. L. Bendezú 
Instituto Tecgraf e Departamento de Engenharia Civil, Rio de Janeiro – RJ, Brasil, 
markini@tecgraf.puc-rio.br 
 
Celso Romanel 
PUC-Rio, Departamento de Engenharia Civil, Rio de Janeiro – RJ, Brasil,romanel@puc-rio.br 
 
RESUMO: Quando explosivos são usados para desmonte de rocha a fogo, a técnica mais comum é 
perfurar a rocha e colocar explosivos no interior dos furos. O resultado da detonação provoca o 
fraturamento dinâmico da rocha pela interação de ondas de choque e da pressurização do gás no 
interior dos furos e das fraturas à medida que estas se desenvolvem. O objetivo deste trabalho é 
investigar pelo método dos elementos finitos estendidos (XFEM) o fraturamento dinâmico de um 
maciço de rocha sã, homogêneo e isotrópico. O modelo constitutivo da zona coesiva foi 
considerado para representar o comportamento mecânico da rocha durante o fraturamento. Nesta 
pesquisa investigou-se apenas o fraturamento dinâmico do maciço rochoso causado pelos altos 
níveis de tensão gerados pelas ondas de choque, desconsiderando os efeitos da pressurização dos 
gases. O desmonte de uma bancada a fogo é estudado numericamente, discutindo-se os aspectos 
que influenciam a solução numérica do problema, como o número e a distribuição de fissuras 
preexistentes. 
 
PALAVRAS-CHAVE: desmonte de rochas, fraturamento por explosão, método dos elementos 
finitos estendidos. 
 
 
1 INTRODUÇÃO 
 
Quando explosivos são usados para desmonte 
de rocha a fogo, uma pequena massa de 
explosivo químico é transformada em um 
grande volume de gás, com geração de altas 
pressões e aumento da temperatura. O resultado 
da detonação provoca o fraturamento dinâmico 
da rocha pela interação de ondas de choque e da 
pressurização do gás no interior dos furos e das 
fraturas à medida que estas se desenvolvem 
(Fig. 1). O desmonte de rocha é afetado pelas 
propriedades do maciço rochoso, condições de 
carregamento e pela geometria do problema, 
tais como a existência de bordas livres e 
descontinuidades. 
 A maioria dos modelos para previsão de 
desmonte em rocha é baseada em formulações 
empíricas, dada à grande complexidade do 
problema. Os pesquisadores ainda hoje 
discutem diferentes hipóteses para explicar os 
mecanismos fundamentais responsáveis pelo 
fraturamento, apesar dos enormes esforços de 
pesquisas experimentais, teóricas e numéricas, 
feitas nos últimas décadas para melhor entender 
o fenômeno. 
 Os rápidos avanços nos métodos de 
modelagem numérica fizeram da simulação 
computacional uma ferramenta promissora para 
estudar os processos dinâmicos de fraturamento 
de rocha. Um dos métodos mais utilizados é o 
método dos elementos finitos, que tipicamente 
acompanha no tempo a evolução das fraturas, 
com atualizações frequentes da malha de 
elementos finitos para representar a nova 
geometria do material recém-fraturado. Esta 
metodologia, além de ser computacionalmente 
demorada e difícil pela necessidade da 
reconstrução constante de malhas, também 
resulta na perda de precisão numérica quando as 
variáveis de interesse (tensões, deslocamentos) 
são mapeadas e interpoladas da malha antiga 
para os pontos de Gauss e pontos nodais da 
malha nova. 
 
 
SBMR 2014 
 
Figura 1. (a) Propagação de onda de choque; (b) 
expansão dos gases no furo e nas fraturas. 
 O método estendido dos elementos finitos 
(Extended Finite Element Method - XFEM) 
permite a incorporação de enriquecimentos 
locais, i.e. de um conjunto de funções de 
interpolação enriquecidas que fornecem valores 
das variáveis de interesse com maior precisão e 
eficiência computacional. Além disso, é 
importante ressaltar que nesta metodologia a 
presença da fratura, e sua propagação no tempo 
através do maciço rochoso, não é 
geometricamente modelada e a malha de 
elementos não precisa ser reatualizada, como 
nas aplicações convencionais do método dos 
elementos finitos neste tipo de problema. 
 O método estendido dos elementos finitos foi 
aqui aplicado para simular o fraturamento 
dinâmico de rochas por explosão, utilizando o 
modelo constitutivo de zona coesiva para 
simular a propagação de fraturas. Os resultados 
numéricos obtidos foram comparados com 
aqueles computados pelo método convencional 
dos elementos finitos (Lima, 2001) e pela 
técnica de eliminação de elementos (Saharam e 
Mitro, 2008), na qual elementos da malha são 
removidos quando o critério de ruptura por 
tração de Rankine é satisfeito. 
 
 
 
 
2 MÉTODO ESTENDIDO DOS 
ELEMENTOS FINITOS (XFEM) 
 
O método estendido dos elementos finitos 
(XFEM), introduzido por Belytschko e Black 
(1999) e Moës et al. (1999), é uma ferramenta 
poderosa para simulação de descontinuidades, 
tais como fraturamento em rocha, incorporando 
em sua formulação o modelo constitutivo da 
zona coesiva (Moës e Belytschko, 2002). 
 Uma variação do XFEM foi apresentada por 
Song et al. (2006), denominada de técnica dos 
nós fantasmas, cuja diferença básica é a 
maneira de enriquecimento das funções de 
interpolação. Enquanto que na formulação 
original do XFEM graus de liberdade adicionais 
são introduzidos para simulação da cinemática 
do fraturamento, a técnica dos nós fantasmas 
utiliza a sobreposição de elementos sem 
necessidade de criar graus de liberdade extras. 
 
2.1 Técnica dos nós fantasmas 
 
Considere um corpo que está sendo fraturado 
(Fig. 2) e a correspondente discretização, 
mostrando um elemento finito interceptado por 
uma fratura. Para ter um conjunto de funções 
de interpolação, a parte dos elementos 
fraturados que pertence ao domínio real Ω0 é 
estendido para o domínio fantasma ΩP. Em 
seguida, o deslocamento no domínio real pode 
ser interpolado usando os graus de liberdade 
associados aos nós do domínio fantasma. A 
aproximação do campo de deslocamento é 
então dada por: 
 
0
0
{ , }
{ , }
( , ) ( ) ( ) ( ( ))
( ) ( ) ( ( ))
P
P
h
I I
I w w
J J
J w w
u X t u t N X H X
u t N X H X
+ −
− +
∈
∈
= +
−
∑
∑
f
f
 (1) 
 
 onde f(X) representa um conjunto de funções 
denominado level set, H(x) é a função 
Heaviside e w0+, w0-, wP+ e wP- são os nós 
pertencentes a 0+Ω , 0−Ω , P+Ω e P−Ω
respectivamente. Como pode ser observado na 
Fig. 2, elementos interceptados por fraturas têm 
nós reais (círculos sólidos) e fantasmas 
(círculos vazados). A descontinuidade do 
campo de deslocamentos causada pela fratura é 
SBMR 2014 
simulada fazendo-se integrações numéricas 
separadas nas áreas 0+Ω e 0−Ω do elemento, 
ambas delimitadas pela fratura e contendo o 
lado definido pelos nós reais. 
 Song et al. (2006) e Remmers et al. (2008) 
demonstraram que os resultados obtidos tem 
pouca dependência da malha de elementos 
finitos se esta for suficientemente refinada. 
 
 
Figura 2. Técnica dos nós fantasmas onde as áreas 
sombreadas 0
+Ω e 0
−Ω do elemento são integradas 
numericamente para representação da descontinuidade do 
campo de deslocamentos causada pela fratura. 
 
 Para a localização da fratura e o 
acompanhamento de sua propagação no tempo 
é utilizado o conjunto level-set, consistindo de 
duas funções Ф e Ψ (Fig. 3), a primeira das 
quais (Ф) descrevendo a trajetória de 
propagação da fratura e a segunda (Ψ) 
consistindo de uma superfície ortogonal que 
localiza a ponta da fratura. 
 
 
Figura 3. (a) Domínio Ω parcialmente cortado por uma 
fratura; (b) função Ф descreve a trajetória de propagação 
da fratura; (c) função Ψ localiza a ponta da fratura. 
 
2.2 Modelo da zona coesiva 
 
Na análise numérica da propagação defraturas 
é necessário utilizar um modelo constitutivo 
específico para a região próxima à ponta da 
fratura para levar em consideração as 
deformações inelásticas que ali ocorrem. 
Barenblatt (1959) e Dugdale (1960) propuseram 
o modelo da zona coesiva, na qual as 
deformações na ponta da fratura antes da 
propagação são avaliadas e a dissipação de 
energia ocorre em uma região finita nas 
proximidades. As tensões aplicadas sobre a 
superfície da fratura decrescem com o aumento 
da abertura da mesma, mas não subitamente 
para zero, sem causar singularidade na ponta da 
fratura, o que não permite, portanto, a utilização 
de um fator de intensidade de tensão (Fig. 4). 
 Várias versões do modelo de zona coesiva 
foram propostas na literatura, sendo as 
principais diferenças entre elas quanto à forma 
da resposta tensão versus deslocamento e as 
constantes usadas para descrição do modelo. 
 O modelo da zona coesiva é ativado quando 
um critério de dano for atingido, quando então a 
zona coesiva é estendida através do elemento, 
simulando a propagação da fratura. Na Fig. 5 o 
conceito esquemático de zona coesiva com 
amolecimento no modo misto é ilustrado 
(Camacho e Dávila, 2002). O critério do início 
do dano é ativado em função da tensão principal 
máxima, de acordo com a Eq. (2). 
 
max maxmax
max max maxmax
0, se 0 
1,
, se 0T
σ σσ
σ σ σ
= <
= =
= ≥
f (2) 
 
 onde Tmax é a resistência à tração da rocha e 
tensões positivas são interpretadas como de 
tração. 
 A evolução do dano descreve a taxa na qual 
a rigidez do material é degradado uma vez que 
o critério de dano seja atingido. Neste trabalho 
foi utilizado o critério de propagação de modo 
misto I-II definido pela seguinte lei de potência, 
considerando η=1, 
 
1sn
C C
n s
GG
G G
ηη
  
+ =   
   
 (3) 
 
 onde Gn, Gs representam a energia de 
deformação nos modos I e II, respectivamente, 
e ��
�
	e ��
�
 os valores necessários que devem ser 
atingidos para causar a propagação da fratura. 
 
 
SBMR 2014 
 
Figura 4. Modelo de zona coesiva de Barenblatt (1959). 
 
 
Figura 5. Resposta no modo misto no modelo de dano de 
Camacho e Dávila (2002). 
 
 
3 SIMULAÇÃO NUMÉRICA 
 
3.1 Detonação simultânea em dois furos 
 
Com auxílio do XFEM este exemplo investiga 
os efeitos da explosão simultânea em dois furos 
de detonação de raio 0.0254m, posicionados 
conforme a Fig. 6, que também apresenta a 
numeração das oito fissuras preexistentes ao 
longo do perímetro de cada furo, com 
comprimento igual a um terço do raio. A malha 
foi formada por elementos quadrilaterais de 4 
nós, sendo 37.040 elementos finitos e 428 
elementos infinitos (Fig. 7). O processamento 
numérico do problema foi feito com o programa 
computacional ABAQUS. 
 
 
Figura 6. Localização dos furos de detonação. 
 
 
Figura 7. Malha de elementos do modelo. 
 
 A rocha (granito) tem massa especifica 
ρ=2800 kg/m3, parâmetros elásticos E=60 GPa 
e ν=0.25, resistência à tração de 3MPa, energia 
de fraturamento no modo I e II de 42.54 e 16.58 
Pa·m, respetivamente. A variação com o tempo 
da pressão dos gases no interior do furo devido 
à detonação do explosivo é mostrada 
esquematicamente na Fig. 8. 
 
 
Figura 8. Variação com o tempo da pressão dos gases de 
explosão nas paredes do furo. 
 
 No instante t = 90.7 µs as fraturas 1 e 13, 
paralelas à face livre do modelo, se encontram 
(Fig. 9), após seguirem trajetórias levemente 
assimétricas devido à influência do padrão de 
discretização da malha de elementos finitos. 
 
furo 1
3
2
1
8
7
4
5
6
furo 2
11
10
9
16
15
12
13
14
face livre
0.4 m
0.4 m
SBMR 2014 
 
Figura 9. Fraturas 1 e 13 se interceptam em t = 90.7 µs. 
 
No instante t = 200 µs as fraturas 2-12 e 8-14 
também se interceptam (Fig. 10-a) enquanto 
que as fraturas 7 e 15 atingem a face livre, 
formando um bloco de rocha geometricamente 
delimitado pelas fraturas 7-2 e 12-15. 
Para o instante t = 467 µs após a detonação, 
quando a pressão no interior do furo atinge seu 
valor máximo P0 = 1000 MPa, as fraturas 
dominantes desviam-se significativamente da 
direção radial, conforme mostra a Fig. 10b. As 
fraturas 6 e 16 não interceptam a face livre, e 
seguem em propagação praticamente paralela à 
face livre. As fraturas restantes atingem os 
contornos do modelo numérico. 
A Fig. 11 mostra o campo de deslocamentos 
verticais calculado no instante t = 538 µs, 
observando-se que a rocha existente entre as 
fraturas 6 e 16 tende a ser arremessada para fora 
do maciço rochoso. Os valores máximos de 
deslocamento computados são da ordem de 
6mm, nas cunhas formadas entre as fraturas 6 e 
7 e entre as fraturas 15 e 16. 
Para o mesmo problema, um comportamento 
similar foi numericamente obtido pelo método 
convencional dos elementos finitos por Lima 
(2001), por meio da atualização de sucessivas 
malhas para acompanhamento da propagação 
das fraturas no tempo e elementos finitos 
singulares quarter-points para cálculo dos 
fatores de intensidade de tensão. A principal 
desvantagem deste tipo de simulação é o tempo 
dispendido para geração de malhas e 
mapeamento de variáveis entre elas, com perda 
gradual de precisão numérica neste processo. 
Para fins comparativos, a Fig. 12 apresenta a 
situação das fraturas computadas por Lima 
(2001) no tempo t = 538 µs, observando-se que 
o padrão do faturamento é o mesmo do 
observado na Fig. 11. 
 
Figura 10. Propagação de fraturas: (a) fraturas 2-12 e 8-
14 se interceptam em t = 200 µs; (b) fraturamento no 
tempo t = 467µs quando a pressão nos furos de detonação 
atinge a pressão de pico. 
 
 
Figura 11. Campo de deslocamentos verticais e 
propagação das fraturas no tempo t = 538 µs (ampliação 
de 10x). 
 
 
Figura 12. Propagação das fraturas no tempo t = 538 µs 
computadas com o método convencional dos elementos 
finitos (Lima, 2001). 
 
SBMR 2014 
Não há indicações claras na literatura sobre o 
número de fraturas radiais que se desenvolvem 
em torno do furo de detonação. Evidências 
experimentais com explosivos de alta energia 
de choque sugiram de 8 a 12 fraturas radiais 
principais (Ghost e Daemen, 1995), enquanto 
que estudos numéricos de Song e Kim (1995) 
tenham observado o desenvolvimento de 10 a 
12 fraturas dominantes. 
Assim, foi novamente realizada a simulação 
numérica do fraturamento da rocha, porém 
desta vez sem consideração da existência prévia 
de fissuras ao longo do perímetro dos furos de 
detonação. A Fig. 13 apresenta os resultados no 
t = 90.7 µs que, quando comparados com 
aqueles calculados anteriormente na Fig. 9, 
permite observar que muitas pequenas fraturas 
se desenvolvem inicialmente mas só algumas 
realmente predominam. A sugestão anterior de 
12 fraturas radiais parece se confirmar, embora 
o número de fraturas no furo 2 seja superior ao 
observado no furo 1. Adicionalmente, estas 
fraturas não estão simetricamente distribuídas 
ao redor do perímetro dos furos e, neste 
instante, 3 delas já se interceptam. 
Embora uma simulação de fraturamento sem 
fissuras preexistentes seja possível, alguns 
inconvenientes foram observados. O tempo da 
simulação é bastante maior quando comparado 
ao modelo com fissuras preexistentes e, neste 
exemplo, o programa computacional apresentou 
problemas de convergência sendo interrompido 
no tempo t = 404 µs (Fig. 14) devido a 
problemas numéricos com a interseção de 
fraturas. Nota-se também que, além dos efeitos 
de malha observados anteriormente para o caso 
de fissuras preexistentes, o problema 
nitidamente deixou de ser simétrico em virtude 
da criação assimétrica das fraturas dominantes. 
 
 
Figura 13. Propagação das fraturasno tempo t = 90.7 µs 
(modelagem XFEM sem fissuras preexistentes). 
 
Figura 14. Propagação de fraturas sem fissuras 
preexistentes no tempo t = 404 µs, na interrupção da 
análise. 
 
3.3 Fraturamento com critério de Rankine 
 
O mesmo problema foi resolvido com a técnica 
de eliminação do elemento (TEE) que, 
incorpora o critério de tração máxima de 
Rankine e o modelo constitutivo de material 
frágil para o granito. 
O carregamento nas paredes do furo 
considera novamente a aplicação de um pulso 
de pressão transiente, mas a pressão do gás nas 
superfícies da fratura, anteriormente 
desconsiderada, é indiretamente levada em 
conta nesta metodologia. As ondas de choque 
são a causa primária do fraturamento da rocha 
enquanto a energia da pressão de gás tem mais 
influência no processo de fragmentação da 
rocha. A técnica de eliminação do elemento 
utiliza o método convencional dos elementos 
finitos com solução explícita no tempo 
(Saharam e Mitri, 2008). 
O critério de ruptura de Rankine, ou da 
tensão de tração máxima, é utilizado (Fig. 15-
a). O elemento no qual o critério é satisfeito é 
removido imediatamente da malha, zerando-se 
as tensões nele atuantes, o que representa, no 
contexto desta metodologia, a criação e 
propagação de fraturas através do maciço 
rochoso (Fig. 15-b). 
O modelo de faturamento de material frágil é 
intrinsicamente incorporado, com o 
fraturamento fundamentalmente controlado pela 
resistência à tração da rocha. Nesta abordagem, 
assume-se que certa quantidade de energia (Gf) 
é absorvida pela formação da superfície da 
fratura, a qual pode ser calculada com base nas 
tensões de tração e deslocamentos na abertura 
da fratura. 
 
SBMR 2014 
 
Figura 15. (a) Critério de ruptura de Rankine; (b) 
Representação da técnica de eliminação do elemento. 
 
Os resultados numéricos obtidos estão 
mostrados na Fig. 16. Para facilitar a 
comparação visual com os resultados anteriores 
(elementos finitos convencional, XFEM), todas 
as figuras representam uma mesma área de 
dimensão 2.1m x 1.45m. 
Pode ser observado que ao redor do furo de 
detonação existe um anel de rocha muito 
fragmentada, sem a criação de fraturas 
dominantes, pelo simples fato que elementos 
finitos são removidos uma vez que seja atingido 
o critério de tração máxima de Rankine. A 
influência da face livre também pode ser 
observada à medida que as “fraturas” se 
propagam. 
O tempo de simulação computacional é 
baixo apesar de utilizar mais elementos do que 
na malha XFEM (solução implícita), porque a 
técnica de eliminação do elemento faz uso de 
solução explícita no tempo, o que lhe confere 
maior rapidez na solução das equações. 
No presente exemplo não foram 
consideradas fissuras preexistentes porque elas 
não teriam influência alguma na solução do 
problema, dado o caráter desta técnica, baseada 
na simples remoção de elementos e não no 
acompanhamento da fratura através de uma 
malha de elementos, como nas metodologias 
anteriores (XFEM, método convencional dos 
elementos finitos). 
 
4 CONCLUSÕES 
 
Neste trabalho o método dos elementos 
estendido (XFEM) e técnica de eliminação de 
elementos (TEE) foram empregados na análise 
do fraturamento de um maciço rochoso por 
explosão (desmonte a fogo). O conceito do 
método XFEM foi brevemente explicado, no 
qual se agregam funções de enriquecimento nos 
elementos finitos interceptados pela fratura para 
melhor representar a descontinuidade no campo 
de deslocamentos. 
 
 
 
 
 
 
Figura 16. Evolução das fraturas e fragmentação do 
maciço rochoso na técnica de eliminação de elementos. 
 
SBMR 2014 
Cada um destes elementos é substituído por 
dois outros elementos com nós fantasmas, 
conservando as mesmas funções de interpolação 
dos demais elementos intactos da malha, 
facilitando a implementação computacional do 
método XFEM nos programas baseados no 
método convencional dos elementos finitos. Os 
resultados numéricos do exemplo analisado 
neste trabalho sugerem que o método XFEM 
produz resultados satisfatórios para problemas 
de desmonte de rocha a fogo (por explosão), 
confirmada pelos resultados anteriormente 
obtidos por Lima (2001) com o método 
convencional dos elementos finitos. 
Adicionalmente, o XFEM, por não representar 
fisicamente a fratura na malha de elementos 
finitos, não necessita de atualizações frequentes 
de malha com a consequente perda de precisão 
dos resultados entre mapeamentos das 
variáveis. 
Também não é necessário que fissuras 
preexistentes sejam consideradas para iniciação 
do fraturamento e, computacionalmente, requer 
menor tempo de solução que no método 
convencional. Os resultados numéricos obtidos 
com o XFEM também foram comparados com 
aqueles determinados com a aplicação da 
técnica de eliminação de elementos, sendo 
qualitativamente condizentes. 
 
AGRADECIMENTOS 
 
O primeiro autor agradece ao Instituto 
Tecgraf/PUC-Rio e ao Conselho Nacional de 
Desenvolvimento Científico e Tecnológico 
(CNPq) por ter possibilitado a realização desta 
pesquisa no curso de doutorado em Engenharia 
Civil da PUC-Rio, na área de Geotecnia. 
 
REFERÊNCIAS BIBLIOGRÁFICAS 
 
Barenblatt, G.I. (1959). Mathematical theory of 
equilibrium cracks in brittle fracture. Advances in 
Applied Mechanics, v.7, p. 55-129, Academic Press, 
New York. 
Belytschko, T. e Black, T. (1999). Elastic crack growth in 
finite elements with minimal remeshing. International 
Journal for Numerical Methods in Engineering; 
45(5), pp. 601–620. 
Camacho, P. e Dávila, C. (2002). Mixed-mode 
decohesion finite elements for the simulation of 
delamination in composite materials. NASA/TM-
2002–211737, pp. 1–37. 
Dugdale, D.S. (1960). Yielding of steel sheets containing 
slits. Journal of Mechanics and Physics of Solids, v.8, 
p. 100-104. 
Ghost, A. e Daemen, J.J.K. (1995). Rock fragmentation 
in bench blasting – A numerical study. Rock 
Mechanics, Balkema, Rotterdam, 553-558. 
Lima, A.D.R. (2001). Análise Numérica do Fraturamento 
Dinâmico de Rochas por Explosão. Departamento de 
Engenharia Civil, Pontifícia Universidade Católica 
do Rio de Janeiro. 110p. Tese (Doutorado). 
Moës, N., Dolbow, J. e Belytschko, T. (1999). A finite 
element method for crack growth without remeshing. 
International Journal for Numerical Methods in 
Engineering; 46(1): pp. 131–150. 
Moës, N. e Belytschko, T. (2002). Extended finite 
element method for cohesive crack growth. 
Engineering Fracture Mechanics; 69, pp. 813-834. 
Remmers, J., Borst, R. e Needleman, A. (2008). The 
simulation of dynamic crack propagation using the 
cohesive segments method. Journal of the Mechanics 
and Physics of Solids, 56, pp. 70-92. 
Rots, J. G. e Blaauwendraad, J. (1989). Crack models for 
concrete: discrete or smeared? fixed, multi-directional 
or rotating? Heron, 34 (1), pp. 1-59. 
Saharam, M.R. e Mitri, H.S. (2008). Numerical procedure 
for dynamic simulation of discrete fractures due to 
blasting. Rock Mech. Rock Eng. 41 (5), pp. 641–670. 
Song, J. H., Areias, P.M.A. e Belytschko, T. (2006). A 
method for dynamic crack and shear band propagation 
with phantom nodes. International Journal for 
Numerical Methods in Engineering, 67, pp. 868–893. 
Song, J. e Kim, K. (1995). Blasting induced fracturing 
and stress field evolution at fracture tips. Rock 
Mechanics, Balkema, Rotterdam, pp. 546-552.