• Nie Znaleziono Wyników

Solution of a load diffusion problem by relaxation methods

N/A
N/A
Protected

Academic year: 2021

Share "Solution of a load diffusion problem by relaxation methods"

Copied!
20
0
0

Pełen tekst

(1)

Note No. 17

- - " - ^ T E S C H O O L /KUNDE, - DELFT

2 8 Jliil 1955

THE COLLEGE OF AERONAUTICS

CRANFIELD

SOLUTION OF A LOAD DIFFUSION PROBLEM

BY RELAXATION METHODS

by

G. VAISEY, M.A., and W. S. HEMP, M.A.

(2)

TECHNISCHE HOGESCHOOL

VLIEGTUIGBOUWKUNDE Kanaalstraat 10 - DELFT

(XlIxiaC-E ITOTE NO. 17 NOVKBER. 1954, T H E C O L L E G E Q F A E R O N A U T I C S C R A W F I E L D S o l u t i o n of a Load D i f f u s i o n Problem b y R e l a x a t i o n i l e t h o d s

-by-lliss Gillian Vaisey, ï.i.A,, and Prof, \ToS, Hemp, II,A,

S U 1.1 1.1 A a Y

The need to generalise the usual assiJLiptions made in the analysis of load diffusion problems has been emphasised by recent experimental work (Ref, 3)» v/hich has shown the import-ance of bending of the edge members. Direct mathematical solution of the plate problems, which arise, is hardly feasible and so in this report a numerical solution using the 'relaxation method' is carried out. Results show the method to be suitable for design purposes, but ccoparison tTith experiment still shows the need for further physical generalisations e These will form the subject of future Y/ork,

(3)

-2-Introduction

The problem of the diffusion of lo?.d from edge members into a panel of stringer reinforced skin has been the subject of numerous theoretical investigations (c.g, Refs 1,2), In the majority of these the assumption is made that the transverse direct strain component can be neglected, \7ith the consequence that the edge members remain straight. In some solutions for example that of Ref, 2, a general, two-dimensional system of stress and strain is assumed to exist in the panel, but in this case mathematical difficulties liiiiit the known solutions to the cases T/here the edge member has infinitely small flexural stiff-ness, Recent experimental evidence (Ref, 3) has shovm however, the importance of edge aembcr bending and the acccmpanying trans-verse direct stresses. It is therefore necessary for purposes of design to hcive available, methods, v/hich will allow for these effects. Since direct mathematical solution does not seem to be feasible, the present report tackles the diffusion problem by

the numerical method of 'relaxation',

Formulation of the Problem

The problem considered in this repoirb is illustrated in Fig, 1, It is the 'classical' diffusion problem of tlie

liters.ture, for the special case of edge members of constant area, The skin is shown as reinforced by stringers of area A and

s

pitch a , It TÓ.11 be assumed that these meinbers are distributed s

over the skin to form a uniform 'orthotropic' plate,

The components of displacement in the plate are v/ritten (u,v) and the corresponding strain components e , e and e are then given by

da 3v dv èu , .\

XX ÖX ' yy dy ' xy ox 3y ^ '

(4)

-3-Estimating these a.s resultants T., ï- and S per unit length we find, -n* EA Et / , \ s \ rp. = — (e + oe ) + e \ ^1 (1-0-2) ^ XX yy^ % ^ '

T„ = ^ \ (e + ere ) V (2)

2 (1-0-2) ^ yy xx^ ^

Et

^ •" 2(1+0-) ®xy

T/here t is the skin thickness, E is Young's modulus and C Boisson's Ratio, The stringers of coiorse only contribute to T. ,

The conditions of equilibrixjm for the plate are

^ T . .O AC 3T„

1_ ciS _ dS £ _ ^ /,s

ÖX dy ~ dx <3y ~ ••••• f\Jj

Substituting from (l) into (2) and from (2) into (3) we find,

^ -^ a t 72 •" f — 72 n"2-i a3^ = °

L s J ÖX V , ay ^ • -y ^ ^^^

••l-a-\ ö^v ^ 2 ^ /1+0-1 a^u / i - o - j a V a V /1+0:1

ax- ay- V - ^ ^^^y

= 0

Our problem is thus reduced to the solution of (4) subject to the appropriate boundary conditions at the edges of the plate,

The end rib at x = 0 v/ill be assumed to be completely flexible in bending» This implies that ( T, \ = 0 , A second

V V x = 0

condition at this edge follovre from the balance of shear input into the rib against the build up of end load within it. This last is given by EA'(e ) • Vfe thus find,

^ x=0

r (l-a-2)A ) ^

{1 + — — T ,i -r— -^ cr -r— = 0 J

) a t ( ax

dy I

^ ^ J \ at x = 0 ,,. (5)

'^^v t /av auV ^ j

(5)

_ ) ,

-The conditions assumed at x = I are those of complete fixity, ¥e shall thus ha>,ve,

u = v = 0 at x = l ,,, ..(6)

On the edge y = b we have first of all the shear input balance and secondly the relationi from the theory of beams, between the transverse loading ( Tp] and the displacement

(v) 6 These yield, -^"^ 'y=b 2

a u _ t (av au

^^2 - 2(14<T)A \dx •" dy/

f av

au\

^

I a7 + a7J /

^, at y = b , , . . , . , ( 7 )

d\ _ t / av au', (

ax^ ^ " (1-0-2)1 ^'y' W j

Since t h e p a n e l and J.oading a r e b o t h s y m m e t r i c a l about Ox, we need o n l y c o n s i d e r t h e r e g i o n y •^ 0 , T h i s r e q u i r e s b o u n d a r y c o n d i t i o n s a t y = 0 , Symmetry demands t h a t b o t h

(v) _^ and ( s ) „ s h o u l d be z e r o and s o ,

v = | ^ = 0 a t y = 0 (8)

Finally certain special conditions must hold at the comers. At x = I v/e assume the edge members fully built-in.

implying,-1 ^ = 0 at X = Z, y = b (9)

At X = 0 T/e have a given end load P in the edge member, We assume too that this member is pinned to the end rib, so that its bending moment v/ill be zero and its sheajr force equal to the end

load in the rib. These conditions yield,

au _ p_ a^v _ „ d \ _ A^ av ax ~ E A ' ^ 2 ~ » , 3 ~ ' " l ay

^^ ^^ (10) at X = 0, y = b

(6)

^^;^-THE REL^ÏXATIOINIAL S O I U T I O N

The specific problem solved

Since equations (4) to (10) are solved (approximately) by Relaxation Methods it is necessejry to select a skin having

specific shape and properties. The quantities chosen v/ere as

follows#-Length, I = 51 ins,

Width befr.TOen 4 of booms, 2b = 45 i n s , (_, 2 Boom a r e a , A = 2 , 2 3 i n s ,

,„ 4

Bocxi moment of i n e r t i a , I = 1.689 i n ,

Slcin t h i c k n e s s , t = 0,0813 i n ,

2

Stringer a r e a , A^^ = 0,0732 i n ,

S t r i n g e r spacing, a = 2 , 0 0 i n ,

s 2

End-rib a r e a . A' = 0,684 i n ,

Poisson's ratio, o- = 0,3

Yoiing's modulus, E = 1 0 ' p , s , i , )

Loading, F/A = 4000 p.s,i. < AE ~ ^^4

P

These are identical v/ith the quantities appropriate to the test specimen used in the experimental v/ork of Ref, 3»

Non-dimensional transformation of the eqi;iations

Since a numerical solution is contenplated it is necessary to put all the equations (4) to (10) into

non-dimen-sional form, Tfith I chosen as representative dimension let u,v,x,y,t,a^ = Z(ü,v,x,y,t,a^) j

A,A',A^ = l2ü,A',A )j

I = i^'i ,

Then a l l barred q u a n t i t i e s are purel;^'" numerical. Those

sub-s t i - t u t i o n sub-s made i n equationsub-s (4) to (10) leave thosub-se equationsub-s

u n a l t e r e d apart from the f a c t t h a t b a r r e d q u a n t i t i e s novir replace

t h e i r o r i g i n a l c o u n t e r p a r t s j i t i s therefore unnecessary to r e

-vjrite the equations. Hereafter v/e s h a l l v/ork i n terms of the

(7)

-6-Finite-difference approximations

The next step, in preparation for a relaxational solution, is to convert the governing equations and boundary conditions into finite-difference approximations. The chosen

(square) net at noda.1 points of v/hich values of u and v v/ere fovmd is shovm in figTJre 2, The mesh-size over most of the field is 2 inches, and in fact a solution v/as first obtained for which the net was liniformly of this size all over. The final

solution, incorporating values of u and v on the finer net (of mesh-size 1 inch), gave more detail in the region of stress-concentration and also provided sane slight indication of the accuracy to be expected from the results. The net is so arranged that net-lines fall on the boundaories y = b and x = Oj

'irregular stars' occur near the other tv/o boundaries. Clearly a rectangular net might ha.ve been used (at a cost in more cumber-seme arithr.aetic) v/hich would ha.ve presented no irregulax stars, 'Je decided against this since the conditions on the boundary X = I and centre-line y = 0 are particularly simple making the treatment of 'irregular stars' there s-traightfomTard,

Finite-difference approximations to the equations (4) to (10) v/ill be set dovm next, for use on any square net of mesh-size h inches. Let h = Zh so that h represents the non-dimensional mesh-size,

The numbering scheme of figure 3 TO-H be used to indicate, by suffix notation, the relative positions of values of the functions at a typical nodal point, 0, of the net and at surrounding nodal points.

The governing equations

The governing equations (4) in non-dimensional form and v/ith coefficients evaluated are

1.40967 ^ + 0,35 ^ -^ 0.65 - ^ = 0 )

ax

dy dxdy

( (4A.)

n xc o v a v _ ^^ a U ^ ;

0,35 —^2 + ~i2 "^ ^'^ ~r~r ~ ° ^

(8)

_ " 7 —

T y p i c a l a p p r o x i m a t i o n s needed i n t h e s e equa-tions a r e t h o s e f o r a u / a x Diic

were u s e d ,

' f d u 'i ,

-a ü / -a x ;:!nd -a v/dxèy -an,d the f o l l -a i / i n g s t -a n d -a r d e x p r e s s i o n s

r2/a2ü"\ ,

-h —rr = u^ + u , - 2u

,r2 /'a2v . .

^ ( 7 7 7 ^5 " ^6 -^ ^7 - ^8

\dxdy/^

F i n i t e " d i f f e r e n c e a p p r o x i m a t i o n s t o e q u a t i o n s (4A.) can now be w r i t t e n dovm and ( a f t e r m u l t i p l y i n g t h r o u g h by 2 f o r c o n -v e n i e n c e ) t-v/o r e s i d u a l s ( F - J and ('P-) , c o r r e s p o n d i n g -v/ith

V •' o "-~ / o t h e f i r s t and second r e s p e c t i v e l y of ( 4 A ) , a r e d e f i n e d b y

fp-'' = 2,8193 (Ü, + Ü J + 0»7 (ü„ + Ü, ) - 7.0387 ü^ ,.

+So,325(v -v.+v^-Vo)/ /

(p-)^ ^ 0-7 (v^+v3) + 2(v2 + v,^) - 5.4 v^ I

+jO,325(ü^-üg+ü_,-üg)^- -^

The r e a s o n f o r t h e s u f f i x e s u aaid v i n ( F - ) and F - ! i s

I u.f ( v; fairly obviousj the first expression in (11) defining /P-)

^ ' o

contains predominantly terns in u, the second, defining terns in v. The same notation v/ill be adopted v/hen V/e define residuals corresponding v/ith the boundary conditions,

The boundary conditions I, Along end-rib and boon All that has been described so faa- is stajidard

technique, details of v/hich can be found, for example, in Refs» 4 and 5<. We nov/ come to the bound.ary conditions. Those along the end-iT.b, equ.ations {d)) •'^^nd along the boom, equations (7)> together v/ith the conditions at the corner v/here the force is applied are somewhat unusual ajid present the crux of the problem relaxationally,

(9)

TtCHNISCHE HOGESCHOOL

VLIEGTUIGDüüVv'KUNDE Kanaalstraat 10 - DELFT

Leaving aside the corner x = O, y = b f o r the p r e s e n t consider the conditions along end-rib and boon, In p a r t i c u l a r l e t us examine equations (5)« I n t h e i r non-dinaensional form and v\dth c o e f f i c i e n t s evaluated these become

1.4097 ^ + 0,3 ^ =

'^ J^ ,... ( ...(5A)

^ . 0.0457(51)! 4 - 4 i =0 )

dy l^ax dy j •-'^

Equations (5A) hold on the vertical boundary x = 0, a piece of which is shov/n in Figure 4.

First derivatives appear in these equaitions and the simplest approximation to use, exer.Tplified by — , is

this.-dx

2 h ( 4 ) f Ü - Ü (12)

Vax^ ^

An approximation for the second derivate has already been stated so equations (5A) could now be put into finite-difference form. Applied at a point 0 on the boundary (figiore 4)> values of u , and V, appear in the resulting expressions and the node 3 is a 'fictitious' node, i.e. lies outside the field. But we need a pair of equations at 0 v/hich do not involve ü and v at fictitious nodes, and the most obvious way to proceed is to employ the governing equations (4A) together v/ith the boundary conditions (5A) eliminating be-tv/een them all vaJ.ues at fictitious nodes. (Cf, Ref, 6 ) , The governing equations, applied at 0 in figure 4, introduce values of ü and v at a further number of fictitious points and it beccmes necessary to use the boundary conditions not only at 0 itself but at other points sudh as 2 and 4 in order to get a pair of equations to be satisfied dt Oo The same method of elimination could be used along the boom

V7here equations (7) hold. This procedure was investigated at first but v/as abandoned (perhaps inadvisedly) because the

(10)

-9-along the boom (v/hei-e a fourth dei'ivative is present in one of the boundary conditions). It v/as also difficult to see ha\7 to get adequate approximations near the comer where the force acts.

The method finally adopted on these boundaries was the simple one of satisfying the boundary conditions only and .avoiding the introduction of values at fictitious nodes by using

end-difference foraiulae vihere necessary in place of central-differences, j The accuracy of such a procedure ought to be investigated moreflilly before firm, reliance is placed on the results of this reportJ

TreatJJig equations (5-^-) in this v/ay xre replace the approxina.tion (l2) by

2'hl — j = 4Ü - 3Ü - Ü , (see figure 4 ) . ....«(13) vax.' .

o

f — I i s approximated s i m i l a r l y and the r e s i d u a l expressions for

V dx '

o

the end-rib boundary can then be defined by

{F-j = 5.6387 ü^ - 1,4097 ü^ - 4.2290 u^+f 0,3(v2-v^)/ indept,of h , ^

(•p-^ = 0.0914 v^ - 0.0229 v^ + (v2+v^)-2.0686v^+|0,0229(Ü2"ü^)jwhere i (14)

h = 1/51 > i t s value for tlie f i n e r n e t , ]

Consider equations (7) nextj i n non-dimensional forra and

v/ith c o e f f i c i e n t s evaluated they became

av au ^

. 2

-^ - 0,01402(51)

^^' L^x öyj (7A)

„ i _ ï . 0.05290(51)^ 1 4 + 0.3 4 i =0

dx i dy dx i

Now i t i s I — ) and j — j v/hich must be approximated by t h e

5-point end formula i l l u s t r a t e d i n ( 1 3 ) . Cen-fcral-differenoes are

used for a l l o t h e r d e r i v a t i v e s . The (standard) approximation

(11)

1 0

-used for — r was

dx''

- ^ { ^ = M v ^ + v ) - ( v + v ^ ^ ) . 6 v ^ .

Residuals may now be defined as

follov/s,-/ F - 1 5 0,0280Ü,- 0.0070{I^2 +(^1+^3)" 2,02ial^- lo.0070(v^-v )follov/s,-/,

M = 0.l058v^- 0.0264v^2^^^i'"S^"^V''l1^" ^•0793^^0,0079(^^-^3)^'i^''^^

' ° - 1 "1 for h = -TT , its value for the first net, /

I Some coefficients in both of (15) need modification when h = -r:^ , its value for the coarse netJ

Finally, in detail, the approximations used at the comer x = 0, y = b v/ill be given. Equations (IO) hold here| in non-dimensional form and v/ith data inserted these bee one

^ = 4/10^

dx

4|= 0 / (10^)

dx I

- ^ - 0,40497(51)^-^=0 1

ax-'

dy J

Corresponding v/ith the first of (IOA) an approximation using the 3-point end-formula v/as adopted, a residual F- being defined ty

(FA

=

l^ - Ü - 3n - - \ ,{2h)

^..(16)

V H'o ' '' ° 10^

The second tv/o equations in (IOA) v/ere used together to find a residua.l expression P-. The 3-point cnd-fomula v/as again employed, for dv/dy, and the approximation taken for a v/ax

2— —2

(derived v/ith a use of d v/aïc = O) v/as

(12)

-11-\7ith h = 1/51 a residual F- at this corner can then be

' V

defined by

/P^'l = 0.8099 V - 0.2025 v^2 + ^v^ " ^q " 1-6075 v^.

..:.e (17)

The boundary conditions, II, Along built-in end and centre-line The methods used for dealing \/ith the irregular stars which occur near the built-in end of the skin (x = Z) and near the centre-line (y = O) v/ere mostly standard and call for no detailed description. The boundary conditions v/ere here applied by using them to eliminate values of u and v at fictitious nodes which occur in equations (ll) when the point 0 lies on either of the net-lines immediately adjacent to tliese boundaries,

The relaxation process

Knowledge of relaxa.tional technique must be assumed here,* but it is perhaps appropriate to add a fev/ remarks about the technique as applied to the particular problem of this report,

The residual expressions (ll), (14), (15)> (16) and (17)> together with those (not given) for use at tlie centres of the

irregular stars discussed in the last section, form the basis of the relaxation solution,

It is possible, and most convenient, to work the problejn in stages, first relaxing u only for a time and then sv/itching over to releix v alone. The significance of the notation F-and F- nov/ becoraes apparent, Vilriile v/orking on ti, for example, v/e concentrate' on the F- i-esidual expressions ei'-iploying incomplete patterns (derived from these) which involve changes in F- residuals only J the terms in v, in cui-ly brackets, remain (temporarily) constant at their values obtained from the cmrrent distribution of V, This is a standard approach which is not invalidated by the particular boundary conditions of the problem \mder discussion,

(13)

TECHNISCHE HOGESCHOOL

VLIEGTUIGBOUWKUNDE Kanaalstraat 10 - DELFT 1 2

-The p r o b l e m i s of t h e k i n d vvhich r e q u i r e s m e t i c u l o u s trea-tment - e s p e c i a l l y on tlie boom and e n d - r i b b o x m d a r i e s ,

I n d i v i d u a l r e s i d u a l s axe e a s i l y red.uced t o a c c e p t a b l e magnitudes b u t l a r g e changes i n ü and v may s t i l l be n e c e s s a r y i n o r d e r

t o malce t h e sun of t h e r e s i d u a l s n e g l i g i b l y s m a l l ,

RESULTS

The complete solution is not recorded for reasons of space, hovrever figixres 5 and 6 give sane indication of the nature of the distributions of u and v over the field. Figure 5 shov/s the variation of u along a fev7 selected lines of constant xj figure 6 the variation of v along certain lines of constant y, u find V are veiy smooth functions and the curves irm directly through a largo nuLiber of plotted points,

Figures 7 shov/ the strains by means of contours | since values of e , e and e must be found by differencing these are iess accurate than u and v V7hich are computed directly, There is also the difficulty that u varies very rapidly near the comer \-jhere the force is applied v/hich makes close estimation of the gradient there impossible. It is perhaps v/orth recording that the preliLiinary solution found for u and v - v/orking on a net of mesh-size 2 inches throughout the field - v/as very little different from the final solution (here presented) v/hich incorpor-ated a finer net in the region of stress-concentration,

Figures 8 and 9 shov/ the variations of strains along bocm and end-rib necar the corner v/here the force is applied,

Conclusions

The numerical results obtained in this report shov/ tjiat it is quite feasible to tack.lc the diffusion problem, even in the complicated case where bending of the edge members is considered, by means of the relaxation method, Ccaiiparison v/ith the experimental evidence of Ref, 3 still shov/s hovvever that further complexity must be introduced into our assuniptions, before full agreement can be

(14)

•13-assumed betv/een p l a t e rjid edge member must be r e p l a c e d b y an e l a s t i c e l e m e n t . This and o t h e r developments must av/ait f u t u r e i n v e s t i g a t i o n , b u t e x p e r i e n c e t o d a t e s u g g e s t s t h a t t h e r e l a x a t i o n method v / i l l b e an adequate t o o l f o r t h e f u r t h e r v/ork,

REFEREICES 1 , Cox, H,L. 2 . I l a n s f i e l d , EcH, 3 . Brov/n, La Verne W. L, A l l e n , D.N.de G, 5 . Pox, L . , and Southv/ell, R.V, 6 , Fox, L. D i f f u s i o n of l o a d s i n t o monocoque s t r u c t u r e s . A.R.C, R. and M. i 8 6 0 , S t r e s s d i s t r i b u t i o n i n p a n e l s bounded b y c o n s t a n t s t r e s s edge members,

R.A.E, Report Structxares l 6 2 , 1954. /in expei'i-nental i n v e s t i g a t i o n i n t o some of t h e ijroblems a s s o c i a t e d v/ith s t r e s s d i f f \ i s i o n i n t h e v i c i n i t y of chord-vvdse c u t - o u t s i n t h e v.ing, and comparison v/ith e x i s t i n g t h e o r i e s ,

C o l l e g e of A e r o n a u t i c s R e p o r t 8 3 , 1954.

Relaxation methods, McGraw Hill, 1954.

Relaxation methods applied to engineering problems. Vila. BihaoTaonic analysis as

applied to the flexure and tension of flat elastic plates.

Ehiil,Trc\ns.Royal Society,A, Vol,239, 1945, pp.419-460.

Iviixed boundary conditions in relaxational treatment of biharmonic problems,

Proc, Royal Society, A, Vol, IG9, 1947,

pp. 535-543.

(15)

FIGS I & 2

STRINGER AREA A , SKIN THICKNESS t

SKETCH OF PROBLEM

FIG.

THE NET USED

(16)

3.4. & 5 •O» 6 2

"_A—a

I 2 i L—2.

GENERAL NUMBERING SCHEME

FIG.3. 6 2 3 7l-A 12,

r

BOUNDARY ( E N D R I B ) ^ 2 2 - 5 - y ) IN INCHES O 4 8 12 16 20 22-5 lOO -200 • 3 0 0 • 4 0 0 -500 - 6 0 0 i/> UI X O Z - 7 0 0 z ^ - 8 0 0 o - 9 0 0 lOOO - I I O O - I 2 0 0 • I 3 0 0 \ \ /

Y

1

\ \ ALONG w = 51" '^^' ^

4

' #

A

A

^

y

^"^ — 1 — .

u ALONG SELECTED LINES OF CONSTANT x

SCHEME IN POSITION ON END-RIB FIG.5.

(17)

I/) lU I

u

to

o

o

• l O O -or^r^ ALONG

V°^

y /

A /

J

1

^

^r

= ^ -- - ^

h

12 16 2 0 24 28 X IN INCHES 32 36 4 0 44 48 51

(18)

FIGS. 7a 7b & 7c

THE INSET DIAGRAM

SHOWS A REGION

NEAR THE CORNER

MAGNIFIED 4 TIMES

CONTOURS OF CONSTANT «xac (x 10*)

CONTOURS OF CONSTANT euvr ( x 10*)

FIG. 7a

FIG. 7b

THE INSET DIAGRAM

SHOWS A REGION

NEAR THE CORNER

— MAGNIFIED 4 TIMES

CONTOURS OF CONSTANT «acu ( x 10^)

(19)

FIGS. 8a & 8b. 4 0 3 0 2 0 H 10 H 0 ?? 0 0 -10 -20

r^

^ * « « ^ ^ _ _ _ ^^^ 1 \ l 2 3 4 5 6 7 8 \ \ \ ^ ac IN INCHeS 1

^__^asu-.-

. — - — "

«ocac AND (Zuu ALONG BOOM NEAR CORNER WHERE FORCE IS APPLIED

FIG.Sa. 0 - 5 0 - 1 0 0 -ISO _5nrtl / / i : ^ \ 3 i ^ \ 5 < ' i 7 B 1 X IN INCHES

ea-u ALONG BOOM NEAR

J

(20)

FIGS, 9a & 9b 4 0 3 0 2 0 IO

i

\ \ « x x \ ~~~. \ ' - ^

1 T f

^ H M O

Cacx AND « u u DOWN END-RIB NEAR CORNER WHERE FORCE IS APPLIED

( 2 2 5 - J ) IN INCHES FIG.9a. O - 5 0 - i n o •ISO 1 2

y

3 4

h —

5 6 7 8

eacu DOWN END-RIB NEAR

CORNER WHERE FORCE IS APPLIED

Cytaty

Powiązane dokumenty