## CHAPTER B - B3 Development of Concept Improvements and Advances in Fundamental Understanding of Fusion Plasmas

3.4.1 Stationary MHD modes in magnetically confined plasmas

Background and Objectives: This activity aims at constructing equilibria and relaxed states of laboratory plasmas of fusion concern (e.g., plasmas of tokamaks and reversed-field-pinches) with flow or/and finite electrical conductivity in connection with improved confinement regimes and investigating their linear and nonlinear stability. Understanding these issues can contribute to improving the current magnetic confinement systems and possibly developing new ones.

Work performed in year 2010 (*in co-operation with the Institutes indicated*):

- We have applied a sufficient criterion for linear stability [Throumoulopoulos and Tasso, Phys. Plasmas 14, 122104 (2007)] to the extended linear Herrnegger-Maschke equilibrium for incompressible flows parallel to the magnetic field of arbitrary shear. i.e. for , where is the Alfvén-Mach-function, labels the magnetic surfaces and is a free parameter associated with the flow shear. It turns out that, unlike the nonlinear equilibrium regime [Throumoulopoulos and Tasso, Phys. Plasmas, 17 032508 (2010)], ITER relevant flows () do not have a stabilising impact in the sense that the aforementioned stability criterion is not satisfied. This result supports a conjecture that equilibrium nonlinearity may activate flow stabilisation effects. For more details see
**Annex 25**. (I*n co-operation with IPP.*) - For field aligned incompressible flows a generalised Grad-Shafranov equation under an integral transformation can be put in a form identical to the respective usual (non flow) equation. This permits smooth extension of known non-flow solutions to the field aligned flow regime. In particular, we have extended a recently obtained Solovév-like equilibrium solution [Srinivasan et al., PPCF 52, 035007 (2010)]. This solution can have an arbitrary number of free parameters which have been exploited to construct configurations with ITER-like geometrical characteristics and confinement figures of merit. i.e. with a lower x-point, safety factor on axis equal to 1.1 and average toroidal beta equal to 0.01. Subsequently, application of the stability condition mentioned in part (i) led to similar conclusions regarding the insensitivity of the condition to the sheared flow. Details are given in
**Annex 25**. (I*n co-operation with IPP.*) - We have extended the generic linear static solution to the Grad-Shafranov equation [Atanasiu et al., Phys. Plasmas Vol.11, 3510 (2004)] to field aligned incompressible flows. The solution is expressed in terms of confluent hypergeometric functions and by using the principle of superposition it can contain an arbitrary number of free parameters. Configurations with either finite or vanishing current density on the boundary can be constructed. In particular, we pursued to construct ITER relevant equilibria by assigning appropriately the free parameters as in milestone (ii) in two ways: either by prescribing an ITER shape plasma boundary by a good number of points (about 30) or by using only few characteristic boundary points, i.e. the inner, outer, upper and lower ones, and imposing conditions not only on the solution but on its first and second derivatives in connection with appropriate shaping. Also, additional conditions associated with desirable values of confinement figures of merit led to a set of algebraic equations with integrals involving some of the unknown parameters. Numerical solutions of this set of equations required more time than anticipated. For this reason, completion of this study is planned for the year 2011 according to the approved work program. (I
*n co-operation with IPP and**MEC-EURATΟM Association.*) - Since application of the sufficient condition for linear stability mentioned in milestone (ii) to the extended generic linear equilibrium solution requires completion of the milestone (iii), this study was also moved to the work programme of year 2011. (I
*n co-operation with IPP and**MEC-EURATΟM Association.*) - We performed a study concerning some critical points on the drift kinetic and gyrokinetic theories in comparison with the Vlasov theory. Motivation was a previous paper [Tasso and Throumoulopoulos, J. Phys. A: Mathem. Theor. 40, F631 (2007)] in which an additional Vlasov constant of motion was found in the vicinity of the magnetic axis of an axisymmetric toroidal equilibrium. Since in the reduced phase space this constant of motion is missed it turns out that the drift kinetic and gyrokinetic current densities parallel to the magnetic field near axis can be remarkably different from the Vlasov ones. Also, unlikely Vlasov, the reduced phase space theories do not distinguish a toroidal parallel current density from a straight one near axis. These discrepancies, being independent of the particular forms of the drift kinetic and gyrokinetic equations, may put some limitation on the implications of computer simulations extensively performed in the framework of gyrokinetic theories, particularly as concerns the physics of energetic particles. For more details see
**Annex 26**. (I*n co-operation with IPP.*) - In view of common interest with IPP in the extension of the HELENA equilibrium code to incompressible flows, familiarisation with the current version of the code has been initiated. (I
*n co-operation with IPP.*)

3.4.2 Turbulence and Anomalous Transport Phenomena

Background and Objectives: This activity addresses several tasks:

In a first part of the activity, we study particle and heat transport in turbulent environments, with the use of realistic models of turbulence and in toroidal topologies, and with particular focus on the edge region, impurities, fast particles, and anomalous transport behaviour. Several tools are used in parallel in this study: (1) We perform test-particle simulations in turbulent electromagnetic field environments. (2) A code for the solution of the Fokker-Planck equation is developed, as an alternative tool to the test-particle simulations and the Langevin equation, and moreover for comparison to the ETS code. (3) Random walk models and fractional diffusion equations are developed as tools to describe anomalous transport, with particular aim to describe in parallel the combined diffusion in position and velocity space. (4) Eulerian-Lagrangian CFD models for the study of near-wall motion of particulate matter and its interaction with stochastic turbulent fluid flow in the presence of deterministic/stochastic magnetic fields.

In a second part of the activity, we explore particular turbulence models that are based on the concept of Self-Organised Criticality (SOC), where the turbulent plasma is considered as a complex system, in which localised instabilities are relaxed in small-scale diffusive events. SOC is usually modelled in the form of Cellular Automata and the sand-pile analogy is made use of, yet it is a particular aim of our approach to use the natural physical variables and that the SOC models have a consistent physical interpretation. In particular, we develop the Extended Cellular Automaton (X-CA) model for the magnetic field, which is fully MHD compatible, as well as a SOC model for turbulence driven by micro instabilities. In parallel, we build the MHD code MYDAS2, as an alternative to the X-CA model and in order to validate results of the latter.

The last part of this activity utilises Hamiltonian methods for obtaining evolution equations of the particle distribution functions under the presence of turbulent or other perturbing modes. These equations describe anomalous transport phenomena more accurately than standard theories, such as the quasilinear theory, since they incorporate the actual underlying particle dynamics and are not based on commonly used statistical assumptions with questionable validity in realistic cases.

Work performed in year 2010 (*in co-operation with the Institutes indicated*):

- Several direct numerical simulations (DNS) and large eddy simulations (LES) of MHD turbulent channel flows laden with ensembles of neutral particles of various properties have been conducted in order to account for situations of particulate impurities in liquid–metals and in less conducting fluids or in weakly ionised gases. For this reason, the Lagrangian particle tracking scheme has been developed further by including additional forces acting on particles. Several subgrid turbulence models and models that account for the effect of subgrid scales on the turbulent transport of particles (Langevin based models) have been developed and tested successfully. Finally, the discrete particle simulation module has also been modified properly in order to include electromagnetic forces on the determination of the trajectories of electrically charged particles. The effects of the changes of the fluid turbulence and its organised structures by the magnetic field on the turbulent diffusion of fluid tracers and the transport and deposition of particulate impurities have been investigated during year 2010 (see
**Annex 27**). As planned, a parametric numerical study has been initiated in the end of the year 2010 in order to continue with the numerical investigation of the turbulent diffusion/dispersion of tracers/charged particles in MHD turbulent flows. (*In co-operation with ULB*.) - No work was performed on this task in the period in subject, since priority was given to the development of the MHD code MYDAS2, which originated from the X-CA SOC model (with common way of calculating the spatial derivatives), and which we expect to yield highly useful information for further developing the X-CA model, once it will be running properly.
- A few trials were made to reduce the stiffness, in order to model also electrons (TEM mode driven turbulence), but they so-far all led to a loss of the SOC state, SOC has a tendency to be inherently stiff. Beyond this, no further work was performed on this task in the period in subject, since priority was given to the development of the Fokker-Planck code, with which we plan to investigate whether heat transport in the SOC model can be described with a Fokker-Planck approach. (in collaboration with Univ. of Bayreuth)
- Test-particle simulations were performed in the electromagnetic fields generated by the EMEDGE3D code for the MHD ballooning mode, and transport coefficients (diffusivities, pinch velocities, both in position and in energy), as well as energy distributions were determined, so-far for helium and electrons. Since the region modelled by the 3D EMEDGE3D code is very narrow in radial extent, the particles had to be re-injected when they leave, in order for the asymptotic regime to be reached. The electric fields turned out to be very strong (super-Dreicer), and particles get accelerated to high energies very fast, so that we need to re-examine the MHD simulations. Since no new simulations with EMEDGE3D could be made in the period in subject (the code is run by Univ. of Marseille), and the code MYDAS2 is still under development, test particle simulations of ions in an analytically given, vacuum-type equilibrium magnetic field were initiated, with a superimposed perturbation that models islands chains as resulting from NTMs. So-far, we find that there is an increased fraction of trapped particles near the region of the islands, with a characteristic shape of the trapped orbits. This work will be continued in the year 2011, with the determination of the relevant transport coefficients, and test-particle simulations in MHD fields will be done again as soon as the respective MHD simulations will be available. (For more details, see
**Annex 28**in collaboration with CEA, Univ. of Marseille, MEdC, Univ. of Craiova.) - The Fokker-Planck solver CHET1 was further tested by applying it to the analytically treatable case of quasi-linear diffusion, and very satisfying coincidence with the analytical results was found. Also, the effects of time-dependent transport parameters ware studied, addressing the question on how far initial transient phases influence the later evolution and asymptotic states. A third application of CHET1, to the modelling of the plasma response to ECRH/ICRH, was only initiated, with basically the elaboration of the relevant transport coefficients. The application of CHET1 to the heat transport as seen in the SOC model for ITG mode driven turbulence could not be addressed yet. Acquaintance with the use of the ETS was improved by participating in a respective code camp. Moreover, given the problems with the boundary conditions in the ETS, the basic solver (integrator) of CHET1 was ported to fit into the standardised programming environment of the ETS, so that it can be used as a new alternative solver module in the ETS, providing newly now the possibility to calculate the derivatives with the pseudospectral (Chebyshev) method in the ETS. The modifications of the solver have been completed, and its first use in the ETS is expected in the year 2011. For more details, see
**Annex 29**. (*In collaboration with Univ. of Bayreuth*.) - Starting from a general conservation law in integral form in state space (including position and velocity), a transport equation was derived that includes a deterministic part (Vlasov type equation), a classical stochastic part (Fokker-Planck type equation), and an anomalous stochastic part (with fractional derivatives), with the fractional part so-far for the case of decoupled dynamics in position and velocity space. We then developed a method to solve fractional differential equations numerically. Thereto, we used the left, right, and symmetric Grunwald Letnikov definition of fractional derivatives, which we implemented in the form of a matrix multiplication, in combination with a 4th order, adaptive step Runge - Kutta solver in time direction (method of lines). Alternatively, we calculated the fractional derivatives through Fourier transform methods, finding good coincidence with the Grunwald Letnikov approach. In this method, the unknown function must be transformed such that the initial conditions appear as a source term. A first parametric study in a set-up with off-axis heating has been made, for varying order of the fractional derivative, and slightly non-local effects have been identified. The new, physically meaningful coupling between position- and momentum space, which was developed for the transport equation in integral form, has not yet been tried to be introduced into the fractional diffusion equation. For more details, see
**Annex 30**. (*In collaboration with MEdC, Univ. of Craiova*.) - In order to further test the MHD code MYDAS, the MHD equations in polar coordinates were linearised, and analytical solutions of wave eigen-modes were derived, which were successfully reproduced with the non-linear, numerical code in the linear, low-amplitude regime. Letting MYDAS evolve into the full non-linear regime, numerical problems appeared due to the formation of small-scale structures, the aliasing caused by the non-linear terms, and also the boundary conditions needed to be modified such as to allow free outflow. To proceed further, the following steps were undertaken: (1) A version of MYDAS in Cartesian coordinates was built, with periodic boundary conditions (Fourier expansion in space), in order to have a code version where the boundary conditions cannot cause numerical problems. (2) In order to be able to treat larger grids with better resolution of small structures, (2a) the standard Runge - Kutta integrator was replaced with the Dormand Prince variant, (2b) the fast Fourier transform was replaced with the FFTW algorithm; both measures together led to a gain in computational time of a factor of roughly 30, and grid-sizes like 512512 are now efficiently treatable. (3) The code was de-aliased by spectral blocking of the high frequency modes prior to the evaluation of the non-linear terms (i.e. by applying the ‘2/3-rule’). In this way, the non-linear evolution could be followed until the formation of shocks in the test case of the Orszag Tang vortex. Turning then back to the case of polar coordinates, the construction of perfectly matched layers as free outflow boundary conditions was started and so-far successfully implemented in the case of a simple, radial advection equation with constant advection velocity. (For more details, see
**Annex 31**.) - The problem of interfacing between the output of the MHD-Turbulence code EMEDGE3D (developed in the Univesity of Marseille) and the module which reconstructs the density fluctuations at the plasma edge has been addressed. The statistics of these fluctuations led to the conclusion that the mildly higher blob density assumption in 1.5.1(viii) is barely valid. Additionally, data from TORPEX and NSTX experiments suggest that the relative blob density increment is in the range of 5 to 95%. Thus, for the study of refraction bu a single blob in this range of density increments, a Snell-Fresnel type of approach is more suitable which then could be extended to a multi-refraction statistical model. This study has been initiated in the second half of 2010 and it is currently in a preliminary phase.
- A Fokker-Planck equation with time-dependent, nonsingular diffusion tensor has been obtained for the case of a toroidal magnetic field when RF waves are also present, with the utilisation of canonical perturbation theory for Hamiltonian systems. The results have been extended in order to include not only steady state but also damped and growing waves. A numerical code has been implemented for the one-dimensional case. The code has been used for considering several cases of different wave spectra and particle distribution functions.
- The relationship between the particle phase space structure and the position and width of the magnetic islands has been studied. It has been shown that this relation is more complex for the case of particles having high energies and/or large Larmor radius, than for particles having low-energy and small Larmor radius.

3.4.3 Discrete kinetic and stochastic models for transport

Background and Objectives: Discrete kinetic simulations may provide an alternative mesoscale type of description of the transport properties of complex physical systems. The main unknown is a particle distribution function, which obeys a suitably chosen kinetic equation, while the bulk quantities of practical interest are obtained by taking moments of the distribution function. The benefit of this description, compared to the more traditional micro- and macroscopic approaches, is due to the fact that physics, at particle level, may be investigated with moderate computational effort. This activity aims to develop kinetic computer codes capable to simulate in a computationally efficient manner non-linear resistive MHD flows and vacuum systems in fusion devices.

Work performed in year 2010 (*in co-operation with the Institutes indicated*):

- The Direct Simulation Monte Carlo (DSMC) algorithm is a powerful computational tool for simulating high speed vacuum flows. Although the algorithm has been extensively implemented in single gas flows the corresponding work in gaseous mixtures is very limited. Here, the UTH “in house” DSMC code is accordingly modified to simulate binary gas flows at the vacuum systems of DT fusion reactors. In particular, almost all code subroutines (streaming, collision, indexing, and sampling) have been rewritten to include molecules of different species. Initial benchmarking indicates that the new binary DSMC code is running properly and it is free of errors. Various binary mixtures have been examined including mixtures consisting of species with about the same as well as quite different molecular masses such as Ne/Ar and He/Xe respectively. It has been found that a binary gas mixture may be simulated with reasonable accuracy by replacing the mixture with a single gas having an equivalent molecular mass only when the ratio of the molecular mass of the heavy over the light species is, roughly speaking, less than 10. Following this result the code will be used to solve binary gas flows of He and hydrogen isotopes through circular tubes.
- Modelling of gas-surface interaction in vacuum flows is one of the main issues in the accurate estimation of quantities of great practical interest such as conductance and pressure drop. All previous work has been based on the Maxwell diffuse-specular scattering law. It is known however that in certain cases there are discrepancies between computations and measurements in vacuum gas flows, which are attributed to the inadequacy of the Maxwell boundary conditions. Therefore, a more advanced and complete model namely the Cercignani-Lampis (CL) model has been developed and successfully implemented for flows through channels of various cross sections. Corresponding experimental work has been performed at the TRANSFLOW facility at KIT. Two different gases (nitrogen and helium) and two different pipe surfaces (rough and elektropolished) have been used. The results have been compared to the computational ones based on the CL boundary conditions and important concluding remarks on the effect of the molecular mass and of the surface characteristics on the parameters characterising the gas-surface interaction under any vacuum conditions have been obtained. Very briefly it seems that computed flow rates based on the Maxwell boundary conditions are in general overestimated about 10-20% (see
**Annex 32**). Comparisons between kinetic results (based on the DSMC algorithm) and corresponding experimental work have been also performed for flows through tubes with abrupt changes in their cross sectional areas. The vacuum piping systems of DT fusion reactors include such pipe elements with sudden enlargements and reductions. No available data exist in such flow configurations. The comparison shows in general good agreement in moderate pressure, which however degrades for low and very low pressures. This is not expected and at this stage is not conclusive if these discrepancies are attributed to computational errors or experimental uncertainties. Therefore still more work is needed in the short future to ensure the efficiency and accuracy of the developed DSMC code for this more complex geometry in the whole range of the Knudsen number (see**Annex 33**). Also, additional work not initially planed has been performed in the advancement of the nonlinear fully deterministic kinetic code. This more advanced version of the code, which includes (a) Wynn-epsilon acceleration, (b) Romberg acceleration, (c) grid refinement and adaptation, (d) full parallelisation and (e) efficient memory handling techniques has been successfully used to solve vacuum gas flows (see**Annex 34**). (*In co-operation with KIT, CEA.*) - The development of a code simulating simple gas piping distribution systems consisting of long pipes in the whole range of the Knudsen number, based on kinetic theory, has been completed. The objective is to compute the flow quantities, i.e. the mass flow rates (or conductance) through each tube and the pressure head at each node of the network, provided that the geometry, i.e., the length and the diameter of each channel of the network is given. The system of equations describing the network consists of the pressure drop equations along each tube and the mass conservation equations at each node of the network. When the network is properly defined the pressure drop equations are reduced to a set of energy balance equations along the closed loops of the network, which along with the mass conservation equations form a closed set, which may be solved for the mass flow rates. Then, the pressure heads at the nodes are estimated through the pressure drop equations. An iterative process between the system of equations and the pressure drop equations yields the converged values of the unknown quantities. The main problem, in the present rarefied conditions, compared to the typical network solvers in the hydrodynamic limit, is that in the later case the pressure drop along a channel is given by explicit algebraic expressions, while far from this limit no such expressions exist. Here, this information is obtained from a data base, which has been developed for this purpose by solving a linearised BGK kinetic equation in a wide range of the Knudsen number and obtaining the corresponding data. Based on the above modelling several pipe networks under vacuum conditions have been solved. It is noted that this algorithm is the first one to simulate in an exact manner gas vacuum systems based on kinetic principals and it is considered as a significant advancement tool in vacuum technology (see
**Annex 35**).

3.4.4 2-D MHD code for plasma temporal evolution in fusion devices including specific modules on magnetic field topology, neutron-pellet interaction and “Development of an alternative source and neutraliser technique aiming at increased efficiency of neutral beam heating systems”

Background and Objectives: The main aim of our investigation is to study a compact magnetic fusion device in open magnetic topology for high neutron flux production. The last few years our work concerns numerical calculations of the spatio-temporal evolution of a high density and temperature plasma in a high external applied magnetic field with mirror-like topology, using a 2-D MHD resistive code. A new idea to improve the initial ‘plasma volume’ and increase the neutron production is based on the effect of the self-guiding propagation of an ultrashort laser beam in the plasma. The laser propagation produces a relatively long filament with an important number of multi-focal spots (or hot spots) due to the non-linear self-guiding laser-plasma interaction. Experimental results show that laser-cluster interaction produce positive and negative ions of hydrogen and deuterium with initial kinetic energies up to 30 keV. The photo-detachment of negative ion beam, using appropriate laser wavelengths, allows efficient generation of high energy neutral beams for plasma heating applications to Tokamaks. The objectives are: (1) Numerical investigations, using the 2-D MHD resistive code in axisymmetric cylindrical geometry, to study (a) the spatio-temporal evolution of the state parameters of the plasma produced by laser filamentation due to the non-linear interaction with the plasma, and (b) the effects of the ratio *B*_{max}/*B*_{min} of the magnetic mirror on the plasma trapping. (2) Investigations on the physical process of pellet interaction with neutrons and estimation of the energy losses. (3) Initiation of studies for neutralisation of high kinetic energy negative ions (H and/or D), employing photo-detachment techniques, for Tokamak application.

Work performed in year 2010 (*in **collaboration with Ecole Polytechnique-Paris, CEA-Saclay-Paris, Bordeaux-France*):

(i) We investigate for the first time the filamentation effects during the propagation of an ultrashort high intensity laser beam with deutereted clusters and its potential application to compact magnetic fusion devices. We consider a double-plasma spot, describing the initial plasma distribution with a plasma density up to 10^{19} cm^{–3} and temperature of 20 keV for the self-guided region and a plasma density of 10^{17} cm^{–3} for the surrounding the multi-plasma spot region. For these conditions of laser-plasma interaction the formation of multi-plasma spots (multi-focal-spots) of high initial plasma density distribution (see the photo in **Annex 36**) was observed. We run the 2-D MHD resistive code with a high external applied magnetic field in mirror-like topology for two different cases. The first corresponds to the double-plasma spot and the second to the single focus case. The main result is that the plasma density in the double-plasma spot case corresponds to a higher density by an order of magnitude compared to the single focus case. For the double-plasma spot case we observe an oscillation of the plasma density producing 2.5 times longer trapping time with high density compare to the single focus case (see **Annex 36**). The results confirm both the efficiency of the proposed self-guiding propagation configuration and the improved neutron production.

The plasma trapping in the mirror-like magnetic configuration is very sensitive to the ratio *B*_{max}/*B*_{min}. We run the 2-D MHD for different mirror-ratio and distance between the coils producing the external high magnetic field. The diameter of the single-turn-double-coil is 1 cm in order to achieve an external magnetic field of the order of 120 Tesla and remains unchanged for all runs. A value of the mirror-ratio up to 3 allows important reduction of the plasma expansion and losses in the axial direction due to high magnetic field. Values for the mirror-ratio higher than 3 enable efficient plasma ‘reflection’ in the axial direction reducing the losses but the losses are very important in the radial direction. Under these conditions the simulations show that the trapping time is less than 100 nanoseconds which decrease the neutron production. By increasing the distance between the coils by a factor of 2.6 the losses in the radial direction decrease. The best values for the mirror-ratio, *B*_{max}/*B*_{min}, of the proposed magnetic fusion device can be varied from 2 to 3.

The kinetic energy losses during inelastic collision of neutrons with protons (or hydrogen or deuterium) are very important due to the equal value of their mass. Experimental data show that energetic neutrons loss their initial energy after 4-5 inelastic collisions with protons. In the case of the ‘bulk’ collision the transfer of energy could be increased by a factor of 2-2.5. The diameter of the pellet considered for the neutron energy transfer evaluation is few mm (2-6 mm). For a neutron flux up to 10^{15}-10^{16} n/cm^{2} the transfer of energy to the pellet is 1/200 less than the corresponding energy transfers due to the local plasma temperature. By increasing the neutron flux and/or the dimension of the pellet the transfer of energy increases but the value remains low compared to the plasma energy transfer.

(ii) Negative Ions: Production and Photo-Neutralisation. We investigate a new source of negative ion (H and/or D) production, from laser –cluster and/or laser-deuterteted gas interaction. The main advantages are the high density of the produced negative ions and the relatively small volume of production. Such a negative ion beam can be post accelerated in a magnetically isolated diode using a pulsed power generator working at 1 MeV and producing currents up 500 A – 1 kA. For the design and details of the proposed configuration see our presentation during the CCNB Sofia Workshop on “Negative ions Alternative Source” (2-4 June 2010, Sofia, Bulgaria) [http/efdasql.ipp.mpg.de/ccnb]. Moreover, we propose a laser “cavity” of approximately 1.5 -2 m to neutralise the pulsed negative ion beam of 1 MeV and 1 kA by photo-detachment. Three different laser systems were considered: the Nd:Yag (1.05 μm), the Ti:Sapphire (800 nm) and the KrF Excimer (248 nm). The energy of the laser systems is up to 20-25 J The efficiencies for the Nd:Yag and the Ti:Sapphire are very close to 99% and for the KrF is 37% for the same dimension of the laser neutraliser. The efficiency of the KrF-Excimer negative ion neutraliser can be improved by increasing the dimension of the neutraliser by a factor of three.

Τελευταία Ενημέρωση (Κυριακή, 01 Απρίλιος 2012 18:35)