Shear stress evaluation on blood cells using computational fluid dynamics

Abstract

BACKGROUND:

Thrombus formation and hemolysis are important factors in developing blood pumps and mechanical heart valve prostheses. These phenomena are induced by flow properties. High shear stress induces platelet and red cell damage. Computational fluid dynamics (CFD) analysis calculates shear stress of fluid and particle pathlines of blood cells.

OBJECTIVE:

We studied blood cell damage in a blood pump by using CFD analysis and proposed a method for estimating blood damage.

METHODS:

We analyzed a pulsatile blood pump that was developed as a totally implantable left ventricular assist system at Tokai University. Shear stress on blood cells throughout pulsatile blood pumps were analyzed using CFD software.

RESULTS:

Based on the assumption that the effect of shear stress on platelets is accumulated along the trace, we proposed a method for estimating blood damage using the damage parameter D. Platelet damage parameter is calculated regardless of the division time 𝛥t which is dependent on the calculation time step. The results of the simulations are in good agreement with Giersiepen’s equation obtained from the experiments.

CONCLUSION:

The history of shear stress on a particle was calculated using CFD analysis. The new damage parameter D yields a value close to that of Giersiepen’s equation with small errors.

Keywords

Computational fluid dynamics (CFD)shear stress thrombus formation hemolysis blood pump

1. Introduction

Thrombus formation and hemolysis are important factors in developing artificial organs such as blood pumps and heart valve prostheses. Hemolysis is a phenomenon in which hemoglobins are lysed from destroyed red blood cells. Thrombus is made by blood-clotting reaction. When peripheral blood vessel flow is filled with thrombi, it causes brain infarction in the brain, economy class syndrome in the leg and renal vein thrombosis in the kidney. Thrombus formation is related to blood vessel wall deformation caused by e.g. the damage of vascular endothelial cells, stagnation or eddy in blood flow, and activation of blood coagulation factors.

As thrombus formation is a complex cascade reaction, there are three systems, namely formation of insoluble fibrin that is formed by coagulation factors, white thrombus made of platelets agglutination and promotion of white cell adhesion caused by complement system. High shear stress of blood flow influences the von Willebrand factor (vWF) and GpIb on a platelet membrane surface. It then increases the concentration of Ca²⁺ in a platelet, and then vWF crosslinks with GpIb and GpIIb/IIIa. Therefore, high shear stress leads to platelets aggregation [1]. Suppressing high shear stress of blood flow is necessary not only for preventing hemolysis but also thrombus formation. To accelerate the development of artificial organs, many repetitions of in vitro and in vivo experiments should be avoided due to cost and time. In addition, the frequency of animal experiments should be suppressed for research ethics. Therefore, it is desirable that in vivo and in vitro experiments are changed to computer simulations as much as possible. Computational fluid dynamics (CFD) analyses have been applied for pulsatile flow blood pumps [2–6] and continuous flow blood pumps [7–18].

In our previous study, we proposed the evaluation function of blood cell damage for developing an intra-cardiac axial flow pump utilizing CFD and estimated hemolysis [19]. This evaluation function was based on the hemolysis evaluation equation results from in vitro experiments by Giersipen et al. [20]. The evaluation function was derived from the damaged function method [21] and the time discretized Giersipen’s equation. However, it was found that results of the evaluation function are markedly varied when discretization time-step is changed. Therefore, in this paper we propose a new equation that can be accommodated against time-step changes.

2. Method

2.1. CFD analysis

Shear stress on blood cells throughout pulsatile blood pumps was analyzed using CFD software. The software package ANSYS-CFX 14.0 (ANSYS Japan, Tokyo, Japan) was used in this study. The structure of the pulsatile pump is shown in Fig. 1. This pump was developed at Tokai University [22]. The pump stroke was 12 mm and the pump stroke volume was 55 ml. The fluid density and viscosity were set 1060 kg/m³ and 0.0047 Pa ⋅ s, respectively. Computational grids were 226104. These grids were used to solve the Navier–Stokes equations that were discretized by finite volume elements. The device Reynolds number of 2343 is computed based on the following equation [23]: $\begin{eqnarray}\displaystyle Re=\frac{4}{{\pi}{\nu}}\frac{SV}{d_{i}}\frac{N}{R} & & \displaystyle\end{eqnarray}$ (1) where 𝜈 is dynamic viscosity, SV is the chamber stroke volume, d _i is the inlet port diameter, N is the beat rate, and R is the ratio of diastolic time to cycle time. SV = 55 × 10⁻⁶ m³, N = 80∕60 s, 𝜈 = 4.43 × 10⁻⁶ m²∕s, d _i = 18 × 10⁻³ m and R = 0.5.

Fig. 1.

The structure of the blood pump and boundary condition: (a) diastole phase, (b) systole phase.

Large eddy simulation (LES) was performed to study pulsatile blood flow, because it is suitable for unsteady flow and is widely used in pulsatile blood flow [24–27]. The pump rate was 80 bpm, and 1 beat cycle was then 0.75 s. Two analyses were conducted. First, the inlet port was an opening at 10 mmHg and the outlet port was a wall in the diastole phase. Second, the inlet port was a wall and the outlet port was an opening at 125 mmHg in the systole phase. The diaphragm was a moving wall in both phases. The wall speed was attempted under the following equation: $\begin{eqnarray}\displaystyle y=0.006\{\cos (8.38t)-1\} & & \displaystyle\end{eqnarray}$ (2) $\begin{eqnarray}\displaystyle \qquad y:\text{pump stroke}∼[\text{m}],\quad t:\text{time}∼[\text{s}]. & & \displaystyle \nonumber\end{eqnarray}$

Figure 1 shows the moving diaphragm and y direction is the horizontal direction. The pump stroke is displayed in the figure. Figure 1(a) shows when y is at the maximum and Fig. 1(b) shows when y is at the minimum.

Using the particle tracking command of CFX, trace data of 10 particles were calculated. The particle diameter and density were 2 μm and 1060 kg/m³ respectively.

2.2. Calculation of shear stress of the flow field

Shear stress in the turbulent flow field equals the sum of Newtonian shearing stress and Reynolds’ stress. In the two-dimensional field, it is expressed as: $\begin{eqnarray}\displaystyle {\tau}={\mu}\frac{du}{dy}+(-{\rho}\overline{u^{\prime }v^{\prime }}). & & \displaystyle\end{eqnarray}$ (3) By the assumption of Boussinesq, Reynolds stress 𝜏_t is described using turbulent viscosity 𝜇_t, by the following equation: $\begin{eqnarray}\displaystyle {\tau}_{t}=-{\rho}\overline{u^{\prime }v^{\prime }}={\mu}_{t}\frac{d\overline{u}}{dy}. & & \displaystyle\end{eqnarray}$ (4) Using the above equations, the stress tensor is described by the following equation: $\begin{eqnarray}\displaystyle {\tau}_{ij}=\left[\begin{array}{@{}ccc@{}}{\tau}_{xx} & {\tau}_{xy} & {\tau}_{xz}\\ {\tau}_{yx} & {\tau}_{yy} & {\tau}_{yz}\\ {\tau}_{zx} & {\tau}_{zy} & {\tau}_{zz}\end{array}\right]=-({\mu}+{\mu}_{t})\left[\begin{array}{@{}ccc@{}}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial u}{\partial x}} & {\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}} & {\displaystyle \frac{\partial u}{\partial z}}+{\displaystyle \frac{\partial w}{\partial x}}\\[10.0pt] {\displaystyle \frac{\partial v}{\partial x}}+{\displaystyle \frac{\partial u}{\partial y}} & {\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial v}{\partial y}} & {\displaystyle \frac{\partial v}{\partial z}}+{\displaystyle \frac{\partial w}{\partial y}}\\[10.0pt] {\displaystyle \frac{\partial w}{\partial x}}+{\displaystyle \frac{\partial u}{\partial z}} & {\displaystyle \frac{\partial w}{\partial y}}+{\displaystyle \frac{\partial v}{\partial z}} & {\displaystyle \frac{\partial w}{\partial z}}+{\displaystyle \frac{\partial w}{\partial z}}\end{array}\right]. & & \displaystyle\end{eqnarray}$ (5) In addition, von Mises stress is expressed using Eq. (5) [17,18,28,29]: $\begin{eqnarray}\displaystyle {\tau}=\left[\frac{1}{3}\{{\tau}_{xx}^{2}+{\tau}_{yy}^{2}+{\tau}_{zz}^{2}-{\tau}_{xx}{\tau}_{yy}-{\tau}_{yy}{\tau}_{zz}-{\tau}_{zz}{\tau}_{xx}+3({\tau}_{xy}^{2}+{\tau}_{yz}^{2}+{\tau}_{zx}^{2})\}\right]^{0.5}. & & \displaystyle\end{eqnarray}$ (6) Each blood cell was assumed to pass along the pathline, and the 𝜏 values were saved as the shear stress history of each blood cell.

3. Results

3.1. Estimation of shear stress acting on particles

Particle trace lines are shown in Fig. 2. Particles flowed around the diaphragm-housing junction of the pump.

Fig. 2.

Particle trace lines in the pulsatile blood pump.

Turbulent viscosity 𝜇_t and Reynolds stress 𝜏_t from CFD results were applied to Eq. (6). Then von Mises stress of 10 particles’ data was averaged. The history of the average shear stress on particles is shown in Fig. 3. Shear stress in the diastole phase was relatively higher than in the systole phase, and the largest value was 0.81 Pa.

4. Discussion

4.1. Estimation of platelet damage by CFD analysis

The platelet damage is shown as the function of exposed shear stress and time. Gieresiepen et al. proposed the following equation [20]: $\begin{eqnarray}\displaystyle \frac{{\Delta}\mathit{LDH}}{\mathit{LDH}}[\%]=3.31\times 10^{-6}\times t_{\text{exp}}^{0.77}\times {\tau}^{3.075}. & & \displaystyle\end{eqnarray}$ (7)

Lactate dehydrogenase (LDH) is the cytoplasm enzyme released from broken platelets. It is proportional to the amount of destroyed platelets. t _exp and 𝜏 are expose time [s] and shear stress on platelets respectively. In Eq. (7), LDL of the denominator is LDL content in all platelets and 𝛥LDL of the numerator is damaged LDL content released by platelets.

Fig. 3.

History of shear stress on particles.

Equation (7) is expressed as a percentage, and is then changed to a ratio expression: $\begin{eqnarray}\displaystyle \frac{{\Delta}\mathit{LDH}}{\mathit{LDH}}=3.31\times 10^{-8}\times t_{\text{exp}}^{0.77}\times {\tau}^{3.075}. & & \displaystyle\end{eqnarray}$ (8)

Based on the assumption that the effect of shear stress on platelets is accumulated along the trace, Eq. (8) was summed along the particle trace. It is assumed that the damage parameter D associated with each particle has a value of zero initially [21,30]. When the particle enters the pump, the parameter increases monotonically due to the accumulation of platelet damage along the particle trace. When D reaches a value of one, all particles are broken. An increase in the damage of a particle over 𝛥t _i is described as: $\begin{eqnarray}\displaystyle d_{p,i}=3.31\times 10^{-8}\times {\Delta}t_{i}^{0.77}\times {\tau}^{3.075} & & \displaystyle\end{eqnarray}$ (9) where 𝛥t _i is time difference between t _i−1 and t _i.

The damage accumulation from time zero to t _i is given by the following equation: $\begin{eqnarray}\displaystyle D_{p,i}=D_{p,i-1}+(1-D_{p,i-1})d_{p,i}. & & \displaystyle\end{eqnarray}$ (10) The image of damage parameter D is shown in Fig. 4.

Fig. 4.

Image of damage accumulation.

The platelet damage index is defined as the mean damage of n particles as: $\begin{eqnarray}\displaystyle E=\frac{1}{n}\mathop{\sum }_{p}^{n}D_{p}. & & \displaystyle\end{eqnarray}$ (11)

Here, we consider the condition that particles are exposed on shear stress 1 Pa and time 1 s. In addition, time difference 𝛥t is substituted under three cases: 1 s, 0.1 s, and 0.01 s. The damages of platelets are different according to time difference 𝛥t, as shown in Table 1. From Table 1, despite of the same exposure shear stress and time, calculation results are changed by time difference 𝛥t. This is considered to be caused by Eq. (8). Equation (8) is differentiated with respect to t _exp: $\begin{eqnarray}\displaystyle \frac{d\left({\displaystyle \frac{{\Delta}\mathit{LDH}}{\mathit{LDH}}}\right)}{dt_{\text{exp}}}=3.31\times 0.77\times 10^{-8}\times t_{\text{exp}}^{-0.23}\times {\tau}^{3.075}. & & \displaystyle\end{eqnarray}$ (12)

Then, Eq. (12) is divided by Eq. (8): $\begin{eqnarray}\displaystyle \left.\frac{d\left({\displaystyle \frac{{\Delta}\mathit{LDH}}{\mathit{LDH}}}\right)}{dt}\right/\left(1-\frac{{\Delta}\mathit{LDH}}{\mathit{LDH}}\right)=\frac{3.31\times 0.77\times 10^{-8}\times t_{\text{exp}}^{-0.23}\times {\tau}^{3.075}}{1-3.31\times 10^{-8}\times t_{\text{exp}}^{0.77}\times {\tau}^{3.075}}. & & \displaystyle\end{eqnarray}$ (13)

Equation (13) shows that the ratio of 𝛥LDL∕LDL is changed by exposure time t _exp. As the result of damage evaluation is changed by time difference 𝛥t, it is difficult to evaluate platelets damage. Then we assumed that platelet damage is proportional to the amount of platelets.

Table 1

Prediction of damages on platelets among different 𝛥t (under 1 Pa shear stress and exposure time 1 s)

	No division	𝛥t = 0.1 s	𝛥t = 0.01 s
Damage of platelets D	3.31 × 10⁻⁶	5.62 × 10⁻⁶	9.55 × 10⁻⁶

As the amount of platelets is N(t), reduction ratio N(t)^′, which is a temporal reduction ratio of the generation amount, is proportional to N(t): $\begin{eqnarray}\displaystyle N^{\prime }(t)=\frac{dN(t)}{dt}=-f({\tau})N(t). & & \displaystyle\end{eqnarray}$ (14) Then, $\begin{eqnarray}\displaystyle N(t)=Ce^{-f\left({\tau}\right)t}. & & \displaystyle\end{eqnarray}$ (15)

When t = 0, the amount of platelets N(t) = 1. Then, we have $\begin{eqnarray}\displaystyle N(t)=e^{-f({\tau})t} & & \displaystyle\end{eqnarray}$ (16) where f(𝜏) should be determined appropriately afterwards.

As Eq. (16) shows the ratio of the remaining platelets exposed to shear stress, the ratio of the remaining platelets at 𝛥t _i is expressed as: $\begin{eqnarray}\displaystyle l_{i}=e^{-f({\tau}){\Delta}t_{i}}. & & \displaystyle\end{eqnarray}$ (17)

Therefore, the amount of remaining platelets can be given as: $\begin{eqnarray}\displaystyle L=\mathop{\prod }_{i=0}li. & & \displaystyle\end{eqnarray}$ (18)

Hence, the amount of destroyed platelets damage is given by the following equation: $\begin{eqnarray}\displaystyle D=1-L=1-\mathop{\prod }_{i=0}li. & & \displaystyle\end{eqnarray}$ (19)

The demonstration of no dependence on 𝛥t length is shown as: $\begin{eqnarray}\displaystyle D & = & \displaystyle 1-\mathop{\prod }_{i=0}li\nonumber\\ \displaystyle & = & \displaystyle 1-\{\exp (f({\tau}){\Delta}t_{1})\times \exp (f({\tau}){\Delta}t_{2})\times \cdots \times \exp (f({\tau}){\Delta}t_{n})\}\nonumber\\ \displaystyle & = & \displaystyle 1-\exp \{f({\tau})({\Delta}t_{1}+{\Delta}t_{2}+\cdots {\Delta}t_{n})\}\nonumber\\ \displaystyle & = & \displaystyle 1-\exp \{f({\tau})\times \text{exposure}\_\text{time}\}\nonumber\end{eqnarray}$ D is independent of 𝛥t.

Fig. 5.

Surviving platelets ratio calculated on the basis of Gieresiepen’s in vitro experiments.

Fig. 6.

Surviving platelets ratio calculated by the new equation that is independent of the length of discretized 𝛥t.

Coefficient f(𝜏) is derived from Eq. (14): $\begin{eqnarray}\displaystyle f({\tau})=-\frac{{\displaystyle \frac{dN(t)}{dt}}}{N(t)}. & & \displaystyle\end{eqnarray}$ (20)

Treating Eq. ((20)) as a time-depending form of f(𝜏), we find f(𝜏). We get f(𝜏) as follows: $\begin{eqnarray}\displaystyle f({\tau})=-\frac{1}{0.7}\ln \{1-3.31\times 10^{-8}\times 0.7^{0.77}\times {\tau}^{3.075}\}. & & \displaystyle\end{eqnarray}$ (21) Then, the amount of remaining platelets is given by the following equation: $\begin{eqnarray}\displaystyle N(t)=(1-3.31\times 10^{-8}\times 0.7^{0.77}\times {\tau}^{3.075})^{\frac{t}{0.7}}\quad (0\leqq t\leqq 0.7∼\text{s},∼{\tau}\leqq 255∼\text{Pa}). & & \displaystyle\end{eqnarray}$ (22)

Calculation results based on Eq. (8) and Eq. (22) are shown in Figs 5 and 6, respectively.

Figure 5 was obtained from in vitro experiments and Fig. 6 was obtained by our new equation, which is independent of the length of discretized 𝛥t _i defining time difference between t _i and t _i−1.

The difference between these figures is from −4.46% to 2.12% in region t ≦0.7 s, 𝜏 ≦ 255 Pa. This result shows that Eq. (22) has enough accuracy for the evaluation of platelets damage under shear stress.

5. Conclusion

The history of shear stress on a particle was calculated using CFD analysis. The new damage parameter D we proposed in this study is hardly affected by the length of discretized 𝛥t _i of computer simulation time step and yields a value close to that of Giersiepen’s equation with small errors.

Footnotes

Conflict of interest

The authors have no conflict of interest to disclose.

References

Ikeda

Murata

and Goto

, Von Willebrand factor-dependent shear-induced platelet aggregation: Basic mechanisms and clinical implications, Ann N Y Acad Sci 811 (1997), 325–336.

Kiris

Kwak

Rogers

and Chang

I.D.

, Computational approach for probing the flow through artificial heart devices, Journal of Biomedical Engineering-Transactions of the ASME 119(4) (1997), 452–460.

Wood

H.G.

Throckmorton

A.L.

Untaroiu

and Song

X.W.

, The medical physics of ventricular assist devices, Reports on Progress in Physics 68(3) (2005), 545–576.

Yang

X.L.

Liu

and Yang

J.M.

, Unsteady flow and diaphragm motion in total artificial heart, Journal of Mechanical Science and Technology 21(11) (2007), 1869–1875.

Medvitz

R.B.

Reddy

Deutsch

Manning

K.B.

and Paterson

E.G.

, Validation of a CFD methodology for positive displacement LVAD analysis using PIV data, Journal of Biomedical Engineering-Transactions of the ASME 131(11) (2009), 111009.

Giridharan

G.A.

Lederer

Berthe

Goubergrits

Hutzenlaub

Slaughter

M.S.

Dowling

R.D.

Spence

P.A.

and Koenig

S.C.

, Flow dynamics of a novel counterpulsation device characterized by CFD and PIV modeling, Medical Engineering & Physics 33(10) (2011), 1193–1202.

Wood

H.G.

Anderson

Allaire

P.E.

McDaniel

J.C.

and Bearnson

, Numerical solution for blood flow in a centrifugal ventricular assist device, ISAO 22(12) (1999), 827–836.

Tamagawa

Akamatsu

and Sano

, Predictions of index of hemolysis in shear blood flow (Effects of contact with wall and particle inertia), JSME International Journal Series B-Fluid and Thermal Engineering 43(2) (2000), 225–232.

Tansley

Vidakovic

and Reizes

, Fluid dynamic characteristics of the VentrAssist rotary blood pump, Artif Organs 24(6) (2000), 483–487.

10.

Anderson

J.B.

Wood

H.G.

Allaire

P.E.

McDaniel

J.C.

Olsen

D.B.

and Bearnson

, Numerical studies of blood shear and washing in a continuous flow ventricular assist device, ASAIO Journal 46(4) (2000), 486–494.

11.

Tsukamoto

Ito

Sawairi

Konishi

Yamane

Nishida

Masuzawa

Tsukiya

Endo

and Taenaka

, Computational fluid dynamics analysis of a centrifugal blood pump with washout holes, Artif Organs 24(8) (2000), 648–652.

12.

Burgreen

G.W.

Antaki

J.F.

Z.J.

and Holmes

A.J.

, Computational fluid dynamics as a development tool for rotary blood pumps, Artif Organs 25(5) (2001), 336–340.

13.

Song

X.W.

Throckmorton

A.L.

Wood

H.G.

Antaki

J.F.

and Olsen

D.B.

, Computational fluid dynamics prediction of blood damage in a centrifugal pump, Artif Organs 27(10) (2003), 938–941.

14.

Song

X.W.

Throckmorton

A.L.

Wood

H.G.

Antaki

J.F.

and Olsen

D.B.

, Quantitative evaluation of blood damage in a centrifugal VAD by computational fluid dynamics, Journal of Fluid Engineering-Transactions of the ASME 126(3) (2004), 410–418.

15.

Chan

W.K.

Wong

Y.W.

and Ding

, Numerical investigations of a centrifugal blood pump, Journal of Mechanics in Medicine and Biology 4(3) (2004), 237–255.

16.

Song

X.W.

Throckmorton

A.L.

Wood

H.G.

Allaire

P.E.

and Olsen

D.B.

, Transient and quasi-steady computational fluid dynamics study of a left ventricular assist device, ASAIO Journal 50(5) (2004), 410–417.

17.

Wood

H.G.

Anderson

Allaire

P.E.

McDaniel

J.C.

and Bearnson

, Numerical solution for blood flow in a centrifugal ventricular assist device, International Journal of Artificial Organs 22(12) (1999), 827–836.

18.

Yano

Sekine

Mitoh

Mitamura

Okamoto

Kim

D.W.

Nishimura

Murabayashi

and Yozu

, An estimation method of hemolysis within an axial flow blood pump by computational fluid dynamics analysis, Artif Organs 27(10) (2003), 920–925.

19.

Mitoh

Yano

Sekine

Mitamura

Okamoto

Kim

D.W.

Yozu

and Kawada

, Computational fluid dynamics analysis of an intra-cardiac axial flow pump, Artif Organs 27(1) (2003), 34–40.

20.

Giersiepen

Wurzinger

L.J.

Opitz

and Reul

, Estimation of shear stress-related blood damage in heart valve prostheses-in vitro comparison of 25 aortic valves, ISAO 13(5) (1990), 300–306.

21.

Bludszuweit

, Model for a general mechanical blood damage prediction, Artif Organs 19(7) (1995), 583–589.

22.

Okamoto

Hashimoto

and Mitamura

, Design of a miniature implantable left ventricular assist device using CAD/CAM technology, J Artif Organs 6 (2003), 162–167.

23.

Medvitz

R.B.

Reddy

Deutsch

Manning

K.B.

and Paterson

E.G.

, Validation of a CFD methodology for positive displacement LVAD analysis using PIV data, ASME Journal of Biomechanical Engineering 131 (2009), 111009.

24.

Mittal

Simmons

S.P.

and Udaykumar

H.S.

, Application of large-eddy simulation to the study of pulsatile flow in a modeled arterial stenosis, ASME Journal of Biomechanical Engineering 123 (2001), 325–332.

25.

Delorme

Anupimdi

Kerlo

A.E.

Shetty

Rodefeid

Chen

and Frankel

, Large eddy simulation of powered fontan hemodynamics, Journal of Biomechanics 46 (2013), 408–422.

26.

Jones

S.C.

Soliropoulos

Healy

T.M.

and Yoganathan

A.P.

, Numerical simulation of flow in mechanical heart valves: Grid resolution and assumption of flow symmetry, ASME Journal of Biomechanical Engineering 125 (2003), 709–718.

27.

Manosh

C.P.

Molla

M.M.

and Roditi

, Large-Eddy simulation of pulsatile blood flow, Medical Engineering & Physics 31 (2009), 153–159.

28.

Bludszuweit

, Three-dimensional numerical prediction of stress loading of blood particles in a centrifugal pump, Artif Organs 19(7) (1995), 590–596.

29.

Bludszuweit

, A theoretical approach to the prediction of haemolysis in centrifugal blood pumps, PhD Thesis, University of Strathclyde, Glasgow, 1994.

30.

Yokobori

, An Interdisciplinary Approach to Fracture and Strength of Solid, Wolters-Noordhoff Scientific Publications Ltd, Groningen, 1968, p. 249.