A PROBABILISTIC METHOD OF DETERMINING
FATIGUE LIVES OF SOME SELECTED STRUCTURAL
COMPONENTS OF AN AIRCRAFT TO MAINTAIN THE
REQUIRED FLIGHT SAFETY LEVEL – AN OUTLINE
ZARYS PROBABILISTYCZNEJ METODY
OKREŚLENIA TRWAŁOŚCI ZMĘCZENIOWEJ
WYBRANYCH ELEMENTÓW KONSTRUKCJI
STATKU POWIETRZNEGO DLA ZACHOWANIA
WYMAGANEGO POZIOMU BEZPIECZEŃSTWA
LOTÓW
Henryk Tomaszek
1, Józef Żurek
2,
Sławomir Klimaszewski
3, Mariusz Zieja
4(1,2,3,4) Air Force Institute of Technology ul. Księcia Bolesława 6, 01-494 Warszawa, Poland E-mails: (1) henryk.tomaszek@itwl.pl, (2) jozef.zurek@itwl.pl,
(3) slawomir.klimaszewski@itwl.pl, (4) mariusz.zieja@itwl.pl
Abstract. The paper has been intended to outline a method of determining fatigue
life of a structural component of an aircraft for some assumed flight safety level. The results gained allowed of finding the density function of time (i.e. flying time) indispensable to exceed the permissible crack length. With this function determined, one could determine reliability of the component to be then used to find fatigue life of this structural component. Two solutions have been given consideration, both depending on the m coefficient in the Paris relationship, i.e. for
m = 2 and m ≠ 2.
Key words: crack, reliability, critical condition, fatigue life
Streszczenie. W artykule przedstawiony zostanie zarys probabilistycznej metody
wyznaczania trwałości zmęczeniowej elementu konstrukcji statku powietrznego wyznaczonej dla przyjętego poziomu bezpieczeństwa lotów. Otrzymane wyniki pozwoliły na wyznaczenie postaci funkcji gęstości czasu (nalotu) niezbędnego do przekroczenia dopuszczalnej długości pęknięcia. Mając wyznaczoną tą funkcję określono niezawodność elementu, którą następnie wykorzystano do wyznaczenia trwałości zmęczeniowej elementu konstrukcji. Rozpatrzono dwa przypadki rozwiązania zależne od współczynnika m występującego we wzorze Parisa tzn. dla
m=2 i m≠2.
1. Introduction
A matter discussed in the paper is a method to determine fatigue life of a structural component of an aircraft. The following assumptions have been made:
the component’s health/maintenance status has been determined with one parameter only, i.e. the length of a crack therein. The actual value of the parameter has been denoted with l;
any change in the crack length may only occur in the course of the system/device being operated;
in the case given consideration the Paris formula takes the following form:
2 2 max m m m m k z l CM dN dl
(1) where: C, m – material constants,Nz – a variable that denotes the number of the component-affecting load cycles due
to the system’s vibration,
Mk – coefficient of the finiteness of the component’s dimensions at the crack
location,
σmax – maximum load defined with equation (2);
the load upon the structure’s component, with the system’s vibration taken into account, is a destructive factor. Let us assume we’ve got a component-affecting-load spectrum, with account taken of vibration. The spectrum allows for the determination of:
- the total number of load cycles Nc in the course of one flight assumed
a standard cycle,
- maximum loads within thresholds in the assumed spectrum amount to
max 1
,
2max, …,
Lmax (the assumed number of thresholds in the spectrum is L),- the number of repetitions of specific threshold values of the loading during one flight (standard load) ni, where:
L i i cn
N
1 ; maximum values of loads within the assumed thresholds are found in the following way: a i sr i i
max (2) where: max i
- maximum value of the cyclic load within the i-th threshold,sr i
- mean value of the cyclic load within the i-th threshold,a i
The following frequencies of the occurrence of loads correspond to values thereof within the thresholds
1max,
2max, …,
Lmax:1 1
P
N
n
c
, 2P
2N
n
c
, …, L c LP
N
n
.Relationship (1) may be expressed against the flying time of the aircraft. Therefore, we assume that:
t
N
z
, (3) where:λ – the occurrence rate of load cycles upon the component, t – flying time of the aircraft.
In the case under consideration
t
1
,where:
∆t – the average duration of the vibration-attributable fatigue-load cycle.
The relationship (1) against the flying time takes the following form:
2 2 max m m m m kl
CM
dt
dl
. (4) Having applied the hitherto made assumptions, one can proceed to determine the relationship that describes the dynamics of the fatigue-crack growth, i.e. of the increase in its length.Let Ul,t denote the probability that at the time t (for the flying time equal to t) the
crack reaches the length l. With the above-shown notation used, the dynamics of the crack length increase can be described with the following difference equation:
t l l L t l l t l l t t l PU PU PU L U, 1 , 2 , ... , 2 1 , (5) where:
Pi – probability that the load
max
i
defined with equation (2) occurs, whereL
i
1
,
2
,
3
,...,
andP
1
P
2
P
3
...
P
L
1
,∆li – crack increment in time t for the load equal to
imax, wherei
1
,
2
,
3
,...,
L
.The increments are to be found on the grounds of the dependence (4). Equation (5) in function notation takes the following form:
L i i i t t lP
u
l
l
t
u
1 ,,
, (6) where:u(l,t) – the probability density function of the crack length, which depends on the
flying time of the aircraft.
The difference equation (6) can be rearranged in the following partial differential equation of the Fokker-Planck type [3]:
2 2
)
,
(
)
(
2
1
)
,
(
)
(
)
,
(
l
t
l
u
t
l
t
l
u
t
t
t
l
u
. (7)A particular solution to equation (7) is the crack-length density function of the following form: ) ( 2 ) ( 2
)
(
2
1
)
,
(
At t B le
t
A
t
l
u
, (8) where:B(t) – an average crack length for the aircraft’s flying time t, A(t) – crack-length variance for the aircraft’s flying time t.
Taking eq (4) into account, two different forms of the solution can be found, depending on the value of the m coefficient:
for m = 2, coefficients B(t) and A(t) are solutions to integrals [3]:
t t C e l dt t t B 0 0 1 ) ( ) (
2 (9)
t t C e C l dt t t A 0 2 2 2 0 1 2 1 ) ( ) (
2 (10) where: ] ) [( max 2 2 2 C
C ;
2 2 CMk C ; 2 2 max 4 max])
)
[(
(
]
)
[(
;]
)
[(
max 2
- the second moment of structure-affecting load;]
)
[(
max 4
- the fourth moment of structure-affecting load;2 2
max
)
])
[(
(
- the square second moment of the structural-component-affecting load. for m ≠ 2 the parameters of distribution take the following forms:
0 2 2 max 2 2 2 0
[(
)
]
]
2
2
[
)
(
t
l
m
CM
t
l
B
m m m m k m
(11) 2 2 0 2 2 max 2 2 2 0 max 2 max 2 [( ) ]) 2 2 ( ] ) [( ] ) [( 2 2 ) ( m m m m m m k m m m m m k CM t l m l CM m t A
(12)2. Finding the probability density function of the flying time until the
boundary condition of the crack length in the structural component
is reached
Using the density function of the crack length (9) dependant on the flying time of the aircraft, one can determine the probability that the actual length of the crack in the aircraft structure’s component exceeds the permissible value within the time interval (0, tN). The relationship is as follows:
d l du
l
t
dl
l
t
Q
,
,
(13) where:ld – the permissible value of the crack length as determined for some assumed risk
of failure to the structural component.
The probability density function of the flying time up to the moment the crack exceeds the permissible value will be determined by the following equation:
Q
t ld
t t f , (14) For m = 2, the component’s unreliability will be determined with the following equation:
dt e e C e l l e C l t l U t Q t t C t C t C d t C d
0 2 2 2 0 2 0 ] ) 1 ( ) ) 1 ( ( [ , 2 2 2 2
(15)where u(ld ,t) is determined with eq (17).
From eq (14) the following is found:
] ) 1 ( ) ) 1 ( ( [ , , 2 2 2 2 2 2 2 0 2 0 t C t C t C d t C d d e e C e l l e C l t l u l t f
(16) where: ) 1 ( )) 1 ( ( 2 2 2 0 2 2 2 2 0 2 2 0 21
))
(
2
1
(
2
1
)
,
(
Ct t C d e C l e l l t C de
e
C
l
t
l
u
(17)The way of finding the probability density function of time of exceeding the permissible condition (16) has been given in [4] pp. 87- 90.
) ( 2 2 2 )) ( ( 2 2 2 2 , , 2 2 2 0 2 2 2 2 0 t A t C m l C t B l C m t C m l t l u l t f m m m d m m m d d
where: ) ( 2 ) ( 2)
(
2
1
)
,
(
At t B l d de
t
A
t
l
u
(19)Relationship (18) determines the probability density function of fatigue life of the selected aircraft’s structural component under operational-conditions spectrum for the Paris formula of the m ≠ 2. Parameters A(t) and B(t) in relationships (18) and (19) have been determined with formulae (11) and (12).
3. Finding fatigue life of the structural component up to the assumed
flight-safety level
The formula for reliability of the aircraft’s structural component can be written down in the following form:
t d dt l t f t R 0 ) , ( 1 ) ( (20)Where probability density function f(t,ld) is determined with the formula (16) for
m = 2 and with the formula (18) for m ≠ 2.
Hence, the unreliability of the component will be given with the following equation:
t d dt l t f t Q 0 ) , ( (21)The integral (21) should be rearranged to a simpler form and the problem reduced to solving an indefinite integral:
f
(
t
,
l
d)
dt
(22) For m = 2, the following change has been made in the integrand:
" 2 " 2 2 2 0 2 0 " 1 " 2 2 2 0 2 0)
1
(
)
)
1
(
(
)
1
(
))
1
(
(
2 2 2 2
t C d t C t C t C de
C
l
l
e
l
e
C
l
e
l
l
(23)Expression “1” is to be replaced with expression “2”, and expression “2” is denoted with z: z e C l l e l t C d t C ) 1 ( ) ) 1 ( ( 2 2 2 2 2 0 2 0
Hence, 2 2 2 2 2 4 0 2 2 2 2 0 2 0 2 2 2 0 2 0 0)
1
(
2
]
)
1
(
[
)]
1
(
[
]
)
1
(
[
2
2 2 2 2 2 2
t C t C d t C t C t C d t Ce
C
l
e
C
l
l
e
l
e
C
l
e
C
l
l
e
l
dt
dz
2 2 2 2 2 4 0 2 2 0 2 2 2 0 2 0 2 2 3 0)
1
(
]
)
1
(
[
2
)
1
](
)
1
(
[
2
2 2 2 2 2 2
t C t C d t C t C t C d t Ce
C
l
e
l
e
l
C
l
e
e
l
e
l
C
l
dt
dz
2 2 2 0 2 2 0 2 0 0 ) 1 ( ] ) 1 ( [ 2 ) 1 ]( ) 1 ( [ 2 2 2 2 2 2 2 t C t C d t C t C t C d t C e l e l e l e e l e l l dt dz
dz e l e l e e l e l l e l dt t C d t C t C t C d t C t C 2 2 2 2 2 2 2 2 0 2 0 0 2 2 2 0 ] ) 1 ( [ 2 ) 1 ]( ) 1 ( [ 2 ) 1 (
Then, the substitution is made in the indefinite integral:
) 1 ( )] 1 ( [( 2 2 2 2 2 2 2 0 2 0 Ct t C t C d t C e e C e l l e C l
] ] ) 1 ( [ 2 ) 1 ]( ) 1 ( [ 2 ) 1 ( [ 2 2 2 2 2 2 2 2 0 2 0 0 2 2 2 0 t C d t C t C t C d t C t C e l e l e e l e l l e l
dz
e
e
C
l
z t C
))
1
(
1
2 2 2 2 0
(24)Therefore, the following is arrived at:
] ) 1 ( )] 1 ( [ ) 1 ( [ 2 2 2 2 2 2 2 2 0 2 2 0 t C t C t C d t C t C e e C e l l e e C l
] ) ) 1 ( ( 2 ) 1 ]( ) 1 ( [ 2 ) 1 ( [ 2 2 2 2 2 2 2 2 0 2 0 0 2 2 2 0 t C d t C t C t C d t C t C e l e l e e l e l l e l
dz
e
e
C
l
z t C
))
1
(
1
2 2 2 2 0
Then we get:
[
2(
2
1
)
[
(
2
1
)]
2 2]
2 0 2 2 0 t C t C d t C t Ce
C
e
l
l
e
e
C
l
] ] ) 1 ( [ 2 ) 1 ]( ) 1 ( [ 2 ) 1 ( [ 2 2 2 2 2 2 2 2 0 2 0 0 2 2 0 t C d t C t C t C d t C t C e l e l e e l e l l e l
e
dz
e
C
l
z t C))
1
(
1
2 2 2 2 0
[
(
1
)
[
(
1
)]
]
2
1
2 2 2 2 2 0 2 0 t C t C d t C t Ce
e
l
l
e
e
l
] ] ) 1 ( [ ) 1 ]( ) 1 ( [ ) 1 ( [ 2 2 2 2 2 2 2 2 0 2 0 0 2 2 2 0 t C d t C t C t C d t C t C e l e l e e l e l l e C l
e
dz
e
C
l
z t C))
1
(
1
2 2 2 2 0
[(
(
1
))
[(
(
1
)]
]
2
1
2 2 2 2 2 0 2 0 t C t C d t C t Ce
e
l
l
e
e
l
e
dz
e
l
e
l
e
e
l
e
l
l
e
C
l
z t C d t C t C t C d t C t C]
]
)
1
(
[
)
1
](
)
1
(
[
)
1
(
[
2 2 2 2 2 2 2 2 0 2 0 0 2 2 2 0
[
(
1
)
(
(
1
))
]
2
1
2 2 2 2 2 0 2 0 t C t C d t C t Ce
e
l
l
e
e
l
dz
e
l
e
l
e
C
l
e
l
e
l
e
e
l
z z d t C t C t C d t C t C t C
1 0 2 2 2 0 2 0 2 0[
(
1
)
]
)
1
(
)
)
1
(
(
)
1
(
1
2 2 2 2 2 2
We assume that:))
1
(
(
)
)
1
(
(
2 2 0 0
t C d d t Ce
l
l
l
e
l
Hence,
[
(
1
)
(
(
1
)]
]
2
1
2 2 2 2 2 0 2 0 t C t C d t C t Ce
e
l
l
e
e
l
dz
e
z
e
e
l
l
e
e
l
z t C t C d t C t C
]
1
]
)
1
(
(
[
)
1
(
1
[
2 2 2 2 2 0 2 0
which means we arrive at:
e
dz
z
z1
2
1
The indefinite integral (24) after rearrangements takes the form:
e
dz
z
z1
2
1
(25)Then, the second substitution has to be made in the integral (25), which should take the form:
w
z
z dz dw 2 1 w dw dz 2 wdw dz2 (26) The dependence (26) is inserted in the integral (25). Hence, the following is arrived at:
e
wdw
e
dw
w
w w21
22
1
2
1
(27)One more substitution:
2
2 2y
w
ydy
wdw
2
dy w y dw 2 2 dy dw (28) Hence, after inserting (28) in (27) the following integral is arrived at:
e dy y 2 2 2 1
(29)) 1 ( ) ) 1 ( ( ; ; 2 2 2 2 2 2 0 2 0 2 2 t C d t C e C l l e l z z w w y
)
1
(
)
)
1
(
(
2 2 2 2 2 0 2 0
t C d t Ce
C
l
l
e
l
w
)
1
(
)
)
1
(
(
2 2 2 2 2 0 0
t C d t Ce
C
l
l
e
l
w
2
2 2y
w
2 22w
y
22w
y
)
1
(
)
)
1
(
(
2
2
2 2 2 2 2 0 0
t C d t Ce
C
l
l
e
l
w
y
(30)Having inserted the results gained in the equation (20) and remembering about a suitable notation of the limits of integration, the following dependence for the reliability is arrived at:
) ( 2 2 2 1 1 ) ( t y y dy e t R
(31)where equation (30) should be substituted for the upper limit of the integral y(t) . The cumulative distribution function for the standard Gaussian (normal) distribution takes the form:
x y dy e x 2 2 2 1 ) (
With the above-shown dependence taken into account, the formula for the reliability of the structure’s component is expressed with the following equation:
)
)
1
(
)
)
1
(
(
2
(
1
)
(
2 2 2 2 2 0 0
t C d t Ce
C
l
l
e
l
t
R
(32)Hence, reliability of the structure’s component will be determined with the following dependence:
) 1 ( ] ) 1 ( [ 2 2 2 2 2 2 0 2 0 22
1
)
(
t C d t C e C l l e l ydy
e
t
Q
(33)Having found (assumed) the level of risk of a failure to the structure’s component, i.e. the level of exceeding the permissible value of the length of a crack in this component, we get: *
)
(
t
Q
Q
(34) Hence,
e dy Q y 2 * 2 2 1 (35)For the assumed value of Q*, the value of the upper limit of the integral (for which the integral on the right side of the equation (35) takes value Q*) is to be found in the standard Gaussian distribution tables.
Hence, the following dependence is arrived at:
)
1
(
)
)
1
(
(
2
2 2 2 2 2 0 0
t C d t Ce
C
l
l
e
l
(36))
1
(
)
)
1
(
(
2
2 2 2 2 2 0 0
t C d t Ce
C
l
l
e
l
We assume that * 2
Hence,)
1
(
)
)
1
(
(
2 2 2 2 2 0 0 *
t C d t Ce
C
l
l
e
l
(37)From (37) we can find time
t
*, for which the equality relation (34) takes place. Time t* will be the searched life of the structure’s component, i.e. it will be the aircraft’s flying time for the assumed risk of exceeding the permissible value of the crack length. We assume thatx
e
C2t
(38) Hence, ) 1 ( ) ) 1 ( ( 2 2 2 0 0 * x C l l x l d
(39)From (39) we can find x. With some specific value of x gained from the dependence (38), we can find t*:
x
e
C2t
x t C2 ln
2 *ln
C
x
t
(40) Formula (40) determines fatigue life of the aircraft structure’s component t* for the assumed risk of exceeding the boundary condition Q*.For m ≠ 2, the unreliability of the structural component can be determined with the following formula:
dt t A t A t B l t B t l u t Q t d d
0 2 ( ) ) ( ' )) ( ( ) ( ' , (41) where:
0 2 2 2 2 02
2
)
(
t
l
m
C
t
l
B
m m
;
2 2 0 2 2 2 2 02
2
2
2
)
(
m m m ml
t
C
m
l
C
m
t
A
;
l
t
u
d,
- function determined with the dependence (19).The integral (41) should be rearranged to a simpler form and the problem reduced to solving an indefinite integral:
f
(
t
,
l
d)
dt
(42) The following change has been made in the integrand for m ≠ 2: " 2 " 2 " 1 " 2
)
(
2
)
)
(
(
)
(
2
))
(
(
t
A
l
t
B
t
A
t
B
l
d
d
(43)Expression “1” is to be replaced with expression “2”, and expression “2” is denoted with z:
z
t
A
l
t
B
d
)
(
2
)
)
(
(
2 (44)The derivative of the relationship (44) is calculated
2 2
))
(
2
(
)
(
'
2
)
)
(
(
)
(
2
)
(
'
)
)
(
(
2
t
A
t
A
l
t
B
t
A
t
B
l
t
B
dt
dz
d
d
Hence,)
(
'
2
)
)
(
(
)
(
2
)
(
'
)
)
(
(
2
)
(
'
2
)
)
(
(
)
(
2
)
(
'
)
)
(
(
2
))
(
2
(
2 2 2t
A
l
t
B
t
A
t
B
l
t
B
t
A
l
t
B
t
A
t
B
l
t
B
t
A
dt
d d d d
The substitution is made in the indefinite integral (42):
dz
t
A
l
t
B
t
A
t
B
l
t
B
t
A
t
A
t
A
t
B
l
t
B
e
t
A
d d d z)]
(
'
2
)
)
(
(
)
(
2
)
(
'
)
)
(
[(
2
))
(
2
(
)
(
2
)
(
'
))
(
(
)
(
'
)
(
2
1
2 2
dz
l
t
B
t
A
l
t
B
t
A
t
B
t
A
t
A
t
A
t
B
l
t
A
t
B
e
t
A
d d d z)
)
(
)](
(
'
)
)
(
(
)
(
)
(
'
2
[
))
(
2
(
)
(
2
)
(
'
))
(
(
)
(
)
(
'
2
)
(
2
1
2
1
2
dz e l t B t A dz l t B t A e t A z d d z
) ) ( ( ) ( 2 2 1 ) ) ( ( )) ( 2 ( ) ( 2 1 2 1
For z l t B t A d 1 ) ) ( ( ) ( 2 ,After some rearrangement we arrive at:
e
dz
z
l
t
f
d1
z2
1
)
,
(
(45)Then, we introduce conversions analogous to those presented for the m=2 case. Having done suitable rearrangements, the following integral is arrived at:
e dy y 2 2 2 1
where y takes value determined with the dependence (46), since
)
(
2
)
)
(
(
;
;
2
2 2 2t
A
l
t
B
z
z
w
w
y
d
) ( 2 ) ) ( ( 2 t A l t B w d)
(
2
)
)
(
(
t
A
l
t
B
w
d2
2 2y
w
2 22w
y
)
(
2
)
)
(
(
2
t
A
l
t
B
y
d (46)Having inserted the results gained in the equation (20) and remembering about a suitable notation of the limits of integration, the following dependence for the reliability is arrived at:
) ( 2 2 2 1 1 ) ( t y y dy e t R
(47)Hence, the unreliability of the structural component is determined with the relationship
) ( ) ) ( ( 2 22
1
)
(
t A l t B y ddy
e
t
Q
(48)Having found (assumed) the level of risk of a failure to the structure’s component, i.e. the level of exceeding the permissible value of the length of a crack in this component, we get: *
)
(
t
Q
Q
(49) Hence,
e dy Q y 2 * 2 2 1 (50)For the assumed value of Q*, the value of the upper limit of the integral is to be found in the standard Gaussian distribution tables. In this way we get value of θ. Hence, we arrive at the equation that allows for the determination of the component’s life for the assumed risk level:
)
(
)
)
(
(
t
A
l
t
B
d
(51)In the relationship (51) we look for such a value of t*, for which the left side of the equation equals the right side thereof. Solving the dependence (51) in this way, we find life of the structural component of an aircraft we have been looking for.
4. Final remarks
A probabilistic method to determine fatigue lives of some selected structural components of an aircraft has been presented for an assumed flight safety level (reliability). For the needs of the deterministic approach the physical part of the study has been based on the Paris formula. Solution to this formula depends on the value of the m coefficient because of the crack growth rate in the component. Therefore, two solutions are accepted:
for m = 2, for m ≠ 2.
Some random operation-induced loading in the form of a load spectrum is a fatigue-provoking destructive factor in the model of the crack growth in a structural component. An assumption has been made in the study that the sequence of load cycles, as far as values thereof are concerned, remains of no effect upon the crack growth rate. All the dependences arrived at enable specific calculations, if we have values of material constants and data on the load spectrum.
References
1. Smirnov N.N, Ickovicz A.A.: Obsłuzivanije i remont aviacjonnej techniki po
sostojaniju. Transport, 1980.
2. Kałmucki W.S.: Prognozirowanije resursov detali maszin i elementov
konstrukcji. Kisziniev,1989.
3. Tomaszek H., Klimaszewski S., Zieja M.: Zarys probabilistycznej metody
wyznaczania trwałości zmęczeniowej elementu konstrukcji z wykorzystaniem wzoru Parisa i funkcji gęstości czasu przekraczania stanu granicznego. 35th
International Scientific Congress on Powertrain and Transport Means – European KONES 2009, Zakopane 2009.
4. Tomaszek H., Żurek J., Jasztal M.: Prognozowanie uszkodzeń zagrażających
bezpieczeństwu lotów statków powietrznych. Wydawnictwo Naukowe JTE –
PIB Radom 2008.
5. Tomaszek H., Ważny M.: „Zarys metody oceny trwałości na zużycie
powierzchniowe elementów konstrukcji z wykorzystaniem rozkładu czasu przekraczania stanu granicznego”. ZEM 2008.
Prof. Henryk Tomaszek DSc. Eng, Scientific Worker of Air Force Institute of Technology. Field of science: aviation transport, reliability, safety. Specialisation: Air Armaments, Machine Building and Operation, transport, safety and reliability,
DSc. Eng. Józef Żurek, professor at Air Force Institute of Technology in Warsaw. Specialisation: Machine Building and Operation, transport, safety and reliability,
Major Mariusz Zieja PhD. Eng, Air Force Institute of Technology, graduated from Military University of Technology in 2000. M.Sc. in Mechatronics specialized in Aircraft Avionics. In 2008 achieved Ph.D. in Mechanical Engineering. He is engaged in development and implementation of IT systems to support aircraft maintenance, safety and reliability management.
!["Ancient Hebrew Marriage Laws with special references to general Semitic Laws and Customs", E. Neufeld, London-New York-Toronto 1944 : [recenzja]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==)